RegularGridInterpolatorTest.php 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790
  1. <?php
  2. namespace MathPHP\Tests\NumericalAnalysis\Interpolation;
  3. use MathPHP\Exception\BadDataException;
  4. use MathPHP\NumericalAnalysis\Interpolation\RegularGridInterpolator;
  5. class RegularGridInterpolatorTest extends \PHPUnit\Framework\TestCase
  6. {
  7. /**
  8. * @test Interpolated regular grid function computes expected values: p(x) = expected
  9. * @dataProvider dataProviderForRegularGridAgrees
  10. * @param array $point
  11. * @param float $expected
  12. * @throws \Exception
  13. */
  14. public function testRegularGridAgrees(array $point, float $expected)
  15. {
  16. // Given
  17. [$points, $values] = $this->getSample4d();
  18. // And
  19. $p = new RegularGridInterpolator($points, $values);
  20. // When
  21. $evaluated = $p($point);
  22. // Then
  23. $this->assertEqualsWithDelta($expected, $evaluated, 0.00001);
  24. }
  25. /**
  26. * Test data based on SciPy test cases
  27. * https://github.com/scipy/scipy/blob/c734dacd61c5962a86ab3cc4bf2891fc94b720a6/scipy/interpolate/tests/test_interpolate.py#L2361
  28. * @return array (x, expected)
  29. */
  30. public function dataProviderForRegularGridAgrees(): array
  31. {
  32. return [
  33. // test_linear_xi3d
  34. [[0.1, 0.1, 1., .9], 1001.1],
  35. [[0.2, 0.1, .45, .8], 846.2],
  36. [[0.5, 0.5, .5, .5], 555.5],
  37. // test_linear_edges
  38. [[0., 0., 0., 0.], 0.],
  39. [[1., 1., 1., 1.], 1111.],
  40. ];
  41. }
  42. /**
  43. * @test Interpolated regular grid function computes expected values: p(x) = expected
  44. * @dataProvider dataProviderForRegularGridNearestAgrees
  45. * @param array $point
  46. * @throws \Exception
  47. */
  48. public function testRegularGridNearestAgrees(array $point, $expected)
  49. {
  50. // Given
  51. [$points, $values] = $this->getSample4d();
  52. // And
  53. $p = new RegularGridInterpolator($points, $values, RegularGridInterpolator::METHOD_NEAREST);
  54. // When
  55. $evaluated = $p($point);
  56. // Then
  57. $this->assertEquals($expected, $evaluated);
  58. }
  59. /**
  60. * @return array (x, expected)
  61. */
  62. public function dataProviderForRegularGridNearestAgrees(): array
  63. {
  64. return [
  65. [[0.1, 0.1, .9, .9], 1100],
  66. [[0.1, 0.1, 0.1, 0.1], 0.],
  67. [[0., 0., 0., 0.], 0.],
  68. [[1., 1., 1., 1.], 1111.],
  69. [[0.1, 0.4, 0.6, 0.9], 1055.]
  70. ];
  71. }
  72. /**
  73. * @test Example from Github issue 382
  74. * https://github.com/markrogoyski/math-php/issues/382
  75. *
  76. * >>> xs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
  77. * >>> ys = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
  78. * >>> zs = [110, 111, 112, 113, 114, 115, 116, 117, 118, 119];
  79. *
  80. * >>> def func(x, y, z):
  81. * ... return 2 * x + 3 * y - z
  82. *
  83. * >>> values = [[[func(x, y, z) for z in zs] for y in ys] for x in xs]
  84. *
  85. * >>> my_interpolating_function = RegularGridInterpolator((xs, ys, zs), values, method='linear')
  86. * >>> my_interpolating_function([2.21, 12.1, 115.9])
  87. * array([-75.18])
  88. */
  89. public function testIssue382ExampleLinear()
  90. {
  91. // Given
  92. $xs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
  93. $ys = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
  94. $zs = [110, 111, 112, 113, 114, 115, 116, 117, 118, 119];
  95. // And
  96. $func = function ($x, $y, $z) {
  97. return 2 * $x + 3 * $y - $z;
  98. };
  99. // And
  100. $data = [];
  101. foreach ($xs as $i => $x) {
  102. foreach ($ys as $j => $y) {
  103. foreach ($zs as $k => $z) {
  104. $data[$i][$j][$k] = $func($x, $y, $z);
  105. }
  106. }
  107. }
  108. // When
  109. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'linear');
  110. $result = $interp([2.21, 12.1, 115.9]);
  111. // Then
  112. $this->assertEqualsWithDelta(-75.18, $result, 0.00001);
  113. }
  114. /**
  115. * @test Example from Github issue 382
  116. * https://github.com/markrogoyski/math-php/issues/382
  117. *
  118. * >>> xs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
  119. * >>> ys = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
  120. * >>> zs = [110, 111, 112, 113, 114, 115, 116, 117, 118, 119];
  121. *
  122. * >>> def func(x, y, z):
  123. * ... return 2 * x + 3 * y - z
  124. *
  125. * >>> values = [[[func(x, y, z) for z in zs] for y in ys] for x in xs]
  126. *
  127. * >>> my_interpolating_function = RegularGridInterpolator((xs, ys, zs), values, method='nearest')
  128. * >>> my_interpolating_function([2.21, 12.1, 115.9])
  129. * array([-76.])
  130. */
  131. public function testIssue382ExampleNearest()
  132. {
  133. // Given
  134. $xs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
  135. $ys = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
  136. $zs = [110, 111, 112, 113, 114, 115, 116, 117, 118, 119];
  137. // And
  138. $func = function ($x, $y, $z) {
  139. return 2 * $x + 3 * $y - $z;
  140. };
  141. // And
  142. $data = [];
  143. foreach ($xs as $i => $x) {
  144. foreach ($ys as $j => $y) {
  145. foreach ($zs as $k => $z) {
  146. $data[$i][$j][$k] = $func($x, $y, $z);
  147. }
  148. }
  149. }
  150. // When
  151. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'nearest');
  152. $result = $interp([2.21, 12.1, 115.9]);
  153. // Then
  154. $this->assertEqualsWithDelta(-76, $result, 0.00001);
  155. }
  156. /**
  157. * @test SciPy documentation example 1
  158. * https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html#scipy.interpolate.RegularGridInterpolator
  159. *
  160. * from scipy.interpolate import RegularGridInterpolator
  161. * def f(x, y, z):
  162. * return 2 * x**3 + 3 * y**2 - z
  163. * x = np.linspace(1, 4, 11)
  164. * y = np.linspace(4, 7, 22)
  165. * z = np.linspace(7, 9, 33)
  166. * data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))
  167. * my_interpolating_function = RegularGridInterpolator((x, y, z), data)
  168. * pts = np.array([[2.1, 6.2, 8.3], [3.3, 5.2, 7.1]])
  169. * my_interpolating_function(pts)
  170. * array([ 125.80469388, 146.30069388])
  171. */
  172. public function testSciPyExample1()
  173. {
  174. // Given
  175. $xs = [1. , 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.7, 4.];
  176. $ys = [
  177. 4. , 4.14285714, 4.28571429, 4.42857143, 4.57142857,
  178. 4.71428571, 4.85714286, 5. , 5.14285714, 5.28571429,
  179. 5.42857143, 5.57142857, 5.71428571, 5.85714286, 6. ,
  180. 6.14285714, 6.28571429, 6.42857143, 6.57142857, 6.71428571,
  181. 6.85714286, 7.
  182. ];
  183. $zs = [
  184. 7. , 7.0625, 7.125 , 7.1875, 7.25 , 7.3125, 7.375 , 7.4375,
  185. 7.5 , 7.5625, 7.625 , 7.6875, 7.75 , 7.8125, 7.875 , 7.9375,
  186. 8. , 8.0625, 8.125 , 8.1875, 8.25 , 8.3125, 8.375 , 8.4375,
  187. 8.5 , 8.5625, 8.625 , 8.6875, 8.75 , 8.8125, 8.875 , 8.9375,
  188. 9.
  189. ];
  190. // And
  191. $func = function ($x, $y, $z) {
  192. return 2 * $x ** 3 + 3 * $y ** 2 - $z;
  193. };
  194. // And
  195. $data = [];
  196. foreach ($xs as $i => $x) {
  197. foreach ($ys as $j => $y) {
  198. foreach ($zs as $k => $z) {
  199. $data[$i][$j][$k] = $func($x, $y, $z);
  200. }
  201. }
  202. }
  203. // When
  204. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'linear');
  205. $result = $interp([2.1, 6.2, 8.3]);
  206. // Then
  207. $this->assertEqualsWithDelta(125.80469388, $result, 0.00001);
  208. }
  209. /**
  210. * @test SciPy documentation example 2
  211. * https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html#scipy.interpolate.RegularGridInterpolator
  212. *
  213. * from scipy.interpolate import RegularGridInterpolator
  214. * def f(x, y, z):
  215. * return 2 * x**3 + 3 * y**2 - z
  216. * x = np.linspace(1, 4, 11)
  217. * y = np.linspace(4, 7, 22)
  218. * z = np.linspace(7, 9, 33)
  219. * data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))
  220. * my_interpolating_function = RegularGridInterpolator((x, y, z), data)
  221. * pts = np.array([[3.3, 5.2, 7.1]])
  222. * my_interpolating_function(pts)
  223. * array([146.30069388])
  224. */
  225. public function testSciPyExample2()
  226. {
  227. // Given
  228. $xs = [1. , 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.7, 4.];
  229. $ys = [
  230. 4. , 4.14285714, 4.28571429, 4.42857143, 4.57142857,
  231. 4.71428571, 4.85714286, 5. , 5.14285714, 5.28571429,
  232. 5.42857143, 5.57142857, 5.71428571, 5.85714286, 6. ,
  233. 6.14285714, 6.28571429, 6.42857143, 6.57142857, 6.71428571,
  234. 6.85714286, 7.
  235. ];
  236. $zs = [
  237. 7. , 7.0625, 7.125 , 7.1875, 7.25 , 7.3125, 7.375 , 7.4375,
  238. 7.5 , 7.5625, 7.625 , 7.6875, 7.75 , 7.8125, 7.875 , 7.9375,
  239. 8. , 8.0625, 8.125 , 8.1875, 8.25 , 8.3125, 8.375 , 8.4375,
  240. 8.5 , 8.5625, 8.625 , 8.6875, 8.75 , 8.8125, 8.875 , 8.9375,
  241. 9.
  242. ];
  243. // And
  244. $func = function ($x, $y, $z) {
  245. return 2 * $x ** 3 + 3 * $y ** 2 - $z;
  246. };
  247. // And
  248. $data = [];
  249. foreach ($xs as $i => $x) {
  250. foreach ($ys as $j => $y) {
  251. foreach ($zs as $k => $z) {
  252. $data[$i][$j][$k] = $func($x, $y, $z);
  253. }
  254. }
  255. }
  256. // When
  257. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'linear');
  258. $result = $interp([3.3, 5.2, 7.1]);
  259. // Then
  260. $this->assertEqualsWithDelta(146.30069388, $result, 0.00001);
  261. }
  262. /**
  263. * @test Interpolated point values are outside the domain of the input data grid. Values outside the domain are extrapolated.
  264. * Linear method.
  265. * This test will hit the condition in the findIndices method where i > gridSize - 2.
  266. *
  267. * https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html#scipy.interpolate.RegularGridInterpolator
  268. *
  269. * from scipy.interpolate import RegularGridInterpolator
  270. * def f(x, y, z):
  271. * return 2 * x**3 + 3 * y**2 - z
  272. * x = np.linspace(1, 4, 11)
  273. * y = np.linspace(4, 7, 22)
  274. * z = np.linspace(7, 9, 33)
  275. * data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))
  276. * my_interpolating_function = RegularGridInterpolator((x, y, z), data, method='linear', bounds_error=False, fill_value=None)
  277. * pts = np.array([[3.3, 7.2, 7.1]]) # 7.2 is outside the bounds of the grid
  278. * my_interpolating_function(pts)
  279. * array([220.48028571])
  280. */
  281. public function testInterpolatedPointValuesOutsideDomainOfInputDataGridAreExtrapolatedLinear()
  282. {
  283. // Given
  284. $xs = [1. , 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.7, 4.];
  285. $ys = [
  286. 4. , 4.14285714, 4.28571429, 4.42857143, 4.57142857,
  287. 4.71428571, 4.85714286, 5. , 5.14285714, 5.28571429,
  288. 5.42857143, 5.57142857, 5.71428571, 5.85714286, 6. ,
  289. 6.14285714, 6.28571429, 6.42857143, 6.57142857, 6.71428571,
  290. 6.85714286, 7.
  291. ];
  292. $zs = [
  293. 7. , 7.0625, 7.125 , 7.1875, 7.25 , 7.3125, 7.375 , 7.4375,
  294. 7.5 , 7.5625, 7.625 , 7.6875, 7.75 , 7.8125, 7.875 , 7.9375,
  295. 8. , 8.0625, 8.125 , 8.1875, 8.25 , 8.3125, 8.375 , 8.4375,
  296. 8.5 , 8.5625, 8.625 , 8.6875, 8.75 , 8.8125, 8.875 , 8.9375,
  297. 9.
  298. ];
  299. // And
  300. $func = function ($x, $y, $z) {
  301. return 2 * $x ** 3 + 3 * $y ** 2 - $z;
  302. };
  303. // And
  304. $data = [];
  305. foreach ($xs as $i => $x) {
  306. foreach ($ys as $j => $y) {
  307. foreach ($zs as $k => $z) {
  308. $data[$i][$j][$k] = $func($x, $y, $z);
  309. }
  310. }
  311. }
  312. // When
  313. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'linear');
  314. $result = $interp([3.3, 7.2, 7.1]); // 7.2 is outside the bounds of the grid
  315. // Then
  316. $this->assertEqualsWithDelta(220.48028571, $result, 0.00001);
  317. }
  318. /**
  319. * @test Interpolated point values are outside the domain of the input data grid. Values outside the domain are extrapolated.
  320. * Nearest method.
  321. * This test will hit the condition in the findIndices method where i > gridSize - 2.
  322. *
  323. * https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html#scipy.interpolate.RegularGridInterpolator
  324. *
  325. * from scipy.interpolate import RegularGridInterpolator
  326. * def f(x, y, z):
  327. * return 2 * x**3 + 3 * y**2 - z
  328. * x = np.linspace(1, 4, 11)
  329. * y = np.linspace(4, 7, 22)
  330. * z = np.linspace(7, 9, 33)
  331. * data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))
  332. * my_interpolating_function = RegularGridInterpolator((x, y, z), data, method='nearest', bounds_error=False, fill_value=None)
  333. * pts = np.array([[3.3, 7.2, 7.1]]) # 7.2 is outside the bounds of the grid
  334. * my_interpolating_function(pts)
  335. * array([220.48028571])
  336. */
  337. public function testInterpolatedPointValuesOutsideDomainOfInputDataGridAreExtrapolatedNearest()
  338. {
  339. // Given
  340. $xs = [1. , 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.7, 4.];
  341. $ys = [
  342. 4. , 4.14285714, 4.28571429, 4.42857143, 4.57142857,
  343. 4.71428571, 4.85714286, 5. , 5.14285714, 5.28571429,
  344. 5.42857143, 5.57142857, 5.71428571, 5.85714286, 6. ,
  345. 6.14285714, 6.28571429, 6.42857143, 6.57142857, 6.71428571,
  346. 6.85714286, 7.
  347. ];
  348. $zs = [
  349. 7. , 7.0625, 7.125 , 7.1875, 7.25 , 7.3125, 7.375 , 7.4375,
  350. 7.5 , 7.5625, 7.625 , 7.6875, 7.75 , 7.8125, 7.875 , 7.9375,
  351. 8. , 8.0625, 8.125 , 8.1875, 8.25 , 8.3125, 8.375 , 8.4375,
  352. 8.5 , 8.5625, 8.625 , 8.6875, 8.75 , 8.8125, 8.875 , 8.9375,
  353. 9.
  354. ];
  355. // And
  356. $func = function ($x, $y, $z) {
  357. return 2 * $x ** 3 + 3 * $y ** 2 - $z;
  358. };
  359. // And
  360. $data = [];
  361. foreach ($xs as $i => $x) {
  362. foreach ($ys as $j => $y) {
  363. foreach ($zs as $k => $z) {
  364. $data[$i][$j][$k] = $func($x, $y, $z);
  365. }
  366. }
  367. }
  368. // When
  369. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'nearest');
  370. $result = $interp([3.3, 7.2, 7.1]); // 7.2 is outside the bounds of the grid
  371. // Then
  372. $this->assertEqualsWithDelta(218.483, $result, 0.00001);
  373. }
  374. /**
  375. * @test Similar values
  376. *
  377. * >>> xs = [1, 2, 3];
  378. * >>> ys = [1, 2, 3];
  379. * >>> zs = [1, 2, 3];
  380. *
  381. * >>> def func(x, y, z):
  382. * ... return 2 * x + 3 * y - z
  383. *
  384. * >>> values = [[[func(x, y, z) for z in zs] for y in ys] for x in xs]
  385. *
  386. * >>> my_interpolating_function = RegularGridInterpolator((xs, ys, zs), values, method='linear')
  387. * >>> my_interpolating_function(point)
  388. *
  389. * @dataProvider dataProviderForSimilarValues
  390. * @param array $point
  391. * @param float $expected
  392. */
  393. public function testSimilarValues(array $point, float $expected)
  394. {
  395. // Given
  396. $xs = [1, 2, 3];
  397. $ys = [1, 2, 3];
  398. $zs = [1, 2, 3];
  399. // And
  400. $func = function ($x, $y, $z) {
  401. return 2 * $x + 3 * $y - $z;
  402. };
  403. // And
  404. $data = [];
  405. foreach ($xs as $i => $x) {
  406. foreach ($ys as $j => $y) {
  407. foreach ($zs as $k => $z) {
  408. $data[$i][$j][$k] = $func($x, $y, $z);
  409. }
  410. }
  411. }
  412. // When
  413. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'linear');
  414. $result = $interp($point);
  415. // Then
  416. $this->assertEqualsWithDelta($expected, $result, 0.00001);
  417. }
  418. /**
  419. * @return array (point, expected)
  420. */
  421. public function dataProviderForSimilarValues(): array
  422. {
  423. return [
  424. [[1, 1, 1], 4.],
  425. [[1, 1, 2], 3.],
  426. [[1, 1, 3], 2.],
  427. [[1, 2, 1], 7.],
  428. [[1, 2, 2], 6.],
  429. [[1, 2, 3], 5.],
  430. [[1, 3, 1], 10.],
  431. [[1, 3, 2], 9.],
  432. [[1, 3, 3], 8.],
  433. [[2, 1, 1], 6.],
  434. [[2, 1, 2], 5.],
  435. [[2, 1, 3], 4.],
  436. [[2, 2, 1], 9.],
  437. [[2, 2, 2], 8.],
  438. [[2, 2, 3], 7.],
  439. [[2, 3, 1], 12.],
  440. [[2, 3, 2], 11.],
  441. [[2, 3, 3], 10.],
  442. [[3, 1, 1], 8.],
  443. [[3, 1, 2], 7.],
  444. [[3, 1, 3], 6.],
  445. [[3, 2, 1], 11.],
  446. [[3, 2, 2], 10.],
  447. [[3, 2, 3], 9.],
  448. [[3, 3, 1], 14.],
  449. [[3, 3, 2], 13.],
  450. [[3, 3, 3], 12.],
  451. ];
  452. }
  453. /**
  454. * @test Points and values not sorted
  455. *
  456. * >>> xs = [1, 2, 3]
  457. * >>> ys = [1, 2, 3]
  458. * >>> zs = [1, 2, 3]
  459. * >>>
  460. * >>> def func(x, y, z):
  461. * ... return 2 * x + 3 * y - z
  462. * ...
  463. * >>>
  464. * >>> values = [[[func(x, y, z) for z in [3, 2, 1]] for y in [3, 2, 1]] for x in [3, 2, 1]]
  465. * >>>
  466. * >>> my_interpolating_function = RegularGridInterpolator((xs, ys, zs), values, method='linear')
  467. * >>> my_interpolating_function([1, 1, 2])
  468. * array([13.])
  469. */
  470. public function testPointsAndValuesNotSorted()
  471. {
  472. // Given
  473. $xs = [1, 2, 3];
  474. $ys = [1, 2, 3];
  475. $zs = [1, 2, 3];
  476. // And
  477. $func = function ($x, $y, $z) {
  478. return 2 * $x + 3 * $y - $z;
  479. };
  480. // And
  481. $data = [];
  482. foreach ([3, 2, 1] as $i => $x) {
  483. foreach ([3, 2, 1] as $j => $y) {
  484. foreach ([3, 2, 1] as $k => $z) {
  485. $data[$i][$j][$k] = $func($x, $y, $z);
  486. }
  487. }
  488. }
  489. // When
  490. $interp = new RegularGridInterpolator([$xs, $ys, $zs], $data, 'linear');
  491. $result = $interp([1, 1, 2]);
  492. // Then
  493. $this->assertEqualsWithDelta(13, $result, 0.00001);
  494. }
  495. /**
  496. * @test Stack Overflow Example
  497. * https://stackoverflow.com/questions/30056577/correct-usage-of-scipy-interpolate-regulargridinterpolator
  498. *
  499. * from scipy.interpolate import RegularGridInterpolator
  500. * >>> x = np.linspace(0, 1, 3)
  501. * >>> x
  502. * array([0. , 0.5, 1. ])
  503. * >>> X, Y, Z = np.meshgrid(x, x, x, indexing='ij')
  504. * >>> X
  505. * array([[[0. , 0. , 0. ],
  506. * [0. , 0. , 0. ],
  507. * [0. , 0. , 0. ]],
  508. *
  509. * [[0.5, 0.5, 0.5],
  510. * [0.5, 0.5, 0.5],
  511. * [0.5, 0.5, 0.5]],
  512. *
  513. * [[1. , 1. , 1. ],
  514. * [1. , 1. , 1. ],
  515. * [1. , 1. , 1. ]]])
  516. * >>> Y
  517. * array([[[0. , 0. , 0. ],
  518. * [0.5, 0.5, 0.5],
  519. * [1. , 1. , 1. ]],
  520. *
  521. * [[0. , 0. , 0. ],
  522. * [0.5, 0.5, 0.5],
  523. * [1. , 1. , 1. ]],
  524. *
  525. * [[0. , 0. , 0. ],
  526. * [0.5, 0.5, 0.5],
  527. * [1. , 1. , 1. ]]])
  528. * >>> Z
  529. * array([[[0. , 0.5, 1. ],
  530. * [0. , 0.5, 1. ],
  531. * [0. , 0.5, 1. ]],
  532. *
  533. * [[0. , 0.5, 1. ],
  534. * [0. , 0.5, 1. ],
  535. * [0. , 0.5, 1. ]],
  536. *
  537. * [[0. , 0.5, 1. ],
  538. * [0. , 0.5, 1. ],
  539. * [0. , 0.5, 1. ]]])
  540. * >>>
  541. * >>>
  542. * >>> vals = np. sin(X) + np.\cos(Y) + np.tan(Z)
  543. * >>> vals
  544. * array([[[1. , 1.54630249, 2.55740772],
  545. * [0.87758256, 1.42388505, 2.43499029],
  546. * [0.54030231, 1.0866048 , 2.09771003]],
  547. *
  548. * [[1.47942554, 2.02572803, 3.03683326],
  549. * [1.3570081 , 1.90331059, 2.91441583],
  550. * [1.01972784, 1.56603033, 2.57713557]],
  551. *
  552. * [[1.84147098, 2.38777347, 3.39887871],
  553. * [1.71905355, 2.26535604, 3.27646127],
  554. * [1.38177329, 1.92807578, 2.93918102]]])
  555. *
  556. * >>> rgi = RegularGridInterpolator(points=[x, x, x], values=vals)
  557. * >>> tst = (0.47, 0.49, 0.53)
  558. * >>> rgi(tst)
  559. * array(1.93765972)
  560. */
  561. public function testStackOverflowExample()
  562. {
  563. // Given
  564. $x = [0., 0.5, 1.];
  565. // And
  566. $vals = [
  567. [
  568. [1. , 1.54630249, 2.55740772],
  569. [0.87758256, 1.42388505, 2.43499029],
  570. [0.54030231, 1.0866048 , 2.09771003]
  571. ],
  572. [
  573. [1.47942554, 2.02572803, 3.03683326],
  574. [1.3570081 , 1.90331059, 2.91441583],
  575. [1.01972784, 1.56603033, 2.57713557]],
  576. [
  577. [1.84147098, 2.38777347, 3.39887871],
  578. [1.71905355, 2.26535604, 3.27646127],
  579. [1.38177329, 1.92807578, 2.93918102]
  580. ]
  581. ];
  582. // And
  583. $tst = [0.47, 0.49, 0.53];
  584. // When
  585. $rgi = new RegularGridInterpolator([$x, $x, $x], $vals, 'linear');
  586. $result = $rgi($tst);
  587. // Then
  588. $this->assertEqualsWithDelta(1.93765972, $result, 0.00001);
  589. }
  590. /**
  591. * @test Bad method (not defined)
  592. * @throws BadDataException
  593. */
  594. public function testBadMethodException()
  595. {
  596. // Given
  597. $invalidMethod = 'methodDoesNotExist';
  598. // Then
  599. $this->expectException(BadDataException::class);
  600. // When
  601. $rgi = new RegularGridInterpolator([0], [0], $invalidMethod);
  602. }
  603. /**
  604. * @test Bad values - there are two point arrays, but values have one dimension
  605. * @throws BadDataException
  606. */
  607. public function testBadValuesException()
  608. {
  609. // Given
  610. $points = [[0], [1]];
  611. $values = [0];
  612. // Then
  613. $this->expectException(BadDataException::class);
  614. // When
  615. $rgi = new RegularGridInterpolator($points, $values);
  616. }
  617. /**
  618. * @test Bad pint dimension - The requested sample points xi have dimension 2, but this RegularGridInterpolator has dimension 1
  619. * @throws BadDataException
  620. */
  621. public function testInvokeBadPointDimensionException()
  622. {
  623. // Given
  624. $interp = new RegularGridInterpolator([0], [0]);
  625. // Then
  626. $this->expectException(BadDataException::class);
  627. // When
  628. $interp([0, 2]);
  629. }
  630. /**
  631. * Test fixture - Create a 4-D grid of 3 points in each dimension
  632. *
  633. * Based off of Python SciPy unit test fixture
  634. * def _get_sample_4d(self):
  635. * # create a 4-D grid of 3 points in each dimension
  636. * points = [(0., .5, 1.)] * 4
  637. * values = np.asarray([0., .5, 1.])
  638. * values0 = values[:, np.newaxis, np.newaxis, np.newaxis]
  639. * values1 = values[np.newaxis, :, np.newaxis, np.newaxis]
  640. * values2 = values[np.newaxis, np.newaxis, :, np.newaxis]
  641. * values3 = values[np.newaxis, np.newaxis, np.newaxis, :]
  642. * values = (values0 + values1 * 10 + values2 * 100 + values3 * 1000)
  643. * return points, values
  644. *
  645. * PHP code to generate data would look like:
  646. * $points = [[0.0, 0.5, 1.0], [0.0, 0.5, 1.0], [0.0, 0.5, 1.0], [0.0, 0.5, 1.0]];
  647. * $values = [];
  648. * $v = [0, 0.5, 1.];
  649. * for ($x = 0; $x < 3; $x++) {
  650. * for ($y = 0; $y < 3; $y++) {
  651. * for ($z = 0; $z < 3; $z++) {
  652. * for ($m = 0; $m < 3; $m++) {
  653. * $values[$x][$y][$z][$m] = $v[$x] + $v[$y] * 10 + $v[$z] * 100 + +$v[$m] * 1000;
  654. * }
  655. * }
  656. * }
  657. * }
  658. * return [$points, $values];
  659. *
  660. * @return array (points, values)
  661. */
  662. private function getSample4d(): array
  663. {
  664. $points = [[0.0, 0.5, 1.0], [0.0, 0.5, 1.0], [0.0, 0.5, 1.0], [0.0, 0.5, 1.0]];
  665. $values = [
  666. [
  667. [
  668. [0.0000e+00, 5.0000e+02, 1.0000e+03],
  669. [5.0000e+01, 5.5000e+02, 1.0500e+03],
  670. [1.0000e+02, 6.0000e+02, 1.1000e+03]
  671. ],
  672. [
  673. [5.0000e+00, 5.0500e+02, 1.0050e+03],
  674. [5.5000e+01, 5.5500e+02, 1.0550e+03],
  675. [1.0500e+02, 6.0500e+02, 1.1050e+03]
  676. ],
  677. [
  678. [1.0000e+01, 5.1000e+02, 1.0100e+03],
  679. [6.0000e+01, 5.6000e+02, 1.0600e+03],
  680. [1.1000e+02, 6.1000e+02, 1.1100e+03]
  681. ]
  682. ],
  683. [
  684. [
  685. [5.0000e-01, 5.0050e+02, 1.0005e+03],
  686. [5.0500e+01, 5.5050e+02, 1.0505e+03],
  687. [1.0050e+02, 6.0050e+02, 1.1005e+03]
  688. ],
  689. [
  690. [5.5000e+00, 5.0550e+02, 1.0055e+03],
  691. [5.5500e+01, 5.5550e+02, 1.0555e+03],
  692. [1.0550e+02, 6.0550e+02, 1.1055e+03]
  693. ],
  694. [
  695. [1.0500e+01, 5.1050e+02, 1.0105e+03],
  696. [6.0500e+01, 5.6050e+02, 1.0605e+03],
  697. [1.1050e+02, 6.1050e+02, 1.1105e+03]
  698. ]
  699. ],
  700. [
  701. [
  702. [1.0000e+00, 5.0100e+02, 1.0010e+03],
  703. [5.1000e+01, 5.5100e+02, 1.0510e+03],
  704. [1.0100e+02, 6.0100e+02, 1.1010e+03]
  705. ],
  706. [
  707. [6.0000e+00, 5.0600e+02, 1.0060e+03],
  708. [5.6000e+01, 5.5600e+02, 1.0560e+03],
  709. [1.0600e+02, 6.0600e+02, 1.1060e+03]
  710. ],
  711. [
  712. [1.1000e+01, 5.1100e+02, 1.0110e+03],
  713. [6.1000e+01, 5.6100e+02, 1.0610e+03],
  714. [1.1100e+02, 6.1100e+02, 1.1110e+03]
  715. ]
  716. ]
  717. ];
  718. return [$points, $values];
  719. }
  720. }