Package nMOLDYN :: Package Core :: Module Mathematics |
[hide private]
[frames] | no frames]

# Source Code for Module nMOLDYN.Core.Mathematics

      1  """This modules implements the mathematics-related classes, functions and procedures.
      2
      3  Classes:
      4      * QVectors: the class that actually performs the q vectors generation.
      5
      6  Functions:
      7      * differentiate      : performs a numerical differentiation of 1D Numeric array
      8      * correlation        : performs the numerical correlation between two 1D Numeric arrays
      9      * convolution        : performs the numerical convolution between two 1D Numeric arrays
     10      * FFT                : performs the FFT of a 1D Numeric array
     11      * invFFT             : performs the inverse FFT of a 1D Numeric array
     12      * gaussianWindow     : performs a Gaussian smoothing of 1D Numeric array.
     13      * factorial          : computes factorial (n) where n is an integer.
     14      * basisVectors       : computes the basis vectors of the simulation cell from a set of values defining its geometry (3 distances and 3 angles).
     15      * randomPointInCircle: returns a vector within a circle of radius |r| and orthogonal to a given direction.
     16      * randomDirection2D  : returns a normalized vector generated from a unit circle orthogonal to a given direction.
     17      * randomPlane2D      : generates a normalized random q-vector on a plane defined by vect1, vect2.
     18      * qVectorGenerator   : sets up and returns a set of q vectors generated from different user-defined parameters.
     19      * sphericalHarmonics : calculates the spherical functions Y from a set of j, m, n Wigner indexes
     20      * preparePP          : sets up the calculation of spherical harmonics.
     21  """
     22
     23  # The python distribution modules
     24  from random import uniform
     25  import sys
     26
     27  # The ScientificPython modules
     28  from Scientific import N
     29  from Scientific.FFT import fft, inverse_fft
     30  from Scientific.Geometry import Tensor, Vector
     31
     32  # The MMTK distribution modules
     33  from MMTK.Random import randomDirection
     34
     35  # The nMOLDYN modules
     36  from nMOLDYN.Core.Error import Error
     37  from nMOLDYN.Core.Logger import LogMessage
     38
     39  ###############################################################################
     40  # Analysis
     41  ###############################################################################
     42
     43  """a2 = array used to perform order 2 numerical differentiation scheme."""
     44  a2 = N.array([[   -3.,   4.,  -1.],
     45                [   -1.,   0.,   1.],
     46                [    1.,  -4.,   3.]])
     47
     48  """a3 = array used to perform order 3 numerical differentiation scheme."""
     49  a3 = N.array([[  -11.,  18.,  -9.,   2.],
     50                [   -2.,  -3.,   6.,  -1.],
     51                [    1.,  -6.,   3.,   2.],
     52                [   -2.,   9., -18.,  11.]])
     53
     54  """a4 = array used to perform order 4 numerical differentiation scheme."""
     55  a4 = N.array([[  -50.,  96., -72.,  32.,  -6.],
     56                [   -6., -20.,  36., -12.,   2.],
     57                [    2., -16.,   0.,  16.,  -2.],
     58                [   -2.,  12., -36.,  20.,   6.],
     59                [    6., -32.,  72., -96.,  50.]])
     60
     61  """a5 = N.array used to perform order 5 numerical differentiation scheme."""
     62  a5 = N.array([[  -274.,  600., -600.,  400., -150.,  24.],
     63                [   -24., -130.,  240., -120.,   40.,  -6.],
     64                [     6.,  -60.,  -40.,  120.,  -30.,   4.],
     65                [    -4.,   30., -120.,   40.,   60.,  -6.],
     66                [     6.,  -40.,  120., -240.,  130.,  24.],
     67                [   -24.,  150., -400.,  600., -600., 274.]])
     68
     69 -def differentiate(inputSeries, order, dx):
     70      """Returns the numerical derivative of order |order| of the signal |inputSeries| using the
     71      differentiation step |dx|.
     72
     73      @param inputSeries: the signal to differentiate.
     74      @index inputSeries: NumPy array.
     75
     76      @param order: an integer in [1,5] specifying the numerical differentiation order.
     77      @index order: integer
     78
     79      @param dx: a float specifying the differentiation step. Assumed to be constant over all the spectrum.
     80      @index dx: float
     81
     82      @return: the differentiated signal.
     83      @rtype: NumPy array
     84
     85      @see: M. Abramowitz, I.A. Stegun; 'Handbook of mathematical functions', Dover, New-York, 1972 p.914.
     86      """
     87
     88      if len(inputSeries) == 0:
     89          raise Error('The inputSeries is empty.')
     90
     91      if dx <= 0.0:
     92          raise Error('The time step for derivation must strictly positive.')
     93
     94      # outputSeries is the output resulting from the differentiation
     95      outputSeries = N.zeros(len(inputSeries),'d')
     96
     97      if order == 1:
     98          # Classical order 1
     99          # The formula used here is the double-sided differentiation where
    100          # x'(t) = (x(t+1) - x(t-1))/2dt
    101          # this formula triggers two special indexes cases 0 and len(x) - 1
    102
    103          # An intermediate factor representing the denominator of the derivation formula
    104          fact = 1.0/(2.0*dx)
    105
    106          # array of (x(ti+1) - x(ti)) with i ranging from 0 to len(x) - 2
    107          gj = inputSeries[1:] - inputSeries[:-1]
    108
    109          # Special case for the first and last elements
    110          outputSeries[0] = N.add.reduce(a2[0,:]*inputSeries[:3])*fact
    111          outputSeries[-1] = N.add.reduce(a2[2,:]*inputSeries[-3:])*fact
    112
    113          # gj[1:] = array of (x(ti+1) - x(ti)) with i ranging from 1 to len(x) - 2
    114          # gj[:-1] = array of (x(ti+1) - x(ti)) with i ranging from 0 to len(x) - 3
    115          # This is the double-sided differentiation formula
    116          outputSeries[1:-1] = (gj[1:]+gj[:-1])*fact
    117
    118      elif order == 2:
    119          # Case of the order 2
    120
    121          # An intermediate factor representing the denominator of the derivation formula
    122          fact = 1.0/(2.0*dx)
    123
    124          # Special case for the first and last elements
    125          outputSeries[0]  = N.add.reduce(a2[0,:]*inputSeries[:3])*fact
    126          outputSeries[-1] = N.add.reduce(a2[2,:]*inputSeries[-3:])*fact
    127
    128          # General case
    129          gj      = N.zeros((len(inputSeries)-2,3),'d')
    130          gj[:,0] = a2[1,0]*inputSeries[:-2]
    131          gj[:,1] = a2[1,1]*inputSeries[1:-1]
    132          gj[:,2] = a2[1,2]*inputSeries[2:]
    133          outputSeries[1:-1] = N.add.reduce(gj,-1)*fact
    134
    135      elif order == 3:
    136          # Case of the order 3
    137
    138          # An intermediate factor representing the denominator of the derivation formula
    139          fact = 1./(6.*dx)
    140
    141          # Special case for the first and last elements
    142          outputSeries[0]  = N.add.reduce(a3[0,:]*inputSeries[:4])*fact
    143          outputSeries[1]  = N.add.reduce(a3[1,:]*inputSeries[:4])*fact
    144          outputSeries[-1] = N.add.reduce(a3[3,:]*inputSeries[-4:])*fact
    145
    146          # General case
    147          gj      = N.zeros((len(inputSeries)-3,4),'d')
    148          gj[:,0] = a3[2,0]*inputSeries[:-3]
    149          gj[:,1] = a3[2,1]*inputSeries[1:-2]
    150          gj[:,2] = a3[2,2]*inputSeries[2:-1]
    151          gj[:,3] = a3[2,3]*inputSeries[3:]
    152          outputSeries[2:-1] = N.add.reduce(gj,-1)*fact
    153
    154      elif order == 4:
    155          # Case of the order 4
    156
    157          # An intermediate factor representing the denominator of the derivation formula
    158          fact = 1./(24.*dx)
    159
    160          # Special case for the first and last elements
    161          outputSeries[0]  = N.add.reduce(a4[0,:]*inputSeries[:5])*fact
    162          outputSeries[1]  = N.add.reduce(a4[1,:]*inputSeries[:5])*fact
    163          outputSeries[-2] = N.add.reduce(a4[3,:]*inputSeries[-5:])*fact
    164          outputSeries[-1] = N.add.reduce(a4[4,:]*inputSeries[-5:])*fact
    165
    166          # General case
    167          gj      = N.zeros((len(inputSeries)-4,5),'d')
    168          gj[:,0] = a4[2,0]*inputSeries[:-4]
    169          gj[:,1] = a4[2,1]*inputSeries[1:-3]
    170          gj[:,2] = a4[2,2]*inputSeries[2:-2]
    171          gj[:,3] = a4[2,3]*inputSeries[3:-1]
    172          gj[:,4] = a4[2,4]*inputSeries[4:]
    173          outputSeries[2:-2] = N.add.reduce(gj,-1)*fact
    174
    175      elif order == 5:
    176          # Case of the order 5
    177
    178          # An intermediate factor representing the denominator of the derivation formula
    179          fact = 1./(120.*dx)
    180
    181          # Special case for the first and last elements
    182          outputSeries[0]  = N.add.reduce(a5[0,:]*inputSeries[:6])*fact
    183          outputSeries[1]  = N.add.reduce(a5[1,:]*inputSeries[:6])*fact
    184          outputSeries[2]  = N.add.reduce(a5[2,:]*inputSeries[:6])*fact
    185          outputSeries[-2] = N.add.reduce(a5[4,:]*inputSeries[-6:])*fact
    186          outputSeries[-1] = N.add.reduce(a5[5,:]*inputSeries[-6:])*fact
    187
    188          # General case
    189          gj      = N.zeros((len(inputSeries)-5,6),'d')
    190          gj[:,0] = a5[3,0]*inputSeries[:-5]
    191          gj[:,1] = a5[3,1]*inputSeries[1:-4]
    192          gj[:,2] = a5[3,2]*inputSeries[2:-3]
    193          gj[:,3] = a5[3,3]*inputSeries[3:-2]
    194          gj[:,4] = a5[3,4]*inputSeries[4:-1]
    195          gj[:,5] = a5[3,5]*inputSeries[5:]
    196          outputSeries[3:-2] = N.add.reduce(gj,-1)*fact
    197
    198      else:
    199          raise Error('Unknown differentiation order.')
    200
    201      return outputSeries
    202
    203  # Just for testing purpose.
    204 -def myautocorrelation(inputSeries):
    205
    206      nTimes = len(inputSeries)
    207      corr = N.zeros((nTimes,), N.Float)
    208      for i in range(nTimes):
    209          s = 0.0
    210          for j in range(nTimes-i):
    211              try:
    212                  s += sum(inputSeries[j]*inputSeries[j+i])
    213              except:
    214                  s += inputSeries[j]*inputSeries[j+i]
    215          corr[i] = s/float(nTimes-i)
    216
    217      return corr
    218
    219 -def correlation(inputSeries1, inputSeries2 = None):
    220      """Returns the numerical correlation between |inputSeries1| and |inputSeries2| multidimensional
    221      NumPy arrays.
    222
    223      @param inputSeries1: the first signal
    224      @type inputSeries1: NumPy array
    225
    226      @param inputSeries2: if not None, the second signal to correlate with |inputSeries1| otherwise the
    227          correlation will be an autocorrelation.
    228      @type inputSeries2: NumPy array
    229
    230      @return: an array (length(|inputSeries1|)) storing the result of the correlation.
    231      @rtype: NumPy array
    232
    233      @note: if |inputSeries1| is a multidimensional array the correlation calculation is performed on
    234      the first dimension.
    235
    236      @note: The correlation is computed using the FCA algorithm.
    237      """
    238
    239      if len(inputSeries1) <= 0:
    240          raise Error('One or both time series are empty.')
    241
    242      # The length of inputSeries1 is stored in inputSeries1Length
    243      inputSeries1Length = len(inputSeries1)
    244
    245      # extendedLength = 2*len(inputSeries1)
    246      extendedLength = 2*inputSeries1Length
    247
    248      # The FCA algorithm:
    249
    250      # 1) computation of the FFT of inputSeries1 zero-padded until extendedLength
    251      # The computation is done along the 0-axis
    252      FFTSeries1 = fft(inputSeries1,extendedLength,0)
    253
    254      if inputSeries2 is None:
    255              # Autocorrelation case
    256          FFTSeries2 = FFTSeries1
    257      else:
    258          # 2) computation of the FFT of inputSeries2 zero-padded until extendedLength
    259          # The computation is  done along the 0-axis
    260          FFTSeries2 = fft(inputSeries2,extendedLength,0)
    261
    262      # 3) Product between FFT(inputSeries1)* and FFT(inputSeries2)
    263      FFTSeries1 = N.conjugate(FFTSeries1)*FFTSeries2
    264
    265      # 4) inverse FFT of the product
    266      # The computation is done along the 0-axis
    267      FFTSeries1 = inverse_fft(FFTSeries1,len(FFTSeries1),0)
    268
    269      # This refers to (1/(N-m))*Sab in the published algorithm.
    270      # This is the correlation function defined for positive indexes only.
    271      if len(FFTSeries1.shape) == 1:
    272          corr = FFTSeries1.real[:inputSeries1Length] / (inputSeries1Length-N.arange(inputSeries1Length))
    273      else:
    274          corr = N.add.reduce(FFTSeries1.real[:inputSeries1Length],1) / (inputSeries1Length-N.arange(inputSeries1Length))
    275
    276      return corr
    277
    278 -def convolution(inputSeries1, inputSeries2):
    279      """Returns the numerical convolution between |inputSeries1| and |inputSeries2| one-dimensional
    280      NumPy arrays.
    281
    282      @param inputSeries1: the first signal
    283      @type inputSeries1: NumPy array
    284
    285      @param inputSeries2: the second signal to convolve with |inputSeries1|.
    286      @type inputSeries2: NumPy array
    287
    288      @return: an array (length(|inputSeries1|)) storing the result of the convolution.
    289      @rtype: NumPy array
    290
    291      @note: if |inputSeries1| is a multidimensional array the convolution calculation is performed on
    292      the first dimension.
    293
    294      @note: the convolution is computed using the convolve function of NumPy package.
    295      """
    296
    297      if (len(inputSeries1) == 0) | (len(inputSeries2) == 0):
    298          raise Error('One or both inputSeries are empty.')
    299
    300      # The convolution between inputSeries1 and inputSeries2.
    301      # A 1D Numeric array of dimension len(inputSeries1)*len(inputSeries2)-1
    302      return N.convolve(inputSeries1,inputSeries2)
    303
    304 -def FFT(inputSeries):
    305      """Returns the FFT of |inputSeries| multidimensional NumPy array.
    306
    307      @param inputSeries: the array on which to computes the FFT.
    308      @type inputSeries: NumPy array
    309
    310      @return: the FFT transformed array.
    311      @rtype: NumPy array
    312
    313      @note: the FFT is computed using the fft function of Scientific.FFT package.
    314      """
    315
    316      if len(inputSeries) == 0:
    317          raise Error('The inputSeries is empty.')
    318
    319      # The FFT of inputSeries
    320      return fft(inputSeries,len(inputSeries),0)
    321
    322 -def invFFT(inputSeries):
    323      """Returns the inverse FFT of |inputSeries| multidimensional NumPy array.
    324
    325      @param inputSeries: the array on which to computes the inverse FFT.
    326      @type inputSeries: NumPy array
    327
    328      @return: the inverse FFT transformed array.
    329      @rtype: NumPy array
    330
    331      @note: the inverse FFT is computed using the inverse_fft function of Scientific.FFT package.
    332      """
    333
    334      if len(inputSeries) == 0:
    335          raise Error('The inputSeries is empty.')
    336
    337      # The inverse FFT of inputSeries
    338      return inverse_fft(inputSeries,len(inputSeries),0)
    339
    340 -def gaussianWindow(inputSeries, alpha):
    341      """Returns a smoothed signal using |inputSeries| input signal and a gaussian kernel of width |alpha|.
    342
    343      @param inputSeries: the signal to smooth.
    344      @type inputSeries: NumPy array
    345
    346      @param alpha: a float specifying the width of the Gaussian.
    347      @type alpha: float
    348
    349      @return: an array (length = 2*len(|inputSeries|) - 1) containing the smoothed signal.
    350      @rtype: NumPy array
    351      """
    352
    353      # outputSeries is an array of length 2*len(series)-1 directionized to zero
    354      # outputSeries is the smoothed version of inputSeries obtained by applying
    355      # a gaussian kernel to intputSeries
    356      outputSeries = N.zeros((2*len(inputSeries) - 2,), N.Float)
    357
    358      # exp(...) is the gaussian window used to smooth the spectrum series
    359      res = inputSeries*N.exp(-0.5*(alpha*N.arange(len(inputSeries))/(len(inputSeries)-1))**2)
    360
    361      # The second half of outputSeries is filled using periodic conditions
    362      outputSeries[:len(inputSeries)] = res
    363      outputSeries[len(inputSeries):] = res[-2:0:-1]
    364
    365      return outputSeries
    366
    367 -def factorial(n):
    368      """Returns n!
    369
    370      @param n: the n of n!.
    371      @type: integer
    372
    373      @return: n!.
    374      @rtype: integer
    375      """
    376      if n < 2:
    377          return 1
    378      else:
    379          # multiply is a ufunc of Numeric module
    380          return N.multiply.reduce(N.arange(2,n+1))
    381
    382 -def preparePP(j, m, n):
    383      """Intermediate function used to setup the calculation of spherical harmonics."""
    384
    385      # c2 is the sqrt[(j+m)! * (j-m)!] * j! of equation 3.55
    386      c2 = N.sqrt((factorial(j+m) * factorial(j-m))) * factorial(j)
    387
    388      c3 = j+m
    389      pp = []
    390      aa = []
    391      for p in range(0,j+1,1):
    392          if (p <= c3) and (p >= m):
    393              # p1 is the j+m-p of equation 3.55
    394              p1 = c3-p
    395              # p2 is the j-p of equation 3.55
    396              p2 = j-p
    397              p3 = p-m
    398
    399              # a1 is the first part of equation 3.55 before the x
    400              a1 = ((-1)**p)*c2/(factorial(p1)*factorial(p2)*factorial(p)*factorial(p3))
    401
    402              pp.append((p1,p2,p3,p))
    403
    404              aa.append(a1)
    405
    406      return N.array(pp), N.array(aa)
    407
    408 -def sphericalCoordinates(x, y, z):
    409      """This function returns the r, theta and phi spherical coordinates from the x, y z cartesian coordinates.
    410
    411      @param x: the cartesian x.
    412      @type x: float
    413
    414      @param y: the cartesian y.
    415      @type y: float
    416
    417      @param z: the cartesian z.
    418      @type z: float
    419
    420      @return: the r, theta and phi spherical coordinates..
    421      @rtype: a list of three floats
    422      """
    423
    424      # The spherical radius is computed
    425      r = N.sqrt(x**2 + y**2 + z**2)
    426
    427      # The spherical theta is computed
    428      theta = N.arccos(z/r)
    429
    430      # The spherical phi is computed
    431      phi = N.arctan2(y,x)
    432
    433      return r, theta, phi
    434
    435 -def changeBasis(pt, op, ip, jp, kp):
    436      """This function return the coordinates
    437
    438      @param pt: the coordinates of the point in the old basis.
    439      @type pt: Scientific Vector
    440
    441      @param op: the coordinates of the new origin in the old basis.
    442      @type op: Scientific Vector
    443
    444      @param ip: the coordinates of the new x axis in the old basis.
    445      @type ip: Scientific Vector
    446
    447      @param jp: the coordinates of the new y axis in the old basis.
    448      @type jp: Scientific Vector
    449
    450      @param kp: the coordinates of the new z axis in the old basis.
    451      @type kp: Scientific Vector
    452
    453      @return: the coordinates of the point in the new basis.
    454      @rtype: Scientific Vector
    455      """
    456
    457      mat = Tensor(N.array([ip,jp,kp]).transpose())
    458      matinv = mat.inverse()
    459
    460      mprim = matinv * (pt - op)
    461
    462      return mprim
    463
    464  ###############################################################################
    465  # Geometry
    466  ###############################################################################
    467
    468 -def basisVectors(parameters):
    469      """Returns the basis vectors for the simulation cell from the six crystallographic parameters.
    470
    471      @param parameters: a list of six floats defining the simulation cell geometry.
    472      @type: parameters: list
    473
    474      @return: a list of three Scientific.Geometry.Vector objects representing respectively a, b and c
    475          basis vectors.
    476      @rtype: list
    477      """
    478
    479      a, b, c, alpha, beta, gamma = parameters
    480      e1 = Vector(a, 0.0, 0.0)
    481      e2 = b*Vector(N.cos(gamma), N.sin(gamma), 0.0)
    482      e3_x = N.cos(beta)
    483      e3_y = (N.cos(alpha) - N.cos(beta)*N.cos(gamma)) / N.sin(gamma)
    484      e3_z = N.sqrt(1.0 - e3_x**2 - e3_y**2)
    485      e3 = c*Vector(e3_x, e3_y, e3_z)
    486
    487      return (e1, e2, e3)
    488
    489 -def randomPointInCircle(r,axis):
    490      """Returns a vector drawn from an uniform distribution within a circle of radius |r| and
    491      orthogonal to vector |axis|.
    492
    493      @param r: float specifying the radius of the circle.
    494      @type r: float
    495
    496      @param axis: the axis orthogonal to the plane where the vectors have to be generated.
    497      @type axis: Scientific.Geometry.Vector object
    498
    499      @return: a vector pointing to a random point of the circle.
    500      @rtype: Scientific.Geometry.Vector object
    501      """
    502      rsq = r*r
    503      while 1:
    504          x = N.array([uniform(-r, r), uniform(-r,r), uniform(-r,r)])
    505          y = Vector(x) - (Vector(x)*axis/axis.length()**2)*axis
    506          if y*y < rsq :
    507              break
    508      return y
    509
    510 -def randomDirection2D(axis):
    511      """Returns a normalized vector drawn from an uniform distribution on the surface of a unit circle on a
    512      plane orthogonal to |axis|.
    513
    514      @param axis: the axis orthogonal to the plane where the vectors have to be generated.
    515      @type axis: Scientific.Geometry.Vector object
    516
    517      @return: A normalized vector defined in a unit disk orthogonal to |axis|
    518      @rtype: Scientific.Geometry.Vector object
    519      """
    520      r = randomPointInCircle(1.0,axis)
    521      return r.normal()
    522
    523 -def randomVector(directions = None):
    524      """Returns a normalized random vector on a plane or in space.
    525
    526      @param directions: if not None, a list of 2 Scientific.Vector that will define the plane on which
    527          the vector should be generated.
    528      @type directions: list of 2 Scientific.Vector or None
    529
    530      @return: a normalized random vector on a plane defined by |directions| or in space (|directions| = None).
    531      @rtype: Scientific.Geometry.Vector object
    532      """
    533
    534      if directions is None:
    535          return randomDirection()
    536
    537      if len(directions) == 1:
    538          return directions[0]
    539
    540      elif len(directions) == 2:
    541          # The first vector of the plane.
    542          v1, v2 = directions
    543
    544          # v1 /\ v2.
    545          axis = (v1.cross(v2))
    546
    547          if axis.length() == 0.0:
    548              raise Error('Your basis vectors are colinear. Can not build a plane out of them.')
    549
    550          return randomDirection2D(axis)
    551