Repopulate MDANSE

MDANSE/DistributedComputing/Slave.py

+"""This modules contains the handlers functions used by Pyro slaves to perform the job.
+
+Functions:
+    * do_run_step: performs the analysis step by step.
+"""
+
+ # Define (or import) all the task handlers.
+def do_run_step(job, step):
+             
+    return job.run_step(step)
+
+from nMOLDYN.DistributedComputing.MasterSlave import startSlaveProcess
+startSlaveProcess(master_host="localhost:%d")
+    
+

MDANSE/Mathematics/Arithmetic.py

+import cmath
+import itertools
+
+import numpy
+
+def factorial(n):
+    """Returns n!
+    
+    @param n: n.
+    @type: integer
+    
+    @return: n!.
+    @rtype: integer
+    """
+    
+    # 0! = 1! = 1
+    if n < 2:
+        return 1
+    else:
+        # multiply is a ufunc of Numeric module
+        return numpy.multiply.reduce(numpy.arange(2,n+1, dtype=numpy.float64))
+
+def pgcd(n):
+    """Computes the pgcd for a set of integers.
+    
+    @param n: n.
+    @type: list of integers
+    
+    @return: pgcd([i1,i2,i3...]).
+    @rtype: integer
+    """   
+    
+    def _pgcd(a,b):
+        while b: a, b = b, a%b
+        return a
+    
+    p = _pgcd(n[0],n[1])
+    
+    for x in n[2:]:
+        p = _pgcd(p, x)
+        
+    return p
+
+def get_weights(props, contents, dim):
+                                        
+    normFactor = None
+
+    weights = {}
+                
+    cartesianProduct = set(itertools.product(props.keys(), repeat=dim))
+    for elements in cartesianProduct:
+    
+        n = numpy.product([contents[el] for el in elements])        
+        p = numpy.product(numpy.array([props[el] for el in elements]),axis=0)
+                
+        fact = n*p
+
+        weights[elements] = numpy.float64(numpy.copy(fact))
+                
+        if normFactor is None:
+            normFactor = numpy.abs(fact)
+        else:
+            normFactor += numpy.abs(fact)
+        
+    for k in weights.keys():
+        weights[k] /= numpy.float64(normFactor)
+       
+    return weights, normFactor
+
+def weight(props,values,contents,dim,key,symmetric=True):
+    weights = get_weights(props,contents,dim)[0]
+    weightedSum = None      
+    matches = dict([(key%k,k) for k in weights.keys()])
+            
+    for k,val in values.iteritems():
+        
+        if not matches.has_key(k):
+            continue
+                    
+        if symmetric:
+            permutations = set(itertools.permutations(matches[k], r=dim))
+            w = sum([weights.pop(p) for p in permutations])
+        else:
+            w = weights.pop(matches[k])
+                                
+        if weightedSum is None:
+            weightedSum = w*val
+        else:
+            weightedSum += w*val                
+                                
+    return weightedSum
+
+class ComplexNumber(complex):
+    
+    def argument(self):
+        return cmath.phase(self)
+        
+    def modulus(self):
+        return numpy.sqrt(self*self.conjugate()).real
+                
\ No newline at end of file

MDANSE/Mathematics/Geometry.py

