Periodic Boundary Conditions

[alvarozamora]

Any plans on implementing periodic boundary conditions?

[sdd]

None so far. I'd potentially be open to it as long as the implementation did not affect the performance when the feature was not in use, which it shouldn't if done correctly.

[alvarozamora]

I've been learning Rust and would like to work on adding this. I've been trying to parse through SciPy's implementation but I can't quite figure out how they implement periodic BCs. Any chance you could point me in the right direction for how to implement this efficiently?

They claim to offset the components by integer multiples of the boxsize, but when I try doing that manually I basically double the run time whereas using their periodic boundary functionality only increases the run time by O(≈30%)

[sdd]

Interesting. Can you share a link to their implementation? And possibly a code example where you see this 30% runtime penalty? I'm interested to take a look at their code to see how they do it in order to try the same thing in Rust.

[alvarozamora]

I actually found a bug in my code. Creating a larger tree by replicating the data 6 times seems to be as fast as their periodic implementation, so I think that's what they do. I think there's a tradeoff between choosing to

1) replicate half of the data per shift which comes with cost of finding points that are in a particular half three times, which may make building the tree slower than 2) but makes queries faster than 2) and

2) replicate the whole box 6 times, which may make building the tree faster than 1) but makes queries slower than 1)

Here is their implementation.

Forgive me for the ipynb, but here's what I did:

#!/usr/bin/env python
# coding: utf-8

# In[1]:

import numpy as np
from scipy.spatial import cKDTree as Tree

# In[2]:

Nd = 10**5; Nr = 10**6;
data = np.random.uniform(size=(Nd, 3))
rand = np.random.uniform(size=(Nr, 3))#/2 + 0.25

# In[3]:

# Time periodic
get_ipython().run_line_magic('timeit', '-n 1 -r 10 tree = Tree(data, leafsize=1, balanced_tree=True, compact_nodes=True, boxsize=1); r, ids = tree.query(rand)')

#2.82 s ± 323 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

# In[4]:

# Time nonperiodic
get_ipython().run_line_magic('timeit', '-n 1 -r 10 tree = Tree(data, leafsize=1, balanced_tree=True, compact_nodes=True); r, ids = tree.query(rand)')

#2.02 s ± 13.8 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

# In[5]:

def offset(data):
   return np.concatenate((data, 
                    data + np.array([1,0,0])[None,:],
                    data - np.array([1,0,0])[None,:],
                    data + np.array([0,1,0])[None,:],
                    data - np.array([0,1,0])[None,:],
                    data + np.array([0,0,1])[None,:],
                    data - np.array([0,0,1])[None,:],
                   ), axis=0)

# In[6]:

# time manual periodic
# this makes a larger tree by shifting data
get_ipython().run_line_magic('timeit', '-n 1 -r 10 tree = Tree(offset(data), leafsize=1, balanced_tree=True, compact_nodes=True); r, ids = tree.query(rand)')

# 2.66 s ± 58.5 ms per loop (mean ± std. dev. of 10 runs, 1 loop each)

# In[7]:

#time manual periodic alternate
# this uses the same tree and shifts queries
# much slower so only 3
get_ipython().run_line_magic('timeit', '-n 1 -r 3 tree = Tree(data, leafsize=1, balanced_tree=True, compact_nodes=True); r, ids = tree.query(offset(rand)); amn = r.reshape(7, Nr).argmin(axis=0); r = r.flatten()[amn]; ids = ids.flatten()[amn]')

# 22.6 s ± 176 ms per loop (mean ± std. dev. of 3 runs, 1 loop each)

[alvarozamora]

going to give implementing this a shot and if successful will submit a PR

[cavemanloverboy]

I implemented a version that uses only arrays. Will share tomorrow.

[sdd]

See https://github.com/sdd/kiddo/pull/35 for the most recent PR for this feature