<jupyter_start><jupyter_text>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.<jupyter_code>import numpy
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
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.<jupyter_output><empty_output><jupyter_text>Now we create a list of `n` random particles, define a function to calculate their interaction via direct summation and run!<jupyter_code>n = 1000
particles = [Particle(m = 1 / n) for i in range(n)]
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
direct_sum(particles)<jupyter_output><empty_output><jupyter_text>There was a noticeable lag there.  How long does this thing take for 1000 particles?<jupyter_code>orig_time = %timeit -o direct_sum(particles)
%load_ext line_profiler
%lprun -f direct_sum direct_sum(particles)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>particle_dtype = numpy.dtype({'names':['x','y','z','m','phi'], 
                             'formats':[numpy.double, 
                                        numpy.double, 
                                        numpy.double, 
                                        numpy.double, 
                                        numpy.double]})
myarray = numpy.ones(3, dtype=particle_dtype)
myarray<jupyter_output><empty_output><jupyter_text>You can access an individual "attribute" like this:<jupyter_code>myarray[0]['x'] = 2.0
myarray<jupyter_output><empty_output><jupyter_text>[Exercise 1](./exercises/04.Direct.Summation.Exercises.ipynbExercise-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.<jupyter_code>from numba import njit
# %load snippets/nbody/create_n.py<jupyter_output><empty_output><jupyter_text>**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.<jupyter_code>parts = create_n_random_particles(1000, .001, 1)
parts[:3]<jupyter_output><empty_output><jupyter_text>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.ipynbExercise-2)Write a JITted function `distance` to calculate the distance between two particles of dtype `particle_dtype`<jupyter_code># %load snippets/nbody/distance.py
distance(parts[0], parts[1])
%%timeit
distance(parts[0], parts[1])<jupyter_output><empty_output><jupyter_text>[Exercise 3](./exercises/04.Direct.Summation.Exercises.ipynbExercise-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.<jupyter_code># %load snippets/nbody/direct_sum.py
numba_time = %timeit -o direct_sum(parts)
orig_time.best / numba_time.best<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10602496833725349, 'fraction_numerical': 0.010132078885471322, 'mean_word_length': 3.6259414225941424, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '142315', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 1717, 'n_tokens_neox': 1654, 'n_words': 542, 'path': 'notebooks/04.Direct.Summation.ipynb', 'repo_name': 'IsabelAverill/Scipy-2017---Numba', 'size': '5527'}