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 Ilya Yaroshenko, Kaleidic Associates Advisory Limited, Symmetry Investments 10 Authors: Ilya Yaroshenko 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.interpolate: Repeat; 23 import mir.math.common; 24 import mir.ndslice.slice; 25 import mir.primitives; 26 import mir.rc.array; 27 import mir.utility: min, max; 28 import std.meta: AliasSeq, staticMap; 29 import std.traits: Unqual; 30 public import mir.interpolate: atInterval; 31 32 @fmamath: 33 34 /// 35 @safe pure @nogc version(mir_test) unittest 36 { 37 import mir.algorithm.iteration: all; 38 import mir.math.common: approxEqual; 39 import mir.ndslice.slice: sliced; 40 import mir.ndslice.allocation: rcslice; 41 import mir.ndslice.topology: vmap; 42 43 static immutable xdata = [-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 44 auto x = xdata.rcslice; 45 static immutable ydata = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 46 auto y = ydata.sliced; 47 48 auto interpolant = spline!double(x, y); // constructs Spline 49 auto xs = x + 0.5; // input X values for cubic spline 50 51 static immutable test_data0 = [ 52 -0.68361541, 7.28568719, 10.490694 , 0.36192032, 53 11.91572713, 16.44546433, 17.66699525, 4.52730869, 54 19.22825394, -2.3242592 ]; 55 /// not-a-knot (default) 56 assert(xs.vmap(interpolant).all!approxEqual(test_data0)); 57 58 static immutable test_data1 = [ 59 10.85298372, 5.26255911, 10.71443229, 0.1824536 , 60 11.94324989, 16.45633939, 17.59185094, 4.86340188, 61 17.8565408 , 2.81856494]; 62 /// natural cubic spline 63 interpolant = spline!double(x, y, SplineBoundaryType.secondDerivative); 64 assert(xs.vmap(interpolant).all!approxEqual(test_data1)); 65 66 static immutable test_data2 = [ 67 9.94191781, 5.4223652 , 10.69666392, 0.1971149 , 11.93868415, 68 16.46378847, 17.56521661, 4.97656997, 17.39645585, 4.54316446]; 69 /// set both boundary second derivatives to 3 70 interpolant = spline!double(x, y, SplineBoundaryType.secondDerivative, 3); 71 assert(xs.vmap(interpolant).all!approxEqual(test_data2)); 72 73 static immutable test_data3 = [ 74 16.45728263, 4.27981687, 10.82295092, 0.09610695, 75 11.95252862, 16.47583126, 17.49964521, 5.26561539, 76 16.21803478, 8.96104974]; 77 /// set both boundary derivatives to 3 78 interpolant = spline!double(x, y, SplineBoundaryType.firstDerivative, 3); 79 assert(xs.vmap(interpolant).all!approxEqual(test_data3)); 80 81 static immutable test_data4 = [ 82 16.45730084, 4.27966112, 10.82337171, 0.09403945, 83 11.96265209, 16.44067375, 17.6374694 , 4.67438921, 84 18.6234452 , -0.05582876]; 85 /// different boundary conditions 86 interpolant = spline!double(x, y, 87 SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3), 88 SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5)); 89 assert(xs.vmap(interpolant).all!approxEqual(test_data4)); 90 91 92 static immutable test_data5 = [ 93 12.37135558, 4.99638066, 10.74362441, 0.16008641, 94 11.94073593, 16.47908148, 17.49841853, 5.26600921, 95 16.21796051, 8.96102894]; 96 // ditto 97 interpolant = spline!double(x, y, 98 SplineBoundaryCondition!double(SplineBoundaryType.secondDerivative, -5), 99 SplineBoundaryCondition!double(SplineBoundaryType.firstDerivative, 3)); 100 assert(xs.vmap(interpolant).all!approxEqual(test_data5)); 101 102 static immutable test_data6 = [ 103 11.40871379, 2.64278898, 9.55774317, 4.84791141, 11.24842121, 104 16.16794267, 18.58060557, 5.2531411 , 17.45509005, 1.86992521]; 105 /// Akima spline 106 interpolant = spline!double(x, y, SplineType.akima); 107 assert(xs.vmap(interpolant).all!approxEqual(test_data6)); 108 109 /// Double Quadratic spline 110 interpolant = spline!double(x, y, SplineType.doubleQuadratic); 111 import mir.interpolate.utility: ParabolaKernel; 112 auto kernel1 = ParabolaKernel!double(x[2], x[3], x[4], y[2], y[3], y[4]); 113 auto kernel2 = ParabolaKernel!double( x[3], x[4], x[5], y[3], y[4], y[5]); 114 // weighted sum of quadratic functions 115 auto c = 0.35; // from [0 .. 1] 116 auto xp = c * x[3] + (1 - c) * x[4]; 117 auto yp = c * kernel1(xp) + (1 - c) * kernel2(xp); 118 assert(interpolant(xp).approxEqual(yp)); 119 // check parabolic extrapolation of the boundary intervals 120 kernel1 = ParabolaKernel!double(x[0], x[1], x[2], y[0], y[1], y[2]); 121 kernel2 = ParabolaKernel!double(x[$ - 3], x[$ - 2], x[$ - 1], y[$ - 3], y[$ - 2], y[$ - 1]); 122 assert(interpolant(x[0] - 23.421).approxEqual(kernel1(x[0] - 23.421))); 123 assert(interpolant(x[$ - 1] + 23.421).approxEqual(kernel2(x[$ - 1] + 23.421))); 124 } 125 126 /// 127 @safe pure version(mir_test) unittest 128 { 129 import mir.rc.array: rcarray; 130 import mir.algorithm.iteration: all; 131 import mir.functional: aliasCall; 132 import mir.math.common: approxEqual; 133 import mir.ndslice.allocation: uninitSlice; 134 import mir.ndslice.slice: sliced; 135 import mir.ndslice.topology: vmap, map; 136 137 auto x = rcarray!(immutable double)(-1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22).asSlice; 138 auto y = rcarray( 139 8.77842512, 140 7.96429686, 141 7.77074363, 142 1.10838032, 143 2.69925191, 144 1.84922654, 145 1.48167283, 146 2.8267636 , 147 0.40200172, 148 7.78980608).asSlice; 149 150 auto interpolant = x.spline!double(y); // default boundary condition is 'not-a-knot' 151 152 auto xs = x + 0.5; 153 154 auto ys = xs.vmap(interpolant); 155 156 auto r = 157 [5.56971848, 158 9.30342403, 159 4.44139761, 160 -0.74740285, 161 3.00994108, 162 1.50750417, 163 1.73144979, 164 2.64860361, 165 0.64413911, 166 10.81768928]; 167 168 assert(all!approxEqual(ys, r)); 169 170 // first derivative 171 auto d1 = xs.vmap(interpolant.aliasCall!"withDerivative").map!"a[1]"; 172 auto r1 = 173 [-4.51501279, 174 2.15715986, 175 -7.28363308, 176 -2.14050449, 177 0.03693092, 178 -0.49618999, 179 0.58109933, 180 -0.52926703, 181 0.7819035 , 182 6.70632693]; 183 assert(all!approxEqual(d1, r1)); 184 185 // second derivative 186 auto d2 = xs.vmap(interpolant.aliasCall!"withTwoDerivatives").map!"a[2]"; 187 auto r2 = 188 [7.07104751, 189 -2.62293241, 190 -0.01468508, 191 5.70609505, 192 -2.02358911, 193 0.72142061, 194 0.25275483, 195 -0.6133589 , 196 1.26894416, 197 2.68067146]; 198 assert(all!approxEqual(d2, r2)); 199 200 // third derivative (6 * a) 201 auto d3 = xs.vmap(interpolant.aliasCall!("opCall", 3)).map!"a[3]"; 202 auto r3 = 203 [-3.23132664, 204 -3.23132664, 205 14.91047457, 206 -3.46891432, 207 1.88520325, 208 -0.16559031, 209 -0.44056064, 210 0.47057577, 211 0.47057577, 212 0.47057577]; 213 assert(all!approxEqual(d3, r3)); 214 } 215 216 /// R -> R: Cubic interpolation 217 version(mir_test) 218 @safe unittest 219 { 220 import mir.algorithm.iteration: all; 221 import mir.math.common: approxEqual; 222 import mir.ndslice; 223 224 static immutable x = [0, 1, 2, 3, 5.00274, 7.00274, 10.0055, 20.0137, 30.0192]; 225 static immutable y = [0.0011, 0.0011, 0.0030, 0.0064, 0.0144, 0.0207, 0.0261, 0.0329, 0.0356,]; 226 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]; 227 228 auto interpolation = spline!double(x.rcslice, y.sliced); 229 230 auto data = 231 [ 0.0011 , 0.003 , 0.0064 , 0.01042622, 0.0144 , 232 0.01786075, 0.0207 , 0.02293441, 0.02467983, 0.0261 , 233 0.02732764, 0.02840225, 0.0293308 , 0.03012914, 0.03081002, 234 0.03138766, 0.03187161, 0.03227637, 0.03261468, 0.0329 , 235 0.03314357, 0.03335896, 0.03355892, 0.03375674, 0.03396413, 236 0.03419436, 0.03446018, 0.03477529, 0.03515072, 0.0356 ]; 237 238 assert(all!approxEqual(xs.sliced.vmap(interpolation), data)); 239 } 240 241 /// R^2 -> R: Bicubic interpolation 242 version(mir_test) 243 unittest 244 { 245 import mir.math.common: approxEqual; 246 import mir.ndslice; 247 alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10); 248 249 ///// set test function //// 250 const double y_x0 = 2; 251 const double y_x1 = -7; 252 const double y_x0x1 = 3; 253 254 // this function should be approximated very well 255 alias f = (x0, x1) => y_x0 * x0 + y_x1 * x1 + y_x0x1 * x0 * x1 - 11; 256 257 ///// set interpolant //// 258 static immutable x0 = [-1.0, 2, 8, 15]; 259 static immutable x1 = [-4.0, 2, 5, 10, 13]; 260 auto grid = cartesian(x0, x1); 261 262 auto interpolant = spline!(double, 2)(x0.rcslice, x1.rcslice, grid.map!f); 263 264 ///// compute test data //// 265 auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 3.23); 266 // auto test_grid = cartesian(x0 + 0, x1 + 0); 267 auto real_data = test_grid.map!f; 268 auto interp_data = test_grid.vmap(interpolant); 269 270 ///// verify result //// 271 assert(all!appreq(interp_data, real_data)); 272 273 //// check derivatives //// 274 auto z0 = 1.23; 275 auto z1 = 3.21; 276 // writeln("-----------------"); 277 // writeln("-----------------"); 278 auto d = interpolant.withDerivative(z0, z1); 279 assert(appreq(interpolant(z0, z1), f(z0, z1))); 280 // writeln("d = ", d); 281 // writeln("interpolant.withTwoDerivatives(z0, z1) = ", interpolant.withTwoDerivatives(z0, z1)); 282 // writeln("-----------------"); 283 // writeln("-----------------"); 284 // writeln("interpolant(z0, z1) = ", interpolant(z0, z1)); 285 // writeln("y_x0 + y_x0x1 * z1 = ", y_x0 + y_x0x1 * z1); 286 // writeln("y_x1 + y_x0x1 * z0 = ", y_x1 + y_x0x1 * z0); 287 // writeln("-----------------"); 288 // writeln("-----------------"); 289 // assert(appreq(d[0][0], f(z0, z1))); 290 // assert(appreq(d[1][0], y_x0 + y_x0x1 * z1)); 291 // assert(appreq(d[0][1], y_x1 + y_x0x1 * z0)); 292 // assert(appreq(d[1][1], y_x0x1)); 293 } 294 295 /// R^3 -> R: Tricubic interpolation 296 version(mir_test) 297 unittest 298 { 299 import mir.math.common: approxEqual; 300 import mir.ndslice; 301 alias appreq = (a, b) => approxEqual(a, b, 10e-10, 10e-10); 302 303 ///// set test function //// 304 enum y_x0 = 2; 305 enum y_x1 = -7; 306 enum y_x2 = 5; 307 enum y_x0x1 = 10; 308 enum y_x0x1x2 = 3; 309 310 // this function should be approximated very well 311 alias f = (x0, x1, x2) => y_x0 * x0 + y_x1 * x1 + y_x2 * x2 312 + y_x0x1 * x0 * x1 + y_x0x1x2 * x0 * x1 * x2 - 11; 313 314 ///// set interpolant //// 315 static immutable x0 = [-1.0, 2, 8, 15]; 316 static immutable x1 = [-4.0, 2, 5, 10, 13]; 317 static immutable x2 = [3, 3.7, 5]; 318 auto grid = cartesian(x0, x1, x2); 319 320 auto interpolant = spline!(double, 3)(x0.rcslice, x1.rcslice, x2.rcslice, grid.map!f); 321 322 ///// compute test data //// 323 auto test_grid = cartesian(x0.sliced + 1.23, x1.sliced + 3.23, x2.sliced - 3); 324 auto real_data = test_grid.map!f; 325 auto interp_data = test_grid.vmap(interpolant); 326 327 ///// verify result //// 328 assert(all!appreq(interp_data, real_data)); 329 330 //// check derivatives //// 331 auto z0 = 1.23; 332 auto z1 = 3.23; 333 auto z2 = -3; 334 auto d = interpolant.withDerivative(z0, z1, z2); 335 assert(appreq(interpolant(z0, z1, z2), f(z0, z1, z2))); 336 assert(appreq(d[0][0][0], f(z0, z1, z2))); 337 338 // writeln("-----------------"); 339 // writeln("-----------------"); 340 // auto d = interpolant.withDerivative(z0, z1); 341 assert(appreq(interpolant(z0, z1, z2), f(z0, z1, z2))); 342 // writeln("interpolant(z0, z1, z2) = ", interpolant(z0, z1, z2)); 343 // writeln("d = ", d); 344 // writeln("interpolant.withTwoDerivatives(z0, z1, z2) = ", interpolant.withTwoDerivatives(z0, z1, z2)); 345 // writeln("-----------------"); 346 // writeln("-----------------"); 347 // writeln("interpolant(z0, z1) = ", interpolant(z0, z1)); 348 // writeln("y_x0 + y_x0x1 * z1 = ", y_x0 + y_x0x1 * z1); 349 // writeln("y_x1 + y_x0x1 * z0 = ", y_x1 + y_x0x1 * z0); 350 // writeln("-----------------"); 351 // writeln("-----------------"); 352 353 // writeln("y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2 = ", y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2); 354 // assert(appreq(d[1][0][0], y_x0 + y_x0x1 * z1 + y_x0x1x2 * z1 * z2)); 355 // writeln("y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2 = ", y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2); 356 // assert(appreq(d[0][1][0], y_x1 + y_x0x1 * z0 + y_x0x1x2 * z0 * z2)); 357 // writeln("y_x0x1 + y_x0x1x2 * z2 = ", y_x0x1 + y_x0x1x2 * z2); 358 // assert(appreq(d[1][1][0], y_x0x1 + y_x0x1x2 * z2)); 359 // writeln("y_x0x1x2 = ", y_x0x1x2); 360 // assert(appreq(d[1][1][1], y_x0x1x2)); 361 } 362 363 364 /// Monotone PCHIP 365 version(mir_test) 366 @safe unittest 367 { 368 import mir.math.common: approxEqual; 369 import mir.algorithm.iteration: all; 370 import mir.ndslice.allocation: rcslice; 371 import mir.ndslice.slice: sliced; 372 import mir.ndslice.topology: vmap; 373 374 static immutable x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 375 static immutable y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 376 auto interpolant = spline!double(x.rcslice, y.sliced, SplineType.monotone); 377 378 auto xs = x[0 .. $ - 1].sliced + 0.5; 379 380 auto ys = xs.vmap(interpolant); 381 382 assert(ys.all!approxEqual([ 383 5.333333333333334, 384 2.500000000000000, 385 10.000000000000000, 386 4.288971807628524, 387 11.202580845771145, 388 16.250000000000000, 389 17.962962962962962, 390 5.558593750000000, 391 17.604662698412699, 392 ])); 393 } 394 395 // Check direction equality 396 version(mir_test) 397 @safe unittest 398 { 399 import mir.math.common: approxEqual; 400 import mir.ndslice.slice: sliced; 401 import mir.ndslice.allocation: rcslice; 402 import mir.ndslice.topology: retro, vmap; 403 404 static immutable points = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 405 static immutable values = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 406 407 auto results = [ 408 5.333333333333334, 409 2.500000000000000, 410 10.000000000000000, 411 4.288971807628524, 412 11.202580845771145, 413 16.250000000000000, 414 17.962962962962962, 415 5.558593750000000, 416 17.604662698412699, 417 ]; 418 auto interpolant = spline!double(points.rcslice, values.sliced, SplineType.monotone); 419 420 auto pointsR = rcslice(-points.retro); 421 auto valuesR = values.retro.rcslice; 422 auto interpolantR = spline!double(pointsR, valuesR, SplineType.monotone); 423 424 version(X86_64) 425 assert(vmap(points[0 .. $ - 1].sliced + 0.5, interpolant) == vmap(pointsR.retro[0 .. $ - 1] - 0.5, interpolantR)); 426 } 427 428 /++ 429 Cubic Spline types. 430 431 The first derivatives are guaranteed to be continuous for all cubic splines. 432 +/ 433 extern(C++, "mir", "interpolate") 434 enum SplineType 435 { 436 /++ 437 Spline with contiguous second derivative. 438 +/ 439 c2, 440 /++ 441 $(HTTPS en.wikipedia.org/wiki/Cubic_Hermite_spline#Cardinal_spline, Cardinal) and Catmull–Rom splines. 442 +/ 443 cardinal, 444 /++ 445 The interpolant preserves monotonicity in the interpolation data and does not overshoot if the data is not smooth. 446 It is also known as $(HTTPS docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.PchipInterpolator.html, PCHIP) 447 in numpy and Matlab. 448 +/ 449 monotone, 450 /++ 451 Weighted sum of two nearbor quadratic functions. 452 It is used in $(HTTPS s3-eu-west-1.amazonaws.com/og-public-downloads/smile-interpolation-extrapolation.pdf, financial analysis). 453 +/ 454 doubleQuadratic, 455 /++ 456 $(HTTPS en.wikipedia.org/wiki/Akima_spline, Akima spline). 457 +/ 458 akima, 459 } 460 461 /++ 462 Constructs multivariate cubic spline in symmetrical form with nodes on rectilinear grid. 463 Result has continues second derivatives throughout the curve / nd-surface. 464 +/ 465 template spline(T, size_t N = 1, X = T) 466 if (isFloatingPoint!T && is(T == Unqual!T) && N <= 6) 467 { 468 /++ 469 Params: 470 grid = immutable `x` values for interpolant 471 values = `f(x)` values for interpolant 472 typeOfBoundaries = $(LREF SplineBoundaryType) for both tails (optional). 473 valueOfBoundaryConditions = value of the boundary type (optional). 474 Constraints: 475 `grid` and `values` must have the same length >= 3 476 Returns: $(LREF Spline) 477 +/ 478 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 479 Repeat!(N, Slice!(RCI!(immutable X))) grid, 480 Slice!(yIterator, N, ykind) values, 481 SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, 482 in T valueOfBoundaryConditions = 0, 483 ) 484 { 485 return spline(grid, values, SplineType.c2, 0, typeOfBoundaries, valueOfBoundaryConditions); 486 } 487 488 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 489 Repeat!(N, Slice!(RCI!(immutable X))) grid, 490 Slice!(yIterator, N, ykind) values, 491 SplineType kind, 492 in T param = 0, 493 SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot, 494 in T valueOfBoundaryConditions = 0, 495 ) 496 { 497 return spline(grid, values, SplineBoundaryCondition!T(typeOfBoundaries, valueOfBoundaryConditions), kind, param); 498 } 499 500 /++ 501 Params: 502 grid = immutable `x` values for interpolant 503 values = `f(x)` values for interpolant 504 boundaries = $(LREF SplineBoundaryCondition) for both tails. 505 kind = $(LREF SplineType) type of cubic spline. 506 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 507 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 508 Constraints: 509 `grid` and `values` must have the same length >= 3 510 Returns: $(LREF Spline) 511 +/ 512 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 513 Repeat!(N, Slice!(RCI!(immutable X))) grid, 514 Slice!(yIterator, N, ykind) values, 515 SplineBoundaryCondition!T boundaries, 516 SplineType kind = SplineType.c2, 517 in T param = 0, 518 ) 519 { 520 return spline(grid, values, boundaries, boundaries, kind, param); 521 } 522 523 /++ 524 Params: 525 grid = immutable `x` values for interpolant 526 values = `f(x)` values for interpolant 527 rBoundary = $(LREF SplineBoundaryCondition) for left tail. 528 lBoundary = $(LREF SplineBoundaryCondition) for right tail. 529 kind = $(LREF SplineType) type of cubic spline. 530 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 531 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 532 Constraints: 533 `grid` and `values` must have the same length >= 3 534 Returns: $(LREF Spline) 535 +/ 536 Spline!(T, N, X) spline(yIterator, SliceKind ykind)( 537 Repeat!(N, Slice!(RCI!(immutable X))) grid, 538 Slice!(yIterator, N, ykind) values, 539 SplineBoundaryCondition!T rBoundary, 540 SplineBoundaryCondition!T lBoundary, 541 SplineType kind = SplineType.c2, 542 in T param = 0, 543 ) 544 { 545 auto ret = typeof(return)(forward!grid); 546 ret._values = values; 547 ret._computeDerivatives(kind, param, rBoundary, lBoundary); 548 return ret; 549 } 550 } 551 552 /++ 553 Cubic Spline Boundary Condition Type. 554 555 See_also: $(LREF SplineBoundaryCondition) $(LREF SplineType) 556 +/ 557 extern(C++, "mir", "interpolate") 558 enum SplineBoundaryType 559 { 560 /++ 561 Not implemented. 562 +/ 563 periodic = -1, 564 /++ 565 Not-a-knot (or cubic) boundary condition. 566 It is an aggresive boundary condition that is used only for C2 splines and is default for all API calls. 567 For other then C2 splines, `notAKnot` is changed internally to 568 a default boundary type for used $(LREF SplineType). 569 +/ 570 notAKnot, 571 /++ 572 Set the first derivative. 573 +/ 574 firstDerivative, 575 /++ 576 Set the second derivative. 577 +/ 578 secondDerivative, 579 /++ 580 Default for Cardinal and Double-Quadratic splines. 581 +/ 582 parabolic, 583 /++ 584 Default for monotone (aka PHCIP ) splines. 585 +/ 586 monotone, 587 /++ 588 Default for Akima splines. 589 +/ 590 akima, 591 } 592 593 /++ 594 Cubic Spline Boundary Condition 595 596 See_also: $(LREF SplineBoundaryType) 597 +/ 598 extern(C++, "mir", "interpolate") 599 struct SplineBoundaryCondition(T) 600 { 601 /// type (default is $(LREF SplineBoundaryType.notAKnot)) 602 SplineBoundaryType type = SplineBoundaryType.notAKnot; 603 /// value (default is 0) 604 T value = 0; 605 } 606 607 /++ 608 Multivariate cubic spline with nodes on rectilinear grid. 609 +/ 610 struct Spline(F, size_t N = 1, X = F) 611 if (N && N <= 6) 612 { 613 import mir.rc.array; 614 615 /// Aligned buffer allocated with `mir.internal.memory`. $(RED For internal use.) 616 Slice!(RCI!(F[2 ^^ N]), N) _data; 617 /// Grid iterators. $(RED For internal use.) 618 Repeat!(N, RCI!(immutable X)) _grid; 619 620 @fmamath extern(D): 621 622 /++ 623 +/ 624 this(Repeat!(N, Slice!(RCI!(immutable X))) grid) @safe @nogc 625 { 626 size_t length = 1; 627 size_t[N] shape; 628 enum msg = "spline interpolant: minimal allowed length for the grid equals 2."; 629 version(D_Exceptions) 630 static immutable exc = new Exception(msg); 631 foreach(i, ref x; grid) 632 { 633 if (x.length < 2) 634 { 635 version(D_Exceptions) throw exc; 636 else assert(0, msg); 637 } 638 length *= shape[i] = x.length; 639 this._grid[i] = x._iterator.move; 640 } 641 import mir.ndslice.allocation: rcslice; 642 this._data = shape.rcslice!(F[2 ^^ N]); 643 } 644 645 package static auto pickDataSubslice(D)(auto scope ref D data, size_t index) @trusted 646 { 647 auto strides = data.strides; 648 foreach (i; Iota!(strides.length)) 649 strides[i] *= DeepElementType!D.length; 650 return Slice!(F*, strides.length, Universal)(data.shape, strides, data._iterator.ptr + index); 651 } 652 653 /++ 654 Assigns function values to the internal memory. 655 $(RED For internal use.) 656 +/ 657 void _values(SliceKind kind, Iterator)(Slice!(Iterator, N, kind) values) scope @property @trusted 658 { 659 assert(values.shape == _data.shape, "'values' should have the same shape as the .gridShape"); 660 pickDataSubslice(_data.lightScope, 0)[] = values; 661 } 662 663 /++ 664 Computes derivatives and stores them in `_data`. 665 `_data` is assumed to be preinitialized with function values filled in `F[2 ^^ N][0]`. 666 Params: 667 lbc = left boundary condition 668 rbc = right boundary condition 669 temp = temporal buffer length points count (optional) 670 671 $(RED For internal use.) 672 +/ 673 void _computeDerivatives(SplineType kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc) scope @trusted nothrow @nogc 674 { 675 import mir.algorithm.iteration: maxLength; 676 auto ml = this._data.maxLength; 677 auto temp = RCArray!F(ml); 678 auto tempSlice = temp[].sliced; 679 _computeDerivativesTemp(kind, param, lbc, rbc, tempSlice); 680 } 681 682 /// ditto 683 pragma(inline, false) 684 void _computeDerivativesTemp(SplineType kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, Slice!(F*) temp) scope @system nothrow @nogc 685 { 686 import mir.algorithm.iteration: maxLength, each; 687 import mir.ndslice.topology: map, byDim, evertPack; 688 689 assert(temp.length >= _data.maxLength); 690 691 static if (N == 1) 692 { 693 splineSlopes!(F, F)(_grid[0].lightConst.sliced(_data._lengths[0]), pickDataSubslice(_data.lightScope, 0), pickDataSubslice(_data.lightScope, 1), temp[0 .. _data._lengths[0]], kind, param, lbc, rbc); 694 } 695 else 696 foreach_reverse(i; Iota!N) 697 { 698 // if (i == _grid.length - 1) 699 _data 700 .lightScope 701 .byDim!i 702 .evertPack 703 .each!((d){ 704 enum L = 2 ^^ (N - 1 - i); 705 foreach(l; Iota!L) 706 { 707 auto y = pickDataSubslice(d, l); 708 auto s = pickDataSubslice(d, L + l); 709 // 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); 710 splineSlopes!(F, F)(_grid[i].sliced(_data._lengths[i]), y, s, temp[0 .. _data._lengths[i]], kind, param, lbc, rbc); 711 // debug{ 712 // (cast(void delegate() @nogc)(){ 713 // writeln("y = ", y); 714 // writeln("s = ", s); 715 // })(); 716 // } 717 } 718 }); 719 } 720 } 721 722 @trusted: 723 724 /// 725 Spline lightConst() const @property { return *cast(Spline*)&this; } 726 /// 727 Spline lightImmutable() immutable @property { return *cast(Spline*)&this; } 728 729 /// 730 Slice!(RCI!(immutable X)) grid(size_t dimension = 0)() scope return const @property 731 if (dimension < N) 732 { 733 return _grid[dimension].lightConst.sliced(_data._lengths[dimension]); 734 } 735 736 /// 737 immutable(X)[] gridScopeView(size_t dimension = 0)() scope return const @property @trusted 738 if (dimension < N) 739 { 740 return _grid[dimension]._iterator[0 .. _data._lengths[dimension]]; 741 } 742 743 /++ 744 Returns: intervals count. 745 +/ 746 size_t intervalCount(size_t dimension = 0)() scope const @property 747 { 748 assert(_data._lengths[dimension] > 1); 749 return _data._lengths[dimension] - 1; 750 } 751 752 /// 753 size_t[N] gridShape() scope const @property 754 { 755 return _data.shape; 756 } 757 758 /// 759 alias withDerivative = opCall!1; 760 /// 761 alias withTwoDerivatives = opCall!2; 762 763 /// 764 enum uint derivativeOrder = 3; 765 766 /// 767 template opCall(uint derivative : 2) 768 { 769 auto opCall(X...)(in X xs) scope const 770 if (X.length == N) 771 // @FUTURE@ 772 // X.length == N || derivative == 0 && X.length && X.length <= N 773 { 774 auto d4 = this.opCall!3(xs); 775 SplineReturnType!(F, N, 3) d3; 776 void fun(size_t d, A, B)(ref A a, ref B b) 777 { 778 static if (d) 779 foreach(i; Iota!3) 780 fun!(d - 1)(a[i], b[i]); 781 else 782 b = a; 783 } 784 fun!N(d4, d3); 785 return d3; 786 } 787 } 788 789 /// 790 template opCall(uint derivative = 0) 791 if (derivative == 0 || derivative == 1 || derivative == 3) 792 { 793 static if (N > 1 && derivative) pragma(msg, "Warning: multivariate cubic spline with derivatives was not tested!!!"); 794 795 /++ 796 `(x)` operator. 797 Complexity: 798 `O(log(points.length))` 799 +/ 800 auto opCall(X...)(in X xs) scope const @trusted 801 if (X.length == N) 802 // @FUTURE@ 803 // X.length == N || derivative == 0 && X.length && X.length <= N 804 { 805 import mir.ndslice.topology: iota; 806 alias Kernel = AliasCall!(SplineKernel!F, "opCall", derivative); 807 enum rp2d = derivative == 3 ? 2 : derivative; 808 809 size_t[N] indices; 810 Kernel[N] kernels; 811 812 foreach(i; Iota!N) 813 { 814 static if (isInterval!(typeof(xs[i]))) 815 { 816 indices[i] = xs[i][1]; 817 auto x = xs[i][0]; 818 } 819 else 820 { 821 alias x = xs[i]; 822 indices[i] = this.findInterval!i(x); 823 } 824 kernels[i] = SplineKernel!F(_grid[i][indices[i]], _grid[i][indices[i] + 1], x); 825 } 826 827 align(64) F[2 ^^ N * 2 ^^ N][2] local; 828 829 void load(sizediff_t i)(const(F[2 ^^ N])* from, F[2 ^^ N]* to) 830 { 831 version(LDC) pragma(inline, true); 832 static if (i == -1) 833 { 834 // copyvec(*from, *to); 835 // not aligned: 836 *to = *from; 837 } 838 else 839 { 840 from += strides[i] * indices[i]; 841 load!(i - 1)(from, to); 842 from += strides[i]; 843 enum s = 2 ^^ (N - 1 - i); 844 to += s; 845 load!(i - 1)(from, to); 846 } 847 } 848 849 immutable strides = _data._lengths.iota.strides; 850 load!(N - 1)(_data.ptr, cast(F[2 ^^ N]*) local[0].ptr); 851 852 // debug{ 853 854 // printf("0local[0] = "); 855 // foreach(ref e; local[0][]) 856 // printf("%f ", e); 857 // printf("\n"); 858 // } 859 860 foreach(i; Iota!N) 861 { 862 enum P = 2 ^^ (N - 1 - i) * 2 ^^ (i * rp2d); 863 enum L = (2 ^^ N) ^^ 2 / (2 ^^ (i * (2 - rp2d))) / 4; 864 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]); 865 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]); 866 // debug 867 // { 868 // printf("0local[1] = "); 869 // foreach(ref e; local[1][0 .. L* 4]) 870 // printf("%f ", e); 871 // printf("\n"); 872 // } 873 local[0][] = F.init; 874 vectorize( 875 kernels[i], 876 local[1][0 * L .. 1 * L], local[1][2 * L .. 3 * L], 877 local[1][1 * L .. 2 * L], local[1][3 * L .. 4 * L], 878 *cast(F[L][2 ^^ rp2d]*) local[0].ptr, 879 ); 880 881 // debug{ 882 883 // printf("1local[0] = "); 884 // foreach(ref e; local[0][0 .. L* 2 ^^ rp2d]) 885 // printf("%f ", e); 886 // printf("\n"); 887 // } 888 // printf("local[0][0]", local[0][0]); 889 static if (i + 1 == N) 890 { 891 return *cast(SplineReturnType!(F, N, 2 ^^ rp2d)*) local[0].ptr; 892 } 893 else 894 { 895 static if (rp2d == 1) 896 { 897 shuffle3!1(local[0][0 .. L], local[0][L .. 2 * L], local[1][0 .. L], local[1][L .. 2 * L]); 898 copyvec(local[1][0 .. 1 * L], local[0][0 .. 1 * L]); 899 copyvec(local[1][L .. 2 * L], local[0][L .. 2 * L]); 900 } 901 else 902 static if (rp2d == 2) 903 { 904 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]); 905 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]); 906 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]); 907 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]); 908 } 909 910 // debug{ 911 912 // printf("2local[0] = "); 913 // foreach(ref e; local[0][0 .. L * 2 ^^ rp2d]) 914 // printf("%f ", e); 915 // printf("\n"); 916 // } 917 918 } 919 } 920 } 921 } 922 } 923 924 /++ 925 Piecewise cubic hermite interpolating polynomial. 926 Params: 927 points = `x` values for interpolant 928 values = `f(x)` values for interpolant 929 slopes = uninitialized ndslice to write slopes into 930 temp = uninitialized temporary ndslice 931 kind = $(LREF SplineType) type of cubic spline. 932 param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types). 933 Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline. 934 lbc = left boundary condition 935 rbc = right boundary condition 936 Constraints: 937 `points`, `values`, and `slopes`, must have the same length > 3; 938 `temp` must have length greater or equal to points less minus one. 939 +/ 940 void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind skind)( 941 Slice!(IP, 1, gkind) points, 942 Slice!(IV, 1, vkind) values, 943 Slice!(IS, 1, skind) slopes, 944 Slice!(T*) temp, 945 SplineType kind, 946 F param, 947 SplineBoundaryCondition!F lbc, 948 SplineBoundaryCondition!F rbc, 949 ) @trusted 950 { 951 import mir.ndslice.topology: diff, zip, slide; 952 953 assert (points.length >= 2); 954 assert (points.length == values.length); 955 assert (points.length == slopes.length); 956 assert (temp.length == points.length); 957 958 auto n = points.length; 959 960 typeof(slopes[0]) first, last; 961 962 auto xd = points.diff; 963 auto yd = values.diff; 964 auto dd = yd / xd; 965 auto dd2 = points.zip(values).slide!(3, "(c[1] - a[1]) / (c[0] - a[0])"); 966 967 with(SplineType) final switch(kind) 968 { 969 case c2: 970 break; 971 case cardinal: 972 if (lbc.type == SplineBoundaryType.notAKnot) 973 lbc.type = SplineBoundaryType.parabolic; 974 if (rbc.type == SplineBoundaryType.notAKnot) 975 rbc.type = SplineBoundaryType.parabolic; 976 break; 977 case monotone: 978 if (lbc.type == SplineBoundaryType.notAKnot) 979 lbc.type = SplineBoundaryType.monotone; 980 if (rbc.type == SplineBoundaryType.notAKnot) 981 rbc.type = SplineBoundaryType.monotone; 982 break; 983 case doubleQuadratic: 984 if (lbc.type == SplineBoundaryType.notAKnot) 985 lbc.type = SplineBoundaryType.parabolic; 986 if (rbc.type == SplineBoundaryType.notAKnot) 987 rbc.type = SplineBoundaryType.parabolic; 988 break; 989 case akima: 990 if (lbc.type == SplineBoundaryType.notAKnot) 991 lbc.type = SplineBoundaryType.akima; 992 if (rbc.type == SplineBoundaryType.notAKnot) 993 rbc.type = SplineBoundaryType.akima; 994 break; 995 } 996 997 if (n <= 3) 998 { 999 if (lbc.type == SplineBoundaryType.notAKnot) 1000 lbc.type = SplineBoundaryType.parabolic; 1001 if (rbc.type == SplineBoundaryType.notAKnot) 1002 rbc.type = SplineBoundaryType.parabolic; 1003 1004 if (n == 2) 1005 { 1006 if (lbc.type == SplineBoundaryType.monotone 1007 || lbc.type == SplineBoundaryType.akima) 1008 lbc.type = SplineBoundaryType.parabolic; 1009 if (rbc.type == SplineBoundaryType.monotone 1010 || rbc.type == SplineBoundaryType.akima) 1011 rbc.type = SplineBoundaryType.parabolic; 1012 } 1013 /// special case 1014 if (rbc.type == SplineBoundaryType.parabolic 1015 && lbc.type == SplineBoundaryType.parabolic) 1016 { 1017 import mir.interpolate.utility; 1018 if (n == 3) 1019 { 1020 auto derivatives = parabolaDerivatives(points[0], points[1], points[2], values[0], values[1], values[2]); 1021 slopes[0] = derivatives[0]; 1022 slopes[1] = derivatives[1]; 1023 slopes[2] = derivatives[2]; 1024 } 1025 else 1026 { 1027 assert(slopes.length == 2); 1028 slopes.back = slopes.front = yd.front / xd.front; 1029 } 1030 return; 1031 } 1032 } 1033 1034 with(SplineBoundaryType) final switch(lbc.type) 1035 { 1036 case periodic: 1037 1038 assert(0); 1039 1040 case notAKnot: 1041 1042 auto dx0 = xd[0]; 1043 auto dx1 = xd[1]; 1044 auto dy0 = yd[0]; 1045 auto dy1 = yd[1]; 1046 auto dd0 = dy0 / dx0; 1047 auto dd1 = dy1 / dx1; 1048 1049 slopes.front = dx1; 1050 first = dx0 + dx1; 1051 temp.front = ((dx0 + 2 * first) * dx1 * dd0 + dx0 ^^ 2 * dd1) / first; 1052 break; 1053 1054 case firstDerivative: 1055 1056 slopes.front = 1; 1057 first = 0; 1058 temp.front = lbc.value; 1059 break; 1060 1061 case secondDerivative: 1062 1063 slopes.front = 2; 1064 first = 1; 1065 temp.front = 3 * dd.front - 0.5 * lbc.value * xd.front; 1066 break; 1067 1068 case parabolic: 1069 1070 slopes.front = 1; 1071 first = 1; 1072 temp.front = 2 * dd.front; 1073 break; 1074 1075 case monotone: 1076 1077 slopes.front = 1; 1078 first = 0; 1079 temp.front = pchipTail(xd[0], xd[1], dd[0], dd[1]); 1080 break; 1081 1082 case akima: 1083 1084 slopes.front = 1; 1085 first = 0; 1086 temp.front = akimaTail(dd[0], dd[1]); 1087 break; 1088 1089 } 1090 1091 with(SplineBoundaryType) final switch(rbc.type) 1092 { 1093 case periodic: 1094 assert(0); 1095 1096 case notAKnot: 1097 1098 auto dx0 = xd[$ - 1]; 1099 auto dx1 = xd[$ - 2]; 1100 auto dy0 = yd[$ - 1]; 1101 auto dy1 = yd[$ - 2]; 1102 auto dd0 = dy0 / dx0; 1103 auto dd1 = dy1 / dx1; 1104 slopes.back = dx1; 1105 last = dx0 + dx1; 1106 temp.back = ((dx0 + 2 * last) * dx1 * dd0 + dx0 ^^ 2 * dd1) / last; 1107 break; 1108 1109 case firstDerivative: 1110 1111 slopes.back = 1; 1112 last = 0; 1113 temp.back = rbc.value; 1114 break; 1115 1116 case secondDerivative: 1117 1118 slopes.back = 2; 1119 last = 1; 1120 temp.back = 3 * dd.back + 0.5 * rbc.value * xd.back; 1121 break; 1122 1123 case parabolic: 1124 1125 slopes.back = 1; 1126 last = 1; 1127 temp.back = 2 * dd.back; 1128 break; 1129 1130 case monotone: 1131 1132 slopes.back = 1; 1133 last = 0; 1134 temp.back = pchipTail(xd[$ - 1], xd[$ - 2], dd[$ - 1], dd[$ - 2]); 1135 break; 1136 1137 case akima: 1138 1139 slopes.back = 1; 1140 last = 0; 1141 temp.back = akimaTail(dd[$ - 1], dd[$ - 2]); 1142 break; 1143 1144 } 1145 1146 with(SplineType) final switch(kind) 1147 { 1148 case c2: 1149 1150 foreach (i; 1 .. n - 1) 1151 { 1152 auto dx0 = xd[i - 1]; 1153 auto dx1 = xd[i - 0]; 1154 auto dy0 = yd[i - 1]; 1155 auto dy1 = yd[i - 0]; 1156 slopes[i] = 2 * (dx0 + dx1); 1157 temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0); 1158 } 1159 break; 1160 1161 case cardinal: 1162 1163 foreach (i; 1 .. n - 1) 1164 { 1165 slopes[i] = 1; 1166 temp[i] = (1 - param) * dd2[i - 1]; 1167 } 1168 break; 1169 1170 case monotone: 1171 { 1172 auto step0 = cast()xd[0]; 1173 auto step1 = cast()xd[1]; 1174 auto diff0 = cast()yd[0]; 1175 auto diff1 = cast()yd[1]; 1176 diff0 /= step0; 1177 diff1 /= step1; 1178 1179 for(size_t i = 1;;) 1180 { 1181 slopes[i] = 1; 1182 if (diff0 && diff1 && copysign(1f, diff0) == copysign(1f, diff1)) 1183 { 1184 auto w0 = step1 * 2 + step0; 1185 auto w1 = step0 * 2 + step1; 1186 temp[i] = (w0 + w1) / (w0 / diff0 + w1 / diff1); 1187 } 1188 else 1189 { 1190 temp[i] = 0; 1191 } 1192 if (++i == n - 1) 1193 { 1194 break; 1195 } 1196 step0 = step1; 1197 diff0 = diff1; 1198 step1 = xd[i]; 1199 diff1 = yd[i]; 1200 diff1 /= step1; 1201 } 1202 } 1203 break; 1204 1205 case doubleQuadratic: 1206 1207 foreach (i; 1 .. n - 1) 1208 { 1209 slopes[i] = 1; 1210 temp[i] = dd[i - 1] + dd[i] - dd2[i - 1]; 1211 } 1212 break; 1213 1214 case akima: 1215 { 1216 auto d3 = dd[1]; 1217 auto d2 = dd[0]; 1218 auto d1 = 2 * d2 - d3; 1219 auto d0 = d1; 1220 foreach (i; 1 .. n - 1) 1221 { 1222 d0 = d1; 1223 d1 = d2; 1224 d2 = d3; 1225 d3 = i == n - 2 ? 2 * d2 - d1 : dd[i + 1]; 1226 slopes[i] = 1; 1227 temp[i] = akimaSlope(d0, d1, d2, d3); 1228 } 1229 break; 1230 } 1231 } 1232 1233 foreach (i; 0 .. n - 1) 1234 { 1235 auto c = i == 0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0; 1236 auto a = i == n - 2 ? last : kind == SplineType.c2 ? xd[i + 1] : 0; 1237 auto w = slopes[i] == 1 ? a : a / slopes[i]; 1238 slopes[i + 1] -= w * c; 1239 temp[i + 1] -= w * temp[i]; 1240 } 1241 1242 slopes.back = temp.back / slopes.back; 1243 1244 foreach_reverse (i; 0 .. n - 1) 1245 { 1246 auto c = i == 0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0; 1247 auto v = temp[i] - c * slopes[i + 1]; 1248 slopes[i] = slopes[i] == 1 ? v : v / slopes[i]; 1249 } 1250 } 1251 1252 private F akimaTail(F)(in F d2, in F d3) 1253 { 1254 auto d1 = 2 * d2 - d3; 1255 auto d0 = 2 * d1 - d2; 1256 return akimaSlope(d0, d1, d2, d3); 1257 } 1258 1259 private F akimaSlope(F)(in F d0, in F d1, in F d2, in F d3) 1260 { 1261 if (d1 == d2) 1262 return d1; 1263 if (d0 == d1 && d2 == d3) 1264 return (d1 + d2) * 0.5f; 1265 if (d0 == d1) 1266 return d1; 1267 if (d2 == d3) 1268 return d2; 1269 auto w0 = fabs(d1 - d0); 1270 auto w1 = fabs(d3 - d2); 1271 auto ws = w0 + w1; 1272 w0 /= ws; 1273 w1 /= ws; 1274 return w0 * d2 + w1 * d1; 1275 } 1276 1277 /// 1278 struct SplineKernel(X) 1279 { 1280 X step = 0; 1281 X w0 = 0; 1282 X w1 = 0; 1283 X wq = 0; 1284 1285 /// 1286 this(X x0, X x1, X x) 1287 { 1288 step = x1 - x0; 1289 auto c0 = x - x0; 1290 auto c1 = x1 - x; 1291 w0 = c0 / step; 1292 w1 = c1 / step; 1293 wq = w0 * w1; 1294 } 1295 1296 /// 1297 template opCall(uint derivative = 0) 1298 if (derivative <= 3) 1299 { 1300 /// 1301 auto opCall(Y)(in Y y0, in Y y1, in Y s0, in Y s1) const 1302 { 1303 auto diff = y1 - y0; 1304 auto z0 = s0 * step - diff; 1305 auto z1 = s1 * step - diff; 1306 auto a0 = z0 * w1; 1307 auto a1 = z1 * w0; 1308 auto pr = a0 - a1; 1309 auto b0 = y0 * w1; 1310 auto b1 = y1 * w0; 1311 auto pl = b0 + b1; 1312 auto y = pl + wq * pr; 1313 static if (derivative) 1314 { 1315 Y[derivative + 1] ret = 0; 1316 ret[0] = y; 1317 auto wd = w1 - w0; 1318 auto zd = z1 + z0; 1319 ret[1] = (diff + (wd * pr - wq * zd)) / step; 1320 static if (derivative > 1) 1321 { 1322 auto astep = zd / (step * step); 1323 ret[2] = -3 * wd * astep + (s1 - s0) / step; 1324 static if (derivative > 2) 1325 ret[3] = 6 * astep / step; 1326 } 1327 return ret; 1328 } 1329 else 1330 { 1331 return y; 1332 } 1333 } 1334 } 1335 1336 /// 1337 alias withDerivative = opCall!1; 1338 /// 1339 alias withTwoDerivatives = opCall!2; 1340 } 1341 1342 package T pchipTail(T)(in T step0, in T step1, in T diff0, in T diff1) 1343 { 1344 import mir.math.common: copysign, fabs; 1345 if (!diff0) 1346 { 1347 return 0; 1348 } 1349 auto slope = ((step0 * 2 + step1) * diff0 - step0 * diff1) / (step0 + step1); 1350 if (copysign(1f, slope) != copysign(1f, diff0)) 1351 { 1352 return 0; 1353 } 1354 if ((copysign(1f, diff0) != copysign(1f, diff1)) && (fabs(slope) > fabs(diff0 * 3))) 1355 { 1356 return diff0 * 3; 1357 } 1358 return slope; 1359 }