+import numpy
+
+from Scientific.Geometry import Vector
+
+def get_basis_vectors_from_cell_parameters(parameters):
+    """Returns the basis vectors for the simulation cell from the six crystallographic parameters.
+    
+    @param parameters: the a, b, c, alpha, bete and gamma of the simulation cell.
+    @type: parameters: list of 6 floats
+    
+    @return: a list of three Scientific.Geometry.Vector objects representing respectively a, b and c 
+        basis vectors.
+    @rtype: list
+    """
+
+    # The simulation cell parameters.
+    a, b, c, alpha, beta, gamma = parameters
+    
+    # By construction the a vector is aligned with the x axis.
+    e1 = Vector(a, 0.0, 0.0)
+
+    # By construction the b vector is in the xy plane.
+    e2 = b*Vector(numpy.cos(gamma), numpy.sin(gamma), 0.0)
+    
+    e3_x = numpy.cos(beta)
+    e3_y = (numpy.cos(alpha) - numpy.cos(beta)*numpy.cos(gamma)) / numpy.sin(gamma)
+    e3_z = numpy.sqrt(1.0 - e3_x**2 - e3_y**2)
+    e3 = c*Vector(e3_x, e3_y, e3_z)
+
+    return (e1, e2, e3)
+
+        
+def center_of_mass(coords, masses=None):
+    """Computes the center of massfor a set of coordinates and masses 
+    :Parameters:
+        #. coords  ((n,3)-numpy.array): the input coordinates
+        #. masses  ((n, )-numpy.array): the input masses
+    :Returns:
+        #. (3,)-numpy.array: the center of mass
+    """
+        
+    return numpy.average(coords,weights=masses,axis=0)
+    
+center = center_of_mass
+    
+def build_cartesian_axes(origin, p1, p2, dtype = numpy.float64):
+    
+    origin = numpy.array(origin, dtype = dtype)
+    p1 = numpy.array(p1, dtype = dtype)
+    p2 = numpy.array(p2, dtype = dtype)
+    
+    
+    v1 = p1 - origin
+    v2 = p2 - origin
+    
+    n1 = (v1 + v2)
+    n1 /= numpy.linalg.norm(n1)
+    
+    n3 = numpy.cross(v1,n1)
+    n3 /= numpy.linalg.norm(n3)
+    
+    n2 = numpy.cross(n3,n1)
+    n2 /= numpy.linalg.norm(n2)
+    
+    return n1,n2,n3
+    
+def generate_sphere_points(n):
+    """Returns list of 3d coordinates of points on a sphere using the
+    Golden Section Spiral algorithm.
+    """
+    
+    points = []
+    
+    inc = numpy.pi * (3 - numpy.sqrt(5))
+    
+    offset = 2 / float(n)
+    
+    for k in range(int(n)):
+        y = k * offset - 1 + (offset / 2)
+        r = numpy.sqrt(1 - y*y)
+        phi = k * inc
+        points.append([numpy.cos(phi)*r, y, numpy.sin(phi)*r])
+    
+    return points    
+
+def random_points_on_sphere(radius=1.0, nPoints=100):
+
+    points = numpy.zeros((3,nPoints),dtype=numpy.float)
+    
+    theta = 2.0*numpy.pi*numpy.random.uniform(nPoints)
+    u = numpy.random.uniform(-1.0,1.0,nPoints)
+    points[0,:] = radius * numpy.sqrt(1 - u**2) * numpy.cos(theta)
+    points[1,:] = radius * numpy.sqrt(1 - u**2) * numpy.sin(theta)
+    points[2,:] = radius * u
+    
+    return points
+
+def random_points_on_disk(axis, radius=1.0, nPoints=100):
+
+    axis = Vector(axis).normal().array
+    
+    points = numpy.random.uniform(-radius,radius,3*nPoints)
+    points = points.reshape((3,nPoints))
+
+    proj = numpy.dot(axis,points)
+    proj = numpy.dot(axis[:,numpy.newaxis],proj[numpy.newaxis,:])
+    
+    points -= proj
+                        
+    return points    
+
+def random_points_on_circle(axis, radius=1.0, nPoints=100):
+
+    axis = Vector(axis).normal().array
+    
+    points = numpy.random.uniform(-radius,radius,3*nPoints)
+    points = points.reshape((3,nPoints))
+
+    proj = numpy.dot(axis,points)
+    proj = numpy.dot(axis[:,numpy.newaxis],proj[numpy.newaxis,:])
+    
+    points -= proj
+    
+    points *= (radius/numpy.sqrt(numpy.sum(points**2,axis=0)))
+                        
+    return points    
+
+
+
+    
+    
+    
+    
\ No newline at end of file

MDANSE/Mathematics/Signal.py

