In [8]:
'''Based on the paper doi=10.1.1.89.7835
the tricubic interpolation of a regular possibly non uniform grid can be seen as a computation of 21 cubic splines.
A cubic spline is a special case of cubic interpolation, and in general these 21 cubic splines perform many redundant
calulations. Here we formulate the full tricubic interpolation where the value of a scalar function defined on a 3d
grid can be reconstructed to allow full C1, and thus langrangian structures to persist.'''
import numpy as np
import dill
dill.settings['recurse'] = True
from numba import jit
@jit(nogil=True,nopython=True)
def bisection(array,value):
'''Given an ``array`` , and given a ``value`` , returns an index j such that ``value`` is between array[j]
and array[j+1]. ``array`` must be monotonic increasing. j=-1 or j=len(array) is returned
to indicate that ``value`` is out of range below and above respectively.'''
n = len(array)
if (value < array[0]):
return -1
res = -1# Then set the output
elif (value > array[n-1]):
return n
#array = np.append(np.append(-np.inf,array),np.inf)
jl = 0# Initialize lower
ju = n-1# and upper limits.
while (ju-jl > 1):# If we are not yet done,
jm=(ju+jl) >> 1# compute a midpoint,
if (value >= array[jm]):
jl=jm# and replace either the lower limit
else:
ju=jm# or the upper limit, as appropriate.
# Repeat until the test condition is satisfied.
if (value == array[0]):
return 0
res = -1# Then set the output
elif (value == array[n-1]):
return n-1
else:
return jl
@jit(nogil=True,nopython=True)
def get_bVec(i,j,k,nx,ny,nz,xvec,yvec,zvec,m):
nx = nx
ny = ny
nz = nz
x0 = ny*i
x1 = x0 + j
x2 = nz*x1
x3 = x2 + k
x4 = m[x3]
x5 = i - 1.0
x6 = xvec[x5]
x7 = i + 1.0
x8 = xvec[x7]
x9 = 1/(x6 - x8)
x10 = 0.166666666666667*x9
x11 = ny*(i - 2.0)
x12 = x11 + j
x13 = nz*x12
x14 = x13 + k
x15 = m[x14]
x16 = ny*x5
x17 = x16 + j
x18 = nz*x17
x19 = x18 + k
x20 = m[x19]
x21 = -8.0*x20
x22 = ny*x7
x23 = x22 + j
x24 = nz*x23
x25 = x24 + k
x26 = m[x25]
x27 = 8.0*x26
x28 = i + 2.0
x29 = ny*x28
x30 = x29 + j
x31 = nz*x30
x32 = x31 + k
x33 = m[x32]
x34 = j - 1.0
x35 = j + 1.0
x36 = 1/(yvec[x34] - yvec[x35])
x37 = 0.166666666666667*x36
x38 = nz*(x1 - 2.0)
x39 = x38 + k
x40 = m[x39]
x41 = nz*(x0 + x34)
x42 = x41 + k
x43 = m[x42]
x44 = -8.0*x43
x45 = nz*(x0 + x35)
x46 = x45 + k
x47 = m[x46]
x48 = 8.0*x47
x49 = j + 2.0
x50 = nz*(x0 + x49)
x51 = x50 + k
x52 = m[x51]
x53 = k - 1.0
x54 = k + 1.0
x55 = 1/(zvec[x53] - zvec[x54])
x56 = 0.166666666666667*x55
x57 = m[x3 - 2]
x58 = m[x2 + x53]
x59 = -8.0*x58
x60 = m[x2 + x54]
x61 = 8.0*x60
x62 = k + 2.0
x63 = m[x2 + x62]
x64 = 0.0277777777777778*x36*x9
x65 = nz*(x11 + x34)
x66 = x65 + k
x67 = m[x66]
x68 = -8.0*x67
x69 = nz*(x16 + x34)
x70 = x69 + k
x71 = m[x70]
x72 = 64.0*x71
x73 = nz*(x22 + x34)
x74 = x73 + k
x75 = m[x74]
x76 = 64.0*x75
x77 = -x76
x78 = nz*(x12 - 2.0)
x79 = x78 + k
x80 = m[x79]
x81 = nz*(x29 + x49)
x82 = x81 + k
x83 = m[x82]
x84 = nz*(x11 + x35)
x85 = x84 + k
x86 = m[x85]
x87 = 8.0*x86
x88 = nz*(x11 + x49)
x89 = x88 + k
x90 = m[x89]
x91 = nz*(x17 - 2.0)
x92 = x91 + k
x93 = m[x92]
x94 = 8.0*x93
x95 = nz*(x16 + x35)
x96 = x95 + k
x97 = m[x96]
x98 = 64.0*x97
x99 = nz*(x16 + x49)
x100 = x99 + k
x101 = m[x100]
x102 = 8.0*x101
x103 = nz*(x23 - 2.0)
x104 = x103 + k
x105 = m[x104]
x106 = 8.0*x105
x107 = nz*(x22 + x35)
x108 = x107 + k
x109 = m[x108]
x110 = 64.0*x109
x111 = nz*(x22 + x49)
x112 = x111 + k
x113 = m[x112]
x114 = 8.0*x113
x115 = -x114
x116 = nz*(x30 - 2.0)
x117 = x116 + k
x118 = m[x117]
x119 = nz*(x29 + x34)
x120 = x119 + k
x121 = m[x120]
x122 = 8.0*x121
x123 = nz*(x29 + x35)
x124 = x123 + k
x125 = m[x124]
x126 = 8.0*x125
x127 = -x126
x128 = 0.0277777777777778*x55*x9
x129 = m[x19 - 2]
x130 = -8.0*x129
x131 = m[x18 + x53]
x132 = 64.0*x131
x133 = m[x18 + x54]
x134 = 64.0*x133
x135 = -x134
x136 = m[x14 - 2]
x137 = m[x31 + x62]
x138 = m[x13 + x53]
x139 = -8.0*x138
x140 = m[x13 + x54]
x141 = 8.0*x140
x142 = m[x13 + x62]
x143 = m[x18 + x62]
x144 = 8.0*x143
x145 = m[x25 - 2]
x146 = 8.0*x145
x147 = m[x24 + x53]
x148 = -64.0*x147
x149 = m[x24 + x54]
x150 = 64.0*x149
x151 = m[x24 + x62]
x152 = -8.0*x151
x153 = m[x32 - 2]
x154 = m[x31 + x53]
x155 = 8.0*x154
x156 = m[x31 + x54]
x157 = 8.0*x156
x158 = -x157
x159 = 0.0277777777777778*x36*x55
x160 = m[x42 - 2]
x161 = -8.0*x160
x162 = m[x41 + x53]
x163 = 64.0*x162
x164 = m[x39 - 2]
x165 = m[x50 + x62]
x166 = m[x38 + x53]
x167 = 8.0*x166
x168 = m[x38 + x54]
x169 = 8.0*x168
x170 = m[x38 + x62]
x171 = m[x41 + x54]
x172 = -64.0*x171
x173 = m[x41 + x62]
x174 = 8.0*x173
x175 = m[x46 - 2]
x176 = 8.0*x175
x177 = m[x45 + x53]
x178 = 64.0*x177
x179 = -x178
x180 = m[x45 + x54]
x181 = 64.0*x180
x182 = m[x45 + x62]
x183 = -8.0*x182
x184 = m[x51 - 2]
x185 = m[x50 + x53]
x186 = 8.0*x185
x187 = m[x50 + x54]
x188 = -8.0*x187
x189 = x55*x9
x190 = m[x62 + x91]
x191 = m[x62 + x69]
x192 = 64.0*x191
x193 = m[x62 + x99]
x194 = 8.0*x193
x195 = m[x116 + x53]
x196 = m[x119 + x53]
x197 = 64.0*x196
x198 = m[x123 + x53]
x199 = 64.0*x198
x200 = m[x53 + x81]
x201 = 8.0*x200
x202 = -8.0*x190 + x192 + x194 - 8.0*x195 + x197 - x199 + x201
x203 = m[x62 + x78]
x204 = m[x89 - 2]
x205 = m[x117 - 2]
x206 = m[x62 + x81]
x207 = m[x53 + x78]
x208 = m[x54 + x78]
x209 = m[x66 - 2]
x210 = m[x53 + x65]
x211 = m[x54 + x65]
x212 = m[x62 + x65]
x213 = 8.0*x212
x214 = m[x85 - 2]
x215 = m[x53 + x84]
x216 = m[x54 + x84]
x217 = m[x62 + x84]
x218 = 8.0*x217
x219 = m[x53 + x88]
x220 = 8.0*x219
x221 = m[x54 + x88]
x222 = 8.0*x221
x223 = m[x62 + x88]
x224 = m[x92 - 2]
x225 = m[x53 + x91]
x226 = m[x54 + x91]
x227 = m[x70 - 2]
x228 = m[x53 + x69]
x229 = m[x54 + x69]
x230 = m[x96 - 2]
x231 = m[x53 + x95]
x232 = m[x54 + x95]
x233 = m[x62 + x95]
x234 = -64.0*x233
x235 = m[x100 - 2]
x236 = 8.0*x235
x237 = m[x53 + x99]
x238 = 64.0*x237
x239 = m[x54 + x99]
x240 = 64.0*x239
x241 = -x240
x242 = m[x104 - 2]
x243 = m[x103 + x53]
x244 = m[x103 + x54]
x245 = m[x103 + x62]
x246 = 8.0*x245
x247 = m[x74 - 2]
x248 = m[x53 + x73]
x249 = m[x54 + x73]
x250 = m[x62 + x73]
x251 = 64.0*x250
x252 = -x251
x253 = m[x108 - 2]
x254 = m[x107 + x53]
x255 = m[x107 + x54]
x256 = m[x107 + x62]
x257 = 64.0*x256
x258 = m[x112 - 2]
x259 = 8.0*x258
x260 = m[x111 + x53]
x261 = 64.0*x260
x262 = m[x111 + x54]
x263 = 64.0*x262
x264 = m[x111 + x62]
x265 = 8.0*x264
x266 = -x265
x267 = m[x116 + x54]
x268 = 8.0*x267
x269 = m[x116 + x62]
x270 = m[x120 - 2]
x271 = 8.0*x270
x272 = m[x119 + x54]
x273 = 64.0*x272
x274 = m[x119 + x62]
x275 = 8.0*x274
x276 = m[x124 - 2]
x277 = 8.0*x276
x278 = m[x123 + x54]
x279 = 64.0*x278
x280 = m[x123 + x62]
x281 = 8.0*x280
x282 = -x281
x283 = m[x82 - 2]
x284 = m[x54 + x81]
x285 = 8.0*x284
x286 = -x285
x287 = -8.0*x133
x288 = 8.0*x149
x289 = -8.0*x171
x290 = 8.0*x180
x291 = 1/(-zvec[x62] + zvec[k])
x292 = 0.166666666666667*x291
x293 = m[x3 + 3]
x294 = 8.0*x4
x295 = -8.0*x63
x296 = -8.0*x226
x297 = 64.0*x229
x298 = 8.0*x211
x299 = 8.0*x216
x300 = -64.0*x232
x301 = 8.0*x239
x302 = 8.0*x244
x303 = 64.0*x249
x304 = -x303
x305 = 64.0*x255
x306 = -8.0*x262
x307 = 8.0*x272
x308 = -8.0*x278
x309 = 0.0277777777777778*x291*x9
x310 = 8.0*x142
x311 = 64.0*x151
x312 = 8.0*x137
x313 = -x312
x314 = m[x32 + 3]
x315 = -8.0*x15
x316 = 64.0*x20
x317 = 64.0*x26
x318 = -x317
x319 = 8.0*x33
x320 = m[x14 + 3]
x321 = 8.0*x131
x322 = -x321
x323 = -64.0*x143
x324 = m[x19 + 3]
x325 = 8.0*x324
x326 = 8.0*x147
x327 = m[x25 + 3]
x328 = 8.0*x327
x329 = 0.0277777777777778*x291*x36
x330 = 64.0*x43
x331 = 8.0*x162
x332 = -x331
x333 = 64.0*x173
x334 = -x333
x335 = m[x42 + 3]
x336 = 8.0*x335
x337 = m[x51 + 3]
x338 = 8.0*x40
x339 = 64.0*x47
x340 = 8.0*x52
x341 = 8.0*x170
x342 = m[x39 + 3]
x343 = 8.0*x177
x344 = 64.0*x182
x345 = m[x46 + 3]
x346 = 8.0*x345
x347 = 8.0*x165
x348 = 0.00462962962962963*x291*x36
x349 = 64.0*x228
x350 = 8.0*x237
x351 = 64.0*x274
x352 = 8.0*x206
x353 = 8.0*x225 + 8.0*x269 - x349 - x350 - x351 - x352
x354 = m[x82 + 3]
x355 = 8.0*x90
x356 = 64.0*x101
x357 = 64.0*x113
x358 = -x357
x359 = -8.0*x118
x360 = 64.0*x121
x361 = -64.0*x125
x362 = 8.0*x83
x363 = 8.0*x210
x364 = m[x66 + 3]
x365 = -8.0*x215
x366 = m[x85 + 3]
x367 = 8.0*x223
x368 = m[x89 + 3]
x369 = m[x92 + 3]
x370 = m[x70 + 3]
x371 = 64.0*x231
x372 = m[x96 + 3]
x373 = 64.0*x193
x374 = m[x100 + 3]
x375 = 8.0*x374
x376 = -8.0*x243
x377 = m[x104 + 3]
x378 = 64.0*x248
x379 = m[x74 + 3]
x380 = -64.0*x254
x381 = m[x108 + 3]
x382 = 8.0*x260
x383 = 64.0*x264
x384 = m[x112 + 3]
x385 = -8.0*x384
x386 = m[x117 + 3]
x387 = 8.0*x196
x388 = m[x120 + 3]
x389 = 8.0*x388
x390 = 8.0*x198
x391 = 64.0*x280
x392 = m[x124 + 3]
x393 = -8.0*x392
x394 = -8.0*x97
x395 = 8.0*x109
x396 = 1/(-yvec[x49] + yvec[j])
x397 = 0.166666666666667*x396
x398 = nz*(x1 + 3.0)
x399 = x398 + k
x400 = m[x399]
x401 = -x340
x402 = -x343
x403 = 0.0277777777777778*x396*x9
x404 = 8.0*x71
x405 = -x404
x406 = -x356
x407 = nz*(x17 + 3.0)
x408 = x407 + k
x409 = m[x408]
x410 = 8.0*x409
x411 = nz*(x30 + 3.0)
x412 = x411 + k
x413 = m[x412]
x414 = nz*(x12 + 3.0)
x415 = x414 + k
x416 = m[x415]
x417 = 8.0*x75
x418 = nz*(x23 + 3.0)
x419 = x418 + k
x420 = m[x419]
x421 = 8.0*x420
x422 = -8.0*x230
x423 = 8.0*x233
x424 = 8.0*x253
x425 = -8.0*x256
x426 = 0.0277777777777778*x396*x55
x427 = -8.0*x184
x428 = 64.0*x185
x429 = m[x399 - 2]
x430 = 8.0*x57
x431 = 64.0*x58
x432 = -x431
x433 = 64.0*x60
x434 = -64.0*x187
x435 = m[x398 + x53]
x436 = 8.0*x435
x437 = m[x398 + x54]
x438 = 8.0*x437
x439 = m[x398 + x62]
x440 = 64.0*x154
x441 = 8.0*x191
x442 = m[x407 + x62]
x443 = 64.0*x200
x444 = m[x411 + x53]
x445 = x373 + x387 - x440 + x441 - 8.0*x442 + x443 - 8.0*x444
x446 = m[x414 + x62]
x447 = m[x412 - 2]
x448 = 8.0*x153
x449 = 64.0*x156
x450 = m[x414 + x53]
x451 = m[x414 + x54]
x452 = 8.0*x227
x453 = m[x408 - 2]
x454 = m[x407 + x53]
x455 = m[x407 + x54]
x456 = 8.0*x247
x457 = 8.0*x250
x458 = m[x419 - 2]
x459 = m[x418 + x53]
x460 = m[x418 + x54]
x461 = m[x418 + x62]
x462 = 8.0*x461
x463 = 8.0*x283
x464 = 64.0*x284
x465 = m[x411 + x54]
x466 = 8.0*x465
x467 = m[x411 + x62]
x468 = -8.0*x232
x469 = 8.0*x255
x470 = 8.0*x229
x471 = 8.0*x455
x472 = 8.0*x249
x473 = 8.0*x460
x474 = -x110
x475 = 8.0*x231
x476 = -x475
x477 = 8.0*x372
x478 = 8.0*x254
x479 = -8.0*x381
x480 = 0.0277777777777778*x291*x396
x481 = 64.0*x165
x482 = -8.0*x439
x483 = 64.0*x4
x484 = m[x399 + 3]
x485 = 64.0*x52
x486 = -x485
x487 = 8.0*x400
x488 = -64.0*x63
x489 = 8.0*x293
x490 = 8.0*x337
x491 = 0.00462962962962963*x291*x396
x492 = 8.0*x228
x493 = 64.0*x206
x494 = -x238 - x275 + 8.0*x454 + 8.0*x467 - x492 - x493
x495 = -64.0*x33
x496 = 64.0*x83
x497 = -8.0*x413
x498 = 64.0*x137
x499 = -8.0*x314
x500 = 8.0*x370
x501 = m[x408 + 3]
x502 = 8.0*x248
x503 = 8.0*x379
x504 = -8.0*x459
x505 = m[x419 + 3]
x506 = 8.0*x354
x507 = m[x412 + 3]
x508 = xvec[i]
x509 = xvec[x28]
x510 = 1/(x508 - x509)
x511 = 0.166666666666667*x510
x512 = ny*(i + 3.0)
x513 = x512 + j
x514 = nz*x513
x515 = x514 + k
x516 = m[x515]
x517 = -x417
x518 = -x326
x519 = 0.0277777777777778*x36*x510
x520 = nz*(x513 - 2.0)
x521 = x520 + k
x522 = m[x521]
x523 = nz*(x34 + x512)
x524 = x523 + k
x525 = m[x524]
x526 = 8.0*x525
x527 = nz*(x35 + x512)
x528 = x527 + k
x529 = m[x528]
x530 = 8.0*x529
x531 = nz*(x49 + x512)
x532 = x531 + k
x533 = m[x532]
x534 = 0.0277777777777778*x510*x55
x535 = -x449
x536 = m[x514 + x54]
x537 = 8.0*x536
x538 = m[x515 - 2]
x539 = m[x514 + x53]
x540 = 8.0*x539
x541 = m[x514 + x62]
x542 = -x344
x543 = m[x531 + x62]
x544 = m[x520 + x53]
x545 = m[x520 + x54]
x546 = m[x520 + x62]
x547 = m[x524 - 2]
x548 = m[x523 + x53]
x549 = m[x523 + x54]
x550 = m[x523 + x62]
x551 = 8.0*x550
x552 = m[x528 - 2]
x553 = m[x527 + x53]
x554 = m[x527 + x54]
x555 = m[x527 + x62]
x556 = -8.0*x555
x557 = m[x532 - 2]
x558 = m[x53 + x531]
x559 = 8.0*x558
x560 = m[x531 + x54]
x561 = -8.0*x560
x562 = -x472
x563 = -8.0*x549
x564 = -x279
x565 = 8.0*x554
x566 = 0.0277777777777778*x291*x510
x567 = -8.0*x541
x568 = m[x515 + 3]
x569 = 8.0*x516
x570 = m[x532 + 3]
x571 = 8.0*x533
x572 = 8.0*x548
x573 = m[x524 + 3]
x574 = 8.0*x553
x575 = m[x528 + 3]
x576 = 8.0*x543
x577 = 0.0277777777777778*x396*x510
x578 = nz*(x513 + 3.0)
x579 = x578 + k
x580 = m[x579]
x581 = m[x53 + x578]
x582 = m[x54 + x578]
x583 = m[x578 + x62]
return np.array([x4, -x10*(x15 + x21 + x27 - x33), -x37*(x40 + x44 + x48 - x52), -x56*(x57 + x59 + x61 - x63),
x64*(x102 + x106 + x110 + x115 - x118 + x122 + x127 + x68 + x72 + x77 + x80 + x83 + x87 - x90 -
x94 - x98),
x128*(x130 + x132 + x135 + x136 + x137 + x139 + x141 - x142 + x144 + x146 + x148 + x150 + x152 -
x153 + x155 + x158),
x159*(x161 + x163 + x164 + x165 - x167 + x169 - x170 + x172 + x174 + x176 + x179 + x181 + x183 -
x184 + x186 + x188),
0.00462962962962963*x189*x36*(x202 + x203 + x204 + x205 + x206 + 8.0*x207 - 8.0*x208 + 8.0*x209 -
64.0*x210 + 64.0*x211 - x213 - 8.0*x214 + 64.0*x215 - 64.0*x216 +
x218 - x220 + x222 - x223 + 8.0*x224 - 64.0*x225 + 64.0*x226 -
64.0*x227 + 512.0*x228 - 512.0*x229 + 64.0*x230 - 512.0*x231 +
512.0*x232 + x234 - x236 + x238 + x241 - 8.0*x242 + 64.0*x243 -
64.0*x244 + x246 + 64.0*x247 - 512.0*x248 + 512.0*x249 + x252 -
64.0*x253 + 512.0*x254 - 512.0*x255 + x257 + x259 - x261 + x263 +
x266 + x268 - x269 - x271 - x273 + x275 + x277 + x279 + x282 - x283 +
x286 - m[x79 - 2]), x60, -x10*(x140 - x156 + x287 + x288),
-x37*(x168 - x187 + x289 + x290), x292*(x293 + x294 + x295 - x58),
x64*(x208 - x221 - x267 + x284 + x296 + x297 - x298 + x299 + x300 + x301 + x302 + x304 + x305 + x306 +
x307 + x308), x309*(x138 - x154 + x310 + x311 + x313 + x314 + x315 + x316 + x318 + x319 - x320 +
x322 + x323 + x325 + x326 - x328), x329*(x166 - x185 + x330 + x332 + x334 +
x336 + x337 - x338 - x339 + x340 +
x341 - x342 + x343 + x344 - x346 -
x347),
x348*x9*(64.0*x105 + 512.0*x109 + 64.0*x190 - 512.0*x191 + x195 - x200 - 8.0*x203 - x207 + 64.0*x212 -
64.0*x217 + x219 + 512.0*x233 - 64.0*x245 + 512.0*x250 - 512.0*x256 + x353 + x354 - x355 +
x356 + x358 + x359 + x360 + x361 + x362 + x363 - 8.0*x364 + x365 + 8.0*x366 + x367 - x368 -
8.0*x369 + 64.0*x370 + x371 - 64.0*x372 - x373 + x375 + x376 + 8.0*x377 + x378 - 64.0*x379 +
x380 + 64.0*x381 + x382 + x383 + x385 - x386 - x387 + x389 + x390 + x391 + x393 - 64.0*x67 +
512.0*x71 - 512.0*x75 + 8.0*x80 + 64.0*x86 - 64.0*x93 - 512.0*x97 + m[x79 + 3]),
x47, -x10*(-x125 + x394 + x395 + x86), x397*(x294 + x400 + x401 - x43), -x56*(x175 - x182 + x290 + x402),
x403*(-x121 + x315 + x316 + x318 + x319 + x355 + x357 - x362 + x405 + x406 + x410 + x413 - x416 + x417 -
x421 + x67), x128*(x214 - x217 - x276 + x280 + x299 + x300 + x305 + x308 + x365 + x371 + x380 +
x390 + x422 + x423 + x424 + x425), -x426*(-x160 + x173 + x289 + x295 + x331 +
x347 + x427 + x428 + x429 + x430 +
x432 + x433 + x434 - x436 + x438 -
x439),
-0.00462962962962963*x189*x396*(64.0*x129 - 512.0*x131 + 512.0*x133 - 8.0*x136 + 64.0*x138 -
64.0*x140 - 64.0*x145 + 512.0*x147 - 512.0*x149 + 8.0*x204 + x209 -
x212 - 64.0*x219 + 64.0*x221 - 64.0*x235 + 512.0*x237 - 512.0*x239 +
64.0*x258 - 512.0*x260 + 512.0*x262 - x270 + x274 - x297 + x298 +
x303 - x307 + x310 + x311 + x313 + x323 + x349 + x352 - x363 - x367 -
x378 - x383 + x445 + x446 + x447 + x448 + x449 + 8.0*x450 - 8.0*x451 -
x452 + 8.0*x453 - 64.0*x454 + 64.0*x455 + x456 - x457 - 8.0*x458 +
64.0*x459 - 64.0*x460 + x462 - x463 - x464 + x466 - x467 -
m[x415 - 2]), x180, -x10*(x216 - x278 + x468 + x469),
x397*(-x171 + x188 + x437 + x61), x292*(-x177 + x183 + x345 + x48), x403*(x134 - x141 - x150 + x157 +
x211 + x222 + x241 + x263 -
x272 + x286 - x451 + x465 -
x470 + x471 + x472 - x473),
x309*(x126 - x198 + x215 + x218 + x234 + x257 + x282 - x366 + x392 + x474 + x476 + x477 + x478 +
x479 - x87 + x98), x480*(x162 + x174 + x186 - x335 - x435 + x44 + x481 + x482 + x483 + x484 +
x486 + x487 + x488 + x489 - x490 + x59), x491*(512.0*x101 - 512.0*x113 +
x122 + x132 + x139 -
64.0*x142 + 512.0*x143 +
x148 + 64.0*x15 -
512.0*x151 + x155 - x192 -
512.0*x193 - x196 -
512.0*x20 - x201 + x210 +
x213 + x220 + 64.0*x223 +
x251 + 512.0*x26 + x261 +
512.0*x264 + 8.0*x320 -
64.0*x324 + 64.0*x327 -
x364 - 8.0*x368 +
64.0*x374 - 64.0*x384 +
x388 - 64.0*x409 + 8.0*x416 +
64.0*x420 + 64.0*x442 + x444 -
8.0*x446 - x450 - 64.0*x461 +
x494 + x495 + x496 + x497 + x498 +
x499 + x500 - 8.0*x501 + x502 - x503 +
x504 + 8.0*x505 + x506 - x507 + x68 +
x72 + x77 - 64.0*x90 +
m[x415 + 3])/(-x6 + x8), x26,
x511*(-x20 + x294 - x319 + x516), -x37*(x105 - x113 + x395 + x517), -x56*(x145 - x151 + x288 + x518),
-x519*(x101 - x330 + x338 + x339 + x359 + x360 + x361 + x362 + x394 + x401 + x404 + x522 - x526 + x530 - x533 - x93),
-x534*(-x129 + x143 + x287 + x295 + x312 + x321 + x430 + x432 + x433 + x440 - x448 + x535 + x537 + x538 - x540 - x541),
x159*(x242 - x245 - x258 + x264 + x302 + x304 + x305 + x306 + x376 + x378 + x380 + x382 + x424 + x425 - x456 + x457),
0.00462962962962963*x36*x510*x55*(-64.0*x160 + 512.0*x162 + 8.0*x164 - 64.0*x166 + 64.0*x168 - 512.0*x171 +
64.0*x175 - 512.0*x177 + 512.0*x180 + x190 - x193 + 64.0*x195 - 512.0*x196 + 512.0*x198 -
8.0*x205 - x224 + x235 - 64.0*x267 + 64.0*x270 + 512.0*x272 - 64.0*x276 - 512.0*x278 +
x296 + x297 + x300 + x301 + x333 - x341 + x347 + x353 + x371 + x391 + x422 + x423 + x427 +
x428 + x434 - x441 - x443 + x452 + x463 + x464 + x542 + x543 - 8.0*x544 + 8.0*x545 - x546 -
8.0*x547 + 64.0*x548 - 64.0*x549 + x551 + 8.0*x552 - 64.0*x553 + 64.0*x554 + x556 - x557 +
x559 + x561 + m[x521 - 2]), x149, x511*(-x133 + x158 + x536 + x61),
-x37*(x244 - x262 + x469 + x562), x292*(-x147 + x152 + x27 + x327), -x519*(x169 + x172 + x181 + x188 - x226 + x239 - x268 +
x273 + x285 + x468 + x470 + x545 - x560 + x563 +
x564 + x565), x566*(x131 + x144 + x155 + x21 -
x324 + x483 + x488 + x489 + x495 + x498 + x499 - x539 + x567 + x568 + x569 + x59), x329*(-x106 + x114 + x243 + x246 + x252 + x257 - x260 + x266 - x377 + x384 + x474 + x478 + x479 - x502 + x503 + x76), x348*x510*(-x102 + 64.0*x118 - 512.0*x121 + 512.0*x125 - x163 + x167 + 64.0*x170 - 512.0*x173 + x178 + 512.0*x182 - x186 + x202 - x225 + x234 + x237 - 64.0*x269 + 512.0*x274 - 512.0*x280 + 64.0*x335 - 8.0*x342 - 64.0*x345 + x369 - x374 + 8.0*x386 - 64.0*x388 + 64.0*x392 - 64.0*x40 + 512.0*x43 - 512.0*x47 + x476 + x477 - x481 + x485 + x490 + x492 + x493 - x496 - x500 - x506 - 8.0*x522 + 64.0*x525 - 64.0*x529 + x544 + 8.0*x546 - 64.0*x550 + 64.0*x555 - x558 + x570 + x571 - x572 + 8.0*x573 + x574 - 8.0*x575 - x576 - x72 + x94 + x98 - m[x521 + 3]), x109, x511*(x127 + x48 + x529 - x97), x397*(x115 + x27 + x420 - x75), -x56*(x253 - x256 + x469 - x478), x577*(x102 + x122 + x21 - x409 + x44 + x483 + x486 + x487 + x495 + x496 + x497 - x525 + x569 - x571 + x580 + x71), -x534*(x176 + x179 + x181 + x183 + x199 - x230 + x233 - x277 + x281 + x468 + x475 + x552 - x555 + x564 + x565 - x574), -x426*(x146 + x148 + x150 + x152 - x247 + x250 - x259 + x261 - x263 + x265 + x458 - x461 + x473 + x502 + x504 + x562), -0.00462962962962963*x396*x510*x55*(x130 + x132 + x135 + x144 - 64.0*x153 + 512.0*x154 - 512.0*x156 + x161 + x163 + x172 + x174 - 64.0*x184 + 512.0*x185 - 512.0*x187 - x191 - x194 - x197 - 512.0*x200 + x227 + x236 + x240 + x271 + x273 + 64.0*x283 + 512.0*x284 + 8.0*x429 - 64.0*x435 + 64.0*x437 + x442 + 64.0*x444 - 8.0*x447 - x453 - 64.0*x465 + x470 - x471 + x481 + x482 + x488 + x494 + x498 + 64.0*x536 + 8.0*x538 - 64.0*x539 - x547 + x550 - 8.0*x557 + 64.0*x558 - 64.0*x560 + x563 + x567 + 64.0*x57 + x572 + x576 - 512.0*x58 - 8.0*x581 + 8.0*x582 - x583 + 512.0*x60 + m[x579 - 2]), x255, x511*(-x232 + x290 + x308 + x554), x397*(-x249 + x288 + x306 + x460), x292*(-x254 + x381 + x395 + x425), x577*(x229 + x287 + x289 + x301 + x307 + x433 + x434 + x438 - x455 + x464 - x466 + x535 + x537 - x549 + x561 + x582), x566*(x231 + x339 + x346 + x361 - x372 + x390 + x391 + x393 + x394 + x402 + x423 + x530 + x542 - x553 + x556 + x575), x480*(x248 - x311 + x317 + x328 + x358 - x379 + x382 + x383 + x385 + x421 + x457 - x459 - x462 + x505 + x517 + x518), x491*(-512.0*x137 - 512.0*x165 + 512.0*x206 + x228 - 64.0*x293 + 64.0*x314 + x316 + x322 + x323 + x325 + 512.0*x33 + x330 + x332 + x334 + x336 + 64.0*x337 + x350 + x351 - 64.0*x354 - x360 - x370 - x375 - x389 - 512.0*x4 - 64.0*x400 + x405 + x406 + x410 + 64.0*x413 - x428 + x431 + x436 + 64.0*x439 + x445 - x454 - 64.0*x467 - 8.0*x484 + x501 + 8.0*x507 - 64.0*x516 + 512.0*x52 + x526 + 64.0*x533 + x540 + 64.0*x541 - 64.0*x543 - x548 - x551 - x559 - 8.0*x568 + 8.0*x570 + x573 - 8.0*x580 + x581 + 8.0*x583 + 512.0*x63 - 512.0*x83 - m[x579 + 3])/(-x508 + x509)])
@jit(nogil=True,nopython=True)
def interp_diff(x,y,z,xvec,yvec,zvec,A_ijk):
u = (x - xvec[xi])/(xvec[xi+1] - xvec[xi])
v = (y - yvec[yi])/(yvec[yi+1] - yvec[yi])
w = (z - zvec[zi])/(zvec[zi+1] - zvec[zi])
#print(u,v,w)
x,y,z = u,v,w
x0 = z**2
x1 = z**3
x2 = y*z
x3 = x0*y
x4 = x1*y
x5 = y**2
x6 = x5*z
x7 = x0*x5
x8 = x1*x5
x9 = y**3
x10 = x9*z
x11 = x0*x9
x12 = x1*x9
x13 = x*z
x14 = x*x0
x15 = x*x1
x16 = x*y
x17 = x16*z
x18 = x0*x16
x19 = x*x5
x20 = x13*x5
x21 = x0*x19
x22 = x*x9
x23 = x0*x22
x24 = x**2
x25 = x24*z
x26 = x0*x24
x27 = x1*x24
x28 = x24*y
x29 = x2*x24
x30 = x0*x28
x31 = x24*x5
x32 = x25*x5
x33 = x0*x31
x34 = x24*x9
x35 = x0*x34
x36 = x**3
x37 = x36*z
x38 = x0*x36
x39 = x1*x36
x40 = x36*y
x41 = x0*x40
x42 = x36*x5
x43 = x37*x5
x44 = x0*x42
x45 = x36*x9
x46 = 2*x
x47 = x46*z
x48 = x0*x46
x49 = x1*x46
x50 = x46*y
x51 = x2*x46
x52 = x0*x50
x53 = x1*x50
x54 = x46*x5
x55 = x47*x5
x56 = x46*x9
x57 = x47*x9
x58 = 3*x24
x59 = x58*z
x60 = x0*x58
x61 = x1*x58
x62 = x58*y
x63 = x3*x58
x64 = x5*x58
x65 = x58*x6
x66 = x58*x7
x67 = x58*x8
x68 = x58*x9
x69 = x11*x58
x70 = 2*y
x71 = x70*z
x72 = x0*x70
x73 = x1*x70
x74 = 3*x5
x75 = x74*z
x76 = x0*x74
x77 = x1*x74
x78 = x*x74
x79 = x14*x74
x80 = x24*x70
x81 = x24*x71
x82 = x36*x70
x83 = x36*x71
x84 = x36*x74
x85 = x38*x74
x86 = 2*z
x87 = 3*x0
x88 = x87*y
x89 = x5*x86
x90 = x86*x9
x91 = x87*x9
x92 = x*x87
x93 = x24*x86
x94 = x36*x86
x95 = x36*x87
x96 = 4*x16
x97 = x96*z
x98 = 6*x19
x99 = 6*x20
x100 = x0*x98
x101 = 6*x28
x102 = 6*x29
x103 = x0*x101
x104 = 9*x31
x105 = x0*x104
x106 = 4*x13
x107 = 6*x14
x108 = 6*x18
x109 = 6*x25
x110 = 9*x26
x111 = x109*x5
x112 = 4*x2
x113 = 6*x3
x114 = 6*x6
x115 = 9*x7
B = np.array([[1, z, x0, x1, y, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x, x13, x14, x15, x16, x17, x18, x1*x16, x19, x20, x21, x1*x19, x22, x13*x9, x23, x1*x22, x24, x25, x26, x27, x28, x29, x30, x1*x28, x31, x32, x33, x1*x31, x34, x25*x9, x35, x1*x34, x36, x37, x38, x39, x40, x2*x36, x41, x1*x40, x42, x43, x44, x1*x42, x45, x37*x9, x0*x45, x1*x45],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, z, x0, x1, y, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x46, x47, x48, x49, x50, x51, x52, x53, x54, x55, x0*x54, x1*x54, x56, x57, x0*x56, x1*x56, x58, x59, x60, x61, x62, x2*x58, x63, x4*x58, x64, x65, x66, x67, x68, x10*x58, x69, x12*x58],
[0, 0, 0, 0, 1, z, x0, x1, x70, x71, x72, x73, x74, x75, x76, x77, 0, 0, 0, 0, x, x13, x14, x15, x50, x51, x52, x53, x78, x13*x74, x79, x15*x74, 0, 0, 0, 0, x24, x25, x26, x27, x80, x81, x0*x80, x1*x80, x64, x65, x66, x67, 0, 0, 0, 0, x36, x37, x38, x39, x82, x83, x0*x82, x1*x82, x84, x37*x74, x85, x39*x74],
[0, 1, x86, x87, 0, y, x71, x88, 0, x5, x89, x76, 0, x9, x90, x91, 0, x, x47, x92, 0, x16, x51, x16*x87, 0, x19, x55, x79, 0, x22, x57, x22*x87, 0, x24, x93, x60, 0, x28, x81, x63, 0, x31, x5*x93, x66, 0, x34, x9*x93, x69, 0, x36, x94, x95, 0, x40, x83, x40*x87, 0, x42, x5*x94, x85, 0, x45, x9*x94, x45*x87],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, z, x0, x1, x70, x71, x72, x73, x74, x75, x76, x77, 0, 0, 0, 0, x46, x47, x48, x49, x96, x97, x0*x96, x1*x96, x98, x99, x100, x1*x98, 0, 0, 0, 0, x58, x59, x60, x61, x101, x102, x103, x1*x101, x104, 9*x32, x105, x1*x104],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, x86, x87, 0, y, x71, x88, 0, x5, x89, x76, 0, x9, x90, x91, 0, x46, x106, x107, 0, x50, x97, x108, 0, x54, x106*x5, x100, 0, x56, x106*x9, 6*x23, 0, x58, x109, x110, 0, x62, x102, 9*x30, 0, x64, x111, x105, 0, x68, x109*x9, 9*x35],
[0, 0, 0, 0, 0, 1, x86, x87, 0, x70, x112, x113, 0, x74, x114, x115, 0, 0, 0, 0, 0, x, x47, x92, 0, x50, x97, x108, 0, x78, x99, 9*x21, 0, 0, 0, 0, 0, x24, x93, x60, 0, x80, x112*x24, x103, 0, x64, x111, x105, 0, 0, 0, 0, 0, x36, x94, x95, 0, x82, x112*x36, 6*x41, 0, x84, 6*x43, 9*x44],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, x86, x87, 0, x70, x112, x113, 0, x74, x114, x115, 0, 0, 0, 0, 0, x46, x106, x107, 0, x96, 8*x17, 12*x18, 0, x98, 12*x20, 18*x21, 0, 0, 0, 0, 0, x58, x109, x110, 0, x101, 12*x29, 18*x30, 0, x104, 18*x32, 27*x33]],
dtype=np.double)
bvec = B.dot(A_ijk)
return bvec
@jit(nogil=True,nopython=True)
def interp_scalar(x,y,z,xvec,yvec,zvec,A_ijk):
'''Interpolate and return f,fx,fy,fz if doDouble is True then return also fxy,fxz, fyz.
Double derivatives do not yet work I think.'''
u = (x - xvec[xi])/(xvec[xi+1] - xvec[xi])
v = (y - yvec[yi])/(yvec[yi+1] - yvec[yi])
w = (z - zvec[zi])/(zvec[zi+1] - zvec[zi])
#print(u,v,w)
x,y,z = u,v,w
x0 = z**2
x1 = z**3
x2 = y*z
x3 = x0*y
x4 = x1*y
x5 = y**2
x6 = x5*z
x7 = x0*x5
x8 = x1*x5
x9 = y**3
x10 = x9*z
x11 = x0*x9
x12 = x1*x9
x13 = x*z
x14 = x*x0
x15 = x*x1
x16 = x*y
x17 = x16*z
x18 = x0*x16
x19 = x*x5
x20 = x13*x5
x21 = x0*x19
x22 = x*x9
x23 = x0*x22
x24 = x**2
x25 = x24*z
x26 = x0*x24
x27 = x1*x24
x28 = x24*y
x29 = x2*x24
x30 = x0*x28
x31 = x24*x5
x32 = x25*x5
x33 = x0*x31
x34 = x24*x9
x35 = x0*x34
x36 = x**3
x37 = x36*z
x38 = x0*x36
x39 = x1*x36
x40 = x36*y
x41 = x0*x40
x42 = x36*x5
x43 = x37*x5
x44 = x0*x42
x45 = x36*x9
x46 = 2*x
x47 = x46*z
x48 = x0*x46
x49 = x1*x46
x50 = x46*y
x51 = x2*x46
x52 = x0*x50
x53 = x1*x50
x54 = x46*x5
x55 = x47*x5
x56 = x46*x9
x57 = x47*x9
x58 = 3*x24
x59 = x58*z
x60 = x0*x58
x61 = x1*x58
x62 = x58*y
x63 = x3*x58
x64 = x5*x58
x65 = x58*x6
x66 = x58*x7
x67 = x58*x8
x68 = x58*x9
x69 = x11*x58
x70 = 2*y
x71 = x70*z
x72 = x0*x70
x73 = x1*x70
x74 = 3*x5
x75 = x74*z
x76 = x0*x74
x77 = x1*x74
x78 = x*x74
x79 = x14*x74
x80 = x24*x70
x81 = x24*x71
x82 = x36*x70
x83 = x36*x71
x84 = x36*x74
x85 = x38*x74
x86 = 2*z
x87 = 3*x0
x88 = x87*y
x89 = x5*x86
x90 = x86*x9
x91 = x87*x9
x92 = x*x87
x93 = x24*x86
x94 = x36*x86
x95 = x36*x87
x96 = 4*x16
x97 = x96*z
x98 = 6*x19
x99 = 6*x20
x100 = x0*x98
x101 = 6*x28
x102 = 6*x29
x103 = x0*x101
x104 = 9*x31
x105 = x0*x104
x106 = 4*x13
x107 = 6*x14
x108 = 6*x18
x109 = 6*x25
x110 = 9*x26
x111 = x109*x5
x112 = 4*x2
x113 = 6*x3
x114 = 6*x6
x115 = 9*x7
B = np.array([[1, z, x0, x1, y, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x, x13, x14, x15, x16, x17, x18, x1*x16, x19, x20, x21, x1*x19, x22, x13*x9, x23, x1*x22, x24, x25, x26, x27, x28, x29, x30, x1*x28, x31, x32, x33, x1*x31, x34, x25*x9, x35, x1*x34, x36, x37, x38, x39, x40, x2*x36, x41, x1*x40, x42, x43, x44, x1*x42, x45, x37*x9, x0*x45, x1*x45],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, z, x0, x1, y, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x46, x47, x48, x49, x50, x51, x52, x53, x54, x55, x0*x54, x1*x54, x56, x57, x0*x56, x1*x56, x58, x59, x60, x61, x62, x2*x58, x63, x4*x58, x64, x65, x66, x67, x68, x10*x58, x69, x12*x58],
[0, 0, 0, 0, 1, z, x0, x1, x70, x71, x72, x73, x74, x75, x76, x77, 0, 0, 0, 0, x, x13, x14, x15, x50, x51, x52, x53, x78, x13*x74, x79, x15*x74, 0, 0, 0, 0, x24, x25, x26, x27, x80, x81, x0*x80, x1*x80, x64, x65, x66, x67, 0, 0, 0, 0, x36, x37, x38, x39, x82, x83, x0*x82, x1*x82, x84, x37*x74, x85, x39*x74],
[0, 1, x86, x87, 0, y, x71, x88, 0, x5, x89, x76, 0, x9, x90, x91, 0, x, x47, x92, 0, x16, x51, x16*x87, 0, x19, x55, x79, 0, x22, x57, x22*x87, 0, x24, x93, x60, 0, x28, x81, x63, 0, x31, x5*x93, x66, 0, x34, x9*x93, x69, 0, x36, x94, x95, 0, x40, x83, x40*x87, 0, x42, x5*x94, x85, 0, x45, x9*x94, x45*x87],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, z, x0, x1, x70, x71, x72, x73, x74, x75, x76, x77, 0, 0, 0, 0, x46, x47, x48, x49, x96, x97, x0*x96, x1*x96, x98, x99, x100, x1*x98, 0, 0, 0, 0, x58, x59, x60, x61, x101, x102, x103, x1*x101, x104, 9*x32, x105, x1*x104],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, x86, x87, 0, y, x71, x88, 0, x5, x89, x76, 0, x9, x90, x91, 0, x46, x106, x107, 0, x50, x97, x108, 0, x54, x106*x5, x100, 0, x56, x106*x9, 6*x23, 0, x58, x109, x110, 0, x62, x102, 9*x30, 0, x64, x111, x105, 0, x68, x109*x9, 9*x35],
[0, 0, 0, 0, 0, 1, x86, x87, 0, x70, x112, x113, 0, x74, x114, x115, 0, 0, 0, 0, 0, x, x47, x92, 0, x50, x97, x108, 0, x78, x99, 9*x21, 0, 0, 0, 0, 0, x24, x93, x60, 0, x80, x112*x24, x103, 0, x64, x111, x105, 0, 0, 0, 0, 0, x36, x94, x95, 0, x82, x112*x36, 6*x41, 0, x84, 6*x43, 9*x44],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, x86, x87, 0, x70, x112, x113, 0, x74, x114, x115, 0, 0, 0, 0, 0, x46, x106, x107, 0, x96, 8*x17, 12*x18, 0, x98, 12*x20, 18*x21, 0, 0, 0, 0, 0, x58, x109, x110, 0, x101, 12*x29, 18*x30, 0, x104, 18*x32, 27*x33]],
dtype=np.double)
bvec = B.dot(A_ijk)
return bvec[0]
@jit(nogil=True,nopython=True)
def computeInterpolant(xi,yi,zi,nx,ny,nz,xvec,yvec,zvec,m):
Binv = np.array([
[1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-3,0,0,-2,0,0,0,0,3,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[2,0,0,1,0,0,0,0,-2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,-3,0,0,0,-2,0,0,0,3,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,2,0,0,0,1,0,0,0,-2,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-3,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,-3,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,3,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[9,0,6,6,0,0,4,0,-9,0,-6,3,0,0,2,0,-9,0,3,-6,0,0,2,0,9,0,-3,-3,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-6,0,-4,-3,0,0,-2,0,6,0,4,-3,0,0,-2,0,6,0,-2,3,0,0,-1,0,-6,0,2,3,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-6,0,-3,-4,0,0,-2,0,6,0,3,-2,0,0,-1,0,6,0,-3,4,0,0,-2,0,-6,0,3,2,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[4,0,2,2,0,0,1,0,-4,0,-2,2,0,0,1,0,-4,0,2,-2,0,0,1,0,4,0,-2,-2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,-3,0,0,0,-2,0,0,0,3,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,2,0,0,0,1,0,0,0,-2,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,-3,0,0,-2,0,0,0,0,3,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,2,0,0,1,0,0,0,0,-2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,-3,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,3,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,-3,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,9,0,0,6,6,0,4,0,-9,0,0,-6,3,0,2,0,-9,0,0,3,-6,0,2,0,9,0,0,-3,-3,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,-6,0,0,-4,-3,0,-2,0,6,0,0,4,-3,0,-2,0,6,0,0,-2,3,0,-1,0,-6,0,0,2,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,-6,0,0,-3,-4,0,-2,0,6,0,0,3,-2,0,-1,0,6,0,0,-3,4,0,-2,0,-6,0,0,3,2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,4,0,0,2,2,0,1,0,-4,0,0,-2,2,0,1,0,-4,0,0,2,-2,0,1,0,4,0,0,-2,-2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-3,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,-3,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[9,6,0,6,0,4,0,0,-9,-6,0,3,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-9,3,0,-6,0,2,0,0,9,-3,0,-3,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-6,-4,0,-3,0,-2,0,0,6,4,0,-3,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,-2,0,3,0,-1,0,0,-6,2,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,-3,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,-3,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,9,0,6,0,6,4,0,0,-9,0,-6,0,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-9,0,3,0,-6,2,0,0,9,0,-3,0,-3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,-6,0,-4,0,-3,-2,0,0,6,0,4,0,-3,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,-2,0,3,-1,0,0,-6,0,2,0,3,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[9,6,6,0,4,0,0,0,0,0,0,0,0,0,0,0,-9,-6,3,0,2,0,0,0,0,0,0,0,0,0,0,0,-9,3,-6,0,2,0,0,0,0,0,0,0,0,0,0,0,9,-3,-3,0,1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,9,0,6,6,4,0,0,0,0,0,0,0,0,0,0,0,-9,0,-6,3,2,0,0,0,0,0,0,0,0,0,0,0,-9,0,3,-6,2,0,0,0,0,0,0,0,0,0,0,0,9,0,-3,-3,1,0,0,0,0,0,0,0,0],
[-27,-18,-18,-18,-12,-12,-12,-8,27,18,18,-9,12,-6,-6,-4,27,18,-9,18,-6,12,-6,-4,-27,-18,9,9,6,6,-3,-2,27,-9,18,18,-6,-6,12,-4,-27,9,-18,9,6,-3,6,-2,-27,9,9,-18,-3,6,6,-2,27,-9,-9,-9,3,3,3,-1],
[18,12,12,9,8,6,6,4,-18,-12,-12,9,-8,6,6,4,-18,-12,6,-9,4,-6,3,2,18,12,-6,-9,-4,-6,3,2,-18,6,-12,-9,4,3,-6,2,18,-6,12,-9,-4,3,-6,2,18,-6,-6,9,2,-3,-3,1,-18,6,6,9,-2,-3,-3,1],
[-6,-4,-3,0,-2,0,0,0,0,0,0,0,0,0,0,0,6,4,-3,0,-2,0,0,0,0,0,0,0,0,0,0,0,6,-2,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,-6,2,3,0,-1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,-6,0,-4,-3,-2,0,0,0,0,0,0,0,0,0,0,0,6,0,4,-3,-2,0,0,0,0,0,0,0,0,0,0,0,6,0,-2,3,-1,0,0,0,0,0,0,0,0,0,0,0,-6,0,2,3,-1,0,0,0,0,0,0,0,0],
[18,12,9,12,6,8,6,4,-18,-12,-9,6,-6,4,3,2,-18,-12,9,-12,6,-8,6,4,18,12,-9,-6,-6,-4,3,2,-18,6,-9,-12,3,4,-6,2,18,-6,9,-6,-3,2,-3,1,18,-6,-9,12,3,-4,-6,2,-18,6,9,6,-3,-2,-3,1],
[-12,-8,-6,-6,-4,-4,-3,-2,12,8,6,-6,4,-4,-3,-2,12,8,-6,6,-4,4,-3,-2,-12,-8,6,6,4,4,-3,-2,12,-4,6,6,-2,-2,3,-1,-12,4,-6,6,2,-2,3,-1,-12,4,6,-6,-2,2,3,-1,12,-4,-6,-6,2,2,3,-1],
[2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-6,-3,0,-4,0,-2,0,0,6,3,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,-3,0,4,0,-2,0,0,-6,3,0,2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[4,2,0,2,0,1,0,0,-4,-2,0,2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,2,0,-2,0,1,0,0,4,-2,0,-2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,-6,0,-3,0,-4,-2,0,0,6,0,3,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,-3,0,4,-2,0,0,-6,0,3,0,2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,4,0,2,0,2,1,0,0,-4,0,-2,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,2,0,-2,1,0,0,4,0,-2,0,-2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[-6,-3,-4,0,-2,0,0,0,0,0,0,0,0,0,0,0,6,3,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,6,-3,4,0,-2,0,0,0,0,0,0,0,0,0,0,0,-6,3,2,0,-1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,-6,0,-3,-4,-2,0,0,0,0,0,0,0,0,0,0,0,6,0,3,-2,-1,0,0,0,0,0,0,0,0,0,0,0,6,0,-3,4,-2,0,0,0,0,0,0,0,0,0,0,0,-6,0,3,2,-1,0,0,0,0,0,0,0,0],
[18,9,12,12,6,6,8,4,-18,-9,-12,6,-6,3,4,2,-18,-9,6,-12,3,-6,4,2,18,9,-6,-6,-3,-3,2,1,-18,9,-12,-12,6,6,-8,4,18,-9,12,-6,-6,3,-4,2,18,-9,-6,12,3,-6,-4,2,-18,9,6,6,-3,-3,-2,1],
[-12,-6,-8,-6,-4,-3,-4,-2,12,6,8,-6,4,-3,-4,-2,12,6,-4,6,-2,3,-2,-1,-12,-6,4,6,2,3,-2,-1,12,-6,8,6,-4,-3,4,-2,-12,6,-8,6,4,-3,4,-2,-12,6,4,-6,-2,3,2,-1,12,-6,-4,-6,2,3,2,-1],
[4,2,2,0,1,0,0,0,0,0,0,0,0,0,0,0,-4,-2,2,0,1,0,0,0,0,0,0,0,0,0,0,0,-4,2,-2,0,1,0,0,0,0,0,0,0,0,0,0,0,4,-2,-2,0,1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,4,0,2,2,1,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,2,1,0,0,0,0,0,0,0,0,0,0,0,-4,0,2,-2,1,0,0,0,0,0,0,0,0,0,0,0,4,0,-2,-2,1,0,0,0,0,0,0,0,0],
[-12,-6,-6,-8,-3,-4,-4,-2,12,6,6,-4,3,-2,-2,-1,12,6,-6,8,-3,4,-4,-2,-12,-6,6,4,3,2,-2,-1,12,-6,6,8,-3,-4,4,-2,-12,6,-6,4,3,-2,2,-1,-12,6,6,-8,-3,4,4,-2,12,-6,-6,-4,3,2,2,-1],
[8,4,4,4,2,2,2,1,-8,-4,-4,4,-2,2,2,1,-8,-4,4,-4,2,-2,2,1,8,4,-4,-4,-2,-2,2,1,-8,4,-4,-4,2,2,-2,1,8,-4,4,-4,-2,2,-2,1,8,-4,-4,4,2,-2,-2,1,-8,4,4,4,-2,-2,-2,1]],dtype=np.double)
b = get_bVec(xi,yi,zi,nx,ny,nz,xvec,yvec,zvec,m)
A_ijk = Binv.dot(b).ravel()
return A_ijk
class TriCubic(object):
def __init__(self,xvec=None,yvec=None,zvec=None,M=None,useCache = True,default=None,filename=None):
'''Object that handles tri cubic interpolation. In general use the following parameters:
``xvec`` - the xaxis (regular)
``yvec`` - the yaxis (regular)
``zvec`` - the zaxis (regular)
``M`` - the array of shape (len(xvec),len(yvec),len(zvec)) to be interpolated
``M`` will be transformed to a flat vector, member ``m``. This can be replaced and clearCache called.
``useCache`` - whether to save interpolants (keep yes for efficiency when interpolants used often).
``default`` - when point falls out of domain specified by xvec,yvec, zvec what value to return,
None (default) means nearest.
``pad`` - whether to pad array (experimental: leave False)'''
if filename is not None:
self.load(filename,useCache=useCache,default=default)
return
self.default = default
self.nx = np.size(xvec)
self.ny = np.size(yvec)
self.nz = np.size(zvec)
self.xvec = xvec.ravel(order='C')
self.yvec = yvec.ravel(order='C')
self.zvec = zvec.ravel(order='C')
self.m = M.ravel(order='C')
self.checkIndexing(M)
self.useCache = useCache
if self.useCache:
self.cache = {}
else:
self.cache = None
#@jit
def copy(self,**kwargs):
'''Return a copy of the TriCubic object by essentially creating a copy of all the data.
``kwargs`` are the same as constructor.'''
return TriCubic(self.xvec.copy(),self.yvec.copy(),self.zvec.copy(),self.getShapedArray().copy(),**kwargs)
#@jit
def load(self,filename,**kwargs):
f = open(filename,"rb")
dataDict = dill.load(f)
f.close()
self.__init__(dataDict['xvec'],dataDict['yvec'],dataDict['zvec'],dataDict['M'],**kwargs)
#@jit
def save(self,filename):
f = open(filename,"wb")
dill.dump({'xvec':self.xvec,'yvec':self.yvec,'zvec':self.zvec,'M':self.getShapedArray()},f)
f.close()
#@jit
def clearCache(self):
'''Clear the cache, which should be done if overwriting the array'''
self.cache = {}
#@jit
def interp(self,x,y,z,doDiff=False):
xi,yi,zi,A_ijk = self.getInterpolant(x,y,z)
if xi == -1:#outside entire array
xi = 0
if xi == self.nx:
xi = self.nx-2
if yi == -1:
yi = 0
if yi == self.ny:
yi = self.ny - 2
if zi == -1:
zi = 0
if zi == self.nz:
zi = self.nz - 2
if A_ijk is np.nan:#outside array or near edges
if self.default is not None:
f = self.default
else:
f = self.m[self.index(xi,yi,zi)]#nearest
if doDiff:
return f,0.,0.,0.,0.,0.,0.,0.
else:
return f
if doDiff:
return interp_diff(x,y,z,self.xvec,self.yvec,self.zvec,A_ijk)
else:
return interp_scalar(x,y,z,self.xvec,self.yvec,self.zvec,A_ijk)
#@jit(nopython=True)
def index(self,i,j,k):
'''Correct indexing of 3-tensor in ravelled vector, vectorized'''
return k + self.nz*(j + self.ny*i)
#@jit
def getShapedArray(self):
'''Return the model in 3d array with proper ij ordering'''
return self.m.reshape(self.nx,self.ny,self.nz,order='C')
#@jit
def getModelCoordinates(self):
X,Y,Z = np.meshgrid(self.xvec,self.yvec,self.zvec,indexing='ij')
return X.ravel(order='C'),Y.ravel(order='C'),Z.ravel(order='C')
#@jit
def checkIndexing(self,M,N=100):
'''Check that ordering of elements is correct'''
N = min(N,np.size(self.m))
idx = 0
while idx < N:
i = np.random.randint(self.nx)
j = np.random.randint(self.ny)
k = np.random.randint(self.nz)
assert self.m[self.index(i,j,k)] == M[i,j,k], "Ordering of indexing is wrong"
idx += 1
assert np.alltrue(self.getShapedArray() == M), "reshape is not right"
return True
#@jit
def getInterpolant(self,x,y,z):
'''vectorized build the interpolant'''
xi = bisection(self.xvec,x)
yi = bisection(self.yvec,y)
zi = bisection(self.zvec,z)
try:
ijk = self.index(xi,yi,zi)#bottom corner of cube
if self.useCache:
if ijk in self.cache.keys():
A_ijk = self.cache[ijk]
else:
A_ijk = computeInterpolant(xi,yi,zi,self.nx,self.ny,self.nz,self.xvec,self.yvec,self.zvec,self.m)
self.cache[ijk] = A_ijk
else:
A_ijk = computeInterpolant(xi,yi,zi,self.nx,self.ny,self.nz,self.xvec,self.yvec,self.zvec,self.m)
except:
A_ijk = np.nan
return xi,yi,zi,A_ijk
def testResult():
from TricubicInterpolation import TriCubic as TCI_current
xvec = np.arange(50)
yvec = np.arange(40)
zvec = np.arange(30)
M = np.random.uniform(size=[50,40,30])
tci_current = TCI_current(xvec,yvec,zvec,M)
tci = TriCubic(xvec,yvec,zvec,M)
for i in range(10):
x = np.random.uniform()*50
y = np.random.uniform()*40
z = np.random.uniform()*30
y1 = tci_current.interp(x,y,z)
y2 = tci.interp(x,y,z)
assert np.abs(y1 - y2) < 1e-10, "TCI's not equivalent good: {}, new: {}".format(y1,y2)
if __name__=='__main__':
np.random.seed(1234)
#generateBinv()
#optimizeBvecFormation()
testResult()
In [ ]:
In [ ]: