1 /++ 2 $(H2 Cubic Spline Interpolation) 3 4 The module provides common C2 splines, monotone (PCHIP) splines, Akima splines and others. 5 6 See_also: $(LREF SplineType), $(REF_ALTTEXT $(TT interp1), interp1, mir, interpolate) 7 8 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0) 9 Copyright: 2020 Ilia Ki, Kaleidic Associates Advisory Limited, Symmetry Investments 10 Authors: Ilia Ki 11 12 Macros: 13 SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, interpolate, $1)$(NBSP) 14 T2=$(TR $(TDNW $(LREF $1)) $(TD $+)) 15 +/ 16 module mir.interpolate.spline; 17 18 import core.lifetime: move; 19 import mir.functional; 20 import mir.internal.utility; 21 import mir.interpolate; 22 import mir.math.common; 23 import mir.ndslice.slice; 24 import mir.primitives; 25 import mir.rc.array; 26 import mir.utility: min, max; 27 import std.meta: AliasSeq, staticMap; 28 import std.traits: Unqual; 29 30 /// 31 public import mir.interpolate: atInterval; 32 33 static immutable msg_min = "spline interpolant: minimal allowed length for the grid equals 2."; 34 static immutable msg_eq = "spline interpolant: X and Y values length should be equal."; 35 static immutable splineConfigurationMsg = "spline configuration: .boundary method requires equal left and right boundaries"; 36 37 version(D_Exceptions) 38 { 39 static immutable exc_min = new Exception(msg_min); 40 static immutable exc_eq = new Exception(msg_eq); 41 static immutable splineConfigurationException = new Exception(splineConfigurationMsg); 42 } 43 44 private template ValueType(T, X) 45 { 46 static if (__traits(compiles, {enum N = T.dimensionCount;})) 47 alias ValueType = typeof(T.init.opCall(Repeat!(T.dimensionCount, X.init))); 48 else 49 alias ValueType = typeof(T.init.opCall(X.init)); 50 } 51 52 @fmamath: 53 54 /// 55 @safe pure @nogc version(mir_test) unittest 56 { 57 import mir.algorithm.iteration: all; 58 import mir.math.common: approxEqual; 59 import mir.ndslice.slice: sliced; 60 import mir.ndslice.allocation: rcslice; 61 import mir.ndslice.topology: vmap; 62 63 static immutable xdata = [-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 64 auto x = xdata.rcslice; 65 static immutable ydata = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 66 auto y = ydata.sliced; 67 68 auto interpolant = spline!double(x, y, SplineConfiguration!double()); // constructs Spline 69 auto xs = x + 0.5; // input X values for cubic spline 70 71 static immutable test_data0 = [ 72 -0.68361541, 7.28568719, 10.490694 , 0.36192032, 73 11.91572713, 16.44546433, 17.66699525, 4.52730869, 74 19.22825394, -2.3242592 ]; 75 /// not-a-knot (default) 76 assert(xs.vmap(interpolant).all!approxEqual(test_data0)); 77 78 static immutable test_data1 = [ 79 10.85298372, 5.26255911, 10.71443229, 0.1824536 , 80 11.94324989, 16.45633939, 17.59185094, 4.86340188, 81 17.8565408 , 2.81856494]; 82 /// natural cubic spline 83 interpolant = spline!double(x, y, SplineBoundaryType.secondDerivative); 84 assert(xs.vmap(interpolant).all!approxEqual(test_data1)); 85 86 static immutable test_data2 = [ 87 9.94191781, 5.4223652 , 10.69666392, 0.1971149 , 11.93868415, 88 16.46378847, 17.56521661, 4.97656997, 17.39645585, 4.54316446]; 89 /// set both boundary second derivatives to 3 90 interpolant = spline!double(x, y, SplineBoundaryType.secondDerivative, 3); 91 assert(xs.vmap(interpolant).all!approxEqual(test_data2)); 92 93 static immutable test_data3 = [ 94 16.45728263, 4.27981687, 10.82295092, 0.09610695, 95 11.95252862, 16.47583126, 17.49964521, 5.26561539, 96 16.21803478, 8.96104974]; 97 /// set both boundary derivatives to 3 98 interpolant = spline!double(x, y, SplineBoundaryType.firstDerivative, 3); 99 assert(xs.vmap(interpolant).all!approxEqual(test_data3)); 100 101 static immutable test_data4 = [ 102 16.45730084, 4.27966112, 10.82337171, 0.09403945, 103 11.96265209, 16.44067375, 17.6374694 , 4.67438921, 104 18.6234452 , -0.05582876]; 105 /// different boundary conditions 106 interpolant = spline!double(x, y, 107 SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3), 108 SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5)); 109 assert(xs.vmap(interpolant).all!approxEqual(test_data4)); 110 111 112 static immutable test_data5 = [ 113 12.37135558, 4.99638066, 10.74362441, 0.16008641, 114 11.94073593, 16.47908148, 17.49841853, 5.26600921, 115 16.21796051, 8.96102894]; 116 // ditto 117 interpolant = spline!double(x, y, 118 SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5), 119 SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3)); 120 assert(xs.vmap(interpolant).all!approxEqual(test_data5)); 121 122 static immutable test_data6 = [ 123 11.40871379, 2.64278898, 9.55774317, 4.84791141, 11.24842121, 124 16.16794267, 18.58060557, 5.2531411 , 17.45509005, 1.86992521]; 125 /// Akima spline 126 interpolant = spline!double(x, y, SplineType.akima); 127 assert(xs.vmap(interpolant).all!approxEqual(test_data6)); 128 129 static immutable test_data7 = [ 130 12.363463972190182, 3.066421497605127, 9.848385331143952, 4.429544487015751, 11.244211106655975, 131 16.255750063889600, 18.140374440374440, 5.291585909796099, 17.722477222271888, 3.181558328118484, 132 ]; 133 /// Modified Akima spline 134 interpolant = spline!double(x, y, SplineType.makima); 135 assert(xs.vmap(interpolant).all!approxEqual(test_data7)); 136 137 /// Double Quadratic spline 138 interpolant = spline!double(x, y, SplineType.doubleQuadratic); 139 import mir.interpolate.utility: ParabolaKernel; 140 auto kernel1 = ParabolaKernel!double(x[2], x[3], x[4], y[2], y[3], y[4]); 141 auto kernel2 = ParabolaKernel!double( x[3], x[4], x[5], y[3], y[4], y[5]); 142 // weighted sum of quadratic functions 143 auto c = 0.35; // from [0 .. 1] 144 auto xp = c * x[3] + (1 - c) * x[4]; 145 auto yp = c * kernel1(xp) + (1 - c) * kernel2(xp); 146 assert(interpolant(xp).approxEqual(yp)); 147 // check parabolic extrapolation of the boundary intervals 148 kernel1 = ParabolaKernel!double(x[0], x[1], x[2], y[0], y[1], y[2]); 149 kernel2 = ParabolaKernel!double(x[$ - 3], x[$ - 2], x[$ - 1], y[$ - 3], y[$ - 2], y[$ - 1]); 150 assert(interpolant(x[0] - 23.421).approxEqual(kernel1(x[0] - 23.421))); 151 assert(interpolant(x[$ - 1] + 23.421).approxEqual(kernel2(x[$ - 1] + 23.421))); 152 } 153 154 /// 155 @safe pure version(mir_test) unittest 156 { 157 import mir.rc.array: rcarray; 158 import mir.algorithm.iteration: all; 159 import mir.functional: aliasCall; 160 import mir.math.common: approxEqual; 161 import mir.ndslice.allocation: uninitSlice; 162 import mir.ndslice.slice: sliced; 163 import mir.ndslice.topology: vmap, map; 164 165 auto x = rcarray!(immutable double)(-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22).asSlice; 166 auto y = rcarray( 167 8.77842512, 168 7.96429686, 169 7.77074363, 170 1.10838032, 171 2.69925191, 172 1.84922654, 173 1.48167283, 174 2.8267636 , 175 0.40200172, 176 7.78980608).asSlice; 177 178 auto interpolant = x.spline!double(y); // default boundary condition is 'not-a-knot' 179 180 auto xs = x + 0.5; 181 182 auto ys = xs.vmap(interpolant); 183 184 auto r = 185 [5.56971848, 186 9.30342403, 187 4.44139761, 188 -0.74740285, 189 3.00994108, 190 1.50750417, 191 1.73144979, 192 2.64860361, 193 0.64413911, 194 10.81768928]; 195 196 assert(all!approxEqual(ys, r)); 197 198 // first derivative 199 auto d1 = xs.vmap(interpolant.aliasCall!"withDerivative").map!"a[1]"; 200 auto r1 = 201 [-4.51501279, 202 2.15715986, 203 -7.28363308, 204 -2.14050449, 205 0.03693092, 206 -0.49618999, 207 0.58109933, 208 -0.52926703, 209 0.7819035 , 210 6.70632693]; 211 assert(all!approxEqual(d1, r1)); 212 213 // second derivative 214 auto d2 = xs.vmap(interpolant.aliasCall!"withTwoDerivatives").map!"a[2]"; 215 auto r2 = 216 [7.07104751, 217 -2.62293241, 218 -0.01468508, 219 5.70609505, 220 -2.02358911, 221 0.72142061, 222 0.25275483, 223 -0.6133589 , 224 1.26894416, 225 2.68067146]; 226 assert(all!approxEqual(d2, r2)); 227 228 // third derivative (6 * a) 229 auto d3 = xs.vmap(interpolant.aliasCall!("opCall", 3)).map!"a[3]"; 230 auto r3 = 231 [-3.23132664, 232 -3.23132664, 233 14.91047457, 234 -3.46891432, 235 1.88520325, 236 -0.16559031, 237 -0.44056064, 238 0.47057577, 239 0.47057577, 240 0.47057577]; 241 assert(all!approxEqual(d3, r3)); 242 } 243 244 /// R -> R: Cubic interpolation 245 version(mir_test) 246 @safe unittest 247 { 248 import mir.test; 249 import mir.algorithm.iteration: all; 250 import mir.math.common: approxEqual; 251 import mir.ndslice; 252 253 static immutable x = [0, 1, 2, 3, 5.00274, 7.00274, 10.0055, 20.0137, 30.0192]; 254 static immutable y = [0.0011, 0.0011, 0.0030, 0.0064, 0.0144, 0.0207, 0.0261, 0.0329, 0.0356,]; 255 auto xs = [1, 2, 3, 4.00274, 5.00274, 6.00274, 7.00274, 8.00548, 9.00548, 10.0055, 11.0055, 12.0082, 13.0082, 14.0082, 15.0082, 16.011, 17.011, 18.011, 19.011, 20.0137, 21.0137, 22.0137, 23.0137, 24.0164, 25.0164, 26.0164, 27.0164, 28.0192, 29.0192, 30.0192]; 256 257 auto interpolation = spline!double(x.rcslice, y.sliced); 258 259 auto data = 260 [ 0.0011 , 0.003 , 0.0064 , 0.01042622, 0.0144 , 261 0.01786075, 0.0207 , 0.02293441, 0.02467983, 0.0261 , 262 0.02732764, 0.02840225, 0.0293308 , 0.03012914, 0.03081002, 263 0.03138766, 0.03187161, 0.03227637, 0.03261468, 0.0329 , 264 0.03314357, 0.03335896, 0.03355892, 0.03375674, 0.03396413, 265 0.03419436, 0.03446018, 0.03477529, 0.03515072, 0.0356 ]; 266 267 assert(all!approxEqual(xs.sliced.vmap(interpolation), data)); 268 } 269 270 /// R^2 -> R: Bicubic interpolation 271 version(mir_test) 272 unittest 273 { 274 import mir.math.common: approxEqual; 275 import mir.ndslice; 276 alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10); 277 278 ///// set test function //// 279 const double y_x0 = 2; 280 const double y_x1 = -7; 281 const double y_x0x1 = 3; 282 283 // this function should be approximated very well 284 alias f = (x0, x1) => y_x0 * x0 + y_x1 * x1 + y_x0x1 * x0 * x1 - 11; 285 286 ///// set interpolant //// 287 static immutable x0 = [-1.0, 2, 8, 15]; 288 static immutable x1 = [-4.0, 2, 5, 10, 13]; 289 auto grid = cartesian(x0, x1); 290 291 auto interpolant = spline!(double, 2)(x0.rcslice, x1.rcslice, grid.map!f); 292 293 ///// compute test data //// 294 auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 3.23); 295 // auto test_grid = cartesian(x0 + 0, x1 + 0); 296 auto real_data = test_grid.map!f; 297 auto interp_data = test_grid.vmap(interpolant); 298 299 ///// verify result //// 300 assert(all!appreq(interp_data, real_data)); 301 302 //// check derivatives //// 303 auto z0 = 1.23; 304 auto z1 = 3.21; 305 // writeln("-----------------"); 306 // writeln("-----------------"); 307 auto d = interpolant.withDerivative(z0, z1); 308 assert(appreq(interpolant(z0, z1), f(z0, z1))); 309 // writeln("d = ", d); 310 // writeln("interpolant.withTwoDerivatives(z0, z1) = ", interpolant.withTwoDerivatives(z0, z1)); 311 // writeln("-----------------"); 312 // writeln("-----------------"); 313 // writeln("interpolant(z0, z1) = ", interpolant(z0, z1)); 314 // writeln("y_x0 + y_x0x1 * z1 = ", y_x0 + y_x0x1 * z1); 315 // writeln("y_x1 + y_x0x1 * z0 = ", y_x1 + y_x0x1 * z0); 316 // writeln("-----------------"); 317 // writeln("-----------------"); 318 // assert(appreq(d[0][0], f(z0, z1))); 319 // assert(appreq(d[1][0], y_x0 + y_x0x1 * z1)); 320 // assert(appreq(d[0][1], y_x1 + y_x0x1 * z0)); 321 // assert(appreq(d[1][1], y_x0x1)); 322 } 323 324 /// R^3 -> R: Tricubic interpolation 325 version(mir_test) 326 unittest 327 { 328 import mir.math.common: approxEqual; 329 import mir.ndslice; 330 alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10); 331 332 ///// set test function //// 333 enum y_x0 = 2; 334 enum y_x1 = -7; 335 enum y_x2 = 5; 336 enum y_x0x1 = 10; 337 enum y_x0x1x2 = 3; 338 339 // this function should be approximated very well 340 alias f = (x0, x1, x2) => y_x0 * x0 + y_x1 * x1 + y_x2 * x2 341 + y_x0x1 * x0 * x1 + y_x0x1x2 * x0 * x1 * x2 - 11; 342 343 ///// set interpolant //// 344 static immutable x0 = [-1.0, 2, 8, 15]; 345 static immutable x1 = [-4.0, 2, 5, 10, 13]; 346 static immutable x2 = [3, 3.7, 5]; 347 auto grid = cartesian(x0, x1, x2); 348 349 auto interpolant = spline!(double, 3)(x0.rcslice, x1.rcslice, x2.rcslice, grid.map!f); 350 assert(interpolant.convexity == [SplineConvexity.none, SplineConvexity.none, SplineConvexity.convex]); 351 352 ///// compute test data //// 353 auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 3.23, x2.sliced - 3); 354 auto real_data = test_grid.map!f; 355 auto interp_data = test_grid.vmap(interpolant); 356 357 ///// verify result //// 358 assert(all!appreq(interp_data, real_data)); 359 360 //// check derivatives //// 361 auto z0 = 1.23; 362 auto z1 = 3.23; 363 auto z2 = -3; 364 auto d = interpolant.withDerivative(z0, z1, z2); 365 assert(appreq(interpolant(z0, z1, z2), f(z0, z1, z2))); 366 assert(appreq(d[0][0][0], f(z0, z1, z2))); 367 368 // writeln("-----------------"); 369 // writeln("-----------------"); 370 // auto d = interpolant.withDerivative(z0, z1); 371 assert(appreq(interpolant(z0, z1, z2), f(z0, z1, z2))); 372 // writeln("interpolant(z0, z1, z2) = ", interpolant(z0, z1, z2)); 373 // writeln("d = ", d); 374 // writeln("interpolant.withTwoDerivatives(z0, z1, z2) = ", interpolant.withTwoDerivatives(z0, z1, z2)); 375 // writeln("-----------------"); 376 // writeln("-----------------"); 377 // writeln("interpolant(z0, z1) = ", interpolant(z0, z1)); 378 // writeln("y_x0 + y_x0x1 * z1 = ", y_x0 + y_x0x1 * z1); 379 // writeln("y_x1 + y_x0x1 * z0 = ", y_x1 + y_x0x1 * z0); 380 // writeln("-----------------"); 381 // writeln("-----------------"); 382 383 // writeln("y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2 = ", y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2); 384 // assert(appreq(d[1][0][0], y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2)); 385 // writeln("y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2 = ", y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2); 386 // assert(appreq(d[0][1][0], y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2)); 387 // writeln("y_x0x1 + y_x0x1x2 * z2 = ", y_x0x1 + y_x0x1x2 * z2); 388 // assert(appreq(d[1][1][0], y_x0x1 + y_x0x1x2 * z2)); 389 // writeln("y_x0x1x2 = ", y_x0x1x2); 390 // assert(appreq(d[1][1][1], y_x0x1x2)); 391 } 392 393 394 /// Monotone PCHIP 395 version(mir_test) 396 @safe unittest 397 { 398 import mir.test; 399 import mir.math.common: approxEqual; 400 import mir.algorithm.iteration: all; 401 import mir.ndslice.allocation: rcslice; 402 import mir.ndslice.slice: sliced; 403 import mir.ndslice.topology: vmap; 404 405 static immutable x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 406 static immutable y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 407 auto interpolant = spline!double(x.rcslice, y.sliced, SplineType.monotone); 408 interpolant.argmin.should == 2; 409 interpolant.withDerivative(2)[1].should == 0; 410 interpolant.argmax.should == 12; 411 interpolant.withDerivative(12)[1].should == 0; 412 413 auto xs = x[0 .. $ - 1].sliced + 0.5; 414 415 auto ys = xs.vmap(interpolant); 416 417 assert(ys.all!approxEqual([ 418 5.333333333333334, 419 2.500000000000000, 420 10.000000000000000, 421 4.288971807628524, 422 11.202580845771145, 423 16.250000000000000, 424 17.962962962962962, 425 5.558593750000000, 426 17.604662698412699, 427 ])); 428 } 429 430 // Check direction equality 431 version(mir_test) 432 @safe unittest 433 { 434 import mir.math.common: approxEqual; 435 import mir.ndslice.slice: sliced; 436 import mir.ndslice.allocation: rcslice; 437 import mir.ndslice.topology: retro, vmap; 438 439 static immutable points = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 440 static immutable values = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 441 442 auto results = [ 443 5.333333333333334, 444 2.500000000000000, 445 10.000000000000000, 446 4.288971807628524, 447 11.202580845771145, 448 16.250000000000000, 449 17.962962962962962, 450 5.558593750000000, 451 17.604662698412699, 452 ]; 453 auto interpolant = spline!double(points.rcslice, values.sliced, SplineType.monotone); 454 455 auto pointsR = rcslice(-points.retro); 456 auto valuesR = values.retro.rcslice; 457 auto interpolantR = spline!double(pointsR, valuesR, SplineType.monotone); 458 459 version(X86_64) 460 assert(vmap(points[0 .. $ - 1].sliced + 0.5, interpolant) == vmap(pointsR.retro[0 .. $ - 1] - 0.5, interpolantR)); 461 } 462 463 /// argmin, +, - tests 464 version(mir_test) 465 unittest 466 { 467 import mir.test; 468 import mir.ndslice.slice: sliced; 469 import mir.ndslice.allocation: rcslice; 470 471 static immutable points = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 472 static immutable values = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 473 auto interpolant = spline!double(points.rcslice, values.sliced); 474 475 auto argmin = 5.898317667706634; 476 interpolant.argmin.shouldApprox == argmin; 477 interpolant(argmin).shouldApprox == -0.8684781710737299; 478 interpolant.opCall!2(argmin)[1].shouldApprox == 0; 479 480 auto argmax = 19.8816211020945; 481 interpolant.argmax.shouldApprox == argmax; 482 interpolant.withDerivative(argmax)[1].shouldApprox == 0; 483 interpolant(argmax).shouldApprox == 19.54377309088673; 484 485 auto zeroInterpolant = interpolant - interpolant; 486 zeroInterpolant.opCall!3(13.3).should == [0.0, 0.0, 0.0, 0.0]; 487 488 static immutable pointsR = [1.0, 3, 4,]; 489 static immutable valuesR = [13.0, 12, 10]; 490 auto interpolantR = spline!double(pointsR.rcslice, valuesR.sliced); 491 492 auto sumInterpolant = interpolant + interpolantR; 493 494 sumInterpolant(2.3).shouldApprox == interpolant(2.3) + interpolantR(2.3); 495 sumInterpolant(3.3).shouldApprox == interpolant(3.3) + interpolantR(3.3); 496 } 497 498 /++ 499 Cubic Spline types. 500 501 The first derivatives are guaranteed to be continuous for all cubic splines. 502 +/ 503 extern(C++, "mir", "interpolate") 504 enum SplineType 505 { 506 /++ 507 Spline with contiguous second derivative. 508 +/ 509 c2, 510 /++ 511 $(HTTPS en.wikipedia.org/wiki/Cubic_Hermite_spline#Cardinal_spline, Cardinal) and Catmull–Rom splines. 512 +/ 513 cardinal, 514 /++ 515 The interpolant preserves monotonicity in the interpolation data and does not overshoot if the data is not smooth. 516 It is also known as $(HTTPS docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.PchipInterpolator.html, PCHIP) 517 in numpy and Matlab. 518 +/ 519 monotone, 520 /++ 521 Weighted sum of two nearbor quadratic functions. 522 It is used in $(HTTPS s3-eu-west-1.amazonaws.com/og-public-downloads/smile-interpolation-extrapolation.pdf, financial analysis). 523 +/ 524 doubleQuadratic, 525 /++ 526 $(HTTPS en.wikipedia.org/wiki/Akima_spline, Akima spline). 527 +/ 528 akima, 529 /++ 530 $(HTTPS https://www.mathworks.com/help/matlab/ref/makima.html, Modified Akima spline). 531 +/ 532 makima, 533 } 534 535 /++ 536 Spline convexity type 537 +/ 538 enum SplineConvexity 539 { 540 /// Neither convex nor concave spline 541 none = 0, 542 /// Concave spline 543 concave = -1, 544 /// Convex spline 545 convex = 1, 546 } 547 548 /++ 549 Constructs multivariate cubic spline in symmetrical form with nodes on rectilinear grid. 550 Result has continues second derivatives throughout the curve / nd-surface. 551 +/ 552 template spline(T, size_t N = 1, X = T) 553 if (isFloatingPoint!T && is(T == Unqual!T) && N <= 6) 554 { 555 /++ 556 Params: 557 grid = immutable `x` values for interpolant 558 values = `f(x)` values for interpolant 559 typeOfBoundaries = $(LREF SplineBoundaryType) for both tails (optional). 560 valueOfBoundaryConditions = value of the boundary type (optional). 561 Constraints: 562 `grid` and `values` must have the same length >= 3 563 Returns: $(LREF Spline) 564 +/ 565 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 566 Repeat!(N, Slice!(RCI!(immutable X))) grid, 567 Slice!(yIterator, N, ykind) values, 568 SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, 569 in T valueOfBoundaryConditions = 0, 570 ) 571 { 572 import core.lifetime: forward; 573 return spline(forward!grid, forward!values, SplineType.c2, 0, typeOfBoundaries, valueOfBoundaryConditions); 574 } 575 576 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 577 Repeat!(N, Slice!(RCI!(immutable X))) grid, 578 Slice!(yIterator, N, ykind) values, 579 SplineType kind, 580 in T param = 0, 581 SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, 582 in T valueOfBoundaryConditions = 0, 583 ) 584 { 585 import core.lifetime: forward; 586 return spline(forward!grid, forward!values, SplineBoundaryCondition!T(typeOfBoundaries, valueOfBoundaryConditions), kind, param); 587 } 588 589 /++ 590 Params: 591 grid = immutable `x` values for interpolant 592 values = `f(x)` values for interpolant 593 boundaries = $(LREF SplineBoundaryCondition) for both tails. 594 kind = $(LREF SplineType) type of cubic spline. 595 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 596 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 597 Constraints: 598 `grid` and `values` must have the same length >= 3 599 Returns: $(LREF Spline) 600 +/ 601 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 602 Repeat!(N, Slice!(RCI!(immutable X))) grid, 603 Slice!(yIterator, N, ykind) values, 604 SplineBoundaryCondition!T boundaries, 605 SplineType kind = SplineType.c2, 606 in T param = 0, 607 ) 608 { 609 import core.lifetime: forward; 610 return spline(forward!grid, forward!values, boundaries, boundaries, kind, param); 611 } 612 613 /++ 614 Params: 615 grid = immutable `x` values for interpolant 616 values = `f(x)` values for interpolant 617 lBoundary = $(LREF SplineBoundaryCondition) for left tail. 618 rBoundary = $(LREF SplineBoundaryCondition) for right tail. 619 kind = $(LREF SplineType) type of cubic spline. 620 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 621 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 622 Constraints: 623 `grid` and `values` must have the same length >= 3 624 Returns: $(LREF Spline) 625 +/ 626 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 627 Repeat!(N, Slice!(RCI!(immutable X))) grid, 628 Slice!(yIterator, N, ykind) values, 629 SplineBoundaryCondition!T lBoundary, 630 SplineBoundaryCondition!T rBoundary, 631 SplineType kind = SplineType.c2, 632 in T param = 0, 633 ) 634 { 635 auto ret = typeof(return)(forward!grid); 636 ret._values = values; 637 ret._computeDerivatives(kind, param, lBoundary, rBoundary); 638 return ret; 639 } 640 641 /++ 642 Params: 643 grid = immutable `x` values for interpolant 644 values = `f(x)` values for interpolant 645 configuration = $(LREF SplineConfiguration) 646 Constraints: 647 `grid` and `values` must have the same length >= 3 648 Returns: $(LREF Spline) 649 +/ 650 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 651 Repeat!(N, Slice!(RCI!(immutable X))) grid, 652 Slice!(yIterator, N, ykind) values, 653 SplineConfiguration!T configuration, 654 ) 655 { 656 auto ret = typeof(return)(forward!grid); 657 ret._values = values; 658 with(configuration) 659 ret._computeDerivatives(kind, param, leftBoundary, rightBoundary); 660 return ret; 661 } 662 } 663 664 /++ 665 Cubic Spline Boundary Condition Type. 666 667 See_also: $(LREF SplineBoundaryCondition) $(LREF SplineType) 668 +/ 669 extern(C++, "mir", "interpolate") 670 enum SplineBoundaryType 671 { 672 /++ 673 Not-a-knot (or cubic) boundary condition. 674 It is an aggresive boundary condition that is used only for C2 splines and is default for all API calls. 675 For other then C2 splines, `notAKnot` is changed internally to 676 a default boundary type for used $(LREF SplineType). 677 +/ 678 notAKnot, 679 /++ 680 Set the first derivative. 681 +/ 682 firstDerivative, 683 /++ 684 Set the second derivative. 685 +/ 686 secondDerivative, 687 /++ 688 Default for Cardinal and Double-Quadratic splines. 689 +/ 690 parabolic, 691 /++ 692 Default for monotone (aka PHCIP ) splines. 693 +/ 694 monotone, 695 /++ 696 Default for Akima splines. 697 +/ 698 akima, 699 /++ 700 Default for Modified Akima splines. 701 +/ 702 makima, 703 /++ 704 Not implemented. 705 +/ 706 periodic = -1, 707 } 708 709 /++ 710 Cubic Spline Boundary Condition 711 712 See_also: $(LREF SplineBoundaryType) 713 +/ 714 extern(C++, "mir", "interpolate") 715 struct SplineBoundaryCondition(T) 716 if (__traits(isFloating, T)) 717 { 718 import mir.serde: serdeOptional, serdeIgnoreDefault; 719 720 /// type (default is $(LREF SplineBoundaryType.notAKnot)) 721 SplineBoundaryType type; 722 723 @serdeOptional @serdeIgnoreDefault: 724 725 /// value (default is 0) 726 T value = 0; 727 } 728 729 /// Spline configuration 730 struct SplineConfiguration(T) 731 if (__traits(isFloating, T)) 732 { 733 import mir.serde: serdeOptional, serdeIgnoreDefault, serdeIgnoreOutIfAggregate, serdeIgnore; 734 735 /// 736 @serdeOptional @serdeIgnoreDefault 737 SplineType kind; 738 /// 739 @serdeOptional @serdeIgnoreOutIfAggregate!"a.symmetric" 740 SplineBoundaryCondition!T leftBoundary; 741 /// 742 @serdeOptional @serdeIgnoreOutIfAggregate!"a.symmetric" 743 SplineBoundaryCondition!T rightBoundary; 744 745 /++ 746 Returns: 747 true of `leftBoundary` equals `rightBoundary`. 748 +/ 749 @serdeIgnore 750 bool symmetric() const @property 751 { 752 return leftBoundary == rightBoundary; 753 } 754 755 /// 756 @serdeOptional 757 void boundary(SplineBoundaryCondition!T boundary) @property 758 { 759 leftBoundary = rightBoundary = boundary; 760 } 761 762 /// 763 @serdeIgnoreOutIfAggregate!"!a.symmetric" 764 SplineBoundaryCondition!T boundary() const @property 765 { 766 assert(!symmetric, splineConfigurationMsg); 767 return leftBoundary; 768 } 769 770 /++ 771 Tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 772 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 773 +/ 774 @serdeOptional @serdeIgnoreDefault 775 T param = 0; 776 } 777 778 /// Spline configuration with two boundaries 779 struct SplineSymmetricConfiguration(T) 780 if (__traits(isFloating, T)) 781 { 782 import mir.serde: serdeOptional, serdeIgnoreDefault; 783 784 @serdeOptional @serdeIgnoreDefault: 785 786 /// 787 SplineType type; 788 /// 789 SplineBoundaryCondition!T boundary; 790 /++ 791 Tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 792 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 793 +/ 794 T param = 0; 795 } 796 797 /++ 798 Multivariate cubic spline with nodes on rectilinear grid. 799 +/ 800 struct Spline(F, size_t N = 1, X = F) 801 if (N && N <= 6) 802 { 803 import mir.rc.array; 804 805 /// Aligned buffer allocated with `mir.internal.memory`. $(RED For internal use.) 806 Slice!(RCI!(F[2 ^^ N]), N) _data; 807 /// Grid iterators. $(RED For internal use.) 808 Repeat!(N, RCI!(immutable X)) _grid; 809 /// 810 SplineConvexity[N] convexity; 811 812 enum uint dimensionCount = N; 813 814 @fmamath extern(D): 815 816 bool opEquals()(auto ref scope const typeof(this) rhs) scope const @trusted pure nothrow @nogc 817 { 818 if (rhs._data != this._data) 819 return false; 820 static foreach (d; 0 .. N) 821 if (gridScopeView!d != rhs.gridScopeView!d) 822 return false; 823 return true; 824 } 825 826 /++ 827 +/ 828 this(Repeat!(N, Slice!(RCI!(immutable X))) grid) @safe @nogc 829 { 830 size_t length = 1; 831 size_t[N] shape; 832 foreach(i, ref x; grid) 833 { 834 if (x.length < 2) 835 { 836 version(D_Exceptions) { import mir.exception : toMutable; throw exc_min.toMutable; } 837 else assert(0, msg_min); 838 } 839 length *= shape[i] = x.length; 840 this._grid[i] = x._iterator.move; 841 } 842 import mir.ndslice.allocation: rcslice; 843 this._data = shape.rcslice!(F[2 ^^ N]); 844 } 845 846 package static auto pickDataSubslice(D)(auto scope ref D data, size_t index) @trusted 847 { 848 auto strides = data.strides; 849 foreach (i; Iota!(strides.length)) 850 strides[i] *= DeepElementType!D.length; 851 return Slice!(F*, strides.length, Universal)(data.shape, strides, data._iterator.ptr + index); 852 } 853 854 static if (N == 1) 855 /++ 856 Note: defined only for 1D splines 857 +/ 858 Spline opBinary(string op)(const Spline rhs) @trusted const 859 if (op == "+" || op == "-") 860 { 861 import mir.ndslice.allocation: rcslice; 862 import core.lifetime: move; 863 import mir.algorithm.setops: unionLength, multiwayUnion; 864 865 auto lgrid = this.gridScopeView; 866 auto rgrid = rhs.gridScopeView; 867 scope const(F)[][2] grids = [lgrid, rgrid]; 868 auto length = grids[].unionLength; 869 grids = [lgrid, rgrid]; 870 871 size_t j; 872 auto grid = RCArray!X(length); 873 auto data = length.rcslice!(F[2]); 874 auto un = grids[].multiwayUnion; 875 876 while (!un.empty) 877 { 878 auto x = un.front; 879 un.popFront; 880 auto ly = this.opCall!1(x); 881 auto ry = rhs.opCall!1(x); 882 data[j] = mixin(`[ly[0] ` ~ op ~ ` ry[0], ly[1] ` ~ op ~ ` ry[1]]`); 883 grid[j++] = x; 884 } 885 886 size_t convexCount; 887 size_t concaveCount; 888 foreach_reverse (i; 0 .. length - 1) 889 { 890 auto xdiff = grid[i + 1] - grid[i]; 891 auto ydiff = data[i + 1][0] - data[i][0]; 892 893 auto convex1 = xdiff * (2 * data[i][1] + data[i + 1][1]) <= 3 * ydiff; 894 auto concave1 = xdiff * (2 * data[i][1] + data[i + 1][1]) >= 3 * ydiff; 895 auto convex2 = xdiff * (data[i][1] + 2 * data[i + 1][1]) >= 3 * ydiff; 896 auto concave2 = xdiff * (data[i][1] + 2 * data[i + 1][1]) <= 3 * ydiff; 897 convexCount += convex1 & convex2; 898 convexCount += concave1 & concave2; 899 } 900 901 Spline ret; 902 ret._data = data.move; 903 ret._grid[0] = RCI!(immutable X)(cast(RCArray!(immutable X))grid) ; 904 905 ret.convexity = 906 // convex kind has priority for the linear spline 907 convexCount == length - 1 ? SplineConvexity.convex : 908 concaveCount == length - 1 ? SplineConvexity.concave : 909 SplineConvexity.none; 910 911 return ret; 912 } 913 914 static if (N == 1) 915 template argminImpl(string pred) 916 if (pred == "a < b" || pred == "a > b") 917 { 918 F argminImpl()() @trusted const @property 919 { 920 import mir.functional: naryFun; 921 922 static if (pred == "a < b") 923 auto min = F.max; 924 else 925 auto min = -F.max; 926 927 F argmin; 928 929 auto grid = gridScopeView; 930 foreach (i, ref y; _data.lightScope.field) 931 { 932 if (naryFun!pred(y[0], min)) 933 { 934 min = y[0]; 935 argmin = grid[i]; 936 } 937 } 938 939 foreach (i; 0 .. grid.length - 1) 940 { 941 auto x = grid[i + 0]; 942 943 auto y = SplineKernel!F( 944 grid[i + 0], 945 grid[i + 1], 946 x, // any point between 947 ).opCall!3( 948 _data[i + 0][0], 949 _data[i + 1][0], 950 _data[i + 0][1], 951 _data[i + 1][1], 952 ); 953 954 // 3 ax^2 + 2 bx + c = y[1] 955 // 6 ax + 2 b= y[2] 956 // 6 a = y[3] 957 958 auto a3 = y[3] * 0.5; 959 auto b2 = y[2] - y[3] * x; 960 auto c = y[1] - (b2 + a3 * x) * x; 961 auto d = b2 * b2 - 4 * a3 * c; 962 if (d < 0) 963 continue; 964 965 import mir.math.common: sqrt; 966 F[2] x12 = [ 967 (-b2 - sqrt(d)) / (2 * a3), 968 (-b2 + sqrt(d)) / (2 * a3), 969 ]; 970 971 foreach (xi; x12) 972 if (grid[i + 0] < xi && xi < grid[i + 1]) 973 { 974 auto yi = SplineKernel!F( 975 grid[i + 0], 976 grid[i + 1], 977 xi, 978 )( 979 _data[i + 0][0], 980 _data[i + 1][0], 981 _data[i + 0][1], 982 _data[i + 1][1], 983 ); 984 985 if (naryFun!pred(yi, min)) 986 { 987 min = yi; 988 argmin = xi; 989 } 990 } 991 } 992 993 return argmin; 994 } 995 } 996 997 static if (N == 1) 998 /++ 999 Returns: spline argmin on the interpolation interval 1000 Note: defined only for 1D splines 1001 +/ 1002 alias argmin = argminImpl!"a < b"; 1003 1004 static if (N == 1) 1005 /++ 1006 Returns: spline argmax on the interpolation interval 1007 Note: defined only for 1D splines 1008 +/ 1009 alias argmax = argminImpl!"a > b"; 1010 1011 /++ 1012 Assigns function values to the internal memory. 1013 $(RED For internal use.) 1014 +/ 1015 void _values(SliceKind kind, Iterator)(Slice!(Iterator, N, kind) values) scope @property @trusted 1016 { 1017 assert(values.shape == _data.shape, "'values' should have the same shape as the .gridShape"); 1018 pickDataSubslice(_data.lightScope, 0)[] = values; 1019 } 1020 1021 /++ 1022 Computes derivatives and stores them in `_data`. 1023 `_data` is assumed to be preinitialized with function values filled in `F[2 ^^ N][0]`. 1024 Params: 1025 lBoundary = left boundary condition 1026 rBoundary = right boundary condition 1027 temp = temporal buffer length points count (optional) 1028 1029 $(RED For internal use.) 1030 +/ 1031 void _computeDerivatives(SplineType kind, F param, SplineBoundaryCondition!F lBoundary, SplineBoundaryCondition!F rBoundary) scope @trusted nothrow @nogc 1032 { 1033 import mir.algorithm.iteration: maxLength; 1034 auto ml = this._data.maxLength; 1035 auto temp = RCArray!F(ml); 1036 auto tempSlice = temp[].sliced; 1037 _computeDerivativesTemp(kind, param, lBoundary, rBoundary, tempSlice); 1038 } 1039 1040 /// ditto 1041 pragma(inline, false) 1042 void _computeDerivativesTemp(SplineType kind, F param, SplineBoundaryCondition!F lBoundary, SplineBoundaryCondition!F rBoundary, Slice!(F*) temp) scope @system nothrow @nogc 1043 { 1044 import mir.algorithm.iteration: maxLength, each; 1045 import mir.ndslice.topology: map, byDim, evertPack; 1046 1047 assert(temp.length >= _data.maxLength); 1048 1049 static if (N == 1) 1050 { 1051 convexity[0] = splineSlopes!(F, F)(_grid[0]._iterator.sliced(_data._lengths[0]), pickDataSubslice(_data.lightScope, 0), pickDataSubslice(_data.lightScope, 1), temp[0 .. _data._lengths[0]], kind, param, lBoundary, rBoundary); 1052 } 1053 else 1054 foreach_reverse(i; Iota!N) 1055 { 1056 // if (i == _grid.length - 1) 1057 _data 1058 .lightScope 1059 .byDim!i 1060 .evertPack 1061 .each!((d){ 1062 enum L = 2 ^^ (N - 1 - i); 1063 foreach(l; Iota!L) 1064 { 1065 auto y = pickDataSubslice(d, l); 1066 auto s = pickDataSubslice(d, L + l); 1067 // debug printf("ptr = %ld, stride = %ld, stride = %ld, d = %ld i = %ld l = %ld\n", d.iterator, d._stride!0, y._stride!0, d.length, i, l); 1068 auto c = splineSlopes!(F, F)(_grid[i]._iterator.sliced(_data._lengths[i]), y, s, temp[0 .. _data._lengths[i]], kind, param, lBoundary, rBoundary); 1069 static if (l) 1070 { 1071 if (convexity[i] != c) 1072 convexity[i] = SplineConvexity.none; 1073 } 1074 else 1075 { 1076 convexity[i] = c; 1077 } 1078 // debug{ 1079 // (cast(void delegate() @nogc)(){ 1080 // writeln("y = ", y); 1081 // writeln("s = ", s); 1082 // })(); 1083 // } 1084 } 1085 }); 1086 } 1087 } 1088 1089 @trusted: 1090 1091 /// 1092 Spline lightConst() const @property { return *cast(Spline*)&this; } 1093 /// 1094 Spline lightImmutable() immutable @property { return *cast(Spline*)&this; } 1095 1096 /// 1097 Slice!(RCI!(immutable X)) grid(size_t dimension = 0)() return scope const @property 1098 if (dimension < N) 1099 { 1100 return _grid[dimension].lightConst.sliced(_data._lengths[dimension]); 1101 } 1102 1103 /// 1104 immutable(X)[] gridScopeView(size_t dimension = 0)() return scope const @property @trusted 1105 if (dimension < N) 1106 { 1107 return _grid[dimension]._iterator[0 .. _data._lengths[dimension]]; 1108 } 1109 1110 /++ 1111 Returns: intervals count. 1112 +/ 1113 size_t intervalCount(size_t dimension = 0)() scope const @property 1114 { 1115 assert(_data._lengths[dimension] > 1); 1116 return _data._lengths[dimension] - 1; 1117 } 1118 1119 /// 1120 size_t[N] gridShape() scope const @property 1121 { 1122 return _data.shape; 1123 } 1124 1125 /// 1126 alias withDerivative = opCall!1; 1127 /// 1128 alias withTwoDerivatives = opCall!2; 1129 1130 /// 1131 enum uint derivativeOrder = 3; 1132 1133 /// 1134 template opCall(uint derivative : 2) 1135 { 1136 auto opCall(X...)(in X xs) scope const 1137 if (X.length == N) 1138 // @FUTURE@ 1139 // X.length == N || derivative == 0 && X.length && X.length <= N 1140 { 1141 auto d4 = this.opCall!3(xs); 1142 SplineReturnType!(F, N, 3) d3; 1143 void fun(size_t d, A, B)(ref A a, ref B b) 1144 { 1145 static if (d) 1146 foreach(i; Iota!3) 1147 fun!(d - 1)(a[i], b[i]); 1148 else 1149 b = a; 1150 } 1151 fun!N(d4, d3); 1152 return d3; 1153 } 1154 } 1155 1156 /// 1157 template opCall(uint derivative = 0) 1158 if (derivative == 0 || derivative == 1 || derivative == 3) 1159 { 1160 static if (N > 1 && derivative) pragma(msg, "Warning: multivariate cubic spline with derivatives was not tested!!!"); 1161 1162 /++ 1163 `(x)` operator. 1164 Complexity: 1165 `O(log(points.length))` 1166 +/ 1167 auto opCall(X...)(in X xs) scope const @trusted 1168 if (X.length == N) 1169 // @FUTURE@ 1170 // X.length == N || derivative == 0 && X.length && X.length <= N 1171 { 1172 import mir.ndslice.topology: iota; 1173 alias Kernel = AliasCall!(SplineKernel!F, "opCall", derivative); 1174 enum rp2d = derivative == 3 ? 2 : derivative; 1175 1176 size_t[N] indices; 1177 Kernel[N] kernels; 1178 1179 foreach(i; Iota!N) 1180 { 1181 static if (isInterval!(typeof(xs[i]))) 1182 { 1183 indices[i] = xs[i][1]; 1184 auto x = xs[i][0]; 1185 } 1186 else 1187 { 1188 alias x = xs[i]; 1189 indices[i] = this.findInterval!i(x); 1190 } 1191 kernels[i] = SplineKernel!F(_grid[i][indices[i]], _grid[i][indices[i] + 1], x); 1192 } 1193 1194 align(64) F[2 ^^ N * 2 ^^ N][2] local; 1195 1196 void load(sizediff_t i)(const(F[2 ^^ N])* from, F[2 ^^ N]* to) 1197 { 1198 version(LDC) pragma(inline, true); 1199 static if (i == -1) 1200 { 1201 // copyvec(*from, *to); 1202 // not aligned: 1203 *to = *from; 1204 } 1205 else 1206 { 1207 from += strides[i] * indices[i]; 1208 load!(i - 1)(from, to); 1209 from += strides[i]; 1210 enum s = 2 ^^ (N - 1 - i); 1211 to += s; 1212 load!(i - 1)(from, to); 1213 } 1214 } 1215 1216 immutable strides = _data._lengths.iota.strides; 1217 load!(N - 1)(_data.ptr, cast(F[2 ^^ N]*) local[0].ptr); 1218 1219 // debug{ 1220 1221 // printf("0local[0] = "); 1222 // foreach(ref e; local[0][]) 1223 // printf("%f ", e); 1224 // printf("\n"); 1225 // } 1226 1227 foreach(i; Iota!N) 1228 { 1229 enum P = 2 ^^ (N - 1 - i) * 2 ^^ (i * rp2d); 1230 enum L = (2 ^^ N) ^^ 2 / (2 ^^ (i * (2 - rp2d))) / 4; 1231 shuffle2!P(local[0][0 * L .. 1 * L], local[0][1 * L .. 2 * L], local[1][0 * L .. 1 * L], local[1][1 * L .. 2 * L]); 1232 shuffle2!P(local[0][2 * L .. 3 * L], local[0][3 * L .. 4 * L], local[1][2 * L .. 3 * L], local[1][3 * L .. 4 * L]); 1233 // debug 1234 // { 1235 // printf("0local[1] = "); 1236 // foreach(ref e; local[1][0 .. L* 4]) 1237 // printf("%f ", e); 1238 // printf("\n"); 1239 // } 1240 local[0][] = F.init; 1241 vectorize( 1242 kernels[i], 1243 local[1][0 * L .. 1 * L], local[1][2 * L .. 3 * L], 1244 local[1][1 * L .. 2 * L], local[1][3 * L .. 4 * L], 1245 *cast(F[L][2 ^^ rp2d]*) local[0].ptr, 1246 ); 1247 1248 // debug{ 1249 1250 // printf("1local[0] = "); 1251 // foreach(ref e; local[0][0 .. L* 2 ^^ rp2d]) 1252 // printf("%f ", e); 1253 // printf("\n"); 1254 // } 1255 // printf("local[0][0]", local[0][0]); 1256 static if (i + 1 == N) 1257 { 1258 return *cast(SplineReturnType!(F, N, 2 ^^ rp2d)*) local[0].ptr; 1259 } 1260 else 1261 { 1262 static if (rp2d == 1) 1263 { 1264 shuffle3!1(local[0][0 .. L], local[0][L .. 2 * L], local[1][0 .. L], local[1][L .. 2 * L]); 1265 copyvec(local[1][0 .. 1 * L], local[0][0 .. 1 * L]); 1266 copyvec(local[1][L .. 2 * L], local[0][L .. 2 * L]); 1267 } 1268 else 1269 static if (rp2d == 2) 1270 { 1271 shuffle3!1(local[0][0 * L .. 1 * L], local[0][1 * L .. 2 * L], local[1][0 * L .. 1 * L], local[1][1 * L .. 2 * L]); 1272 shuffle3!1(local[0][2 * L .. 3 * L], local[0][3 * L .. 4 * L], local[1][2 * L .. 3 * L], local[1][3 * L .. 4 * L]); 1273 shuffle3!2(local[1][0 * L .. 1 * L], local[1][2 * L .. 3 * L], local[0][0 * L .. 1 * L], local[0][2 * L .. 3 * L]); 1274 shuffle3!2(local[1][1 * L .. 2 * L], local[1][3 * L .. 4 * L], local[0][1 * L .. 2 * L], local[0][3 * L .. 4 * L]); 1275 } 1276 1277 // debug{ 1278 1279 // printf("2local[0] = "); 1280 // foreach(ref e; local[0][0 .. L * 2 ^^ rp2d]) 1281 // printf("%f ", e); 1282 // printf("\n"); 1283 // } 1284 1285 } 1286 } 1287 } 1288 } 1289 } 1290 1291 /++ 1292 Piecewise cubic hermite interpolating polynomial. 1293 Params: 1294 points = `x` values for interpolant 1295 values = `f(x)` values for interpolant 1296 slopes = uninitialized ndslice to write slopes into 1297 temp = uninitialized temporary ndslice 1298 kind = $(LREF SplineType) type of cubic spline. 1299 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 1300 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 1301 lBoundary = left boundary condition 1302 rBoundary = right boundary condition 1303 Constraints: 1304 `points`, `values`, and `slopes`, must have the same length > 3; 1305 `temp` must have length greater or equal to points less minus one. 1306 Returs: 1307 Spline convexity type 1308 +/ 1309 SplineConvexity splineSlopes(F, T, P, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind skind)( 1310 Slice!(P*, 1, gkind) points, 1311 Slice!(IV, 1, vkind) values, 1312 Slice!(IS, 1, skind) slopes, 1313 Slice!(T*) temp, 1314 SplineType kind, 1315 F param, 1316 SplineBoundaryCondition!F lBoundary, 1317 SplineBoundaryCondition!F rBoundary, 1318 ) @trusted 1319 { 1320 import mir.ndslice.topology: diff, zip, slide; 1321 1322 assert (points.length >= 2); 1323 assert (points.length == values.length); 1324 assert (points.length == slopes.length); 1325 assert (temp.length == points.length); 1326 1327 auto n = points.length; 1328 1329 typeof(slopes[0]) first, last; 1330 1331 auto xd = points.diff; 1332 auto yd = values.diff; 1333 auto dd = yd / xd; 1334 auto dd2 = points.zip(values).slide!(3, "(c[1] - a[1]) / (c[0] - a[0])"); 1335 1336 with(SplineType) final switch(kind) 1337 { 1338 case c2: 1339 break; 1340 case cardinal: 1341 if (lBoundary.type == SplineBoundaryType.notAKnot) 1342 lBoundary.type = SplineBoundaryType.parabolic; 1343 if (rBoundary.type == SplineBoundaryType.notAKnot) 1344 rBoundary.type = SplineBoundaryType.parabolic; 1345 break; 1346 case monotone: 1347 if (lBoundary.type == SplineBoundaryType.notAKnot) 1348 lBoundary.type = SplineBoundaryType.monotone; 1349 if (rBoundary.type == SplineBoundaryType.notAKnot) 1350 rBoundary.type = SplineBoundaryType.monotone; 1351 break; 1352 case doubleQuadratic: 1353 if (lBoundary.type == SplineBoundaryType.notAKnot) 1354 lBoundary.type = SplineBoundaryType.parabolic; 1355 if (rBoundary.type == SplineBoundaryType.notAKnot) 1356 rBoundary.type = SplineBoundaryType.parabolic; 1357 break; 1358 case akima: 1359 if (lBoundary.type == SplineBoundaryType.notAKnot) 1360 lBoundary.type = SplineBoundaryType.akima; 1361 if (rBoundary.type == SplineBoundaryType.notAKnot) 1362 rBoundary.type = SplineBoundaryType.akima; 1363 break; 1364 case makima: 1365 if (lBoundary.type == SplineBoundaryType.notAKnot) 1366 lBoundary.type = SplineBoundaryType.makima; 1367 if (rBoundary.type == SplineBoundaryType.notAKnot) 1368 rBoundary.type = SplineBoundaryType.makima; 1369 break; 1370 } 1371 1372 if (n <= 3) 1373 { 1374 if (lBoundary.type == SplineBoundaryType.notAKnot) 1375 lBoundary.type = SplineBoundaryType.parabolic; 1376 if (rBoundary.type == SplineBoundaryType.notAKnot) 1377 rBoundary.type = SplineBoundaryType.parabolic; 1378 1379 if (n == 2) 1380 { 1381 if (lBoundary.type == SplineBoundaryType.monotone 1382 || lBoundary.type == SplineBoundaryType.makima 1383 || lBoundary.type == SplineBoundaryType.akima) 1384 lBoundary.type = SplineBoundaryType.parabolic; 1385 if (rBoundary.type == SplineBoundaryType.monotone 1386 || rBoundary.type == SplineBoundaryType.makima 1387 || rBoundary.type == SplineBoundaryType.akima) 1388 rBoundary.type = SplineBoundaryType.parabolic; 1389 } 1390 /// special case 1391 if (rBoundary.type == SplineBoundaryType.parabolic 1392 && lBoundary.type == SplineBoundaryType.parabolic) 1393 { 1394 import mir.interpolate.utility; 1395 if (n == 3) 1396 { 1397 auto derivatives = parabolaDerivatives(points[0], points[1], points[2], values[0], values[1], values[2]); 1398 slopes[0] = derivatives[0]; 1399 slopes[1] = derivatives[1]; 1400 slopes[2] = derivatives[2]; 1401 } 1402 else 1403 { 1404 assert(slopes.length == 2); 1405 slopes.back = slopes.front = yd.front / xd.front; 1406 } 1407 return slopes.back >= slopes.front ? SplineConvexity.convex : SplineConvexity.concave; 1408 } 1409 } 1410 1411 with(SplineBoundaryType) final switch(lBoundary.type) 1412 { 1413 case periodic: 1414 1415 assert(0); 1416 1417 case notAKnot: 1418 1419 auto dx0 = xd[0]; 1420 auto dx1 = xd[1]; 1421 auto dy0 = yd[0]; 1422 auto dy1 = yd[1]; 1423 auto dd0 = dy0 / dx0; 1424 auto dd1 = dy1 / dx1; 1425 1426 slopes.front = dx1; 1427 first = dx0 + dx1; 1428 temp.front = ((dx0 + 2 * first) * dx1 * dd0 + dx0 ^^ 2 * dd1) / first; 1429 break; 1430 1431 case firstDerivative: 1432 1433 slopes.front = 1; 1434 first = 0; 1435 temp.front = lBoundary.value; 1436 break; 1437 1438 case secondDerivative: 1439 1440 slopes.front = 2; 1441 first = 1; 1442 temp.front = 3 * dd.front - 0.5 * lBoundary.value * xd.front; 1443 break; 1444 1445 case parabolic: 1446 1447 slopes.front = 1; 1448 first = 1; 1449 temp.front = 2 * dd.front; 1450 break; 1451 1452 case monotone: 1453 1454 slopes.front = 1; 1455 first = 0; 1456 temp.front = pchipTail(xd[0], xd[1], dd[0], dd[1]); 1457 break; 1458 1459 case akima: 1460 1461 slopes.front = 1; 1462 first = 0; 1463 temp.front = akimaTail(dd[0], dd[1]); 1464 break; 1465 1466 case makima: 1467 1468 slopes.front = 1; 1469 first = 0; 1470 temp.front = makimaTail(dd[0], dd[1]); 1471 break; 1472 1473 } 1474 1475 with(SplineBoundaryType) final switch(rBoundary.type) 1476 { 1477 case periodic: 1478 assert(0); 1479 1480 case notAKnot: 1481 1482 auto dx0 = xd[$ - 1]; 1483 auto dx1 = xd[$ - 2]; 1484 auto dy0 = yd[$ - 1]; 1485 auto dy1 = yd[$ - 2]; 1486 auto dd0 = dy0 / dx0; 1487 auto dd1 = dy1 / dx1; 1488 slopes.back = dx1; 1489 last = dx0 + dx1; 1490 temp.back = ((dx0 + 2 * last) * dx1 * dd0 + dx0 ^^ 2 * dd1) / last; 1491 break; 1492 1493 case firstDerivative: 1494 1495 slopes.back = 1; 1496 last = 0; 1497 temp.back = rBoundary.value; 1498 break; 1499 1500 case secondDerivative: 1501 1502 slopes.back = 2; 1503 last = 1; 1504 temp.back = 3 * dd.back + 0.5 * rBoundary.value * xd.back; 1505 break; 1506 1507 case parabolic: 1508 1509 slopes.back = 1; 1510 last = 1; 1511 temp.back = 2 * dd.back; 1512 break; 1513 1514 case monotone: 1515 1516 slopes.back = 1; 1517 last = 0; 1518 temp.back = pchipTail(xd[$ - 1], xd[$ - 2], dd[$ - 1], dd[$ - 2]); 1519 break; 1520 1521 case akima: 1522 1523 slopes.back = 1; 1524 last = 0; 1525 temp.back = akimaTail(dd[$ - 1], dd[$ - 2]); 1526 break; 1527 1528 case makima: 1529 1530 slopes.back = 1; 1531 last = 0; 1532 temp.back = makimaTail(dd[$ - 1], dd[$ - 2]); 1533 break; 1534 1535 } 1536 1537 if (n > 2) with(SplineType) final switch(kind) 1538 { 1539 case c2: 1540 1541 foreach (i; 1 .. n - 1) 1542 { 1543 auto dx0 = xd[i - 1]; 1544 auto dx1 = xd[i - 0]; 1545 auto dy0 = yd[i - 1]; 1546 auto dy1 = yd[i - 0]; 1547 slopes[i] = 2 * (dx0 + dx1); 1548 temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0); 1549 } 1550 break; 1551 1552 case cardinal: 1553 1554 foreach (i; 1 .. n - 1) 1555 { 1556 slopes[i] = 1; 1557 temp[i] = (1 - param) * dd2[i - 1]; 1558 } 1559 break; 1560 1561 case monotone: 1562 { 1563 auto step0 = cast()xd[0]; 1564 auto step1 = cast()xd[1]; 1565 auto diff0 = cast()yd[0]; 1566 auto diff1 = cast()yd[1]; 1567 diff0 /= step0; 1568 diff1 /= step1; 1569 1570 for(size_t i = 1;;) 1571 { 1572 slopes[i] = 1; 1573 if (diff0 && diff1 && copysign(1f, diff0) == copysign(1f, diff1)) 1574 { 1575 auto w0 = step1 * 2 + step0; 1576 auto w1 = step0 * 2 + step1; 1577 temp[i] = (w0 + w1) / (w0 / diff0 + w1 / diff1); 1578 } 1579 else 1580 { 1581 temp[i] = 0; 1582 } 1583 if (++i == n - 1) 1584 { 1585 break; 1586 } 1587 step0 = step1; 1588 diff0 = diff1; 1589 step1 = xd[i]; 1590 diff1 = yd[i]; 1591 diff1 /= step1; 1592 } 1593 } 1594 break; 1595 1596 case doubleQuadratic: 1597 1598 foreach (i; 1 .. n - 1) 1599 { 1600 slopes[i] = 1; 1601 temp[i] = dd[i - 1] + dd[i] - dd2[i - 1]; 1602 } 1603 break; 1604 1605 case akima: 1606 { 1607 auto d3 = dd[1]; 1608 auto d2 = dd[0]; 1609 auto d1 = 2 * d2 - d3; 1610 auto d0 = d1; 1611 foreach (i; 1 .. n - 1) 1612 { 1613 d0 = d1; 1614 d1 = d2; 1615 d2 = d3; 1616 d3 = i == n - 2 ? 2 * d2 - d1 : dd[i + 1]; 1617 slopes[i] = 1; 1618 temp[i] = akimaSlope(d0, d1, d2, d3); 1619 } 1620 break; 1621 } 1622 1623 case makima: 1624 { 1625 auto d3 = dd[1]; 1626 auto d2 = dd[0]; 1627 auto d1 = 2 * d2 - d3; 1628 auto d0 = d1; 1629 foreach (i; 1 .. n - 1) 1630 { 1631 d0 = d1; 1632 d1 = d2; 1633 d2 = d3; 1634 d3 = i == n - 2 ? 2 * d2 - d1 : dd[i + 1]; 1635 slopes[i] = 1; 1636 temp[i] = makimaSlope(d0, d1, d2, d3); 1637 } 1638 break; 1639 } 1640 } 1641 1642 foreach (i; 0 .. n - 1) 1643 { 1644 auto c = i == 0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0; 1645 auto a = i == n - 2 ? last : kind == SplineType.c2 ? xd[i + 1] : 0; 1646 auto w = slopes[i] == 1 ? a : a / slopes[i]; 1647 slopes[i + 1] -= w * c; 1648 temp[i + 1] -= w * temp[i]; 1649 } 1650 1651 slopes.back = temp.back / slopes.back; 1652 1653 size_t convexCount; 1654 size_t concaveCount; 1655 foreach_reverse (i; 0 .. n - 1) 1656 { 1657 auto c = i == 0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0; 1658 auto v = temp[i] - c * slopes[i + 1]; 1659 slopes[i] = slopes[i] == 1 ? v : v / slopes[i]; 1660 1661 auto xdiff = points[i + 1] - points[i]; 1662 auto ydiff = values[i + 1] - values[i]; 1663 1664 auto convex1 = xdiff * (2 * slopes[i] + slopes[i + 1]) <= 3 * ydiff; 1665 auto concave1 = xdiff * (2 * slopes[i] + slopes[i + 1]) >= 3 * ydiff; 1666 auto convex2 = xdiff * (slopes[i] + 2 * slopes[i + 1]) >= 3 * ydiff; 1667 auto concave2 = xdiff * (slopes[i] + 2 * slopes[i + 1]) <= 3 * ydiff; 1668 convexCount += convex1 & convex2; 1669 convexCount += concave1 & concave2; 1670 } 1671 1672 return 1673 // convex kind has priority for the linear spline 1674 convexCount == n - 1 ? SplineConvexity.convex : 1675 concaveCount == n - 1 ? SplineConvexity.concave : 1676 SplineConvexity.none; 1677 } 1678 1679 private F akimaTail(F)(in F d2, in F d3) 1680 { 1681 auto d1 = 2 * d2 - d3; 1682 auto d0 = 2 * d1 - d2; 1683 return akimaSlope(d0, d1, d2, d3); 1684 } 1685 1686 private F akimaSlope(F)(in F d0, in F d1, in F d2, in F d3) 1687 { 1688 if (d1 == d2) 1689 return d1; 1690 if (d0 == d1 && d2 == d3) 1691 return (d1 + d2) * 0.5f; 1692 if (d0 == d1) 1693 return d1; 1694 if (d2 == d3) 1695 return d2; 1696 auto w0 = fabs(d1 - d0); 1697 auto w1 = fabs(d3 - d2); 1698 auto ws = w0 + w1; 1699 w0 /= ws; 1700 w1 /= ws; 1701 return w0 * d2 + w1 * d1; 1702 } 1703 1704 private F makimaTail(F)(in F d2, in F d3) 1705 { 1706 auto d1 = 2 * d2 - d3; 1707 auto d0 = 2 * d1 - d2; 1708 return makimaSlope(d0, d1, d2, d3); 1709 } 1710 1711 private F makimaSlope(F)(in F d0, in F d1, in F d2, in F d3) 1712 { 1713 if (d1 == d2) 1714 return d1; 1715 auto w0 = fabs(d1 - d0) + fabs(d1 + d0) * 0.5f; 1716 auto w1 = fabs(d3 - d2) + fabs(d3 + d2) * 0.5f; 1717 auto ws = w0 + w1; 1718 w0 /= ws; 1719 w1 /= ws; 1720 return w0 * d2 + w1 * d1; 1721 } 1722 1723 /// 1724 struct SplineKernel(X) 1725 { 1726 X step = 0; 1727 X w0 = 0; 1728 X w1 = 0; 1729 X wq = 0; 1730 1731 @fmamath: 1732 1733 /// 1734 this(X x0, X x1, X x) 1735 { 1736 step = x1 - x0; 1737 auto c0 = x - x0; 1738 auto c1 = x1 - x; 1739 w0 = c0 / step; 1740 w1 = c1 / step; 1741 wq = w0 * w1; 1742 } 1743 1744 /// 1745 template opCall(uint derivative = 0) 1746 if (derivative <= 3) 1747 { 1748 /// 1749 auto opCall(Y)(const Y y0, const Y y1, const Y s0, const Y s1) const 1750 { 1751 auto diff = y1 - y0; 1752 auto z0 = s0 * step - diff; 1753 auto z1 = s1 * step - diff; 1754 auto a0 = z0 * w1; 1755 auto a1 = z1 * w0; 1756 auto pr = a0 - a1; 1757 auto b0 = y0 * w1; 1758 auto b1 = y1 * w0; 1759 auto pl = b0 + b1; 1760 auto y = pl + wq * pr; 1761 static if (derivative) 1762 { 1763 Y[derivative + 1] ret = 0; 1764 ret[0] = y; 1765 auto wd = w1 - w0; 1766 auto zd = z1 + z0; 1767 ret[1] = (diff + (wd * pr - wq * zd)) / step; 1768 static if (derivative > 1) 1769 { 1770 auto astep = zd / (step * step); 1771 ret[2] = -3 * wd * astep + (s1 - s0) / step; 1772 static if (derivative > 2) 1773 ret[3] = 6 * astep / step; 1774 } 1775 return ret; 1776 } 1777 else 1778 { 1779 return y; 1780 } 1781 } 1782 1783 static if (derivative) 1784 auto opCall(Y, size_t N)(scope ref const Y[N] y0, scope ref const Y[N] y1, scope ref const Y[N] s0, scope ref const Y[N] s1) 1785 { 1786 Y[N][derivative + 1] ret = void; 1787 Y[derivative + 1][N] temp = void; 1788 1789 static foreach(i; 0 .. N) 1790 temp[i] = this.opCall!derivative(y0[i], y1[i], s0[i], s1[i]); 1791 static foreach(i; 0 .. N) 1792 static foreach(d; 0 .. derivative + 1) 1793 ret[d][i] = temp[i][d]; 1794 return ret; 1795 } 1796 } 1797 1798 /// 1799 alias withDerivative = opCall!1; 1800 /// 1801 alias withTwoDerivatives = opCall!2; 1802 } 1803 1804 package T pchipTail(T)(in T step0, in T step1, in T diff0, in T diff1) 1805 { 1806 import mir.math.common: copysign, fabs; 1807 if (!diff0) 1808 { 1809 return 0; 1810 } 1811 auto slope = ((step0 * 2 + step1) * diff0 - step0 * diff1) / (step0 + step1); 1812 if (copysign(1f, slope) != copysign(1f, diff0)) 1813 { 1814 return 0; 1815 } 1816 if ((copysign(1f, diff0) != copysign(1f, diff1)) && (fabs(slope) > fabs(diff0 * 3))) 1817 { 1818 return diff0 * 3; 1819 } 1820 return slope; 1821 } 1822 1823 /++ 1824 Spline interpolator used for non-rectiliner trapezoid-like greeds. 1825 Params: 1826 grid = rc-array of interpolation grid 1827 data = rc-array of interpolator-like structures 1828 typeOfBoundaries = $(LREF SplineBoundaryType) for both tails (optional). 1829 valueOfBoundaryConditions = value of the boundary type (optional). 1830 Constraints: 1831 `grid` and `values` must have the same length >= 3 1832 Returns: $(LREF Spline) 1833 +/ 1834 MetaSpline!(T, X) metaSpline(F, X, T)( 1835 RCArray!(immutable X) grid, 1836 RCArray!(const T) data, 1837 SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, 1838 const F valueOfBoundaryConditions = 0, 1839 ) 1840 { 1841 import core.lifetime: move; 1842 return metaSpline!(F, X, T)(grid.move, data.move, SplineType.c2, 0, typeOfBoundaries, valueOfBoundaryConditions); 1843 } 1844 1845 /++ 1846 Spline interpolator used for non-rectiliner trapezoid-like greeds. 1847 Params: 1848 grid = rc-array of interpolation grid 1849 data = rc-array of interpolator-like structures 1850 kind = $(LREF SplineType) type of cubic spline. 1851 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 1852 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 1853 typeOfBoundaries = $(LREF SplineBoundaryType) for both tails (optional). 1854 valueOfBoundaryConditions = value of the boundary type (optional). 1855 Constraints: 1856 `grid` and `values` must have the same length >= 3 1857 Returns: $(LREF Spline) 1858 +/ 1859 MetaSpline!(T, X) metaSpline(F, X, T)( 1860 RCArray!(immutable X) grid, 1861 RCArray!(const T) data, 1862 SplineType kind, 1863 const F param = 0, 1864 SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, 1865 const F valueOfBoundaryConditions = 0, 1866 ) 1867 { 1868 import core.lifetime: move; 1869 return metaSpline!(F, X, T)(grid.move, data.move, SplineBoundaryCondition!F(typeOfBoundaries, valueOfBoundaryConditions), kind, param); 1870 } 1871 1872 /++ 1873 Spline interpolator used for non-rectiliner trapezoid-like greeds. 1874 Params: 1875 grid = rc-array of interpolation grid 1876 data = rc-array of interpolator-like structures 1877 boundaries = $(LREF SplineBoundaryCondition) for both tails. 1878 kind = $(LREF SplineType) type of cubic spline. 1879 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 1880 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 1881 Constraints: 1882 `grid` and `values` must have the same length >= 3 1883 Returns: $(LREF Spline) 1884 +/ 1885 MetaSpline!(T, X) metaSpline(F, X, T)( 1886 RCArray!(immutable X) grid, 1887 RCArray!(const T) data, 1888 SplineBoundaryCondition!F boundaries, 1889 SplineType kind = SplineType.c2, 1890 const F param = 0, 1891 ) 1892 { 1893 import core.lifetime: move; 1894 return metaSpline!(F, X, T)(grid.move, data.move, SplineConfiguration!F(kind, boundaries, boundaries, param)); 1895 } 1896 1897 /++ 1898 Spline interpolator used for non-rectiliner trapezoid-like greeds. 1899 Params: 1900 grid = rc-array of interpolation grid 1901 data = rc-array of interpolator-like structures 1902 lBoundary = $(LREF SplineBoundaryCondition) for left tail. 1903 rBoundary = $(LREF SplineBoundaryCondition) for right tail. 1904 kind = $(LREF SplineType) type of cubic spline. 1905 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 1906 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 1907 Constraints: 1908 `grid` and `values` must have the same length >= 3 1909 Returns: $(LREF Spline) 1910 +/ 1911 MetaSpline!(T, X) metaSpline(F, X, T)( 1912 RCArray!(immutable X) grid, 1913 RCArray!(const T) data, 1914 SplineBoundaryCondition!F lBoundary, 1915 SplineBoundaryCondition!F rBoundary, 1916 SplineType kind = SplineType.c2, 1917 const F param = 0, 1918 ) 1919 { 1920 import core.lifetime: move; 1921 return metaSpline(grid.move, data.move, SplineConfiguration!F(kind, lBoundary, rBoundary, param)); 1922 } 1923 1924 /++ 1925 Spline interpolator used for non-rectiliner trapezoid-like greeds. 1926 Params: 1927 grid = rc-array of interpolation grid 1928 data = rc-array of interpolator-like structures 1929 configuration = $(LREF SplineConfiguration) 1930 Constraints: 1931 `grid` and `values` must have the same length >= 3 1932 Returns: $(LREF Spline) 1933 +/ 1934 MetaSpline!(T, X) metaSpline(F, X, T)( 1935 RCArray!(immutable X) grid, 1936 RCArray!(const T) data, 1937 SplineConfiguration!F configuration, 1938 ) 1939 { 1940 import core.lifetime: move; 1941 return MetaSpline!(T, X)(grid.move, data.move, configuration); 1942 } 1943 1944 /// ditto 1945 struct MetaSpline(T, X) 1946 { 1947 import mir.interpolate.utility: DeepType; 1948 // alias ElementInterpolator = Linear!(F, N, X); 1949 alias F = ValueType!(T, X); 1950 /// 1951 private Repeat!(3, Spline!F) splines; 1952 /// 1953 RCArray!(const T) data; 1954 // 1955 private RCArray!F _temp; 1956 /// 1957 SplineConfiguration!F configuration; 1958 1959 /// 1960 this( 1961 RCArray!(immutable X) grid, 1962 RCArray!(const T) data, 1963 SplineConfiguration!F configuration, 1964 ) 1965 { 1966 import core.lifetime: move; 1967 1968 if (grid.length < 2) 1969 { 1970 version(D_Exceptions) { import mir.exception : toMutable; throw exc_min.toMutable; } 1971 else assert(0, msg_min); 1972 } 1973 1974 if (grid.length != data.length) 1975 { 1976 version(D_Exceptions) { import mir.exception : toMutable; throw exc_eq.toMutable; } 1977 else assert(0, msg_eq); 1978 } 1979 1980 this.data = data.move; 1981 this._temp = grid.length; 1982 this.splines[0] = grid.asSlice; 1983 this.splines[1] = grid.asSlice; 1984 this.splines[2] = grid.moveToSlice; 1985 this.configuration = configuration; 1986 } 1987 1988 /// 1989 bool opEquals()(auto ref scope const typeof(this) rhs) scope const @trusted pure nothrow @nogc 1990 { 1991 return this.gridScopeView == rhs.gridScopeView 1992 && this.data == rhs.data 1993 && this.configuration == rhs.configuration; 1994 } 1995 1996 /// 1997 MetaLinear lightConst()() const @property { return *cast(MetaLinear*)&this; } 1998 1999 /// 2000 immutable(X)[] gridScopeView(size_t dimension = 0)() return scope const @property @trusted 2001 if (dimension == 0) 2002 { 2003 return splines[0].gridScopeView; 2004 } 2005 2006 /++ 2007 Returns: intervals count. 2008 +/ 2009 size_t intervalCount(size_t dimension = 0)() scope const @property 2010 if (dimension == 0) 2011 { 2012 assert(data.length > 1); 2013 return data.length - 1; 2014 } 2015 2016 /// 2017 enum uint derivativeOrder = 3; 2018 2019 static if (__traits(compiles, {enum N = T.dimensionCount;})) 2020 /// 2021 enum uint dimensionCount = T.dimensionCount + 1; 2022 2023 /// 2024 template opCall(uint derivative = 0) 2025 // if (derivative <= derivativeOrder) 2026 { 2027 /++ 2028 `(x)` operator. 2029 Complexity: 2030 `O(log(grid.length))` 2031 +/ 2032 auto opCall(X...)(const X xs) scope const @trusted 2033 // if (X.length == dimensionCount) 2034 { 2035 F[2][][derivative + 1] mutable; 2036 2037 static foreach (o; 0 .. derivative + 1) 2038 { 2039 mutable[o] = cast(F[2][]) splines[o]._data.lightScope.field; 2040 assert(mutable[o].length == data.length); 2041 } 2042 2043 static if (!derivative) 2044 { 2045 foreach (i, ref d; data) 2046 mutable[0][i][0] = d(xs[1 .. $]); 2047 } 2048 else 2049 { 2050 foreach (i, ref d; data) 2051 { 2052 auto node = d.opCall!derivative(xs[1 .. $]); 2053 static foreach (o; 0 .. derivative + 1) 2054 mutable[o][i][0] = node[o]; 2055 } 2056 } 2057 2058 static foreach (o; 0 .. derivative + 1) 2059 { 2060 (*cast(Spline!F*)&splines[o])._computeDerivativesTemp( 2061 configuration.kind, 2062 configuration.param, 2063 configuration.leftBoundary, 2064 configuration.rightBoundary, 2065 (cast(F[])_temp[]).sliced); 2066 } 2067 2068 static if (!derivative) 2069 { 2070 return splines[0](xs[0]); 2071 } 2072 else 2073 { 2074 typeof(splines[0].opCall!derivative(xs[0]))[derivative + 1] ret = void; 2075 static foreach (o; 0 .. derivative + 1) 2076 {{ 2077 auto s = splines[o].opCall!derivative(xs[0]); 2078 static foreach (r; 0 .. derivative + 1) 2079 ret[r][o] = s[r]; 2080 2081 }} 2082 return ret; 2083 } 2084 } 2085 } 2086 2087 /// 2088 alias withDerivative = opCall!1; 2089 /// 2090 alias withTwoDerivatives = opCall!2; 2091 } 2092 2093 /// 2D trapezoid-like (not rectilinear) linear interpolation 2094 version(mir_test) 2095 unittest 2096 { 2097 auto x = [ 2098 [0.0, 1, 2, 3, 5], 2099 [-4.0, 3, 4], 2100 [0.0, 10], 2101 ]; 2102 auto y = [ 2103 [4.0, 0, 9, 23, 40], 2104 [9.0, 0, 3], 2105 [4.0, 40], 2106 ]; 2107 2108 auto g = [7.0, 10, 15]; 2109 2110 import mir.rc.array: RCArray; 2111 import mir.ndslice.allocation: rcslice; 2112 2113 auto d = RCArray!(Spline!double)(3); 2114 2115 foreach (i; 0 .. x.length) 2116 d[i] = spline!double(x[i].rcslice!(immutable double), y[i].rcslice!(const double)); 2117 2118 auto trapezoidInterpolator = metaSpline!double(g.rcarray!(immutable double), d.lightConst); 2119 2120 auto val = trapezoidInterpolator(9.0, 1.8); 2121 } 2122 2123 version(mir_test) 2124 unittest 2125 { 2126 import mir.math.common: approxEqual; 2127 import mir.ndslice; 2128 alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10); 2129 2130 //// set test function //// 2131 enum y_x0 = 2; 2132 enum y_x1 = -7; 2133 enum y_x0x1 = 3; 2134 2135 // this function should be approximated very well 2136 alias f = (x0, x1) => y_x0 * x0 + y_x1 * x1 + y_x0x1 * x0 * x1 - 11; 2137 2138 ///// set interpolant //// 2139 static immutable x0 = [-1.0, 2, 8, 15]; 2140 static immutable x1 = [-4.0, 2, 5, 10, 13]; 2141 2142 auto grid = cartesian(x0, x1) 2143 .map!f 2144 .rcslice 2145 .lightConst; 2146 2147 auto int0 = spline!double(x1.rcslice!(immutable double), grid[0]); 2148 auto int1 = spline!double(x1.rcslice!(immutable double), grid[1]); 2149 auto int2 = spline!double(x1.rcslice!(immutable double), grid[2]); 2150 auto int3 = spline!double(x1.rcslice!(immutable double), grid[3]); 2151 2152 auto interpolant = metaSpline(x0.rcarray!(immutable double), rcarray(int0, int1, int2, int3).lightConst, SplineConfiguration!double.init); 2153 assert(interpolant == interpolant); 2154 2155 ///// compute test data //// 2156 auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 3.23); 2157 // auto test_grid = cartesian(x0 + 0, x1 + 0); 2158 auto real_data = test_grid.map!f; 2159 auto interp_data = test_grid.vmap(interpolant); 2160 2161 ///// verify result //// 2162 assert(all!appreq(interp_data, real_data)); 2163 2164 auto z0 = 1.23; 2165 auto z1 = 3.21; 2166 auto d = interpolant.withDerivative(z0, z1); 2167 assert(appreq(d[0][0], f(z0, z1))); 2168 assert(appreq(d[1][0], y_x0 + y_x0x1 * z1)); 2169 assert(appreq(d[0][1], y_x1 + y_x0x1 * z0)); 2170 assert(appreq(d[1][1], y_x0x1)); 2171 }