1 /++ 2 $(H1 Interpolation Algorithms) 3 4 $(BOOKTABLE $(H2 Interpolation modules), 5 $(TR $(TH Module) $(TH Interpolation kind)) 6 $(T2M constant, Constant Interpolant) 7 $(T2M generic, Generic Piecewise Interpolant) 8 $(T2M linear, Linear Interpolant) 9 $(T2M polynomial, Lagrange Barycentric Interpolant) 10 $(T2M spline, Piecewise Cubic Hermite Interpolant Spline: C2 (with contiguous second derivative), cardinal, monotone (aka PCHIP), double-quadratic, Akima) 11 ) 12 ] 13 Copyright: 2020 Ilya Yaroshenko, Kaleidic Associates Advisory Limited, Symmetry Investments 14 Authors: Ilya Yaroshenko 15 16 Macros: 17 SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir,interpolate, $1)$(NBSP) 18 T2M=$(TR $(TDNW $(MREF mir,interpolate,$1)) $(TD $+)) 19 T2=$(TR $(TDNW $(LREF $1)) $(TD $+)) 20 +/ 21 module mir.interpolate; 22 23 import mir.functional: RefTuple; 24 import mir.math.common: optmath; 25 import mir.ndslice.slice: Slice, Contiguous; 26 import mir.primitives; 27 import mir.qualifier; 28 import std.traits: isInstanceOf; 29 30 @optmath: 31 32 package ref iter(alias s)() { return s._iterator; }; 33 package alias GridVector(It) = Slice!It; 34 35 package enum bool isInterval(T) = isInstanceOf!(RefTuple, T); 36 package enum bool isInterval(alias T) = isInstanceOf!(RefTuple, T); 37 38 package template Repeat(ulong n, L...) 39 { 40 static if (n) 41 { 42 static import std.meta; 43 alias Repeat = std.meta.Repeat!(n, L); 44 } 45 else 46 alias Repeat = L[0 .. 0]; 47 } 48 49 /++ 50 Interval index for x value given. 51 +/ 52 template findInterval(size_t dimension = 0) 53 { 54 /++ 55 Interval index for x value given. 56 Params: 57 interpolant = interpolant 58 x = X value 59 +/ 60 size_t findInterval(Interpolant, X)(auto ref const Interpolant interpolant, in X x) @trusted 61 { 62 import mir.ndslice.slice: sliced; 63 import mir.ndslice.sorting: transitionIndex; 64 static if (dimension) 65 { 66 immutable sizediff_t len = interpolant.intervalCount!dimension - 1; 67 auto grid = interpolant.gridScopeView!dimension[1 .. $][0 .. len]; 68 } 69 else 70 { 71 immutable sizediff_t len = interpolant.intervalCount - 1; 72 auto grid = interpolant.gridScopeView[1 .. $][0 .. len]; 73 } 74 assert(len >= 0); 75 return grid.transitionIndex!"a <= b"(x); 76 } 77 } 78 79 /// 80 version(mir_test) unittest 81 { 82 import mir.ndslice.allocation: rcslice; 83 import mir.ndslice.topology: as; 84 import mir.ndslice.slice: sliced; 85 import mir.interpolate.linear; 86 87 static immutable x = [0.0, 1, 2]; 88 static immutable y = [10.0, 2, 4]; 89 auto interpolation = linear!double(x.rcslice, y.as!(const double).rcslice); 90 assert(interpolation.findInterval(1.0) == 1); 91 } 92 93 /++ 94 Lazy interpolation shell with linear complexity. 95 96 Params: 97 range = sorted range 98 interpolant = interpolant structure with `.gridScopeView` method. 99 Complexity: 100 `O(range.length + interpolant.gridScopeView.length)` to evaluate all elements. 101 Returns: 102 Lazy input range. 103 See_also: 104 $(SUBREF linear, linear), 105 $(SUBREF spline, spline). 106 +/ 107 auto interp1(Range, Interpolant)(Range range, Interpolant interpolant, size_t interval = 0) 108 { 109 return Interp1!(Range, Interpolant)(range, interpolant, interval); 110 } 111 112 /// ditto 113 struct Interp1(Range, Interpolant) 114 { 115 /// Sorted range (descending). $(RED For internal use.) 116 private Range _range; 117 /// Interpolant structure. $(RED For internal use.) 118 private Interpolant _interpolant; 119 /// Current interpolation interval. $(RED For internal use.) 120 private size_t _interval; 121 122 /// 123 auto lightScope() 124 { 125 return Interp1!(LightScopeOf!Range, LightScopeOf!Interpolant)(.lightScope(_range), .lightScope(_interpolant), _interval); 126 } 127 128 static if (hasLength!Range) 129 /// Length (optional) 130 size_t length()() const @property { return _range.length; } 131 /// Save primitive (optional) 132 auto save()() @property 133 { 134 auto ret = this; 135 ret._range = _range.save; 136 return ret; 137 } 138 /// Input range primitives 139 bool empty () const @property { return _range.empty; } 140 /// ditto 141 void popFront() { _range.popFront; } 142 /// ditto 143 auto front() @property 144 145 { 146 assert(!empty); 147 auto x = _range.front; 148 return (x) @trusted { 149 auto points = _interpolant.gridScopeView; 150 sizediff_t len = _interpolant.intervalCount - 1; 151 assert(len >= 0); 152 while (x > points[_interval + 1] && _interval < len) 153 _interval++; 154 return _interpolant(x.atInterval(_interval)); 155 } (x); 156 } 157 } 158 159 /++ 160 PCHIP interpolation. 161 +/ 162 version(mir_test) 163 @safe unittest 164 { 165 import mir.math.common: approxEqual; 166 import mir.ndslice.slice: sliced; 167 import mir.ndslice.allocation: rcslice; 168 import mir.interpolate: interp1; 169 import mir.interpolate.spline; 170 171 static immutable x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22]; 172 static immutable y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6]; 173 auto interpolation = spline!double(x.rcslice, y.sliced, SplineType.monotone); 174 175 auto xs = x[0 .. $ - 1].sliced + 0.5; 176 177 auto ys = xs.interp1(interpolation); 178 } 179 180 @safe pure @nogc version(mir_test) unittest 181 { 182 import mir.interpolate.linear; 183 import mir.ndslice; 184 import mir.math.common: approxEqual; 185 186 static immutable x = [0, 1, 2, 3, 5.00274, 7.00274, 10.0055, 20.0137, 30.0192]; 187 static immutable y = [0.0011, 0.0011, 0.0030, 0.0064, 0.0144, 0.0207, 0.0261, 0.0329, 0.0356,]; 188 static immutable 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]; 189 190 auto interpolation = linear!double(x.rcslice, y.as!(const double).rcslice); 191 192 static immutable data = [0.0011, 0.0030, 0.0064, 0.0104, 0.0144, 0.0176, 0.0207, 0.0225, 0.0243, 0.0261, 0.0268, 0.0274, 0.0281, 0.0288, 0.0295, 0.0302, 0.0309, 0.0316, 0.0322, 0.0329, 0.0332, 0.0335, 0.0337, 0.0340, 0.0342, 0.0345, 0.0348, 0.0350, 0.0353, 0.0356]; 193 194 () @trusted { 195 assert(all!((a, b) => approxEqual(a, b, 1e-4, 1e-4))(xs.interp1(interpolation), data)); 196 }(); 197 } 198 199 /++ 200 Optimization utility that can be used with interpolants if 201 x should be extrapolated at interval given. 202 203 By default interpolants uses binary search to find appropriate interval, 204 it has `O(log(.gridScopeView.length))` complexity. 205 If an interval is given, interpolant uses it instead of binary search. 206 +/ 207 RefTuple!(T, size_t) atInterval(T)(in T value, size_t intervalIndex) 208 { 209 return typeof(return)(value, intervalIndex); 210 } 211 212 /// 213 version(mir_test) unittest 214 { 215 import mir.ndslice.allocation; 216 import mir.ndslice.slice; 217 import mir.interpolate.spline; 218 static immutable x = [0.0, 1, 2]; 219 static immutable y = [3.0, 4, -10]; 220 auto interpolant = spline!double(x.rcslice, y.sliced); 221 assert(interpolant(1.3) != interpolant(1.3.atInterval(0))); 222 assert(interpolant(1.3) == interpolant(1.3.atInterval(1))); 223 } 224 225 version(D_AVX2) 226 enum _avx = true; 227 else 228 version(D_AVX) 229 enum _avx = true; 230 else 231 enum _avx = false; 232 233 version(LDC) 234 enum LDC = true; 235 else 236 enum LDC = false; 237 238 version(X86_64) 239 enum x86_64 = true; 240 else 241 enum x86_64 = false; 242 243 auto copyvec(F, size_t N)(ref const F[N] from, ref F[N] to) 244 { 245 import mir.internal.utility; 246 247 static if (LDC && F.mant_dig != 64 && is(__vector(F[N]))) 248 { 249 alias V = __vector(F[N]); // @FUTURE@ vector support 250 *cast(V*) to.ptr = *cast(V*) from.ptr; 251 } 252 else 253 static if (F.sizeof <= double.sizeof && F[N].sizeof >= (double[2]).sizeof && x86_64) 254 { 255 import mir.utility; 256 enum S = _avx ? 32u : 16u; 257 enum M = min(S, F[N].sizeof) / F.sizeof; 258 alias V = __vector(F[M]); // @FUTURE@ vector support 259 enum C = N / M; 260 foreach(i; Iota!C) 261 { 262 *cast(V*)(to.ptr + i * M) = *cast(V*)(from.ptr + i * M); 263 } 264 } 265 else 266 { 267 to = from; 268 } 269 } 270 271 package template SplineReturnType(F, size_t N, size_t P) 272 { 273 static if (P <= 1 || N == 0) 274 alias SplineReturnType = F; 275 else 276 alias SplineReturnType = .SplineReturnType!(F, N - 1, P)[P]; 277 } 278 279 template generateShuffles3(size_t N, size_t P) 280 { 281 enum size_t[N][2] generateShuffles3 = 282 () 283 { 284 auto ret = new size_t[](2 * N); 285 size_t[2] j; 286 bool f = 1; 287 foreach(i; 0 .. N) 288 { 289 ret[i * 2] = i; 290 ret[i * 2 + 1] = i + N; 291 } 292 return [ret[0 .. $ / 2], ret[$ / 2 .. $]]; 293 }(); 294 } 295 296 297 void shuffle3(size_t P, F, size_t N)(ref F[N] a, ref F[N] b, ref F[N] c, ref F[N] d) 298 if (P <= N && N) 299 { 300 static if (P == 0 || N == 1) 301 { 302 copyvec(a, c); 303 copyvec(b, d); 304 } 305 else 306 version(LDC) 307 { 308 enum masks = generateShuffles3!(N, P); 309 import std.meta: aliasSeqOf; 310 import ldc.simd: shufflevector; 311 alias V = __vector(F[N]); // @FUTURE@ vector support 312 auto as = shufflevector!(V, aliasSeqOf!(masks[0]))(*cast(V*)a.ptr, *cast(V*)b.ptr); 313 auto bs = shufflevector!(V, aliasSeqOf!(masks[1]))(*cast(V*)a.ptr, *cast(V*)b.ptr); 314 *cast(V*)c.ptr = as; 315 *cast(V*)d.ptr = bs; 316 } 317 else 318 { 319 foreach(i; Iota!(N / P)) 320 { 321 enum j = 2 * i * P; 322 static if (j < N) 323 { 324 copyvec(a[i * P .. i * P + P], c[j .. j + P]); 325 static assert(j + 2 * P <= c.length); 326 copyvec(b[i * P .. i * P + P], c[j + P .. j + 2 * P]); 327 } 328 else 329 { 330 copyvec(a[i * P .. i * P + P], d[j - N .. j - N + P]); 331 copyvec(b[i * P .. i * P + P], d[j - N + P .. j - N + 2 * P]); 332 } 333 } 334 } 335 } 336 337 void shuffle2(size_t P, F, size_t N)(ref F[N] a, ref F[N] b, ref F[N] c, ref F[N] d) 338 if (P <= N && N) 339 { 340 static if (P == 0 || N == 1) 341 { 342 copyvec(a, c); 343 copyvec(b, d); 344 } 345 else 346 version(LDC) 347 { 348 enum masks = generateShuffles2!(N, P); 349 import std.meta: aliasSeqOf; 350 import ldc.simd: shufflevector; 351 alias V = __vector(F[N]); // @FUTURE@ vector support 352 auto as = shufflevector!(V, aliasSeqOf!(masks[0]))(*cast(V*)a.ptr, *cast(V*)b.ptr); 353 auto bs = shufflevector!(V, aliasSeqOf!(masks[1]))(*cast(V*)a.ptr, *cast(V*)b.ptr); 354 *cast(V*)c.ptr = as; 355 *cast(V*)d.ptr = bs; 356 } 357 else 358 { 359 foreach(i; Iota!(N / P)) 360 { 361 enum j = 2 * i * P; 362 static if (j < N) 363 copyvec(a[j .. j + P], c[i * P .. i * P + P]); 364 else 365 copyvec(b[j - N .. j - N + P], c[i * P .. i * P + P]); 366 } 367 foreach(i; Iota!(N / P)) 368 { 369 enum j = 2 * i * P + P; 370 static if (j < N) 371 copyvec(a[j .. j + P], d[i * P .. i * P + P]); 372 else 373 copyvec(b[j - N .. j - N + P], d[i * P .. i * P + P]); 374 } 375 } 376 } 377 378 void shuffle1(size_t P, F, size_t N)(ref F[N] a, ref F[N] b, ref F[N] c, ref F[N] d) 379 if (P <= N && N) 380 { 381 static if (P == 0 || N == 1) 382 { 383 copyvec(a, c); 384 copyvec(b, d); 385 } 386 else 387 version(LDC) 388 { 389 enum masks = generateShuffles1!(N, P); 390 import std.meta: aliasSeqOf; 391 import ldc.simd: shufflevector; 392 alias V = __vector(F[N]); // @FUTURE@ vector support 393 auto as = shufflevector!(V, aliasSeqOf!(masks[0]))(*cast(V*)a.ptr, *cast(V*)b.ptr); 394 auto bs = shufflevector!(V, aliasSeqOf!(masks[1]))(*cast(V*)a.ptr, *cast(V*)b.ptr); 395 *cast(V*)c.ptr = as; 396 *cast(V*)d.ptr = bs; 397 } 398 else 399 { 400 foreach(i; Iota!(N / P)) 401 { 402 enum j = i * P; 403 static if (i % 2 == 0) 404 copyvec(a[j .. j + P], c[i * P .. i * P + P]); 405 else 406 copyvec(b[j - P .. j], c[i * P .. i * P + P]); 407 } 408 foreach(i; Iota!(N / P)) 409 { 410 enum j = i * P + P; 411 static if (i % 2 == 0) 412 copyvec(a[j .. j + P], d[i * P .. i * P + P]); 413 else 414 copyvec(b[j - P .. j], d[i * P .. i * P + P]); 415 } 416 } 417 } 418 419 template generateShuffles2(size_t N, size_t P) 420 { 421 enum size_t[N][2] generateShuffles2 = 422 () 423 { 424 auto ret = new size_t[][](2, N); 425 size_t[2] j; 426 bool f = 1; 427 foreach(i; 0 .. 2 * N) 428 { 429 if (i % P == 0) 430 f = !f; 431 ret[f][j[f]++] = i; 432 } 433 return ret; 434 }(); 435 } 436 437 438 template generateShuffles1(size_t N, size_t P) 439 { 440 enum size_t[N][2] generateShuffles1 = 441 () 442 { 443 auto ret = new size_t[][](2, N); 444 foreach(i; 0 .. N) 445 { 446 ret[0][i] = (i / P + 1) % 2 ? i : i + N - P; 447 ret[1][i] = ret[0][i] + P; 448 } 449 return ret; 450 }(); 451 } 452 453 unittest 454 { 455 assert(generateShuffles1!(4, 1) == [[0, 4, 2, 6], [1, 5, 3, 7]]); 456 assert(generateShuffles1!(4, 2) == [[0, 1, 4, 5], [2, 3, 6, 7]]); 457 assert(generateShuffles1!(4, 4) == [[0, 1, 2, 3], [4, 5, 6, 7]]); 458 } 459 460 unittest 461 { 462 assert(generateShuffles2!(4, 1) == [[0, 2, 4, 6], [1, 3, 5, 7]]); 463 assert(generateShuffles2!(4, 2) == [[0, 1, 4, 5], [2, 3, 6, 7]]); 464 assert(generateShuffles2!(4, 4) == [[0, 1, 2, 3], [4, 5, 6, 7]]); 465 } 466 467 unittest 468 { 469 enum ai = [0, 1, 2, 3]; 470 enum bi = [4, 5, 6, 7]; 471 align(32) 472 double[4] a = [0, 1, 2, 3], b = [4, 5, 6, 7], c, d; 473 shuffle3!1(a, b, c, d); 474 assert([c, d] == [[0.0, 4, 1, 5], [2.0, 6, 3, 7]]); 475 } 476 477 unittest 478 { 479 enum ai = [0, 1, 2, 3]; 480 enum bi = [4, 5, 6, 7]; 481 align(32) 482 double[4] a = [0, 1, 2, 3], b = [4, 5, 6, 7], c, d; 483 shuffle2!1(a, b, c, d); 484 assert([c, d] == [[0.0, 2, 4, 6], [1.0, 3, 5, 7]]); 485 shuffle2!2(a, b, c, d); 486 assert([c, d] == [[0.0, 1, 4, 5], [2.0, 3, 6, 7]]); 487 // shuffle2!4(a, b, c, d); 488 // assert([c, d] == [[0.0, 1, 2, 3], [4.0, 5, 6, 7]]); 489 } 490 491 unittest 492 { 493 enum ai = [0, 1, 2, 3]; 494 enum bi = [4, 5, 6, 7]; 495 align(32) 496 double[4] a = [0, 1, 2, 3], b = [4, 5, 6, 7], c, d; 497 shuffle1!1(a, b, c, d); 498 assert([c, d] == [[0, 4, 2, 6], [1, 5, 3, 7]]); 499 shuffle1!2(a, b, c, d); 500 assert([c, d] == [[0, 1, 4, 5], [2, 3, 6, 7]]); 501 // shuffle1!4(a, b, c, d); 502 // assert([c, d] == [[0, 1, 2, 3], [4, 5, 6, 7]]); 503 } 504 505 import mir.internal.utility; 506 507 auto vectorize(Kernel, F, size_t N, size_t R)(ref Kernel kernel, ref F[N] a0, ref F[N] b0, ref F[N] a1, ref F[N] b1, ref F[N][R] c) 508 { 509 // static if (LDC && F.mant_dig != 64) 510 // { 511 // alias V = __vector(F[N]); // @FUTURE@ vector support 512 // *cast(V[R]*) c.ptr = kernel( 513 // *cast(V*)a0.ptr, *cast(V*)b0.ptr, 514 // *cast(V*)a1.ptr, *cast(V*)b1.ptr); 515 // } 516 // else 517 // static if (F.sizeof <= double.sizeof && F[N].sizeof >= (double[2]).sizeof) 518 // { 519 // import mir.utility; 520 // enum S = _avx ? 32u : 16u; 521 // enum M = min(S, F[N].sizeof) / F.sizeof; 522 // alias V = __vector(F[M]); // @FUTURE@ vector support 523 // enum C = N / M; 524 // foreach(i; Iota!C) 525 // { 526 // auto r = kernel( 527 // *cast(V*)(a0.ptr + i * M), *cast(V*)(b0.ptr + i * M), 528 // *cast(V*)(a1.ptr + i * M), *cast(V*)(b1.ptr + i * M), 529 // ); 530 // static if (R == 1) 531 // *cast(V*)(c[0].ptr + i * M) = r; 532 // else 533 // foreach(j; Iota!R) 534 // *cast(V*)(c[j].ptr + i * M) = r[j]; 535 // } 536 // } 537 // else 538 // { 539 foreach(i; Iota!N) 540 { 541 auto r = kernel(a0[i], b0[i], a1[i], b1[i]); 542 static if (R == 1) 543 c[0][i] = r; 544 else 545 foreach(j; Iota!R) 546 c[j][i] = r[j]; 547 } 548 // } 549 } 550 551 auto vectorize(Kernel, F, size_t N, size_t R)(ref Kernel kernel, ref F[N] a, ref F[N] b, ref F[N][R] c) 552 { 553 // static if (LDC && F.mant_dig != 64 && is(__vector(F[N]))) 554 // { 555 // alias V = __vector(F[N]); // @FUTURE@ vector support 556 // *cast(V[R]*) c.ptr = kernel(*cast(V*)a.ptr, *cast(V*)b.ptr); 557 // } 558 // else 559 // static if (F.sizeof <= double.sizeof && F[N].sizeof >= (double[2]).sizeof && x86_64) 560 // { 561 // import mir.utility; 562 // enum S = _avx ? 32u : 16u; 563 // enum M = min(S, F[N].sizeof) / F.sizeof; 564 // alias V = __vector(F[M]); // @FUTURE@ vector support 565 // enum C = N / M; 566 // foreach(i; Iota!C) 567 // { 568 // auto r = kernel( 569 // *cast(V*)(a.ptr + i * M), 570 // *cast(V*)(b.ptr + i * M), 571 // ); 572 // static if (R == 1) 573 // *cast(V*)(c[0].ptr + i * M) = r; 574 // else 575 // foreach(j; Iota!R) 576 // *cast(V*)(c[j].ptr + i * M) = r[j]; 577 // } 578 // } 579 // else 580 // { 581 F[N][R] _c;//Temporary array in case "c" overlaps "a" and/or "b". 582 foreach(i; Iota!N) 583 { 584 auto r = kernel(a[i], b[i]); 585 static if (R == 1) 586 _c[0][i] = r; 587 else 588 foreach(j; Iota!R) 589 _c[j][i] = r[j]; 590 } 591 c = _c; 592 // } 593 } 594 595 // version(unittest) 596 // template _test_fun(size_t R) 597 // { 598 // auto _test_fun(T)(T a0, T b0, T a1, T b1) 599 // { 600 // static if (R == 4) 601 // { 602 // return cast(T[R])[a0, b0, a1, b1]; 603 // } 604 // else 605 // static if (R == 2) 606 // { 607 // return cast(T[R])[a0 + b0, a1 + b1]; 608 // } 609 // else 610 // return a0 + b0 + a1 + b1; 611 // } 612 // } 613 614 // unittest 615 // { 616 // import std.meta; 617 618 // foreach(F; AliasSeq!(float, double)) 619 // foreach(N; AliasSeq!(1, 2, 4, 8, 16)) 620 // { 621 // align(F[N].sizeof) F[N] a0 = 4; 622 // align(F[N].sizeof) F[N] b0 = 30; 623 // align(F[N].sizeof) F[N] a1 = 200; 624 // align(F[N].sizeof) F[N] b1 = 1000; 625 // align(F[N].sizeof) F[N][1] c1; 626 // align(F[N].sizeof) F[N][2] c2; 627 // align(F[N].sizeof) F[N][4] c4; 628 // vectorize!(_test_fun!(1))(a0, b0, a1, b1, c1); 629 // vectorize!(_test_fun!(2))(a0, b0, a1, b1, c2); 630 // vectorize!(_test_fun!(4))(a0, b0, a1, b1, c4); 631 // } 632 // }