1 module mir.optim.fit_spline; 2 3 import mir.optim.least_squares; 4 import mir.interpolate.spline; 5 6 /// 7 struct FitSplineResult(T) 8 { 9 /// 10 LeastSquaresResult!T leastSquaresResult; 11 /// 12 Spline!T spline; 13 } 14 15 /++ 16 Params: 17 settings = LMA settings 18 points = points to fit 19 x = fixed X values of the spline 20 l = lower bounds for spline(X) values 21 u = upper bounds for spline(X) values 22 lambda = coefficient for the integral of the square of the second derivative 23 configuration = spline configuration (optional) 24 Returns: $(FitSplineResult) 25 +/ 26 FitSplineResult!T fitSpline(alias d = "a - b", T)( 27 scope const ref LeastSquaresSettings!T settings, 28 scope const T[2][] points, 29 scope const T[] x, 30 scope const T[] l, 31 scope const T[] u, 32 const T lambda = 0, 33 SplineConfiguration!T configuration = SplineConfiguration!T(), 34 ) @nogc @trusted pure 35 if ((is(T == float) || is(T == double))) 36 in (lambda >= 0) 37 { 38 pragma(inline, false); 39 40 import mir.functional: naryFun; 41 import mir.math.common: sqrt, fabs; 42 import mir.ndslice.slice: sliced, Slice; 43 import mir.rc.array; 44 45 if (points.length < x.length && lambda == 0) 46 { 47 static immutable exc = new Exception("fitSpline: points.length has to be greater or equal x.length when lambda is 0.0"); 48 throw exc; 49 } 50 51 FitSplineResult!T ret; 52 53 ret.spline = x.rcarray!(immutable T).moveToSlice.Spline!T; 54 55 auto y = x.length.RCArray!T; 56 y[][] = 0; 57 58 scope f = delegate(scope Slice!(const (T)*) splineY, scope Slice!(T*) y) 59 { 60 assert(y.length == points.length + !lambda); 61 ret.spline._values = splineY; 62 with(configuration) 63 ret.spline._computeDerivatives(kind, param, leftBoundary, rightBoundary); 64 foreach (i, ref point; points) 65 y[i] = naryFun!d(ret.spline(point[0]), point[1]); 66 67 T integral = 0; 68 if (lambda) 69 { 70 T ld = ret.spline.withTwoDerivatives(x[0])[1]; 71 foreach (i; 1 .. x.length) 72 { 73 T rd = ret.spline.withTwoDerivatives(x[i])[1]; 74 integral += (rd * rd + rd * ld + ld * ld) * (x[i] - x[i - 1]); 75 ld = rd; 76 } 77 assert(integral >= 0); 78 } 79 y[$ - 1] = sqrt(integral * lambda / 3); 80 }; 81 82 ret.leastSquaresResult = optimize!(f)(settings, points.length + !lambda, y[].sliced, l[].sliced, u[].sliced); 83 84 return ret; 85 } 86 87 @safe pure 88 unittest 89 { 90 import mir.test; 91 92 LeastSquaresSettings!double settings; 93 94 auto x = [-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 95 96 auto y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 97 98 auto l = new double[x.length]; 99 l[] = -double.infinity; 100 101 auto u = new double[x.length]; 102 u[] = +double.infinity; 103 import mir.stdio; 104 105 double[2][] points = [ 106 [x[0] + 0.5, -0.68361541], 107 [x[1] + 0.5, 7.28568719], 108 [x[2] + 0.5, 10.490694 ], 109 [x[3] + 0.5, 0.36192032], 110 [x[4] + 0.5, 11.91572713], 111 [x[5] + 0.5, 16.44546433], 112 [x[6] + 0.5, 17.66699525], 113 [x[7] + 0.5, 4.52730869], 114 [x[8] + 0.5, 19.22825394], 115 [x[9] + 0.5, -2.3242592 ], 116 ]; 117 118 auto result = settings.fitSpline(points, x, l, u, 0); 119 120 foreach (i; 0 .. x.length) 121 result.spline(x[i]).shouldApprox == y[i]; 122 123 result = settings.fitSpline(points, x, l, u, 1); 124 125 // this case sensetive for numeric noise 126 y = [ 127 0.1971683531479794, 128 5.936895050720581, 129 7.451651002121712, 130 5.122509287945581, 131 11.908292461047825, 132 13.701350302891292, 133 16.97948422229589, 134 7.868130112291985, 135 16.20637990062554, 136 19.58302823176968, 137 ]; 138 139 foreach (i; 0 .. x.length) 140 result.spline(x[i]).shouldApprox == y[i]; 141 }