Skip to content

Commit

Permalink
Refactor: optimized cubic spline (#4081)
Browse files Browse the repository at this point in the history
* redesign eval to enable potential vectorization

* fix interface

* refactored version init commit

* clean up

* clean up

* fix & a few tests

* theoretical error bound test

* clean up

* better unit test for cross-check

* systematic cross-check

* minor change

* clean up

* clean up

* minor internal change

* add explanation

---------

Co-authored-by: Mohan Chen <[email protected]>
  • Loading branch information
jinzx10 and mohanchen committed May 5, 2024
1 parent 378d0e8 commit 48f2b5d
Show file tree
Hide file tree
Showing 16 changed files with 1,487 additions and 1,872 deletions.
711 changes: 454 additions & 257 deletions source/module_base/cubic_spline.cpp

Large diffs are not rendered by default.

749 changes: 525 additions & 224 deletions source/module_base/cubic_spline.h

Large diffs are not rendered by default.

1,755 changes: 373 additions & 1,382 deletions source/module_base/test/cubic_spline_test.cpp

Large diffs are not rendered by default.

7 changes: 7 additions & 0 deletions source/module_base/test/data/cos_periodic.dat
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
periodic periodic
0.00000000000000000e+00 1.00000000000000006e-01 1.00000000000000000e+00 1.50000000000000000e+00 3.00000000000000000e+00 6.28318530717958623e+00
1.00000000000000000e+00 9.95004165278025821e-01 5.40302305868139765e-01 7.07372016677029064e-02 -9.89992496600445415e-01 1.00000000000000000e+00
0.00000000000000000e+00 1.74532925199432948e-01 3.49065850398865896e-01 5.23598775598298816e-01 6.98131700797731791e-01 8.72664625997164767e-01 1.04719755119659763e+00 1.22173047639603061e+00 1.39626340159546358e+00 1.57079632679489656e+00 1.74532925199432953e+00 1.91986217719376251e+00 2.09439510239319526e+00 2.26892802759262846e+00 2.44346095279206121e+00 2.61799387799149441e+00 2.79252680319092716e+00 2.96705972839035992e+00 3.14159265358979312e+00 3.31612557878922587e+00 3.49065850398865907e+00 3.66519142918809182e+00 3.83972435438752502e+00 4.01425727958695777e+00 4.18879020478639053e+00 4.36332312998582328e+00 4.53785605518525692e+00 4.71238898038468967e+00 4.88692190558412243e+00 5.06145483078355518e+00 5.23598775598298882e+00 5.41052068118242158e+00 5.58505360638185433e+00 5.75958653158128708e+00 5.93411945678071984e+00 6.10865238198015348e+00 6.28318530717958623e+00
1.00000000000000000e+00 9.84381299886306516e-01 9.38016522020507670e-01 8.63713424115361450e-01 7.64165387715657451e-01 6.42065794366185383e-01 5.00118352188328785e-01 3.42056178626561569e-01 1.73503016854116870e-01 2.70592234453568648e-04 -1.72386293431320586e-01 -3.39925662651143046e-01 -4.97852121855354546e-01 -6.41670277474296391e-01 -7.66884735938309081e-01 -8.69000103677733948e-01 -9.43520987122911547e-01 -9.85951992704182878e-01 -9.92542393974230608e-01 -9.64850207988183883e-01 -9.06738697819773454e-01 -8.22080502627916765e-01 -7.14748261571530152e-01 -5.88614613809530951e-01 -4.47552198500836496e-01 -2.95433654804362900e-01 -1.36131621879026887e-01 2.64812611162533740e-02 1.88532355022561493e-01 3.46149020680980524e-01 4.95458618932594463e-01 6.32588510618484978e-01 7.53666056579735066e-01 8.54818617657429058e-01 9.32173554692648398e-01 9.81858228526478083e-01 1.00000000000000000e+00
6.18473891098361445e-03 -1.80469237892337431e-01 -3.48259990143118536e-01 -5.00618802880412384e-01 -6.37545676104219083e-01 -7.59040609814538803e-01 -8.64447219749984619e-01 -9.41247940672875094e-01 -9.84665089228650259e-01 -9.95326120656199964e-01 -9.78883265201603403e-01 -9.36683701867835450e-01 -8.68727430654895882e-01 -7.75014451562784812e-01 -6.55544764591502682e-01 -5.10318369741048716e-01 -3.39335267011423580e-01 -1.42595456402626830e-01 6.41233977052624166e-02 2.49507580803381257e-01 4.12703175963197166e-01 5.53710183184709215e-01 6.72528602467918124e-01 7.69158433812823450e-01 8.43599677219425303e-01 8.95852332687724018e-01 9.25916400217719260e-01 9.33791879809410919e-01 9.19478771462799216e-01 8.82977075177883819e-01 8.24286790954665616e-01 7.43407918793143718e-01 6.40340458693318570e-01 5.15084410655190172e-01 3.67639774678758080e-01 1.98006550764022293e-01 6.18473891098361445e-03
-1.16262364588415146e+00 -1.00557944472069671e+00 -9.17160944338298312e-01 -8.28742443955899910e-01 -7.40323943573501397e-01 -6.51905443191102996e-01 -5.35672607329632511e-01 -3.44398827159106469e-01 -1.53125046988580316e-01 2.04230892878108500e-02 1.67998154850598658e-01 3.15573220413386479e-01 4.63148285976174190e-01 6.10723351538962067e-01 7.58298417101749722e-01 9.05873482664537821e-01 1.05344854822732525e+00 1.20102361379011291e+00 1.12573875011125457e+00 9.98607505889287705e-01 8.71476261667320506e-01 7.44345017445353752e-01 6.17213773223386442e-01 4.90082529001419687e-01 3.62951284779452932e-01 2.35820040557485955e-01 1.08688796335518534e-01 -1.84424478864484431e-02 -1.45573692108415198e-01 -2.72704936330382175e-01 -3.99836180552349374e-01 -5.26967424774316351e-01 -6.54098668996283328e-01 -7.81229913218249861e-01 -9.08361157440217282e-01 -1.03549240166218426e+00 -1.16262364588415146e+00
7 changes: 7 additions & 0 deletions source/module_base/test/data/exp_first_deriv.dat
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
first_deriv:1.10517091807564771e+00 first_deriv:2.71828182845904509e+00
1.00000000000000006e-01 1.29154966501488389e-01 1.66810053720005874e-01 2.15443469003188337e-01 2.78255940220712428e-01 3.59381366380462752e-01 4.64158883361277863e-01 5.99484250318940926e-01 7.74263682681126997e-01 1.00000000000000000e+00
1.10517091807564771e+00 1.13786644168535322e+00 1.18152981679905356e+00 1.24041185922219421e+00 1.32082420599103290e+00 1.43244298301307027e+00 1.59067568226275657e+00 1.82117928549843455e+00 2.16899447087117636e+00 2.71828182845904509e+00
1.00000000000000006e-01 1.25000000000000000e-01 1.50000000000000022e-01 1.75000000000000017e-01 2.00000000000000011e-01 2.25000000000000006e-01 2.50000000000000000e-01 2.75000000000000022e-01 3.00000000000000044e-01 3.25000000000000011e-01 3.49999999999999978e-01 3.75000000000000000e-01 4.00000000000000022e-01 4.25000000000000044e-01 4.50000000000000067e-01 4.74999999999999978e-01 5.00000000000000000e-01 5.25000000000000022e-01 5.50000000000000044e-01 5.75000000000000067e-01 5.99999999999999978e-01 6.25000000000000000e-01 6.50000000000000022e-01 6.75000000000000044e-01 7.00000000000000067e-01 7.24999999999999978e-01 7.50000000000000000e-01 7.75000000000000022e-01 8.00000000000000044e-01 8.25000000000000067e-01 8.49999999999999978e-01 8.75000000000000000e-01 9.00000000000000022e-01 9.25000000000000044e-01 9.50000000000000067e-01 9.74999999999999978e-01 1.00000000000000000e+00
1.10517091807564771e+00 1.13314845305620615e+00 1.16183423781052575e+00 1.19124621006218145e+00 1.22140274895196677e+00 1.25232269983173161e+00 1.28402537693915098e+00 1.31653067826587877e+00 1.34985870568413047e+00 1.38403052288246253e+00 1.41906754700129101e+00 1.45499127838897557e+00 1.49182434073154035e+00 1.52959016783122848e+00 1.56831221140395694e+00 1.60801396407294228e+00 1.64872023299870918e+00 1.69045738865219253e+00 1.73325189263525137e+00 1.77713020654974518e+00 1.82211879200369298e+00 1.86824485646834759e+00 1.91553841190619289e+00 1.96403012931547272e+00 2.01375067969443089e+00 2.06473073404130947e+00 2.11700096335435450e+00 2.17059203866204742e+00 2.22553592214226725e+00 2.28186930422749290e+00 2.33962995740218638e+00 2.39885565415081015e+00 2.45958416695782534e+00 2.52185326830769485e+00 2.58570073068487982e+00 2.65116432657384271e+00 2.71828182845904465e+00
1.10517091807564771e+00 1.13314857626511545e+00 1.16183448755634799e+00 1.19124518701385762e+00 1.22140376069562406e+00 1.25232046400808739e+00 1.28402663663421768e+00 1.31653040155273815e+00 1.34985411865053440e+00 1.38403404177963485e+00 1.41907067229022621e+00 1.45497988810568324e+00 1.49182144407227635e+00 1.52960106867554990e+00 1.56831876191550346e+00 1.60798584383446830e+00 1.64869208090593933e+00 1.69045678205178129e+00 1.73327994727199441e+00 1.77716157656657869e+00 1.82210170576561081e+00 1.86818792472837214e+00 1.91550098362087740e+00 1.96404088244312680e+00 2.01380762119512058e+00 2.06480119987685740e+00 2.11702161848833859e+00 2.17046900023482614e+00 2.22529349474582894e+00 2.28162888863530977e+00 2.33947518190326731e+00 2.39883237454970333e+00 2.45970046657461650e+00 2.52207945797800726e+00 2.58596934875987605e+00 2.65137013892022200e+00 2.71828182845904553e+00
1.10510282004738825e+00 1.13310983511001884e+00 1.16189187653084680e+00 1.19124256407428741e+00 1.22144333046701914e+00 1.25229505919741335e+00 1.28419875089301683e+00 1.31610244258862052e+00 1.35006277753420334e+00 1.38433107279383738e+00 1.41859936805347142e+00 1.45490086593012635e+00 1.49242361139733326e+00 1.52994635686453995e+00 1.56746910233174686e+00 1.60708020137141450e+00 1.64941876434625834e+00 1.69175732732110240e+00 1.73409589029594624e+00 1.77643445327079030e+00 1.81891195991557297e+00 1.86798555710533210e+00 1.91705915429509122e+00 1.96613275148485056e+00 2.01520634867460968e+00 2.06427994586436858e+00 2.11335354305412793e+00 2.16276179287056580e+00 2.22319776800966951e+00 2.28363374314877277e+00 2.34406971828787603e+00 2.40450569342697928e+00 2.46494166856608299e+00 2.52537764370518625e+00 2.58581361884428951e+00 2.64624959398339277e+00 2.70668556912249603e+00
64 changes: 64 additions & 0 deletions source/module_base/test/data/gen_ref.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
'''
This script generates reference data for the cross-check in cubic_spline_test.cpp.
To add more tests, append a new Case object to the cases list.
'''
import numpy as np
from scipy.interpolate import CubicSpline

class Case:
def __init__(self, f, x, bc, fname):
self.x = x
self.y = f(x)
self.bc = bc
self.fname = fname
if bc == 'periodic':
self.y[-1] = self.y[0]

cases = [
Case(np.sin, np.logspace(-1, 0, 10), 'not-a-knot',
'sin_not_a_knot.dat'),
Case(np.cos, [0, 0.1, 1, 1.5, 3, 2*np.pi], 'periodic',
'cos_periodic.dat'),
Case(np.exp, np.logspace(-1, 0, 10), ((1, np.exp(0.1)), (1, np.exp(1))),
'exp_first_deriv.dat'),
Case(np.log, np.logspace(0, 1, 10), ((2, -1), (2, -0.01)),
'log_second_deriv.dat'),
Case(np.sqrt, np.logspace(0, 1, 10), ((1, 0.5), 'not-a-knot'),
'sqrt_mix_bc.dat'),
Case(lambda x: np.ones_like(x), [0, 1], 'periodic',
'two_points_periodic.dat'),
Case(np.arccos, [0, 0.5], ((1, -1), (1, -2/np.sqrt(3))),
'two_points_first_deriv.dat'),
Case(np.sqrt, [1, 4], ((2, 0.5), (2, 0.25)),
'two_points_second_deriv.dat'),
Case(np.sin, [0, 0.3, 1], 'not-a-knot',
'three_points_not_a_knot.dat'),
]

n_interp = 37

for case in cases:
cubspl = CubicSpline(case.x, case.y, bc_type=case.bc)
x_interp = np.linspace(case.x[0], case.x[-1], n_interp)
y_interp = cubspl(x_interp)
dy_interp = cubspl(x_interp, 1)
d2y_interp = cubspl(x_interp, 2)

with open(case.fname, 'w') as f:
if case.bc == 'periodic' or case.bc == 'not-a-knot':
f.write('{bc} {bc}\n'.format(bc=case.bc))
else:
for b in case.bc:
if b == 'not-a-knot':
f.write('not-a-knot ')
else:
tag = 'first_deriv' if b[0] == 1 else 'second_deriv'
f.write('{tag}:{val:<24.17e} '.format(tag=tag, val=b[1]))
f.write('\n')

for data in [case.x, case.y, x_interp, y_interp, dy_interp, d2y_interp]:
for elem in data:
f.write('% 24.17e '%(elem))
f.write('\n')

7 changes: 7 additions & 0 deletions source/module_base/test/data/log_second_deriv.dat
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
second_deriv:-1.00000000000000000e+00 second_deriv:-1.00000000000000002e-02
1.00000000000000000e+00 1.29154966501488389e+00 1.66810053720005880e+00 2.15443469003188381e+00 2.78255940220712450e+00 3.59381366380462763e+00 4.64158883361277841e+00 5.99484250318940859e+00 7.74263682681126930e+00 1.00000000000000000e+01
0.00000000000000000e+00 2.55842788110449526e-01 5.11685576220899052e-01 7.67528364331348634e-01 1.02337115244179810e+00 1.27921394055224780e+00 1.53505672866269705e+00 1.79089951677314629e+00 2.04674230488359621e+00 2.30258509299404590e+00
1.00000000000000000e+00 1.25000000000000000e+00 1.50000000000000000e+00 1.75000000000000000e+00 2.00000000000000000e+00 2.25000000000000000e+00 2.50000000000000000e+00 2.75000000000000000e+00 3.00000000000000000e+00 3.25000000000000000e+00 3.50000000000000000e+00 3.75000000000000000e+00 4.00000000000000000e+00 4.25000000000000000e+00 4.50000000000000000e+00 4.75000000000000000e+00 5.00000000000000000e+00 5.25000000000000000e+00 5.50000000000000000e+00 5.75000000000000000e+00 6.00000000000000000e+00 6.25000000000000000e+00 6.50000000000000000e+00 6.75000000000000000e+00 7.00000000000000000e+00 7.25000000000000000e+00 7.50000000000000000e+00 7.75000000000000000e+00 8.00000000000000000e+00 8.25000000000000000e+00 8.50000000000000000e+00 8.75000000000000000e+00 9.00000000000000000e+00 9.25000000000000000e+00 9.50000000000000000e+00 9.75000000000000000e+00 1.00000000000000000e+01
0.00000000000000000e+00 2.23201685990528503e-01 4.05482161721055334e-01 5.59647381157334922e-01 6.93212288502453533e-01 8.10938183318088091e-01 9.16348065656249267e-01 1.01160453727708877e+00 1.09865227971226798e+00 1.17872100366669685e+00 1.25277767969836296e+00 1.32176777395878364e+00 1.38635002079692837e+00 1.44697637249515854e+00 1.50409420887274625e+00 1.55814764821316754e+00 1.60947600084362952e+00 1.65829393447507156e+00 1.70480885095700785e+00 1.74922815213895189e+00 1.79175923970966466e+00 1.83259005049654222e+00 1.87183533226151133e+00 1.90959263384970646e+00 1.94595950410626362e+00 1.98103349187631750e+00 2.01491214600500346e+00 2.04769301510038781e+00 2.07946352550178926e+00 2.11027403528676105e+00 2.14016641953235709e+00 2.16918255331563126e+00 2.19736431171363611e+00 2.22475356980342642e+00 2.25139220266205564e+00 2.27732208536657810e+00 2.30258509299404635e+00
1.00254261613056794e+00 7.98334999625206598e-01 6.66752209256807338e-01 5.71952135021786212e-01 4.99597454958138398e-01 4.44699888470842797e-01 3.99955939076224509e-01 3.63472602732268946e-01 3.33587746439520605e-01 3.07606422584049077e-01 2.85491363057420178e-01 2.66837284058103641e-01 2.50118943859901521e-01 2.35190122938787882e-01 2.22050821294762640e-01 2.10610784319237826e-01 2.00154304624132307e-01 1.90527432327081331e-01 1.81730167428084816e-01 1.73762509927142816e-01 1.66624366318174705e-01 1.60087152540269423e-01 1.53940134142904944e-01 1.48183311126081213e-01 1.42816683489798313e-01 1.37840251234056188e-01 1.33254014358854894e-01 1.29057876274537020e-01 1.25134123654711077e-01 1.21377871343099950e-01 1.17789119339703638e-01 1.14367867644522142e-01 1.11114116257555462e-01 1.08027865178803598e-01 1.05109114408266549e-01 1.02357863945944316e-01 9.97741137918368992e-02
-1.00000000000000377e+00 -6.33660932042887648e-01 -4.49235132334255161e-01 -3.25782694882304436e-01 -2.53054745626878019e-01 -1.95497023679798904e-01 -1.62454571477147702e-01 -1.29412119274496529e-01 -1.11657824079571477e-01 -9.61927667642009099e-02 -8.07277094488303704e-02 -7.04523993469851478e-02 -6.32943222386314758e-02 -5.61362451302777968e-02 -4.89781680219241247e-02 -4.34851335765309294e-02 -4.01667039843129556e-02 -3.68482743920949818e-02 -3.35298447998770149e-02 -3.02114152076590411e-02 -2.69292458727027513e-02 -2.53684643505395734e-02 -2.38076828283763989e-02 -2.22469013062132209e-02 -2.06861197840500430e-02 -1.91253382618868650e-02 -1.75645567397236871e-02 -1.60300110957334006e-02 -1.53600098628741376e-02 -1.46900086300148763e-02 -1.40200073971556132e-02 -1.33500061642963519e-02 -1.26800049314370888e-02 -1.20100036985778258e-02 -1.13400024657185645e-02 -1.06700012328593032e-02 -1.00000000000000401e-02
7 changes: 7 additions & 0 deletions source/module_base/test/data/sin_not_a_knot.dat
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
not-a-knot not-a-knot
1.00000000000000006e-01 1.29154966501488389e-01 1.66810053720005874e-01 2.15443469003188337e-01 2.78255940220712428e-01 3.59381366380462752e-01 4.64158883361277863e-01 5.99484250318940926e-01 7.74263682681126997e-01 1.00000000000000000e+00
9.98334166468281548e-02 1.28796193418703991e-01 1.66037531159675150e-01 2.13780666055298912e-01 2.74679090976688467e-01 3.51695188663471381e-01 4.47670834718957189e-01 5.64216731736942645e-01 6.99189845133651677e-01 8.41470984807896505e-01
1.00000000000000006e-01 1.25000000000000000e-01 1.50000000000000022e-01 1.75000000000000017e-01 2.00000000000000011e-01 2.25000000000000006e-01 2.50000000000000000e-01 2.75000000000000022e-01 3.00000000000000044e-01 3.25000000000000011e-01 3.49999999999999978e-01 3.75000000000000000e-01 4.00000000000000022e-01 4.25000000000000044e-01 4.50000000000000067e-01 4.74999999999999978e-01 5.00000000000000000e-01 5.25000000000000022e-01 5.50000000000000044e-01 5.75000000000000067e-01 5.99999999999999978e-01 6.25000000000000000e-01 6.50000000000000022e-01 6.75000000000000044e-01 7.00000000000000067e-01 7.24999999999999978e-01 7.50000000000000000e-01 7.75000000000000022e-01 8.00000000000000044e-01 8.25000000000000067e-01 8.49999999999999978e-01 8.75000000000000000e-01 9.00000000000000022e-01 9.25000000000000044e-01 9.50000000000000067e-01 9.74999999999999978e-01 1.00000000000000000e+00
9.98334166468281548e-02 1.24674734514712698e-01 1.49438130358688637e-01 1.74108137455688206e-01 1.98669331470162691e-01 2.23106356655036864e-01 2.47403944178155422e-01 2.71546935540866008e-01 2.95520208213219726e-01 3.19308806225840347e-01 3.42897834817652569e-01 3.66272425041630545e-01 3.89418057463319811e-01 4.12320464705872358e-01 4.34965384966057578e-01 4.57338563261367681e-01 4.79425963790102938e-01 5.01213811410824883e-01 5.22688346176908625e-01 5.43835808141729826e-01 5.64642437361874472e-01 5.85094862720930697e-01 6.05181175115716008e-01 6.24889809005897301e-01 6.44209198851141918e-01 6.63127779111117310e-01 6.81633984245490487e-01 6.99716248713929012e-01 7.17363006976099671e-01 7.34562693491669916e-01 7.51303742720306977e-01 7.67574589121678086e-01 7.83363667155450361e-01 7.98659411281291143e-01 8.13450255958867552e-01 8.27724635647847040e-01 8.41470984807896505e-01
9.95004923184505952e-01 9.92197390241733368e-01 9.88771161271814858e-01 9.84726484485526354e-01 9.80066691602079154e-01 9.74793416196813478e-01 9.68912581684626351e-01 9.62425723364175201e-01 9.55336713363020729e-01 9.47651829867646733e-01 9.39371159698343727e-01 9.30499643073893545e-01 9.21055871973040285e-01 9.11041628742959020e-01 9.00456913383649749e-01 8.89303613354681888e-01 8.77596695966640139e-01 8.65339380713611495e-01 8.52531667595595843e-01 8.39173556612593186e-01 8.25265066443230833e-01 8.10851858679044746e-01 7.95976029301547761e-01 7.80637578310739877e-01 7.64836505706621206e-01 7.48572811489191747e-01 7.31846495658451279e-01 7.14657558214399802e-01 6.97005999157037537e-01 6.78891818486364262e-01 6.60315016202380312e-01 6.41275592305085462e-01 6.21773546794479715e-01 6.01808879670563179e-01 5.81381590933335857e-01 5.60491680582797747e-01 5.39139148618948627e-01
-9.99273971679864453e-02 -1.24675238253821624e-01 -1.49423079339714465e-01 -1.74110306895194078e-01 -1.98673123780586686e-01 -2.23112904322202821e-01 -2.47353856652764159e-01 -2.71594808983325553e-01 -2.95479606336378597e-01 -3.19311073293541137e-01 -3.43142540250703731e-01 -3.66341401449573223e-01 -3.89160286618692008e-01 -4.11979171787810738e-01 -4.34798056956929524e-01 -4.57268738221928150e-01 -4.79284652821407420e-01 -5.01300567420886689e-01 -5.23316482020365958e-01 -5.45332396619845228e-01 -5.67275878301223169e-01 -5.85780742833658996e-01 -6.04285607366094824e-01 -6.22790471898530651e-01 -6.41295336430966478e-01 -6.59800200963402195e-01 -6.78305065495838022e-01 -6.96809930028275515e-01 -7.15314794560709233e-01 -7.33819659093142951e-01 -7.52324523625576669e-01 -7.70829388158010387e-01 -7.89334252690444216e-01 -8.07839117222877934e-01 -8.26343981755311652e-01 -8.44848846287745370e-01 -8.63353710820179088e-01
Loading

0 comments on commit 48f2b5d

Please sign in to comment.