+import math
+
+import numpy
+
+from MDANSE.Core.Error import Error
+
+class SignalError(Error):
+    pass
+
+INTERPOLATION_ORDER = []
+
+INTERPOLATION_ORDER.append(numpy.array([[-3., 4.,-1.],
+                                        [-1., 0., 1.],
+                                        [ 1.,-4., 3.]], dtype=numpy.float64))
+
+
+INTERPOLATION_ORDER.append(numpy.array([[-3., 4.,-1.],
+                                        [-1., 0., 1.],
+                                        [ 1.,-4., 3.]], dtype=numpy.float64))
+
+INTERPOLATION_ORDER.append(numpy.array([[-11., 18., -9., 2.],
+                                        [ -2., -3.,  6.,-1.],
+                                        [  1., -6.,  3., 2.],
+                                        [ -2.,  9.,-18.,11.]], dtype=numpy.float64))
+
+INTERPOLATION_ORDER.append(numpy.array([[-50., 96.,-72., 32.,-6.],
+                                        [ -6.,-20., 36.,-12., 2.],
+                                        [  2.,-16.,  0., 16.,-2.],
+                                        [ -2., 12.,-36., 20., 6.],
+                                        [  6.,-32., 72.,-96.,50.]], dtype=numpy.float64))
+
+INTERPOLATION_ORDER.append(numpy.array([[-274., 600.,-600., 400.,-150., 24.],
+                                        [ -24.,-130., 240.,-120.,  40., -6.],
+                                        [   6., -60., -40., 120., -30.,  4.],
+                                        [  -4.,  30.,-120.,  40.,  60., -6.],
+                                        [   6., -40., 120.,-240., 130., 24.],
+                                        [ -24., 150.,-400., 600.,-600.,274.]], dtype=numpy.float64))
+
+
+def correlation(x, y=None, axis=0, sumOverAxis=None, average=None):
+    """Returns the numerical correlation between two signals.
+
+    @param inputSeries1: the first signal.
+    @type inputSeries1: NumPy array   
+    
+    @param inputSeries2: if not None, the second signal otherwise the correlation will be an autocorrelation.
+    @type inputSeries2: NumPy array or None
+    
+    @return: the result of the numerical correlation.
+    @rtype: NumPy array
+
+    @note: if |inputSeries1| is a multidimensional array the correlation calculation is performed on
+    the first dimension.
+
+    @note: The correlation is computed using the FCA algorithm.
+    """
+    
+    x = numpy.array(x)
+    
+    n = x.shape[axis]
+            
+    X = numpy.fft.fft(x, 2*n,axis=axis)
+
+    if y is not None:
+        y = numpy.array(y)
+        Y = numpy.fft.fft(y, 2*n,axis=axis)
+    else:
+        Y = X
+    
+    s = [slice(None)]*x.ndim
+
+    s[axis] = slice(0,x.shape[axis],1)
+    
+    corr = numpy.real(numpy.fft.ifft(numpy.conjugate(X)*Y,axis=axis)[s])
+            
+    norm = n - numpy.arange(n)
+        
+    s = [numpy.newaxis]*x.ndim
+    s[axis] = slice(None)
+        
+    corr = corr/norm[s]
+                    
+    if sumOverAxis is not None:
+        corr = numpy.sum(corr,axis=sumOverAxis)
+    elif average is not None:
+        corr = numpy.average(corr,axis=average)
+            
+    return corr
+
+def normalize(x, axis=0):
+
+    s = [slice(None)]*x.ndim
+    s[axis] = slice(0,1,1)
+    
+    nx = x/x[s]
+    return nx
+    
+
+def differentiate(a, dt=1.0, order=0):
+        
+    try:
+        coefs = INTERPOLATION_ORDER[order]
+    except TypeError:
+        raise SignalError("Invalid type for differentiation order")
+    except IndexError:
+        raise SignalError("Invalid differentiation order")
+        
+    # outputSeries is the output resulting from the differentiation
+    ts = numpy.zeros(a.shape, dtype=numpy.float)
+        
+    if order == 0:
+                        
+        ts[0] = numpy.add.reduce(coefs[0,:]*a[:3])
+        ts[-1] = numpy.add.reduce(coefs[2,:]*a[-3:])
+
+        gj = a[1:] - a[:-1]
+        ts[1:-1] = (gj[1:]+gj[:-1])
+
+    # Case of the order 2
+    elif order == 1:
+
+        ts[0]  = numpy.add.reduce(coefs[0,:]*a[:3])
+        ts[-1] = numpy.add.reduce(coefs[2,:]*a[-3:])
+
+        gj      = numpy.zeros((a.size-2,3), dtype=numpy.float)
+        gj[:,0] = coefs[1,0]*a[:-2]
+        gj[:,1] = coefs[1,1]*a[1:-1]
+        gj[:,2] = coefs[1,2]*a[2:]
+        ts[1:-1] = numpy.add.reduce(gj,-1)
+
+    # Case of the order 3
+    elif order == 2:
+
+        # Special case for the first and last elements
+        ts[0]  = numpy.add.reduce(coefs[0,:]*a[:4])
+        ts[1]  = numpy.add.reduce(coefs[1,:]*a[:4])
+        ts[-1] = numpy.add.reduce(coefs[3,:]*a[-4:])
+
+        # General case
+        gj      = numpy.zeros((a.size-3,4), dtype=numpy.float)
+        gj[:,0] = coefs[2,0]*a[:-3]
+        gj[:,1] = coefs[2,1]*a[1:-2]
+        gj[:,2] = coefs[2,2]*a[2:-1]
+        gj[:,3] = coefs[2,3]*a[3:]
+        ts[2:-1] = numpy.add.reduce(gj,-1)
+
+    # Case of the order 4
+    elif order == 3:
+
+        # Special case for the first and last elements
+        ts[0]  = numpy.add.reduce(coefs[0,:]*a[:5])
+        ts[1]  = numpy.add.reduce(coefs[1,:]*a[:5])
+        ts[-2] = numpy.add.reduce(coefs[3,:]*a[-5:])
+        ts[-1] = numpy.add.reduce(coefs[4,:]*a[-5:])
+
+        # General case
+        gj      = numpy.zeros((a.size-4,5), dtype=numpy.float)
+        gj[:,0] = coefs[2,0]*a[:-4]
+        gj[:,1] = coefs[2,1]*a[1:-3]
+        gj[:,2] = coefs[2,2]*a[2:-2]
+        gj[:,3] = coefs[2,3]*a[3:-1]
+        gj[:,4] = coefs[2,4]*a[4:]
+        ts[2:-2] = numpy.add.reduce(gj,-1)
+
+    # Case of the order 5
+    elif order == 4:
+
+        # Special case for the first and last elements
+        ts[0]  = numpy.add.reduce(coefs[0,:]*a[:6])
+        ts[1]  = numpy.add.reduce(coefs[1,:]*a[:6])
+        ts[2]  = numpy.add.reduce(coefs[2,:]*a[:6])
+        ts[-2] = numpy.add.reduce(coefs[4,:]*a[-6:])
+        ts[-1] = numpy.add.reduce(coefs[5,:]*a[-6:])
+
+        # General case
+        gj      = numpy.zeros((a.size-5,6), dtype=numpy.float)
+        gj[:,0] = coefs[3,0]*a[:-5]
+        gj[:,1] = coefs[3,1]*a[1:-4]
+        gj[:,2] = coefs[3,2]*a[2:-3]
+        gj[:,3] = coefs[3,3]*a[3:-2]
+        gj[:,4] = coefs[3,4]*a[4:-1]
+        gj[:,5] = coefs[3,5]*a[5:]
+        ts[3:-2] = numpy.add.reduce(gj,-1)
+
+
+    fact = 1.0/(dt*math.factorial(max(order+1,2)))
+    
+    ts *= fact
+        
+    return ts        
+
+def symmetrize(signal, axis=0):
+    """Return a symmetrized version of an input signal
+
+    :Parameters:
+        #. signal (numpy.array): the input signal
+        #. axis (int): the axis along which the signal should be symmetrized
+    :Returns:
+        #. numpy.array: the symmetrized signal
+    """
+    
+    s = [slice(None)]*signal.ndim
+    s[axis] = slice(-1,0,-1)
+            
+    signal = numpy.concatenate((signal[s],signal),axis=axis)
+    
+    return signal
+
+def get_spectrum(signal,window=None,timeStep=1.0,axis=0):
+                            
+    signal = symmetrize(signal,axis)
+
+    if window is None:
+        window = numpy.ones(signal.shape[axis])
+
+    s = [numpy.newaxis]*signal.ndim
+    s[axis] = slice(None)
+    
+    return numpy.fft.fftshift(numpy.abs(numpy.fft.fft(signal*window[s],axis=axis)),axes=axis)*timeStep    
+    
+    

