## N-Body problems

Many physical problems require the evaluation of all pairwise interactions of a large number of particles, so-called N-body problems. These problems arise in molecular dynamics, astrodynamics and electromagnetics among others. 

Their pairwise interactions can be expressed as:

\begin{equation}
f_i = \sum_{j=1}^n{P \left(\boldsymbol{x}_i, \boldsymbol{x}_j \right)w_j} \ \ \ \text{for } i=1,2,...,n 
\end{equation}

*  where subscripts $i$,  $j$ respectively denote *target* and *source*
*  $f_i$ can be a *potential* (or *force*) at target point $i$
*  $w_j$ is the *source weight* 
*  $\boldsymbol{x}_i, \boldsymbol{x}_j$ are the *spatial positions* of particles 
*  $P \left(\boldsymbol{x}_i, \boldsymbol{x}_j \right)$ is the *interaction kernel*. 

In order to evalute the potential $f_i$ at a target point $i$, we have to loop over each source particle $j$. Since there are $n$ target points $i$, this 'brute-force' approach costs $\mathcal{O} \left(n^2 \right)$ operations. 

One possible approach in this kind of problem is to define a few classes, say `Point` and `Particle` and then loop over the objects and perform the necessary point-to-point calculations.  


```python
import numpy
```


```python
class Point():
    """    
    Arguments:
        domain: the domain of random generated coordinates x,y,z, 
                default=1.0
    
    Attributes:
        x, y, z: coordinates of the point
    """
    def __init__(self, domain=1.0):
        self.x = domain * numpy.random.random()
        self.y = domain * numpy.random.random()
        self.z = domain * numpy.random.random()
            
    def distance(self, other):
        return ((self.x - other.x)**2 + 
                (self.y - other.y)**2 + 
                (self.z - other.z)**2)**.5
```


```python
class Particle(Point):
    """    
    Attributes:
        m: mass of the particle
        phi: the potential of the particle
    """
    
    def __init__(self, domain=1.0, m=1.0):
        Point.__init__(self, domain)
        self.m = m
        self.phi = 0.
```

Now we create a list of `n` random particles, define a function to calculate their interaction via direct summation and run!


```python
n = 1000
particles = [Particle(m = 1 / n) for i in range(n)]
```


```python
def direct_sum(particles):
    """
    Calculate the potential at each particle
    using direct summation method.

    Arguments:
        particles: the list of particles

    """
    for i, target in enumerate(particles):
        for source in (particles[:i] + particles[i+1:]):
            r = target.distance(source)
            target.phi += source.m / r
```


```python
direct_sum(particles)
```

There was a noticeable lag there.  How long does this thing take for 1000 particles?


```python
orig_time = %timeit -o direct_sum(particles)
```

    1 loop, best of 3: 3 s per loop



```python
%load_ext line_profiler
```


```python
%lprun -f direct_sum direct_sum(particles)
```

## How do we use Numba on this problem?

Problem: Numba doesn't support jitting native Python classes.  There is a `jit_class` structure in Numba but it's still in early development.

But it's nice to have attributes for literate programming.

Solution: NumPy custom dtypes.


```python
particle_dtype = numpy.dtype({'names':['x','y','z','m','phi'], 
                             'formats':[numpy.double, 
                                        numpy.double, 
                                        numpy.double, 
                                        numpy.double, 
                                        numpy.double]})
```


```python
myarray = numpy.ones(3, dtype=particle_dtype)
```


```python
myarray
```




    array([( 1.,  1.,  1.,  1.,  1.), ( 1.,  1.,  1.,  1.,  1.),
           ( 1.,  1.,  1.,  1.,  1.)], 
          dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8'), ('m', '<f8'), ('phi', '<f8')])



You can access an individual "attribute" like this:


```python
myarray[0]['x'] = 2.0
```


```python
myarray
```




    array([( 2.,  1.,  1.,  1.,  1.), ( 1.,  1.,  1.,  1.,  1.),
           ( 1.,  1.,  1.,  1.,  1.)], 
          dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8'), ('m', '<f8'), ('phi', '<f8')])



## [Exercise 1](./exercises/04.Direct.Summation.Exercises.ipynb#Exercise-1)

Write a function `create_n_random_particles` that takes the arguments `n` (number of particles), `m` (mass of every particle) and a domain within to generate a random number (as in the class above).
It should create an array with `n` elements and `dtype=particle_dtype` and then return that array.

For each particle, the mass should be initialized to the value of `m` and the potential `phi` initialized to zero.


```python
from numba import njit
```


```python
# %load snippets/nbody/create_n.py
```

**Note**: You can use "attribute" access on dtypes but there's a caveat.  If you need to debug this function without the decorator, you have to change them back to array access form.  


```python
parts = create_n_random_particles(1000, .001, 1)
```


```python
parts[:3]
```

We don't have a `distance` method anymore, so we need to write a function to take care of that.

## [Exercise 2](./exercises/04.Direct.Summation.Exercises.ipynb#Exercise-2)
Write a JITted function `distance` to calculate the distance between two particles of dtype `particle_dtype`


```python
# %load snippets/nbody/distance.py
```


```python
distance(parts[0], parts[1])
```


```python
%%timeit
distance(parts[0], parts[1])
```

## [Exercise 3](./exercises/04.Direct.Summation.Exercises.ipynb#Exercise-3)
Modify the `direct_sum` function above to instead work a NumPy array of particles.  Loop over each element in the array and calculate its total potential.


```python
# %load snippets/nbody/direct_sum.py
```


```python
numba_time = %timeit -o direct_sum(parts)
```


```python
orig_time.best / numba_time.best
```

257ddf542b83d8ff5ac0022989e668a878faf64f
ipynb
Jupyter Notebook
numba_tutorial_scipy2016-master/notebooks/04.Direct.Summation.ipynb
sunny2309/scipy_conf_notebooks
30a85d5137db95e01461ad21519bc1bdf294044b
MIT
2021-01-09T15:57:26.000Z
2021-11-29T01:44:21.000Z
numba_tutorial_scipy2016-master/notebooks/04.Direct.Summation.ipynb
sunny2309/scipy_conf_notebooks
30a85d5137db95e01461ad21519bc1bdf294044b
MIT
2019-11-15T02:00:26.000Z
2021-01-06T04:26:40.000Z
numba_tutorial_scipy2016-master/notebooks/04.Direct.Summation.ipynb
sunny2309/scipy_conf_notebooks
30a85d5137db95e01461ad21519bc1bdf294044b
MIT
Qwen/Qwen-72B
1. YES
2. YES