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