SINTEF / Splipy

Spline modelling made easy.
GNU General Public License v3.0
102 stars 18 forks source link

Greville points appear outside BSpline domain #44

Closed VikingScientist closed 6 years ago

VikingScientist commented 6 years ago
from splipy import *                                                           
import numpy as np                                                                

b = BSplineBasis(4, [ 0.        , 0.        , 0.        , 0.        , 0.01092385, 0.02200375,\
                      0.03316167, 0.04435861, 0.05526295, 0.0663331 , 0.07748615, 0.08868065,\
                      0.09989588, 0.11080936, 0.12188408, 0.13303942, 0.14423507, 0.15545087,\
                      0.16636463, 0.1774395 , 0.1885949 , 0.19979059, 0.2       , 0.2       ,\  
                      0.2       , 0.2     ]);                                     
print(b.greville())                                                               
print(b.knots[-1])                                                                

Prints the following:

[0.0, 0.0036412833333333335, 0.010975866666666667, 0.022029756666666667, 0.03317467666666666, 0.04426107666666667, 0.05531822, 0.06636073333333334, 0.07749996666666667, 0.08868756, 0.09979529666666669, 0.11086310666666667, 0.12191095333333334, 0.13305285666666666, 0.14424178666666665, 0.15535019, 0.16641833333333333, 0.17746634333333333, 0.18860833, 0.19612849666666668, 0.19993019666666667, 0.20000000000000004]

0.2

Notice the last point greville point is outside the BSpline basis domain [0, 0.2]. This produces all sorts of problems with for instance rebuild calls and other interpolation schemes.

VikingScientist commented 6 years ago

I guess this is a stupid floating point precision error. The last entry is 0.2+4e-17. One solution would be to snap to nearest knot as defined by state.knot_tolerance. The question then becomes where do we put this snapping mechanic in order to produce the least amount of overhead, maintenance and maximum readability?