EHN: Implement cartesian_nearest_index

[oyamad]

Maybe used in #640. (Implement _find_closest_state_index part in #640 (comment).)

Given a point x, find the index of the point nearest to x among those in the Cartesian product generated nodes.

If nodes has n arrays of grid points (so the product is n dimensional) and there are m grid points in each dimension, this function does binary search for each dimension, so that the complexity is O(n log(m)), while if one actually constructs the Cartesian product and does linear search, then the complexity is O(m^n).

[coveralls]

Coverage increased (+0.03%) to 95.283% when pulling 7be79c0 on cartesian_nearest_index into fcde8c9 on master.

[jstac]

array([1, 4])

[jstac]

Thanks @oyamad !! This is clever.

In the following example I map (0.1, 0.1, 0.1) to the nearest grid point in R^3, which should be (0.0, 0.0, 0.0). I get a different result. Am I misunderstanding how to use the function?

[ins] In [14]: nodes = (np.linspace(0, 1, 5), np.linspace(0, 2, 3), np.linspace(0, 1, 4
          ...: ))

[ins] In [15]: nodes
Out[15]: 
(array([0.  , 0.25, 0.5 , 0.75, 1.  ]),
 array([0., 1., 2.]),
 array([0.        , 0.33333333, 0.66666667, 1.        ]))

[ins] In [16]: prod = qe.cartesian(nodes)

[ins] In [17]: qe.cartesian_nearest_index((0.1, 0.1, 0.1), prod)
Out[17]: 22

[ins] In [18]: prod[22]
Out[18]: array([0.25      , 2.        , 0.66666667])

[oyamad]

@jstac Thanks for your review!

You are supposed to pass nodes, not prod.

In [4]: qe.cartesian_nearest_index((0.1, 0.1, 0.1), nodes)
Out[4]: 0

In your call qe.cartesian_nearest_index((0.1, 0.1, 0.1), prod), it searches over the 60-dimensional Cartesian product that would be generated by the 60 arrays in prod, but only the first 3 dimensions (where len(x) = 3) are looked up. I will add lines to raise an error when the dimensions do not match as in this case.

[jstac]

Many thanks @oyamad , this looks good to me, and is greatly appreciated.

[jstac]

@oyamad or @mmcky please merge when ready.