[Feat] Ability to interpolate an array of coordinates

[jrobsontull]

Related Issue

#323

Summary

Ccp4Map.grid.interpolate_values expects the points to be organized on a regular grid. We often want to interpolate a numpy array of 3D coordinates and return the value at each point. I've added an example of what we would like to see added as an API, as well as testing and some documentation. I'm open to any ideas and I am definitely not settled on the naming of this function - the naming will become quite confusing as is...

Example Usage

import timeit
import gemmi
import numpy as np

dmap = gemmi.read_ccp4_map('map.mrc')
grid_arr = dmap.grid.array

coords = np.array([
    [138.462, 128.261, 122.029],
    [139.255, 127.177, 121.461],
    [138.643, 126.674, 120.158],
    [139.349, 126.461, 119.172],
    [139.386, 126.027, 122.46],
    [140.307, 126.324, 123.632],
    [140.408, 125.133, 124.571],
    [141.103, 124.009, 123.954],
    [141.374, 122.867, 124.571],
    [141.022, 122.659, 125.829]
], np.float32)

def gemmi_interp_value():
    return [dmap.grid.interpolate_value(gemmi.Position(*coord)) for coord in coords]

def gemmi_interp_coordinate_array():
    return dmap.grid.interpolate_coordinate_array(coords)

print('** Mean Interpolation Comparison **')
python_loop = gemmi_interp_value()
c_loop = gemmi_interp_coordinate_array()
print('gemmi_interp_value mean:', np.mean(python_loop)) 
print('interpolate_coordinate_array mean:', np.mean(c_loop))

print('\n** Speed Comparison **')
print(f'gemmi_interp_value: {timeit.timeit(gemmi_interp_value, number=100000):.3f} s')
print(f'gemmi_interp_coordinate_array: {timeit.timeit(gemmi_new, number=100000):.3f} s')

** Mean Interpolation Comparison **
gemmi_interp_value mean: 0.1443351425230503
interpolate_coordinate_array mean: 0.14433515

** Speed Comparison **
gemmi_interp_value: 2.333 s
gemmi_interp_coordinate_array: 0.110 s

Testing

I built gemmi with:

python3 -m venv pyenv
source pyenv/bin/activate
python3 -m pip install wheel
python3 -m pip wheel -v .
python3 -m pip install gemmi-*.whl

Unit tests are passing with:

python3 -m unittest discover -v -s tests/

Docs are building with:

python3 -m pip install sphinx
cd docs
sphinx-build -M doctest . _build -n -E

[jrobsontull]

A lot of the changes in this file are formatting changes - do we have auto formatting in this repo?

[wojdyr]

nope

[wojdyr]

Thanks for this PR. I've had several questions regarding grid interpolation over the last few years, and I've been meaning to re-think the API and address various requests for a long time. Let's get back to it next week. I'll be in a meeting most of this week, and I'm in the middle of moving python bindings from pybind11 to nanobind – I need to finish that first.

[wojdyr]

Sorry it takes so long. 5 weeks ago, when I wrote I'll get back to it next week, I was too optimistic. Switching from pybind11 to nanobind is still an ongoing process, in a separate branch.

For the record, I remembered that a few years ago I added another grid interpolating function for a friendly project:
7921286 But I didn't quite understand why it was designed in this way, and then it wasn't used, so I removed it. It's probably not relevant at this point.

I'll get back to this PR in a couple weeks.

[jrobsontull]

Thanks @wojdyr for coming back to this. We rely heavily on grid interpolation to do tasks and would welcome any improvements to this area!

[wojdyr]

Thanks again, it's now documented as interpolate_position_array.

BTW, my colleagues use some of the software you provide. Recently, I've been hearing about coordgenlibs.