MDANSE/MolecularDynamics/Analysis.py

+'''
+MDANSE : Molecular Dynamics Analysis for Neutron Scattering Experiments
+------------------------------------------------------------------------------------------
+Copyright (C)
+2015- Eric C. Pellegrini Institut Laue-Langevin
+BP 156
+6, rue Jules Horowitz
+38042 Grenoble Cedex 9
+France
+pellegrini[at]ill.fr
+
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with this library; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ 
+Created on Mar 23, 2015
+
+@author: pellegrini
+'''
+
+import numpy
+
+from MDANSE.Core.Error import Error
+from MDANSE.Mathematics.Geometry import center_of_mass
+from MDANSE.Mathematics.Signal import correlation
+
+class AnalysisError(Error):
+    pass
+                                              
+def mean_square_deviation(coords1, coords2, masses=None, root=False):
+    """Computes the mean square deviation between two sets of coordinates 
+    :Parameters:
+        #. coords1 ((n,3)-numpy.array): the first set of coordinates
+        #. coords2 ((n,3)-numpy.array): the second set of coordinates
+        #. masses  ((n, )-numpy.array): the input masses
+        #. root    bool: if True, return the square root of the mean square deviation (the so-called RMSD)
+    :Returns:
+        #. float: the radius of gyration
+    """
+
+    if coords1.shape != coords2.shape:
+        raise AnalysisError("The input coordinates shapes do not match")
+
+    if masses is None:
+        masses = numpy.ones((coords1.shape[0]),dtype=numpy.float64)
+
+    rmsd = numpy.sum(numpy.sum((coords1-coords2)**2,axis=1)*masses)/numpy.sum(masses)
+
+    if root:
+        rmsd = numpy.sqrt(rmsd)
+    
+    return rmsd
+
+def mean_square_displacement(coordinates):
+    
+    dsq = numpy.add.reduce(coordinates**2,1)
+
+    # sum_dsq1 is the cumulative sum of dsq
+    sum_dsq1 = numpy.add.accumulate(dsq)
+
+    # sum_dsq1 is the reversed cumulative sum of dsq
+    sum_dsq2 = numpy.add.accumulate(dsq[::-1])
+
+    # sumsq refers to SUMSQ in the published algorithm
+    sumsq = 2.0*sum_dsq1[-1]
+
+    # this line refers to the instruction SUMSQ <-- SUMSQ - DSQ(m-1) - DSQ(N - m) of the published algorithm
+    # In this case, msd is an array because the instruction is computed for each m ranging from 0 to len(traj) - 1
+    # So, this single instruction is performing the loop in the published algorithm
+    Saabb  = sumsq - numpy.concatenate(([0.0], sum_dsq1[:-1])) - numpy.concatenate(([0.0], sum_dsq2[:-1]))
+
+    # Saabb refers to SAA+BB/(N-m) in the published algorithm
+    # Sab refers to SAB(m)/(N-m) in the published algorithm
+    Saabb = Saabb / (len(dsq) - numpy.arange(len(dsq)))
+    Sab   = 2.0*correlation(coordinates, axis=0, reduce=1)
+
+    # The atomic MSD.
+    msd = Saabb - Sab
+                            
+    return msd
+
+def mean_square_fluctuation(coords, root=False):
+    
+    msf = numpy.average(numpy.sum((coords-numpy.average(coords,axis=0))**2,axis=1))
+    
+    if root:
+        msf = numpy.sqrt(msf)
+        
+    return msf    
+
+def radius_of_gyration(coords, masses=None, root=False):
+    """Computes the radius of gyration for a set of coordinates
+    :Parameters:
+        #. coords  ((n,3)-numpy.array): the set of  coordinates
+        #. masses  ((n, )-numpy.array): the input masses
+        #. root    bool: if True, return the square root of the radius of gyration
+    :Returns:
+        #. float: the radius of gyration
+    """
+    
+    if masses is None:
+        masses = numpy.ones((coords.shape[0]),dtype=numpy.float64)