source: trunk/GSASIImath.py @ 1056

Last change on this file since 1056 was 1056, checked in by vondreele, 8 years ago

restore multiprocessing MCSA

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 98.1 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-09-17 19:45:04 +0000 (Tue, 17 Sep 2013) $
5# $Author: vondreele $
6# $Revision: 1056 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1056 2013-09-17 19:45:04Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIImath: computation module*
12================================
13
14Routines for least-squares minimization and other stuff
15
16'''
17import sys
18import os
19import os.path as ospath
20import random as rn
21import numpy as np
22import numpy.linalg as nl
23import numpy.ma as ma
24import cPickle
25import time
26import math
27import copy
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 1056 $")
30import GSASIIElem as G2el
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import GSASIIpwd as G2pwd
34import numpy.fft as fft
35import pypowder as pyd
36
37sind = lambda x: np.sin(x*np.pi/180.)
38cosd = lambda x: np.cos(x*np.pi/180.)
39tand = lambda x: np.tan(x*np.pi/180.)
40asind = lambda x: 180.*np.arcsin(x)/np.pi
41acosd = lambda x: 180.*np.arccos(x)/np.pi
42atand = lambda x: 180.*np.arctan(x)/np.pi
43atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
44
45def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,Print=False):
46   
47    """
48    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
49    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`
50   
51    ::
52   
53                    Nobs
54        x = arg min(sum(func(y)**2,axis=0))
55                    y=0
56
57    :param function func: callable method or function
58        should take at least one (possibly length N vector) argument and
59        returns M floating point numbers.
60    :param np.ndarray x0: The starting estimate for the minimization of length N
61    :param function Hess: callable method or function
62        A required function or method to compute the weighted vector and Hessian for func.
63        It must be a symmetric NxN array
64    :param tuple args: Any extra arguments to func are placed in this tuple.
65    :param float ftol: Relative error desired in the sum of squares.
66    :param float xtol: Relative error desired in the approximate solution.
67    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
68        until other limits are met (ftol, xtol)
69    :param bool Print: True for printing results (residuals & times) by cycle
70
71    :returns: (x,cov_x,infodict) where
72
73      * x : ndarray
74        The solution (or the result of the last iteration for an unsuccessful
75        call).
76      * cov_x : ndarray
77        Uses the fjac and ipvt optional outputs to construct an
78        estimate of the jacobian around the solution.  ``None`` if a
79        singular matrix encountered (indicates very flat curvature in
80        some direction).  This matrix must be multiplied by the
81        residual standard deviation to get the covariance of the
82        parameter estimates -- see curve_fit.
83      * infodict : dict
84        a dictionary of optional outputs with the keys:
85
86         * 'fvec' : the function evaluated at the output
87         * 'num cyc':
88         * 'nfev':
89         * 'lamMax':
90         * 'psing':
91           
92    """
93               
94    x0 = np.array(x0, ndmin=1)      #might be redundant?
95    n = len(x0)
96    if type(args) != type(()):
97        args = (args,)
98       
99    icycle = 0
100    One = np.ones((n,n))
101    lam = 0.001
102    lamMax = lam
103    nfev = 0
104    if Print:
105        print ' Hessian refinement on %d variables:'%(n)
106    while icycle < maxcyc:
107        time0 = time.time()
108        M = func(x0,*args)
109        nfev += 1
110        chisq0 = np.sum(M**2)
111        Yvec,Amat = Hess(x0,*args)
112        Adiag = np.sqrt(np.diag(Amat))
113        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
114        if np.any(psing):                #hard singularity in matrix
115            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
116        Anorm = np.outer(Adiag,Adiag)
117        Yvec /= Adiag
118        Amat /= Anorm
119        while True:
120            Lam = np.eye(Amat.shape[0])*lam
121            Amatlam = Amat*(One+Lam)
122            try:
123                Xvec = nl.solve(Amatlam,Yvec)
124            except nl.LinAlgError:
125                print 'ouch #1'
126                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
127                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
128            Xvec /= Adiag
129            M2 = func(x0+Xvec,*args)
130            nfev += 1
131            chisq1 = np.sum(M2**2)
132            if chisq1 > chisq0:
133                lam *= 10.
134            else:
135                x0 += Xvec
136                lam /= 10.
137                break
138            if lam > 10.e5:
139                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
140                break
141        lamMax = max(lamMax,lam)
142        if Print:
143            print ' Cycle: %d, Time: %.2fs, Chi**2: %.3g, Lambda: %.3g'%(icycle,time.time()-time0,chisq1,lam)
144        if (chisq0-chisq1)/chisq0 < ftol:
145            break
146        icycle += 1
147    M = func(x0,*args)
148    nfev += 1
149    Yvec,Amat = Hess(x0,*args)
150    Amatlam = Amat          #*(One+Lam)/Anorm
151    try:
152        Bmat = nl.inv(Amatlam)          #*(One+Lam)/Anorm
153        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
154    except nl.LinAlgError:
155        print 'ouch #2 linear algebra error in LS'
156        psing = []
157        if maxcyc:
158            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
159        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
160   
161def getVCov(varyNames,varyList,covMatrix):
162    '''obtain variance-covariance terms for a set of variables. NB: the varyList
163    and covMatrix were saved by the last least squares refinement so they must match.
164   
165    :param list varyNames: variable names to find v-cov matric for
166    :param list varyList: full list of all variables in v-cov matrix
167    :param nparray covMatrix: full variance-covariance matrix from the last
168     least squares refinement
169   
170    :returns: nparray vcov: variance-covariance matrix for the variables given
171     in varyNames
172   
173    '''
174    vcov = np.zeros((len(varyNames),len(varyNames)))
175    for i1,name1 in enumerate(varyNames):
176        for i2,name2 in enumerate(varyNames):
177            try:
178                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
179            except ValueError:
180                vcov[i1][i2] = 0.0
181    return vcov
182
183def FindAtomIndexByIDs(atomData,IDs,Draw=True):
184    '''finds the set of atom array indices for a list of atom IDs. Will search
185    either the Atom table or the drawAtom table.
186   
187    :param list atomData: Atom or drawAtom table containting coordinates, etc.
188    :param list IDs: atom IDs to be found
189    :param bool Draw: True if drawAtom table to be searched; False if Atom table
190      is searched
191   
192    :returns: list indx: atom (or drawAtom) indices
193   
194    '''
195    indx = []
196    for i,atom in enumerate(atomData):
197        if Draw and atom[-3] in IDs:
198            indx.append(i)
199        elif atom[-1] in IDs:
200            indx.append(i)
201    return indx
202
203def FillAtomLookUp(atomData):
204    '''create a dictionary of atom indexes with atom IDs as keys
205   
206    :param list atomData: Atom table to be used
207   
208    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
209   
210    '''
211    atomLookUp = {}
212    for iatm,atom in enumerate(atomData):
213        atomLookUp[atom[-1]] = iatm
214    return atomLookUp
215
216def GetAtomsById(atomData,atomLookUp,IdList):
217    '''gets a list of atoms from Atom table that match a set of atom IDs
218   
219    :param list atomData: Atom table to be used
220    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
221    :param list IdList: atom IDs to be found
222   
223    :returns: list atoms: list of atoms found
224   
225    '''
226    atoms = []
227    for id in IdList:
228        atoms.append(atomData[atomLookUp[id]])
229    return atoms
230   
231def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
232    '''gets atom parameters for atoms using atom IDs
233   
234    :param list atomData: Atom table to be used
235    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
236    :param list IdList: atom IDs to be found
237    :param int itemLoc: pointer to desired 1st item in an atom table entry
238    :param int numItems: number of items to be retrieved
239   
240    :returns: type name: description
241   
242    '''
243    Items = []
244    if not isinstance(IdList,list):
245        IdList = [IdList,]
246    for id in IdList:
247        if numItems == 1:
248            Items.append(atomData[atomLookUp[id]][itemLoc])
249        else:
250            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
251    return Items
252   
253def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
254    '''default doc string
255   
256    :param type name: description
257   
258    :returns: type name: description
259   
260    '''
261    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
262    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
263    XYZ = []
264    for ind in indx:
265        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
266        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
267        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
268    return XYZ
269
270def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
271    '''default doc string
272   
273    :param type name: description
274   
275    :returns: type name: description
276   
277    '''
278    for atom in atomData:
279        XYZ = np.inner(Amat,atom[cx:cx+3])
280        if atom[cia] == 'A':
281            UIJ = atom[cia+2:cia+8]
282               
283def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
284    '''default doc string
285   
286    :param type name: description
287   
288    :returns: type name: description
289   
290    '''
291    TLStype,TLS = rbObj['ThermalMotion'][:2]
292    Tmat = np.zeros((3,3))
293    Lmat = np.zeros((3,3))
294    Smat = np.zeros((3,3))
295    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
296        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
297    if 'T' in TLStype:
298        Tmat = G2lat.U6toUij(TLS[:6])
299    if 'L' in TLStype:
300        Lmat = G2lat.U6toUij(TLS[6:12])
301    if 'S' in TLStype:
302        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
303    XYZ = np.inner(Amat,xyz)
304    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
305    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
306    beta = np.inner(np.inner(g,Umat),g)
307    return G2lat.UijtoU6(beta)*gvec   
308       
309def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
310    '''default doc string
311   
312    :param type name: description
313   
314    :returns: type name: description
315   
316    '''
317    cx,ct,cs,cia = atPtrs
318    TLStype,TLS = rbObj['ThermalMotion'][:2]
319    Tmat = np.zeros((3,3))
320    Lmat = np.zeros((3,3))
321    Smat = np.zeros((3,3))
322    G,g = G2lat.A2Gmat(Amat)
323    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
324    if 'T' in TLStype:
325        Tmat = G2lat.U6toUij(TLS[:6])
326    if 'L' in TLStype:
327        Lmat = G2lat.U6toUij(TLS[6:12])
328    if 'S' in TLStype:
329        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
330    for atom in atomData:
331        XYZ = np.inner(Amat,atom[cx:cx+3])
332        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
333        if 'U' in TSLtype:
334            atom[cia+1] = TLS[0]
335            atom[cia] = 'I'
336        else:
337            atom[cia] = 'A'
338            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
339            beta = np.inner(np.inner(g,Umat),g)
340            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
341
342def GetXYZDist(xyz,XYZ,Amat):
343    '''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
344        and are in crystal coordinates; Amat is crystal to Cart matrix
345   
346    :param type name: description
347   
348    :returns: type name: description
349   
350    '''
351    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
352
353def getAtomXYZ(atoms,cx):
354    '''default doc string
355   
356    :param type name: description
357   
358    :returns: type name: description
359   
360    '''
361    XYZ = []
362    for atom in atoms:
363        XYZ.append(atom[cx:cx+3])
364    return np.array(XYZ)
365
366def RotateRBXYZ(Bmat,Cart,oriQ):
367    '''rotate & transform cartesian coordinates to crystallographic ones
368    no translation applied. To be used for numerical derivatives
369   
370    :param type name: description
371   
372    :returns: type name: description
373   
374    '''
375    ''' returns crystal coordinates for atoms described by RBObj
376    '''
377    XYZ = np.zeros_like(Cart)
378    for i,xyz in enumerate(Cart):
379        X = prodQVQ(oriQ,xyz)
380        XYZ[i] = np.inner(Bmat,X)
381    return XYZ
382
383def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
384    '''default doc string
385   
386    :param type name: description
387   
388    :returns: type name: description
389   
390    '''
391    ''' returns crystal coordinates for atoms described by RBObj
392    '''
393    RBRes = RBData[RBType][RBObj['RBId']]
394    if RBType == 'Vector':
395        vecs = RBRes['rbVect']
396        mags = RBRes['VectMag']
397        Cart = np.zeros_like(vecs[0])
398        for vec,mag in zip(vecs,mags):
399            Cart += vec*mag
400    elif RBType == 'Residue':
401        Cart = np.array(RBRes['rbXYZ'])
402        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
403            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
404            for ride in seq[3]:
405                Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
406    XYZ = np.zeros_like(Cart)
407    for i,xyz in enumerate(Cart):
408        X = prodQVQ(RBObj['Orient'][0],xyz)
409        XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
410    return XYZ,Cart
411
412def UpdateMCSAxyz(Bmat,MCSA):
413    '''default doc string
414   
415    :param type name: description
416   
417    :returns: type name: description
418   
419    '''
420    xyz = []
421    atTypes = []
422    iatm = 0
423    for model in MCSA['Models'][1:]:        #skip the MD model
424        if model['Type'] == 'Atom':
425            xyz.append(model['Pos'][0])
426            atTypes.append(model['atType'])
427            iatm += 1
428        else:
429            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
430            Pos = np.array(model['Pos'][0])
431            Ori = np.array(model['Ori'][0])
432            Qori = AVdeg2Q(Ori[0],Ori[1:])
433            if model['Type'] == 'Vector':
434                vecs = RBRes['rbVect']
435                mags = RBRes['VectMag']
436                Cart = np.zeros_like(vecs[0])
437                for vec,mag in zip(vecs,mags):
438                    Cart += vec*mag
439            elif model['Type'] == 'Residue':
440                Cart = np.array(RBRes['rbXYZ'])
441                for itor,seq in enumerate(RBRes['rbSeq']):
442                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
443                    for ride in seq[3]:
444                        Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
445            if model['MolCent'][1]:
446                Cart -= model['MolCent'][0]
447            for i,x in enumerate(Cart):
448                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
449                atType = RBRes['rbTypes'][i]
450                atTypes.append(atType)
451                iatm += 1
452    return np.array(xyz),atTypes
453   
454def SetMolCent(model,RBData):
455    '''default doc string
456   
457    :param type name: description
458   
459    :returns: type name: description
460   
461    '''
462    rideList = []
463    RBRes = RBData[model['Type']][model['RBId']]
464    if model['Type'] == 'Vector':
465        vecs = RBRes['rbVect']
466        mags = RBRes['VectMag']
467        Cart = np.zeros_like(vecs[0])
468        for vec,mag in zip(vecs,mags):
469            Cart += vec*mag
470    elif model['Type'] == 'Residue':
471        Cart = np.array(RBRes['rbXYZ'])
472        for seq in RBRes['rbSeq']:
473            rideList += seq[3]
474    centList = set(range(len(Cart)))-set(rideList)
475    cent = np.zeros(3)
476    for i in centList:
477        cent += Cart[i]
478    model['MolCent'][0] = cent/len(centList)   
479   
480def UpdateRBUIJ(Bmat,Cart,RBObj):
481    '''default doc string
482   
483    :param type name: description
484   
485    :returns: type name: description
486   
487    '''
488    ''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
489    '''
490    TLStype,TLS = RBObj['ThermalMotion'][:2]
491    T = np.zeros(6)
492    L = np.zeros(6)
493    S = np.zeros(8)
494    if 'T' in TLStype:
495        T = TLS[:6]
496    if 'L' in TLStype:
497        L = np.array(TLS[6:12])*(np.pi/180.)**2
498    if 'S' in TLStype:
499        S = np.array(TLS[12:])*(np.pi/180.)
500    g = nl.inv(np.inner(Bmat,Bmat))
501    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
502        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
503    Uout = []
504    Q = RBObj['Orient'][0]
505    for X in Cart:
506        X = prodQVQ(Q,X)
507        if 'U' in TLStype:
508            Uout.append(['I',TLS[0],0,0,0,0,0,0])
509        elif not 'N' in TLStype:
510            U = [0,0,0,0,0,0]
511            U[0] = T[0]+L[1]*X[2]**2+L[2]*X[1]**2-2.0*L[5]*X[1]*X[2]+2.0*(S[2]*X[2]-S[4]*X[1])
512            U[1] = T[1]+L[0]*X[2]**2+L[2]*X[0]**2-2.0*L[4]*X[0]*X[2]+2.0*(S[5]*X[0]-S[0]*X[2])
513            U[2] = T[2]+L[1]*X[0]**2+L[0]*X[1]**2-2.0*L[3]*X[1]*X[0]+2.0*(S[1]*X[1]-S[3]*X[0])
514            U[3] = T[3]+L[4]*X[1]*X[2]+L[5]*X[0]*X[2]-L[3]*X[2]**2-L[2]*X[0]*X[1]+  \
515                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
516            U[4] = T[4]+L[3]*X[1]*X[2]+L[5]*X[0]*X[1]-L[4]*X[1]**2-L[1]*X[0]*X[2]+  \
517                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
518            U[5] = T[5]+L[3]*X[0]*X[2]+L[4]*X[0]*X[1]-L[5]*X[0]**2-L[0]*X[2]*X[1]+  \
519                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
520            Umat = G2lat.U6toUij(U)
521            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
522            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
523        else:
524            Uout.append(['N',])
525    return Uout
526   
527def GetSHCoeff(pId,parmDict,SHkeys):
528    '''default doc string
529   
530    :param type name: description
531   
532    :returns: type name: description
533   
534    '''
535    SHCoeff = {}
536    for shkey in SHkeys:
537        shname = str(pId)+'::'+shkey
538        SHCoeff[shkey] = parmDict[shname]
539    return SHCoeff
540       
541def getMass(generalData):
542    '''default doc string
543   
544    :param type name: description
545   
546    :returns: type name: description
547   
548    '''
549    'Computes mass of unit cell contents'
550    mass = 0.
551    for i,elem in enumerate(generalData['AtomTypes']):
552        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
553    return mass   
554
555def getDensity(generalData):
556    '''default doc string
557   
558    :param type name: description
559   
560    :returns: type name: description
561   
562    '''
563    mass = getMass(generalData)
564    Volume = generalData['Cell'][7]
565    density = mass/(0.6022137*Volume)
566    return density,Volume/mass
567   
568def getWave(Parms):
569    '''default doc string
570   
571    :param type name: description
572   
573    :returns: type name: description
574   
575    '''
576    try:
577        return Parms['Lam'][1]
578    except KeyError:
579        return Parms['Lam1'][1]
580   
581################################################################################
582##### distance, angle, planes, torsion stuff stuff
583################################################################################
584
585def getSyXYZ(XYZ,ops,SGData):
586    '''default doc string
587   
588    :param type name: description
589   
590    :returns: type name: description
591   
592    '''
593    XYZout = np.zeros_like(XYZ)
594    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
595        if op == '1':
596            XYZout[i] = xyz
597        else:
598            oprs = op.split('+')
599            unit = [0,0,0]
600            if oprs[1]:
601                unit = np.array(list(eval(oprs[1])))
602            syop =int(oprs[0])
603            inv = syop/abs(syop)
604            syop *= inv
605            cent = syop/100
606            syop %= 100
607            syop -= 1
608            M,T = SGData['SGOps'][syop]
609            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
610    return XYZout
611   
612def getRestDist(XYZ,Amat):
613    '''default doc string
614   
615    :param type name: description
616   
617    :returns: type name: description
618   
619    '''
620    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
621   
622def getRestDeriv(Func,XYZ,Amat,ops,SGData):
623    '''default doc string
624   
625    :param type name: description
626   
627    :returns: type name: description
628   
629    '''
630    deriv = np.zeros((len(XYZ),3))
631    dx = 0.00001
632    for j,xyz in enumerate(XYZ):
633        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
634            XYZ[j] -= x
635            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
636            XYZ[j] += 2*x
637            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
638            XYZ[j] -= x
639            deriv[j][i] = (d1-d2)/(2*dx)
640    return deriv.flatten()
641
642def getRestAngle(XYZ,Amat):
643    '''default doc string
644   
645    :param type name: description
646   
647    :returns: type name: description
648   
649    '''
650   
651    def calcVec(Ox,Tx,Amat):
652        return np.inner(Amat,(Tx-Ox))
653
654    VecA = calcVec(XYZ[1],XYZ[0],Amat)
655    VecA /= np.sqrt(np.sum(VecA**2))
656    VecB = calcVec(XYZ[1],XYZ[2],Amat)
657    VecB /= np.sqrt(np.sum(VecB**2))
658    edge = VecB-VecA
659    edge = np.sum(edge**2)
660    angle = (2.-edge)/2.
661    angle = max(angle,-1.)
662    return acosd(angle)
663   
664def getRestPlane(XYZ,Amat):
665    '''default doc string
666   
667    :param type name: description
668   
669    :returns: type name: description
670   
671    '''
672    sumXYZ = np.zeros(3)
673    for xyz in XYZ:
674        sumXYZ += xyz
675    sumXYZ /= len(XYZ)
676    XYZ = np.array(XYZ)-sumXYZ
677    XYZ = np.inner(Amat,XYZ).T
678    Zmat = np.zeros((3,3))
679    for i,xyz in enumerate(XYZ):
680        Zmat += np.outer(xyz.T,xyz)
681    Evec,Emat = nl.eig(Zmat)
682    Evec = np.sqrt(Evec)/(len(XYZ)-3)
683    Order = np.argsort(Evec)
684    return Evec[Order[0]]
685   
686def getRestChiral(XYZ,Amat):   
687    '''default doc string
688   
689    :param type name: description
690   
691    :returns: type name: description
692   
693    '''
694    VecA = np.empty((3,3))   
695    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
696    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
697    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
698    return nl.det(VecA)
699   
700def getRestTorsion(XYZ,Amat):
701    '''default doc string
702   
703    :param type name: description
704   
705    :returns: type name: description
706   
707    '''
708    VecA = np.empty((3,3))
709    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
710    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
711    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
712    D = nl.det(VecA)
713    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
714    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
715    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
716    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
717    Ang = 1.0
718    if abs(P12) < 1.0 and abs(P23) < 1.0:
719        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
720    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
721    return TOR
722   
723def calcTorsionEnergy(TOR,Coeff=[]):
724    '''default doc string
725   
726    :param type name: description
727   
728    :returns: type name: description
729   
730    '''
731    sum = 0.
732    if len(Coeff):
733        cof = np.reshape(Coeff,(3,3)).T
734        delt = TOR-cof[1]
735        delt = np.where(delt<-180.,delt+360.,delt)
736        delt = np.where(delt>180.,delt-360.,delt)
737        term = -cof[2]*delt**2
738        val = cof[0]*np.exp(term/1000.0)
739        pMax = cof[0][np.argmin(val)]
740        Eval = np.sum(val)
741        sum = Eval-pMax
742    return sum,Eval
743
744def getTorsionDeriv(XYZ,Amat,Coeff):
745    '''default doc string
746   
747    :param type name: description
748   
749    :returns: type name: description
750   
751    '''
752    deriv = np.zeros((len(XYZ),3))
753    dx = 0.00001
754    for j,xyz in enumerate(XYZ):
755        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
756            XYZ[j] -= x
757            tor = getRestTorsion(XYZ,Amat)
758            p,d1 = calcTorsionEnergy(tor,Coeff)
759            XYZ[j] += 2*x
760            tor = getRestTorsion(XYZ,Amat)
761            p,d2 = calcTorsionEnergy(tor,Coeff)           
762            XYZ[j] -= x
763            deriv[j][i] = (d2-d1)/(2*dx)
764    return deriv.flatten()
765
766def getRestRama(XYZ,Amat):
767    '''Computes a pair of torsion angles in a 5 atom string
768   
769    :param nparray XYZ: crystallographic coordinates of 5 atoms
770    :param nparray Amat: crystal to cartesian transformation matrix
771   
772    :returns: list (phi,psi) two torsion angles in degrees
773   
774    '''
775    phi = getRestTorsion(XYZ[:5],Amat)
776    psi = getRestTorsion(XYZ[1:],Amat)
777    return phi,psi
778   
779def calcRamaEnergy(phi,psi,Coeff=[]):
780    '''Computes pseudo potential energy from a pair of torsion angles and a
781    numerical description of the potential energy surface. Used to create
782    penalty function in LS refinement:     
783    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
784
785    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
786   
787    :param float phi: first torsion angle (:math:`\\phi`)
788    :param float psi: second torsion angle (:math:`\\psi`)
789    :param list Coeff: pseudo potential coefficients
790   
791    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
792      sum is used for penalty function.
793   
794    '''
795    sum = 0.
796    Eval = 0.
797    if len(Coeff):
798        cof = Coeff.T
799        dPhi = phi-cof[1]
800        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
801        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
802        dPsi = psi-cof[2]
803        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
804        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
805        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
806        val = cof[0]*np.exp(val/1000.)
807        pMax = cof[0][np.argmin(val)]
808        Eval = np.sum(val)
809        sum = Eval-pMax
810    return sum,Eval
811
812def getRamaDeriv(XYZ,Amat,Coeff):
813    '''Computes numerical derivatives of torsion angle pair pseudo potential
814    with respect of crystallographic atom coordinates of the 5 atom sequence
815   
816    :param nparray XYZ: crystallographic coordinates of 5 atoms
817    :param nparray Amat: crystal to cartesian transformation matrix
818    :param list Coeff: pseudo potential coefficients
819   
820    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
821     crystallographic xyz coordinates.
822   
823    '''
824    deriv = np.zeros((len(XYZ),3))
825    dx = 0.00001
826    for j,xyz in enumerate(XYZ):
827        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
828            XYZ[j] -= x
829            phi,psi = getRestRama(XYZ,Amat)
830            p,d1 = calcRamaEnergy(phi,psi,Coeff)
831            XYZ[j] += 2*x
832            phi,psi = getRestRama(XYZ,Amat)
833            p,d2 = calcRamaEnergy(phi,psi,Coeff)
834            XYZ[j] -= x
835            deriv[j][i] = (d2-d1)/(2*dx)
836    return deriv.flatten()
837
838def getRestPolefig(ODFln,SamSym,Grid):
839    '''default doc string
840   
841    :param type name: description
842   
843    :returns: type name: description
844   
845    '''
846    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
847    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
848    R = np.where(R <= 1.,2.*atand(R),0.0)
849    Z = np.zeros_like(R)
850    Z = G2lat.polfcal(ODFln,SamSym,R,P)
851    Z = np.reshape(Z,(Grid,Grid))
852    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
853
854def getRestPolefigDerv(HKL,Grid,SHCoeff):
855    '''default doc string
856   
857    :param type name: description
858   
859    :returns: type name: description
860   
861    '''
862    pass
863       
864def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
865    '''default doc string
866   
867    :param type name: description
868   
869    :returns: type name: description
870   
871    '''
872    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
873        TxT = inv*(np.inner(M,Tx)+T)+C+U
874        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
875       
876    inv = Top/abs(Top)
877    cent = abs(Top)/100
878    op = abs(Top)%100-1
879    M,T = SGData['SGOps'][op]
880    C = SGData['SGCen'][cent]
881    dx = .00001
882    deriv = np.zeros(6)
883    for i in [0,1,2]:
884        Oxyz[i] -= dx
885        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
886        Oxyz[i] += 2*dx
887        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
888        Oxyz[i] -= dx
889        Txyz[i] -= dx
890        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
891        Txyz[i] += 2*dx
892        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
893        Txyz[i] -= dx
894    return deriv
895   
896def getAngSig(VA,VB,Amat,SGData,covData={}):
897    '''default doc string
898   
899    :param type name: description
900   
901    :returns: type name: description
902   
903    '''
904    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
905        TxT = inv*(np.inner(M,Tx)+T)+C
906        TxT = G2spc.MoveToUnitCell(TxT)+U
907        return np.inner(Amat,(TxT-Ox))
908       
909    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
910        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
911        VecA /= np.sqrt(np.sum(VecA**2))
912        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
913        VecB /= np.sqrt(np.sum(VecB**2))
914        edge = VecB-VecA
915        edge = np.sum(edge**2)
916        angle = (2.-edge)/2.
917        angle = max(angle,-1.)
918        return acosd(angle)
919       
920    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
921    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
922    invA = invB = 1
923    invA = TopA/abs(TopA)
924    invB = TopB/abs(TopB)
925    centA = abs(TopA)/100
926    centB = abs(TopB)/100
927    opA = abs(TopA)%100-1
928    opB = abs(TopB)%100-1
929    MA,TA = SGData['SGOps'][opA]
930    MB,TB = SGData['SGOps'][opB]
931    CA = SGData['SGCen'][centA]
932    CB = SGData['SGCen'][centB]
933    if 'covMatrix' in covData:
934        covMatrix = covData['covMatrix']
935        varyList = covData['varyList']
936        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
937        dx = .00001
938        dadx = np.zeros(9)
939        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
940        for i in [0,1,2]:
941            OxA[i] -= dx
942            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
943            OxA[i] += 2*dx
944            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
945            OxA[i] -= dx
946           
947            TxA[i] -= dx
948            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
949            TxA[i] += 2*dx
950            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
951            TxA[i] -= dx
952           
953            TxB[i] -= dx
954            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
955            TxB[i] += 2*dx
956            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
957            TxB[i] -= dx
958           
959        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
960        if sigAng < 0.01:
961            sigAng = 0.0
962        return Ang,sigAng
963    else:
964        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
965
966def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
967    '''default doc string
968   
969    :param type name: description
970   
971    :returns: type name: description
972   
973    '''
974    def calcDist(Atoms,SyOps,Amat):
975        XYZ = []
976        for i,atom in enumerate(Atoms):
977            Inv,M,T,C,U = SyOps[i]
978            XYZ.append(np.array(atom[1:4]))
979            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
980            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
981        V1 = XYZ[1]-XYZ[0]
982        return np.sqrt(np.sum(V1**2))
983       
984    Inv = []
985    SyOps = []
986    names = []
987    for i,atom in enumerate(Oatoms):
988        names += atom[-1]
989        Op,unit = Atoms[i][-1]
990        inv = Op/abs(Op)
991        m,t = SGData['SGOps'][abs(Op)%100-1]
992        c = SGData['SGCen'][abs(Op)/100]
993        SyOps.append([inv,m,t,c,unit])
994    Dist = calcDist(Oatoms,SyOps,Amat)
995   
996    sig = -0.001
997    if 'covMatrix' in covData:
998        parmNames = []
999        dx = .00001
1000        dadx = np.zeros(6)
1001        for i in range(6):
1002            ia = i/3
1003            ix = i%3
1004            Oatoms[ia][ix+1] += dx
1005            a0 = calcDist(Oatoms,SyOps,Amat)
1006            Oatoms[ia][ix+1] -= 2*dx
1007            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1008        covMatrix = covData['covMatrix']
1009        varyList = covData['varyList']
1010        DistVcov = getVCov(names,varyList,covMatrix)
1011        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1012        if sig < 0.001:
1013            sig = -0.001
1014   
1015    return Dist,sig
1016
1017def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
1018    '''default doc string
1019   
1020    :param type name: description
1021   
1022    :returns: type name: description
1023   
1024    '''
1025
1026    def calcAngle(Atoms,SyOps,Amat):
1027        XYZ = []
1028        for i,atom in enumerate(Atoms):
1029            Inv,M,T,C,U = SyOps[i]
1030            XYZ.append(np.array(atom[1:4]))
1031            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1032            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1033        V1 = XYZ[1]-XYZ[0]
1034        V1 /= np.sqrt(np.sum(V1**2))
1035        V2 = XYZ[1]-XYZ[2]
1036        V2 /= np.sqrt(np.sum(V2**2))
1037        V3 = V2-V1
1038        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1039        return acosd(cang)
1040
1041    Inv = []
1042    SyOps = []
1043    names = []
1044    for i,atom in enumerate(Oatoms):
1045        names += atom[-1]
1046        Op,unit = Atoms[i][-1]
1047        inv = Op/abs(Op)
1048        m,t = SGData['SGOps'][abs(Op)%100-1]
1049        c = SGData['SGCen'][abs(Op)/100]
1050        SyOps.append([inv,m,t,c,unit])
1051    Angle = calcAngle(Oatoms,SyOps,Amat)
1052   
1053    sig = -0.01
1054    if 'covMatrix' in covData:
1055        parmNames = []
1056        dx = .00001
1057        dadx = np.zeros(9)
1058        for i in range(9):
1059            ia = i/3
1060            ix = i%3
1061            Oatoms[ia][ix+1] += dx
1062            a0 = calcAngle(Oatoms,SyOps,Amat)
1063            Oatoms[ia][ix+1] -= 2*dx
1064            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1065        covMatrix = covData['covMatrix']
1066        varyList = covData['varyList']
1067        AngVcov = getVCov(names,varyList,covMatrix)
1068        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1069        if sig < 0.01:
1070            sig = -0.01
1071   
1072    return Angle,sig
1073
1074def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1075    '''default doc string
1076   
1077    :param type name: description
1078   
1079    :returns: type name: description
1080   
1081    '''
1082
1083    def calcTorsion(Atoms,SyOps,Amat):
1084       
1085        XYZ = []
1086        for i,atom in enumerate(Atoms):
1087            Inv,M,T,C,U = SyOps[i]
1088            XYZ.append(np.array(atom[1:4]))
1089            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1090            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1091        V1 = XYZ[1]-XYZ[0]
1092        V2 = XYZ[2]-XYZ[1]
1093        V3 = XYZ[3]-XYZ[2]
1094        V1 /= np.sqrt(np.sum(V1**2))
1095        V2 /= np.sqrt(np.sum(V2**2))
1096        V3 /= np.sqrt(np.sum(V3**2))
1097        M = np.array([V1,V2,V3])
1098        D = nl.det(M)
1099        Ang = 1.0
1100        P12 = np.dot(V1,V2)
1101        P13 = np.dot(V1,V3)
1102        P23 = np.dot(V2,V3)
1103        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1104        return Tors
1105           
1106    Inv = []
1107    SyOps = []
1108    names = []
1109    for i,atom in enumerate(Oatoms):
1110        names += atom[-1]
1111        Op,unit = Atoms[i][-1]
1112        inv = Op/abs(Op)
1113        m,t = SGData['SGOps'][abs(Op)%100-1]
1114        c = SGData['SGCen'][abs(Op)/100]
1115        SyOps.append([inv,m,t,c,unit])
1116    Tors = calcTorsion(Oatoms,SyOps,Amat)
1117   
1118    sig = -0.01
1119    if 'covMatrix' in covData:
1120        parmNames = []
1121        dx = .00001
1122        dadx = np.zeros(12)
1123        for i in range(12):
1124            ia = i/3
1125            ix = i%3
1126            Oatoms[ia][ix+1] -= dx
1127            a0 = calcTorsion(Oatoms,SyOps,Amat)
1128            Oatoms[ia][ix+1] += 2*dx
1129            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1130            Oatoms[ia][ix+1] -= dx           
1131        covMatrix = covData['covMatrix']
1132        varyList = covData['varyList']
1133        TorVcov = getVCov(names,varyList,covMatrix)
1134        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1135        if sig < 0.01:
1136            sig = -0.01
1137   
1138    return Tors,sig
1139       
1140def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1141    '''default doc string
1142   
1143    :param type name: description
1144   
1145    :returns: type name: description
1146   
1147    '''
1148
1149    def calcDist(Atoms,SyOps,Amat):
1150        XYZ = []
1151        for i,atom in enumerate(Atoms):
1152            Inv,M,T,C,U = SyOps[i]
1153            XYZ.append(np.array(atom[1:4]))
1154            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1155            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1156        V1 = XYZ[1]-XYZ[0]
1157        return np.sqrt(np.sum(V1**2))
1158       
1159    def calcAngle(Atoms,SyOps,Amat):
1160        XYZ = []
1161        for i,atom in enumerate(Atoms):
1162            Inv,M,T,C,U = SyOps[i]
1163            XYZ.append(np.array(atom[1:4]))
1164            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1165            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1166        V1 = XYZ[1]-XYZ[0]
1167        V1 /= np.sqrt(np.sum(V1**2))
1168        V2 = XYZ[1]-XYZ[2]
1169        V2 /= np.sqrt(np.sum(V2**2))
1170        V3 = V2-V1
1171        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1172        return acosd(cang)
1173
1174    def calcTorsion(Atoms,SyOps,Amat):
1175       
1176        XYZ = []
1177        for i,atom in enumerate(Atoms):
1178            Inv,M,T,C,U = SyOps[i]
1179            XYZ.append(np.array(atom[1:4]))
1180            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1181            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1182        V1 = XYZ[1]-XYZ[0]
1183        V2 = XYZ[2]-XYZ[1]
1184        V3 = XYZ[3]-XYZ[2]
1185        V1 /= np.sqrt(np.sum(V1**2))
1186        V2 /= np.sqrt(np.sum(V2**2))
1187        V3 /= np.sqrt(np.sum(V3**2))
1188        M = np.array([V1,V2,V3])
1189        D = nl.det(M)
1190        Ang = 1.0
1191        P12 = np.dot(V1,V2)
1192        P13 = np.dot(V1,V3)
1193        P23 = np.dot(V2,V3)
1194        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1195        return Tors
1196           
1197    Inv = []
1198    SyOps = []
1199    names = []
1200    for i,atom in enumerate(Oatoms):
1201        names += atom[-1]
1202        Op,unit = Atoms[i][-1]
1203        inv = Op/abs(Op)
1204        m,t = SGData['SGOps'][abs(Op)%100-1]
1205        c = SGData['SGCen'][abs(Op)/100]
1206        SyOps.append([inv,m,t,c,unit])
1207    M = len(Oatoms)
1208    if M == 2:
1209        Val = calcDist(Oatoms,SyOps,Amat)
1210    elif M == 3:
1211        Val = calcAngle(Oatoms,SyOps,Amat)
1212    else:
1213        Val = calcTorsion(Oatoms,SyOps,Amat)
1214   
1215    sigVals = [-0.001,-0.01,-0.01]
1216    sig = sigVals[M-3]
1217    if 'covMatrix' in covData:
1218        parmNames = []
1219        dx = .00001
1220        N = M*3
1221        dadx = np.zeros(N)
1222        for i in range(N):
1223            ia = i/3
1224            ix = i%3
1225            Oatoms[ia][ix+1] += dx
1226            if M == 2:
1227                a0 = calcDist(Oatoms,SyOps,Amat)
1228            elif M == 3:
1229                a0 = calcAngle(Oatoms,SyOps,Amat)
1230            else:
1231                a0 = calcTorsion(Oatoms,SyOps,Amat)
1232            Oatoms[ia][ix+1] -= 2*dx
1233            if M == 2:
1234                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1235            elif M == 3:
1236                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1237            else:
1238                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1239        covMatrix = covData['covMatrix']
1240        varyList = covData['varyList']
1241        Vcov = getVCov(names,varyList,covMatrix)
1242        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1243        if sig < sigVals[M-3]:
1244            sig = sigVals[M-3]
1245   
1246    return Val,sig
1247       
1248def ValEsd(value,esd=0,nTZ=False):
1249    '''Format a floating point number with a given level of precision or
1250    with in crystallographic format with a "esd", as value(esd). If esd is
1251    negative the number is formatted with the level of significant figures
1252    appropriate if abs(esd) were the esd, but the esd is not included.
1253    if the esd is zero, approximately 6 significant figures are printed.
1254    nTZ=True causes "extra" zeros to be removed after the decimal place.
1255    for example:
1256
1257      * "1.235(3)" for value=1.2346 & esd=0.003
1258      * "1.235(3)e4" for value=12346. & esd=30
1259      * "1.235(3)e6" for value=0.12346e7 & esd=3000
1260      * "1.235" for value=1.2346 & esd=-0.003
1261      * "1.240" for value=1.2395 & esd=-0.003
1262      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
1263      * "1.23460" for value=1.2346 & esd=0.0
1264
1265    :param float value: number to be formatted
1266    :param float esd: uncertainty or if esd < 0, specifies level of
1267      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
1268    :param bool nTZ: True to remove trailing zeros (default is False)
1269    :returns: value(esd) or value as a string
1270
1271    '''
1272    # Note: this routine is Python 3 compatible -- I think
1273    if math.isnan(value): # invalid value, bail out
1274        return '?'
1275    if math.isnan(esd): # invalid esd, treat as zero
1276        esd = 0
1277        esdoff = 5
1278    elif esd != 0:
1279        # transform the esd to a one or two digit integer
1280        l = math.log10(abs(esd)) % 1
1281        # cut off of 19 1.9==>(19) but 1.95==>(2) N.B. log10(1.95) = 0.2900...
1282        if l < 0.290034611362518: l+= 1.       
1283        intesd = int(round(10**l)) # esd as integer
1284        # determine the number of digits offset for the esd
1285        esdoff = int(round(math.log10(intesd*1./abs(esd))))
1286    else:
1287        esdoff = 5
1288    valoff = 0
1289    if abs(value) < abs(esdoff): # value is effectively zero
1290        pass
1291    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
1292        # where the digit offset is to the left of the decimal place or where too many
1293        # digits are needed
1294        if abs(value) > 1:
1295            valoff = int(math.log10(abs(value)))
1296        elif abs(value) > 0:
1297            valoff = int(math.log10(abs(value))-0.9999999)
1298        else:
1299            valoff = 0
1300    if esd != 0:
1301        if valoff+esdoff < 0:
1302            valoff = esdoff = 0
1303        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1304    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1305        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1306    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1307        if abs(value) > 0:           
1308            extra = -math.log10(abs(value))
1309        else:
1310            extra = 0
1311        if extra > 0: extra += 1
1312        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1313    if esd > 0:
1314        out += ("({:d})").format(intesd)  # add the esd
1315    elif nTZ and '.' in out:
1316        out = out.rstrip('0')  # strip zeros to right of decimal
1317        out = out.rstrip('.')  # and decimal place when not needed
1318    if valoff != 0:
1319        out += ("e{:d}").format(valoff) # add an exponent, when needed
1320    return out
1321
1322################################################################################
1323##### Fourier & charge flip stuff
1324################################################################################
1325
1326def adjHKLmax(SGData,Hmax):
1327    '''default doc string
1328   
1329    :param type name: description
1330   
1331    :returns: type name: description
1332   
1333    '''
1334    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1335        Hmax[0] = ((Hmax[0]+3)/6)*6
1336        Hmax[1] = ((Hmax[1]+3)/6)*6
1337        Hmax[2] = ((Hmax[2]+1)/4)*4
1338    else:
1339        Hmax[0] = ((Hmax[0]+2)/4)*4
1340        Hmax[1] = ((Hmax[1]+2)/4)*4
1341        Hmax[2] = ((Hmax[2]+1)/4)*4
1342
1343def OmitMap(data,reflData):
1344    '''default doc string
1345   
1346    :param type name: description
1347   
1348    :returns: type name: description
1349   
1350    '''
1351    generalData = data['General']
1352    if not generalData['Map']['MapType']:
1353        print '**** ERROR - Fourier map not defined'
1354        return
1355    mapData = generalData['Map']
1356    dmin = mapData['Resolution']
1357    SGData = generalData['SGData']
1358    cell = generalData['Cell'][1:8]       
1359    A = G2lat.cell2A(cell[:6])
1360    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1361    adjHKLmax(SGData,Hmax)
1362    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1363    time0 = time.time()
1364    for ref in reflData:
1365        if ref[4] >= dmin:
1366            Fosq,Fcsq,ph = ref[8:11]
1367            for i,hkl in enumerate(ref[11]):
1368                hkl = np.asarray(hkl,dtype='i')
1369                dp = 360.*ref[12][i]
1370                a = cosd(ph+dp)
1371                b = sind(ph+dp)
1372                phasep = complex(a,b)
1373                phasem = complex(a,-b)
1374                F = np.sqrt(Fosq)
1375                h,k,l = hkl+Hmax
1376                Fhkl[h,k,l] = F*phasep
1377                h,k,l = -hkl+Hmax
1378                Fhkl[h,k,l] = F*phasem
1379    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1380    print 'NB: this is just an Fobs map for now - under development'
1381    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1382    print rho.shape
1383    mapData['rho'] = np.real(rho)
1384    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1385    return mapData
1386
1387def FourierMap(data,reflData):
1388    '''default doc string
1389   
1390    :param type name: description
1391   
1392    :returns: type name: description
1393   
1394    '''
1395    generalData = data['General']
1396    if not generalData['Map']['MapType']:
1397        print '**** ERROR - Fourier map not defined'
1398        return
1399    mapData = generalData['Map']
1400    dmin = mapData['Resolution']
1401    SGData = generalData['SGData']
1402    cell = generalData['Cell'][1:8]       
1403    A = G2lat.cell2A(cell[:6])
1404    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1405    adjHKLmax(SGData,Hmax)
1406    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1407#    Fhkl[0,0,0] = generalData['F000X']
1408    time0 = time.time()
1409    for ref in reflData:
1410        if ref[4] >= dmin:
1411            Fosq,Fcsq,ph = ref[8:11]
1412            for i,hkl in enumerate(ref[11]):
1413                hkl = np.asarray(hkl,dtype='i')
1414                dp = 360.*ref[12][i]
1415                a = cosd(ph+dp)
1416                b = sind(ph+dp)
1417                phasep = complex(a,b)
1418                phasem = complex(a,-b)
1419                if 'Fobs' in mapData['MapType']:
1420                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1421                    h,k,l = hkl+Hmax
1422                    Fhkl[h,k,l] = F*phasep
1423                    h,k,l = -hkl+Hmax
1424                    Fhkl[h,k,l] = F*phasem
1425                elif 'Fcalc' in mapData['MapType']:
1426                    F = np.sqrt(Fcsq)
1427                    h,k,l = hkl+Hmax
1428                    Fhkl[h,k,l] = F*phasep
1429                    h,k,l = -hkl+Hmax
1430                    Fhkl[h,k,l] = F*phasem
1431                elif 'delt-F' in mapData['MapType']:
1432                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1433                    h,k,l = hkl+Hmax
1434                    Fhkl[h,k,l] = dF*phasep
1435                    h,k,l = -hkl+Hmax
1436                    Fhkl[h,k,l] = dF*phasem
1437                elif '2*Fo-Fc' in mapData['MapType']:
1438                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1439                    h,k,l = hkl+Hmax
1440                    Fhkl[h,k,l] = F*phasep
1441                    h,k,l = -hkl+Hmax
1442                    Fhkl[h,k,l] = F*phasem
1443                elif 'Patterson' in mapData['MapType']:
1444                    h,k,l = hkl+Hmax
1445                    Fhkl[h,k,l] = complex(Fosq,0.)
1446                    h,k,l = -hkl+Hmax
1447                    Fhkl[h,k,l] = complex(Fosq,0.)
1448    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1449    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1450    mapData['rho'] = np.real(rho)
1451    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1452    return mapData
1453   
1454# map printing for testing purposes
1455def printRho(SGLaue,rho,rhoMax):                         
1456    '''default doc string
1457   
1458    :param type name: description
1459   
1460    :returns: type name: description
1461   
1462    '''
1463    dim = len(rho.shape)
1464    if dim == 2:
1465        ix,jy = rho.shape
1466        for j in range(jy):
1467            line = ''
1468            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1469                line += (jy-j)*'  '
1470            for i in range(ix):
1471                r = int(100*rho[i,j]/rhoMax)
1472                line += '%4d'%(r)
1473            print line+'\n'
1474    else:
1475        ix,jy,kz = rho.shape
1476        for k in range(kz):
1477            print 'k = ',k
1478            for j in range(jy):
1479                line = ''
1480                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1481                    line += (jy-j)*'  '
1482                for i in range(ix):
1483                    r = int(100*rho[i,j,k]/rhoMax)
1484                    line += '%4d'%(r)
1485                print line+'\n'
1486## keep this
1487               
1488def findOffset(SGData,A,Fhkl):   
1489    '''default doc string
1490   
1491    :param type name: description
1492   
1493    :returns: type name: description
1494   
1495    '''
1496    if SGData['SpGrp'] == 'P 1':
1497        return [0,0,0]   
1498    hklShape = Fhkl.shape
1499    hklHalf = np.array(hklShape)/2
1500    sortHKL = np.argsort(Fhkl.flatten())
1501    Fdict = {}
1502    for hkl in sortHKL:
1503        HKL = np.unravel_index(hkl,hklShape)
1504        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1505        if F == 0.:
1506            break
1507        Fdict['%.6f'%(np.absolute(F))] = hkl
1508    Flist = np.flipud(np.sort(Fdict.keys()))
1509    F = str(1.e6)
1510    i = 0
1511    DH = []
1512    Dphi = []
1513    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1514    for F in Flist:
1515        hkl = np.unravel_index(Fdict[F],hklShape)
1516        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1517        Uniq = np.array(Uniq,dtype='i')
1518        Phi = np.array(Phi)
1519        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1520        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1521        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1522        ang0 = np.angle(Fh0,deg=True)/360.
1523        for H,phi in zip(Uniq,Phi)[1:]:
1524            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1525            dH = H-hkl
1526            dang = ang-ang0
1527            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
1528                continue
1529            DH.append(dH)
1530            Dphi.append((dang+.5) % 1.0)
1531        if i > 20 or len(DH) > 30:
1532            break
1533        i += 1
1534    DH = np.array(DH)
1535    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1536    Dphi = np.array(Dphi)
1537    steps = np.array(hklShape)
1538    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1539    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1540    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1541    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1542    hist,bins = np.histogram(Mmap,bins=1000)
1543#    for i,item in enumerate(hist[:10]):
1544#        print item,bins[i]
1545    chisq = np.min(Mmap)
1546    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1547    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1548#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1549    return DX
1550   
1551def ChargeFlip(data,reflData,pgbar):
1552    '''default doc string
1553   
1554    :param type name: description
1555   
1556    :returns: type name: description
1557   
1558    '''
1559    generalData = data['General']
1560    mapData = generalData['Map']
1561    flipData = generalData['Flip']
1562    FFtable = {}
1563    if 'None' not in flipData['Norm element']:
1564        normElem = flipData['Norm element'].upper()
1565        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1566        for ff in FFs:
1567            if ff['Symbol'] == normElem:
1568                FFtable.update(ff)
1569    dmin = flipData['Resolution']
1570    SGData = generalData['SGData']
1571    cell = generalData['Cell'][1:8]       
1572    A = G2lat.cell2A(cell[:6])
1573    Vol = cell[6]
1574    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1575    adjHKLmax(SGData,Hmax)
1576    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1577    time0 = time.time()
1578    for ref in reflData:
1579        dsp = ref[4]
1580        if dsp >= dmin:
1581            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1582            if FFtable:
1583                SQ = 0.25/dsp**2
1584                ff *= G2el.ScatFac(FFtable,SQ)[0]
1585            if ref[8] > 0.:         #use only +ve Fobs**2
1586                E = np.sqrt(ref[8])/ff
1587            else:
1588                E = 0.
1589            ph = ref[10]
1590            ph = rn.uniform(0.,360.)
1591            for i,hkl in enumerate(ref[11]):
1592                hkl = np.asarray(hkl,dtype='i')
1593                dp = 360.*ref[12][i]
1594                a = cosd(ph+dp)
1595                b = sind(ph+dp)
1596                phasep = complex(a,b)
1597                phasem = complex(a,-b)
1598                h,k,l = hkl+Hmax
1599                Ehkl[h,k,l] = E*phasep
1600                h,k,l = -hkl+Hmax       #Friedel pair refl.
1601                Ehkl[h,k,l] = E*phasem
1602#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1603    CEhkl = copy.copy(Ehkl)
1604    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1605    Emask = ma.getmask(MEhkl)
1606    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1607    Ncyc = 0
1608    old = np.seterr(all='raise')
1609    while True:       
1610        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1611        CEsig = np.std(CErho)
1612        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1613        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
1614        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1615        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1616        phase = CFhkl/np.absolute(CFhkl)
1617        CEhkl = np.absolute(Ehkl)*phase
1618        Ncyc += 1
1619        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1620        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1621        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1622        if Rcf < 5.:
1623            break
1624        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1625        if not GoOn or Ncyc > 10000:
1626            break
1627    np.seterr(**old)
1628    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1629    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1630    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1631    roll = findOffset(SGData,A,CEhkl)
1632       
1633    mapData['Rcf'] = Rcf
1634    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1635    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1636    mapData['rollMap'] = [0,0,0]
1637    return mapData
1638   
1639def SearchMap(data):
1640    '''Does a search of a density map for peaks meeting the criterion of peak
1641    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
1642    mapData is data['General']['mapData']; the map is also in mapData.
1643
1644    :param data: the phase data structure
1645
1646    :returns: (peaks,mags,dzeros) where
1647
1648        * peaks : ndarray
1649          x,y,z positions of the peaks found in the map
1650        * mags : ndarray
1651          the magnitudes of the peaks
1652        * dzeros : ndarray
1653          the distance of the peaks from  the unit cell origin
1654
1655    '''       
1656    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1657   
1658    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
1659   
1660    def noDuplicate(xyz,peaks,Amat):
1661        XYZ = np.inner(Amat,xyz)
1662        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1663            print ' Peak',xyz,' <0.5A from another peak'
1664            return False
1665        return True
1666                           
1667    def fixSpecialPos(xyz,SGData,Amat):
1668        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1669        X = []
1670        xyzs = [equiv[0] for equiv in equivs]
1671        for x in xyzs:
1672            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1673                X.append(x)
1674        if len(X) > 1:
1675            return np.average(X,axis=0)
1676        else:
1677            return xyz
1678       
1679    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1680        Mag,x0,y0,z0,sig = parms
1681        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
1682#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
1683        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
1684       
1685    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1686        Mag,x0,y0,z0,sig = parms
1687        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1688        return M
1689       
1690    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1691        Mag,x0,y0,z0,sig = parms
1692        dMdv = np.zeros(([5,]+list(rX.shape)))
1693        delt = .01
1694        for i in range(5):
1695            parms[i] -= delt
1696            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1697            parms[i] += 2.*delt
1698            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1699            parms[i] -= delt
1700            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1701        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1702        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1703        dMdv = np.reshape(dMdv,(5,rX.size))
1704        Hess = np.inner(dMdv,dMdv)
1705       
1706        return Vec,Hess
1707       
1708    generalData = data['General']
1709    phaseName = generalData['Name']
1710    SGData = generalData['SGData']
1711    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1712    drawingData = data['Drawing']
1713    peaks = []
1714    mags = []
1715    dzeros = []
1716    try:
1717        mapData = generalData['Map']
1718        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1719        rho = copy.copy(mapData['rho'])     #don't mess up original
1720        mapHalf = np.array(rho.shape)/2
1721        res = mapData['Resolution']
1722        incre = np.array(rho.shape,dtype=np.float)
1723        step = max(1.0,1./res)+1
1724        steps = np.array(3*[step,])
1725    except KeyError:
1726        print '**** ERROR - Fourier map not defined'
1727        return peaks,mags
1728    rhoMask = ma.array(rho,mask=(rho<contLevel))
1729    indices = (-1,0,1)
1730    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1731    for roll in rolls:
1732        if np.any(roll):
1733            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1734    indx = np.transpose(rhoMask.nonzero())
1735    peaks = indx/incre
1736    mags = rhoMask[rhoMask.nonzero()]
1737    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1738        rho = rollMap(rho,ind)
1739        rMM = mapHalf-steps
1740        rMP = mapHalf+steps+1
1741        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1742        peakInt = np.sum(rhoPeak)*res**3
1743        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1744        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
1745        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1746            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1747        x1 = result[0]
1748        if not np.any(x1 < 0):
1749            mag = x1[0]
1750            peak = (np.array(x1[1:4])-ind)/incre
1751        peak = fixSpecialPos(peak,SGData,Amat)
1752        rho = rollMap(rho,-ind)       
1753    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1754    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1755   
1756def sortArray(data,pos,reverse=False):
1757    '''data is a list of items
1758    sort by pos in list; reverse if True
1759    '''
1760    T = []
1761    for i,M in enumerate(data):
1762        T.append((M[pos],i))
1763    D = dict(zip(T,data))
1764    T.sort()
1765    if reverse:
1766        T.reverse()
1767    X = []
1768    for key in T:
1769        X.append(D[key])
1770    return X
1771
1772def PeaksEquiv(data,Ind):
1773    '''Find the equivalent map peaks for those selected. Works on the
1774    contents of data['Map Peaks'].
1775
1776    :param data: the phase data structure
1777    :param list Ind: list of selected peak indices
1778    :returns: augmented list of peaks including those related by symmetry to the
1779      ones in Ind
1780
1781    '''       
1782    def Duplicate(xyz,peaks,Amat):
1783        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1784            return True
1785        return False
1786                           
1787    generalData = data['General']
1788    cell = generalData['Cell'][1:7]
1789    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1790    A = G2lat.cell2A(cell)
1791    SGData = generalData['SGData']
1792    mapPeaks = data['Map Peaks']
1793    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1794    Indx = {}
1795    for ind in Ind:
1796        xyz = np.array(mapPeaks[ind][1:4])
1797        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1798#        for x in xyzs: print x
1799        for jnd,xyz in enumerate(XYZ):       
1800            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1801    Ind = []
1802    for ind in Indx:
1803        if Indx[ind]:
1804            Ind.append(ind)
1805    return Ind
1806               
1807def PeaksUnique(data,Ind):
1808    '''Finds the symmetry unique set of peaks from those selected. Works on the
1809    contents of data['Map Peaks'].
1810
1811    :param data: the phase data structure
1812    :param list Ind: list of selected peak indices
1813    :returns: the list of symmetry unique peaks from among those given in Ind
1814
1815    '''       
1816#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1817
1818    def noDuplicate(xyz,peaks,Amat):
1819        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1820            return False
1821        return True
1822                           
1823    generalData = data['General']
1824    cell = generalData['Cell'][1:7]
1825    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1826    A = G2lat.cell2A(cell)
1827    SGData = generalData['SGData']
1828    mapPeaks = data['Map Peaks']
1829    Indx = {}
1830    XYZ = {}
1831    for ind in Ind:
1832        XYZ[ind] = np.array(mapPeaks[ind][1:4])
1833        Indx[ind] = True
1834    for ind in Ind:
1835        if Indx[ind]:
1836            xyz = XYZ[ind]
1837            for jnd in Ind:
1838                if ind != jnd and Indx[jnd]:                       
1839                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1840                    xyzs = np.array([equiv[0] for equiv in Equiv])
1841                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1842    Ind = []
1843    for ind in Indx:
1844        if Indx[ind]:
1845            Ind.append(ind)
1846    return Ind
1847   
1848################################################################################
1849##### single peak fitting profile fxn stuff
1850################################################################################
1851
1852def getCWsig(ins,pos):
1853    '''default doc string
1854   
1855    :param type name: description
1856   
1857    :returns: type name: description
1858   
1859    '''
1860    tp = tand(pos/2.0)
1861    return ins['U']*tp**2+ins['V']*tp+ins['W']
1862   
1863def getCWsigDeriv(pos):
1864    '''default doc string
1865   
1866    :param type name: description
1867   
1868    :returns: type name: description
1869   
1870    '''
1871    tp = tand(pos/2.0)
1872    return tp**2,tp,1.0
1873   
1874def getCWgam(ins,pos):
1875    '''default doc string
1876   
1877    :param type name: description
1878   
1879    :returns: type name: description
1880   
1881    '''
1882    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
1883   
1884def getCWgamDeriv(pos):
1885    '''default doc string
1886   
1887    :param type name: description
1888   
1889    :returns: type name: description
1890   
1891    '''
1892    return 1./cosd(pos/2.0),tand(pos/2.0)
1893   
1894def getTOFsig(ins,dsp):
1895    '''default doc string
1896   
1897    :param type name: description
1898   
1899    :returns: type name: description
1900   
1901    '''
1902    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-q']*dsp
1903   
1904def getTOFsigDeriv(dsp):
1905    '''default doc string
1906   
1907    :param type name: description
1908   
1909    :returns: type name: description
1910   
1911    '''
1912    return 1.0,dsp**2,dsp
1913   
1914def getTOFgamma(ins,dsp):
1915    '''default doc string
1916   
1917    :param type name: description
1918   
1919    :returns: type name: description
1920   
1921    '''
1922    return ins['X']*dsp+ins['Y']*dsp**2
1923   
1924def getTOFgammaDeriv(dsp):
1925    '''default doc string
1926   
1927    :param type name: description
1928   
1929    :returns: type name: description
1930   
1931    '''
1932    return dsp,dsp**2
1933   
1934def getTOFbeta(ins,dsp):
1935    '''default doc string
1936   
1937    :param type name: description
1938   
1939    :returns: type name: description
1940   
1941    '''
1942    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp
1943   
1944def getTOFbetaDeriv(dsp):
1945    '''default doc string
1946   
1947    :param type name: description
1948   
1949    :returns: type name: description
1950   
1951    '''
1952    return 1.0,1./dsp**4,1./dsp
1953   
1954def getTOFalpha(ins,dsp):
1955    '''default doc string
1956   
1957    :param type name: description
1958   
1959    :returns: type name: description
1960   
1961    '''
1962    return ins['alpha']/dsp
1963   
1964def getTOFalphaDeriv(dsp):
1965    '''default doc string
1966   
1967    :param type name: description
1968   
1969    :returns: type name: description
1970   
1971    '''
1972    return 1./dsp
1973   
1974def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1975    '''default doc string
1976   
1977    :param type name: description
1978   
1979    :returns: type name: description
1980   
1981    '''
1982    ind = 0
1983    if useFit:
1984        ind = 1
1985    ins = {}
1986    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
1987        for x in ['U','V','W','X','Y']:
1988            ins[x] = Parms[x][ind]
1989        if ifQ:                              #qplot - convert back to 2-theta
1990            pos = 2.0*asind(pos*wave/(4*math.pi))
1991        sig = getCWsig(ins,pos)
1992        gam = getCWgam(ins,pos)           
1993        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
1994    else:
1995        if ifQ:
1996            dsp = 2.*np.pi/pos
1997            pos = Parms['difC']*dsp
1998        else:
1999            dsp = pos/Parms['difC'][1]
2000        if 'Pdabc' in Parms2:
2001            for x in ['sig-0','sig-1','sig-q','X','Y']:
2002                ins[x] = Parms[x][ind]
2003            Pdabc = Parms2['Pdabc'].T
2004            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2005            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2006        else:
2007            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y']:
2008                ins[x] = Parms[x][ind]
2009            alp = getTOFalpha(ins,dsp)
2010            bet = getTOFbeta(ins,dsp)
2011        sig = getTOFsig(ins,dsp)
2012        gam = getTOFgamma(ins,dsp)
2013        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2014    return XY
2015   
2016################################################################################
2017##### MC/SA stuff
2018################################################################################
2019
2020#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2021# Original Author: Travis Oliphant 2002
2022# Bug-fixes in 2006 by Tim Leslie
2023
2024
2025import numpy
2026from numpy import asarray, tan, exp, ones, squeeze, sign, \
2027     all, log, sqrt, pi, shape, array, minimum, where
2028from numpy import random
2029
2030#__all__ = ['anneal']
2031
2032_double_min = numpy.finfo(float).min
2033_double_max = numpy.finfo(float).max
2034class base_schedule(object):
2035    def __init__(self):
2036        self.dwell = 20
2037        self.learn_rate = 0.5
2038        self.lower = -10
2039        self.upper = 10
2040        self.Ninit = 50
2041        self.accepted = 0
2042        self.tests = 0
2043        self.feval = 0
2044        self.k = 0
2045        self.T = None
2046
2047    def init(self, **options):
2048        self.__dict__.update(options)
2049        self.lower = asarray(self.lower)
2050        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2051        self.upper = asarray(self.upper)
2052        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2053        self.k = 0
2054        self.accepted = 0
2055        self.feval = 0
2056        self.tests = 0
2057
2058    def getstart_temp(self, best_state):
2059        """ Find a matching starting temperature and starting parameters vector
2060        i.e. find x0 such that func(x0) = T0.
2061
2062        Parameters
2063        ----------
2064        best_state : _state
2065            A _state object to store the function value and x0 found.
2066
2067        returns
2068        -------
2069        x0 : array
2070            The starting parameters vector.
2071        """
2072
2073        assert(not self.dims is None)
2074        lrange = self.lower
2075        urange = self.upper
2076        fmax = _double_min
2077        fmin = _double_max
2078        for _ in range(self.Ninit):
2079            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2080            fval = self.func(x0, *self.args)
2081            self.feval += 1
2082            if fval > fmax:
2083                fmax = fval
2084            if fval < fmin:
2085                fmin = fval
2086                best_state.cost = fval
2087                best_state.x = array(x0)
2088
2089        self.T0 = (fmax-fmin)*1.5
2090        return best_state.x
2091
2092    def accept_test(self, dE):
2093        T = self.T
2094        self.tests += 1
2095        if dE < 0:
2096            self.accepted += 1
2097            return 1
2098        p = exp(-dE*1.0/self.boltzmann/T)
2099        if (p > random.uniform(0.0, 1.0)):
2100            self.accepted += 1
2101            return 1
2102        return 0
2103
2104    def update_guess(self, x0):
2105        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2106
2107    def update_temp(self, x0):
2108        pass
2109
2110
2111#  A schedule due to Lester Ingber modified to use bounds - OK
2112class fast_sa(base_schedule):
2113    def init(self, **options):
2114        self.__dict__.update(options)
2115        if self.m is None:
2116            self.m = 1.0
2117        if self.n is None:
2118            self.n = 1.0
2119        self.c = self.m * exp(-self.n * self.quench)
2120
2121    def update_guess(self, x0):
2122        x0 = asarray(x0)
2123        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2124        T = self.T
2125        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2126        xnew = xc*(self.upper - self.lower)+self.lower
2127        return xnew
2128#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2129#        xc = y*(self.upper - self.lower)
2130#        xnew = x0 + xc
2131#        return xnew
2132
2133    def update_temp(self):
2134        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2135        self.k += 1
2136        return
2137
2138class cauchy_sa(base_schedule):     #modified to use bounds - not good
2139    def update_guess(self, x0):
2140        x0 = asarray(x0)
2141        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2142        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2143        xnew = xc*(self.upper - self.lower)+self.lower
2144        return xnew
2145#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2146#        xc = self.learn_rate * self.T * tan(numbers)
2147#        xnew = x0 + xc
2148#        return xnew
2149
2150    def update_temp(self):
2151        self.T = self.T0/(1+self.k)
2152        self.k += 1
2153        return
2154
2155class boltzmann_sa(base_schedule):
2156#    def update_guess(self, x0):
2157#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2158#        x0 = asarray(x0)
2159#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2160#
2161#        xnew = x0 + xc*std*self.learn_rate
2162#        return xnew
2163
2164    def update_temp(self):
2165        self.k += 1
2166        self.T = self.T0 / log(self.k+1.0)
2167        return
2168
2169class log_sa(base_schedule):        #OK
2170
2171    def init(self,**options):
2172        self.__dict__.update(options)
2173       
2174#    def update_guess(self,x0):
2175#        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2176       
2177    def update_temp(self):
2178        self.k += 1
2179        self.T = self.T0*self.slope**self.k
2180       
2181class _state(object):
2182    def __init__(self):
2183        self.x = None
2184        self.cost = None
2185
2186# TODO:
2187#     allow for general annealing temperature profile
2188#     in that case use update given by alpha and omega and
2189#     variation of all previous updates and temperature?
2190
2191# Simulated annealing
2192
2193def anneal(func, x0, args=(), schedule='fast', full_output=0,
2194           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2195           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2196           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=True,dlg=None):
2197    """Minimize a function using simulated annealing.
2198
2199    Schedule is a schedule class implementing the annealing schedule.
2200    Available ones are 'fast', 'cauchy', 'boltzmann'
2201
2202    :param callable func: f(x, \*args)
2203        Function to be optimized.
2204    :param ndarray x0:
2205        Initial guess.
2206    :param tuple args:
2207        Extra parameters to `func`.
2208    :param base_schedule schedule:
2209        Annealing schedule to use (a class).
2210    :param bool full_output:
2211        Whether to return optional outputs.
2212    :param float T0:
2213        Initial Temperature (estimated as 1.2 times the largest
2214        cost-function deviation over random points in the range).
2215    :param float Tf:
2216        Final goal temperature.
2217    :param int maxeval:
2218        Maximum function evaluations.
2219    :param int maxaccept:
2220        Maximum changes to accept.
2221    :param int maxiter:
2222        Maximum cooling iterations.
2223    :param float learn_rate:
2224        Scale constant for adjusting guesses.
2225    :param float boltzmann:
2226        Boltzmann constant in acceptance test
2227        (increase for less stringent test at each temperature).
2228    :param float feps:
2229        Stopping relative error tolerance for the function value in
2230        last four coolings.
2231    :param float quench,m,n:
2232        Parameters to alter fast_sa schedule.
2233    :param float/ndarray lower,upper:
2234        Lower and upper bounds on `x`.
2235    :param int dwell:
2236        The number of times to search the space at each temperature.
2237    :param float slope:
2238        Parameter for log schedule
2239    :param bool ranStart=True:
2240        False for fixed point start
2241
2242    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2243
2244     * xmin (ndarray): Point giving smallest value found.
2245     * Jmin (float): Minimum value of function found.
2246     * T (float): Final temperature.
2247     * feval (int): Number of function evaluations.
2248     * iters (int): Number of cooling iterations.
2249     * accept (int): Number of tests accepted.
2250     * retval (int): Flag indicating stopping condition:
2251
2252              *  0: Points no longer changing
2253              *  1: Cooled to final temperature
2254              *  2: Maximum function evaluations
2255              *  3: Maximum cooling iterations reached
2256              *  4: Maximum accepted query locations reached
2257              *  5: Final point not the minimum amongst encountered points
2258
2259    *Notes*:
2260    Simulated annealing is a random algorithm which uses no derivative
2261    information from the function being optimized. In practice it has
2262    been more useful in discrete optimization than continuous
2263    optimization, as there are usually better algorithms for continuous
2264    optimization problems.
2265
2266    Some experimentation by trying the difference temperature
2267    schedules and altering their parameters is likely required to
2268    obtain good performance.
2269
2270    The randomness in the algorithm comes from random sampling in numpy.
2271    To obtain the same results you can call numpy.random.seed with the
2272    same seed immediately before calling scipy.optimize.anneal.
2273
2274    We give a brief description of how the three temperature schedules
2275    generate new points and vary their temperature. Temperatures are
2276    only updated with iterations in the outer loop. The inner loop is
2277    over xrange(dwell), and new points are generated for
2278    every iteration in the inner loop. (Though whether the proposed
2279    new points are accepted is probabilistic.)
2280
2281    For readability, let d denote the dimension of the inputs to func.
2282    Also, let x_old denote the previous state, and k denote the
2283    iteration number of the outer loop. All other variables not
2284    defined below are input variables to scipy.optimize.anneal itself.
2285
2286    In the 'fast' schedule the updates are ::
2287
2288        u ~ Uniform(0, 1, size=d)
2289        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2290        xc = y * (upper - lower)
2291        x_new = x_old + xc
2292
2293        c = n * exp(-n * quench)
2294        T_new = T0 * exp(-c * k**quench)
2295
2296
2297    In the 'cauchy' schedule the updates are ::
2298
2299        u ~ Uniform(-pi/2, pi/2, size=d)
2300        xc = learn_rate * T * tan(u)
2301        x_new = x_old + xc
2302
2303        T_new = T0 / (1+k)
2304
2305    In the 'boltzmann' schedule the updates are ::
2306
2307        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2308        y ~ Normal(0, std, size=d)
2309        x_new = x_old + learn_rate * y
2310
2311        T_new = T0 / log(1+k)
2312
2313    """
2314    x0 = asarray(x0)
2315    lower = asarray(lower)
2316    upper = asarray(upper)
2317
2318    schedule = eval(schedule+'_sa()')
2319    #   initialize the schedule
2320    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2321                  learn_rate=learn_rate, lower=lower, upper=upper,
2322                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2323
2324    current_state, last_state, best_state = _state(), _state(), _state()
2325    if T0 is None:
2326        x0 = schedule.getstart_temp(best_state)
2327    else:
2328        if ranStart:
2329            x0 = random.uniform(size=len(x0))*(upper-lower) + lower #comment to avoid random start
2330        best_state.x = None
2331        best_state.cost = numpy.Inf
2332
2333    last_state.x = asarray(x0).copy()
2334    fval = func(x0,*args)
2335    schedule.feval += 1
2336    last_state.cost = fval
2337    if last_state.cost < best_state.cost:
2338        best_state.cost = fval
2339        best_state.x = asarray(x0).copy()
2340    schedule.T = schedule.T0
2341    fqueue = [100, 300, 500, 700]
2342    iters = 1
2343    keepGoing = True
2344    bestn = 0
2345    while keepGoing:
2346        retval = 0
2347        for n in xrange(dwell):
2348            current_state.x = schedule.update_guess(last_state.x)
2349            current_state.cost = func(current_state.x,*args)
2350            schedule.feval += 1
2351
2352            dE = current_state.cost - last_state.cost
2353            if schedule.accept_test(dE):
2354                last_state.x = current_state.x.copy()
2355                last_state.cost = current_state.cost
2356                if last_state.cost < best_state.cost:
2357                    best_state.x = last_state.x.copy()
2358                    best_state.cost = last_state.cost
2359                    bestn = n
2360        if dlg:
2361            GoOn = dlg.Update(min(100.,best_state.cost*100),
2362                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2363                    'Best trial:',bestn,  \
2364                    'MC/SA Residual:',best_state.cost*100,'%', \
2365                    ))[0]
2366            if not GoOn:
2367                best_state.x = last_state.x.copy()
2368                best_state.cost = last_state.cost
2369                retval = 5
2370        schedule.update_temp()
2371        iters += 1
2372        # Stopping conditions
2373        # 0) last saved values of f from each cooling step
2374        #     are all very similar (effectively cooled)
2375        # 1) Tf is set and we are below it
2376        # 2) maxeval is set and we are past it
2377        # 3) maxiter is set and we are past it
2378        # 4) maxaccept is set and we are past it
2379        # 5) user canceled run via progress bar
2380
2381        fqueue.append(squeeze(last_state.cost))
2382        fqueue.pop(0)
2383        af = asarray(fqueue)*1.0
2384        if retval == 5:
2385            print ' User terminated run; incomplete MC/SA'
2386            keepGoing = False
2387            break
2388        if all(abs((af-af[0])/af[0]) < feps):
2389            retval = 0
2390            if abs(af[-1]-best_state.cost) > feps*10:
2391                retval = 5
2392#                print "Warning: Cooled to %f at %s but this is not" \
2393#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2394#                      + " the smallest point found."
2395            break
2396        if (Tf is not None) and (schedule.T < Tf):
2397            retval = 1
2398            break
2399        if (maxeval is not None) and (schedule.feval > maxeval):
2400            retval = 2
2401            break
2402        if (iters > maxiter):
2403            print "Warning: Maximum number of iterations exceeded."
2404            retval = 3
2405            break
2406        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2407            retval = 4
2408            break
2409
2410    if full_output:
2411        return best_state.x, best_state.cost, schedule.T, \
2412               schedule.feval, iters, schedule.accepted, retval
2413    else:
2414        return best_state.x, retval
2415
2416def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q):
2417    outlist = []
2418    for n in range(iCyc):
2419        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2420        outlist.append(result[0])
2421        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2422    out_q.put(outlist)
2423
2424def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2425    import multiprocessing as mp
2426   
2427    nprocs = mp.cpu_count()
2428    out_q = mp.Queue()
2429    procs = []
2430    iCyc = np.zeros(nprocs)
2431    for i in range(nCyc):
2432        iCyc[i%nprocs] += 1
2433    for i in range(nprocs):
2434        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q))
2435        procs.append(p)
2436        p.start()
2437    resultlist = []
2438    for i in range(nprocs):
2439        resultlist += out_q.get()
2440    for p in procs:
2441        p.join()
2442    return resultlist
2443
2444def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2445    '''default doc string
2446   
2447    :param type name: description
2448   
2449    :returns: type name: description
2450   
2451    '''
2452   
2453    twopi = 2.0*np.pi
2454    global tsum
2455    tsum = 0.
2456   
2457    def getMDparms(item,pfx,parmDict,varyList):
2458        parmDict[pfx+'MDaxis'] = item['axis']
2459        parmDict[pfx+'MDval'] = item['Coef'][0]
2460        if item['Coef'][1]:
2461            varyList += [pfx+'MDval',]
2462            limits = item['Coef'][2]
2463            lower.append(limits[0])
2464            upper.append(limits[1])
2465                       
2466    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2467        parmDict[pfx+'Atype'] = item['atType']
2468        aTypes |= set([item['atType'],]) 
2469        pstr = ['Ax','Ay','Az']
2470        XYZ = [0,0,0]
2471        for i in range(3):
2472            name = pfx+pstr[i]
2473            parmDict[name] = item['Pos'][0][i]
2474            XYZ[i] = parmDict[name]
2475            if item['Pos'][1][i]:
2476                varyList += [name,]
2477                limits = item['Pos'][2][i]
2478                lower.append(limits[0])
2479                upper.append(limits[1])
2480        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2481           
2482    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2483        parmDict[mfx+'MolCent'] = item['MolCent']
2484        parmDict[mfx+'RBId'] = item['RBId']
2485        pstr = ['Px','Py','Pz']
2486        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
2487        for i in range(3):
2488            name = mfx+pstr[i]
2489            parmDict[name] = item['Pos'][0][i]
2490            if item['Pos'][1][i]:
2491                varyList += [name,]
2492                limits = item['Pos'][2][i]
2493                lower.append(limits[0])
2494                upper.append(limits[1])
2495        AV = item['Ori'][0]
2496        A = AV[0]
2497        V = AV[1:]
2498        for i in range(4):
2499            name = mfx+ostr[i]
2500            if i:
2501                parmDict[name] = V[i-1]
2502            else:
2503                parmDict[name] = A
2504            if item['Ovar'] == 'AV':
2505                varyList += [name,]
2506                limits = item['Ori'][2][i]
2507                lower.append(limits[0])
2508                upper.append(limits[1])
2509            elif item['Ovar'] == 'A' and not i:
2510                varyList += [name,]
2511                limits = item['Ori'][2][i]
2512                lower.append(limits[0])
2513                upper.append(limits[1])
2514        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2515            for i in range(len(item['Tor'][0])):
2516                name = mfx+'Tor'+str(i)
2517                parmDict[name] = item['Tor'][0][i]
2518                if item['Tor'][1][i]:
2519                    varyList += [name,]
2520                    limits = item['Tor'][2][i]
2521                    lower.append(limits[0])
2522                    upper.append(limits[1])
2523        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2524        aTypes |= set(atypes)
2525        atNo += len(atypes)
2526        return atNo
2527               
2528    def GetAtomM(Xdata,SGData):
2529        Mdata = []
2530        for xyz in Xdata:
2531            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2532        return np.array(Mdata)
2533       
2534    def GetAtomTX(RBdata,parmDict):
2535        'Needs a doc string'
2536        Bmat = parmDict['Bmat']
2537        atNo = parmDict['atNo']
2538        nfixAt = parmDict['nfixAt']
2539        Tdata = atNo*[' ',]
2540        Xdata = np.zeros((3,atNo))
2541        keys = {':Atype':Tdata,':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2542        for iatm in range(nfixAt):
2543            for key in keys:
2544                parm = ':'+str(iatm)+key
2545                if parm in parmDict:
2546                    if key == ':Atype':
2547                        keys[key][iatm] = aTypes.index(parmDict[parm])
2548                    else:
2549                        keys[key][iatm] = parmDict[parm]
2550        iatm = nfixAt
2551        for iObj in range(parmDict['nObj']):
2552            pfx = str(iObj)+':'
2553            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2554                if parmDict[pfx+'Type'] == 'Vector':
2555                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2556                    vecs = RBRes['rbVect']
2557                    mags = RBRes['VectMag']
2558                    Cart = np.zeros_like(vecs[0])
2559                    for vec,mag in zip(vecs,mags):
2560                        Cart += vec*mag
2561                elif parmDict[pfx+'Type'] == 'Residue':
2562                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2563                    Cart = np.array(RBRes['rbXYZ'])
2564                    for itor,seq in enumerate(RBRes['rbSeq']):
2565                        tName = pfx+'Tor'+str(itor)
2566                        QuatA = AVdeg2Q(parmDict[tName],Cart[seq[0]]-Cart[seq[1]])
2567                        for ride in seq[3]:
2568                            Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
2569                if parmDict[pfx+'MolCent'][1]:
2570                    Cart -= parmDict[pfx+'MolCent'][0]
2571                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2572                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2573                for i,x in enumerate(Cart):
2574                    X = np.inner(Bmat,prodQVQ(Qori,x))+Pos
2575                    for j in range(3):
2576                        Xdata[j][iatm] = X[j]
2577                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
2578                    iatm += 1
2579            elif parmDict[pfx+'Type'] == 'Atom':
2580                atNo = parmDict[pfx+'atNo']
2581                for key in keys:
2582                    parm = pfx+key[1:]              #remove extra ':'
2583                    if parm in parmDict:
2584                        if key == ':Atype':
2585                            keys[key][atNo] = aTypes.index(parmDict[parm])
2586                        else:
2587                            keys[key][atNo] = parmDict[parm]
2588                iatm += 1
2589            else:
2590                continue        #skips March Dollase
2591        return Tdata,Xdata.T
2592       
2593    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
2594        allX = np.inner(Xdata,SGM)+SGT
2595        allT = np.repeat(Tdata,allX.shape[1])
2596        allM = np.repeat(Mdata,allX.shape[1])
2597        allX = np.reshape(allX,(-1,3))
2598        return allT,allM,allX
2599
2600    def getAllX(Xdata,SGM,SGT):
2601        allX = np.inner(Xdata,SGM)+SGT
2602        allX = np.reshape(allX,(-1,3))
2603        return allX
2604       
2605    def normQuaternions(RBdata,parmDict,varyList,values):
2606        for iObj in range(parmDict['nObj']):
2607            pfx = str(iObj)+':'
2608            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2609                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2610                A,V = Q2AVdeg(Qori)
2611                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
2612                    if i:
2613                        parmDict[pfx+name] = V[i-1]
2614                    else:
2615                        parmDict[pfx+name] = A
2616       
2617    def mcsaCalc(values,refList,rcov,ifInv,allFF,RBdata,varyList,parmDict):
2618        ''' Compute structure factors for all h,k,l for phase
2619        input:
2620            refList: [ref] where each ref = h,k,l,m,d,...
2621            rcov:   array[nref,nref] covariance terms between Fo^2 values
2622            ifInv:  bool True if centrosymmetric
2623            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
2624            RBdata: [dict] rigid body dictionary
2625            varyList: [list] names of varied parameters in MC/SA (not used here)           
2626            ParmDict: [dict] problem parameters
2627        puts result F^2 in each ref[5] in refList
2628        returns:
2629            delt-F*rcov*delt-F/sum(Fo^2)^2
2630        '''       
2631        global tsum
2632        t0 = time.time()
2633        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
2634        Xdata = GetAtomTX(RBdata,parmDict)[1]                   #get new atom coords from RB
2635        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
2636        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
2637        MDaxis = parmDict['0:MDaxis']
2638        Gmat = parmDict['Gmat']
2639        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
2640        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
2641        if ifInv:
2642            refList[5] = Aterm
2643        else:
2644            Bterm = refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
2645            refList[5] = Aterm+Bterm
2646        sumFcsq = np.sum(refList[5])
2647        scale = parmDict['sumFosq']/sumFcsq
2648        refList[5] *= scale
2649        refList[6] = refList[4]-refList[5]
2650        M = np.inner(refList[6],np.inner(rcov,refList[6]))
2651        tsum += (time.time()-t0)
2652#        print M,parmDict['sumFosq'],np.sum(refList[6]**2),np.sum(refList[4]**2)
2653#        print np.sum(refList[6]**2)/np.sum(refList[4]**2)
2654        return M/np.sum(refList[4]**2)
2655
2656    sq8ln2 = np.sqrt(8*np.log(2))
2657    sq2pi = np.sqrt(2*np.pi)
2658    sq4pi = np.sqrt(4*np.pi)
2659    generalData = data['General']
2660    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2661    Gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])[0]
2662    SGData = generalData['SGData']
2663    SGM = [SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))]
2664    SGT = [SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))]
2665    fixAtoms = data['Atoms']                       #if any
2666    cx,ct,cs = generalData['AtomPtrs'][:3]
2667    aTypes = set([])
2668    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2669    varyList = []
2670    atNo = 0
2671    for atm in fixAtoms:
2672        pfx = ':'+str(atNo)+':'
2673        parmDict[pfx+'Atype'] = atm[ct]
2674        aTypes |= set([atm[ct],])
2675        pstr = ['Ax','Ay','Az']
2676        parmDict[pfx+'Amul'] = atm[cs+1]
2677        for i in range(3):
2678            name = pfx+pstr[i]
2679            parmDict[name] = atm[cx+i]
2680        atNo += 1
2681    parmDict['nfixAt'] = len(fixAtoms)       
2682    MCSA = generalData['MCSA controls']
2683    reflName = MCSA['Data source']
2684    phaseName = generalData['Name']
2685    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2686    upper = []
2687    lower = []
2688    for i,item in enumerate(MCSAObjs):
2689        mfx = str(i)+':'
2690        parmDict[mfx+'Type'] = item['Type']
2691        if item['Type'] == 'MD':
2692            getMDparms(item,mfx,parmDict,varyList)
2693        elif item['Type'] == 'Atom':
2694            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2695            parmDict[mfx+'atNo'] = atNo
2696            atNo += 1
2697        elif item['Type'] in ['Residue','Vector']:
2698            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2699    parmDict['atNo'] = atNo                 #total no. of atoms
2700    parmDict['nObj'] = len(MCSAObjs)
2701    aTypes = list(aTypes)
2702    Tdata,Xdata = GetAtomTX(RBdata,parmDict)
2703    Mdata = GetAtomM(Xdata,SGData)
2704    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
2705    FFtables = G2el.GetFFtable(aTypes)
2706    refs = []
2707    allFF = []
2708    sumFosq = 0
2709    if 'PWDR' in reflName:
2710        for ref in reflData:
2711            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2712            if d >= MCSA['dmin']:
2713                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2714                SQ = 0.25/d**2
2715                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2716                refs.append([h,k,l,m,f*m,pos,sig])
2717                sumFosq += f*m
2718        nRef = len(refs)
2719        rcov = np.zeros((nRef,nRef))
2720        for iref,refI in enumerate(refs):
2721            rcov[iref][iref] = 1./(sq4pi*refI[6])
2722            for jref,refJ in enumerate(refs[:iref]):
2723                t1 = refI[6]**2+refJ[6]**2
2724                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2725                if t2 > 10.:
2726                    rcov[iref][jref] = 0.
2727                else:
2728                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2729        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2730        Rdiag = np.sqrt(np.diag(rcov))
2731        Rnorm = np.outer(Rdiag,Rdiag)
2732        rcov /= Rnorm
2733    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
2734        covMatrix = covData['covMatrix']
2735        vList = covData['varyList']
2736        for iref,refI in enumerate(reflData):
2737            h,k,l,m,d,v,f,s = refI
2738            if d >= MCSA['dmin'] and v:       #skip unrefined ones
2739                SQ = 0.25/d**2
2740                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2741                refs.append([h,k,l,m,f*m,iref,0.])
2742                sumFosq += f*m
2743        nRef = len(refs)
2744        pfx = str(data['pId'])+'::PWLref:'
2745        rcov = np.zeros((nRef,nRef))       
2746        for iref,refI in enumerate(refs):
2747            I = refI[5]
2748            nameI = pfx+str(I)
2749            if nameI in vList:
2750                Iindx = vList.index(nameI)
2751                rcov[iref][iref] = covMatrix[Iindx][Iindx]
2752                for jref,refJ in enumerate(refs[:iref]):
2753                    J = refJ[5]
2754                    nameJ = pfx+str(J)
2755                    try:
2756                        rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2757                    except ValueError:
2758                        rcov[iref][jref] = rcov[iref][jref-1]
2759            else:
2760                rcov[iref] = rcov[iref-1]
2761                rcov[iref][iref] = rcov[iref-1][iref-1]
2762        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2763        Rdiag = np.sqrt(np.diag(rcov))
2764        Rnorm = np.outer(Rdiag,Rdiag)
2765        rcov /= Rnorm
2766    elif 'HKLF' in reflName:
2767        for ref in reflData:
2768            [h,k,l,m,d],f = ref[:5],ref[6]
2769            if d >= MCSA['dmin']:
2770                SQ = 0.25/d**2
2771                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2772                refs.append([h,k,l,m,f*m,0.,0.])
2773                sumFosq += f*m
2774        nRef = len(refs)
2775        rcov = np.identity(len(refs))
2776    allFF = np.array(allFF).T
2777    refs = np.array(refs).T
2778    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
2779    print ' Number of parameters varied: %d'%(len(varyList))
2780    parmDict['sumFosq'] = sumFosq
2781    x0 = [parmDict[val] for val in varyList]
2782    ifInv = SGData['SGInv']
2783    results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,allFF,RBdata,varyList,parmDict),
2784        schedule=MCSA['Algorithm'], full_output=True,
2785        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2786        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
2787        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2788        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',True),dlg=pgbar)
2789    M = mcsaCalc(results[0],refs,rcov,ifInv,allFF,RBdata,varyList,parmDict)
2790#    for ref in refs.T:
2791#        print ' %4d %4d %4d %10.3f %10.3f %10.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[4],ref[5],ref[6])
2792#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
2793    Result = [False,False,results[1],results[2],]+list(results[0])
2794    Result.append(varyList)
2795    return Result,tsum
2796
2797       
2798################################################################################
2799##### Quaternion stuff
2800################################################################################
2801
2802def prodQQ(QA,QB):
2803    ''' Grassman quaternion product
2804        QA,QB quaternions; q=r+ai+bj+ck
2805    '''
2806    D = np.zeros(4)
2807    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2808    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2809    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2810    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2811   
2812#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
2813#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
2814   
2815    return D
2816   
2817def normQ(QA):
2818    ''' get length of quaternion & normalize it
2819        q=r+ai+bj+ck
2820    '''
2821    n = np.sqrt(np.sum(np.array(QA)**2))
2822    return QA/n
2823   
2824def invQ(Q):
2825    '''
2826        get inverse of quaternion
2827        q=r+ai+bj+ck; q* = r-ai-bj-ck
2828    '''
2829    return Q*np.array([1,-1,-1,-1])
2830   
2831def prodQVQ(Q,V):
2832    """
2833    compute the quaternion vector rotation qvq-1 = v'
2834    q=r+ai+bj+ck
2835    """
2836    VP = np.zeros(3)
2837    T2 = Q[0]*Q[1]
2838    T3 = Q[0]*Q[2]
2839    T4 = Q[0]*Q[3]
2840    T5 = -Q[1]*Q[1]
2841    T6 = Q[1]*Q[2]
2842    T7 = Q[1]*Q[3]
2843    T8 = -Q[2]*Q[2]
2844    T9 = Q[2]*Q[3]
2845    T10 = -Q[3]*Q[3]
2846    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
2847    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
2848    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
2849    return VP   
2850   
2851def Q2Mat(Q):
2852    ''' make rotation matrix from quaternion
2853        q=r+ai+bj+ck
2854    '''
2855    QN = normQ(Q)
2856    aa = QN[0]**2
2857    ab = QN[0]*QN[1]
2858    ac = QN[0]*QN[2]
2859    ad = QN[0]*QN[3]
2860    bb = QN[1]**2
2861    bc = QN[1]*QN[2]
2862    bd = QN[1]*QN[3]
2863    cc = QN[2]**2
2864    cd = QN[2]*QN[3]
2865    dd = QN[3]**2
2866    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2867        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2868        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2869    return np.array(M)
2870   
2871def AV2Q(A,V):
2872    ''' convert angle (radians) & vector to quaternion
2873        q=r+ai+bj+ck
2874    '''
2875    Q = np.zeros(4)
2876    d = np.sqrt(np.sum(np.array(V)**2))
2877    if d:
2878        V /= d
2879        p = A/2.
2880        Q[0] = np.cos(p)
2881        Q[1:4] = V*np.sin(p)
2882    else:
2883        Q[3] = 1.
2884    return Q
2885   
2886def AVdeg2Q(A,V):
2887    ''' convert angle (degrees) & vector to quaternion
2888        q=r+ai+bj+ck
2889    '''
2890    Q = np.zeros(4)
2891    d = np.sqrt(np.sum(np.array(V)**2))
2892    if d:
2893        V /= d
2894        p = A/2.
2895        Q[0] = cosd(p)
2896        Q[1:4] = V*sind(p)
2897    else:
2898        Q[3] = 1.
2899    return Q
2900   
2901def Q2AVdeg(Q):
2902    ''' convert quaternion to angle (degrees 0-360) & normalized vector
2903        q=r+ai+bj+ck
2904    '''
2905    A = 2.*acosd(Q[0])
2906    V = np.array(Q[1:])
2907    V /= sind(A/2.)
2908    return A,V
2909   
2910def Q2AV(Q):
2911    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
2912        q=r+ai+bj+ck
2913    '''
2914    A = 2.*np.arccos(Q[0])
2915    V = np.array(Q[1:])
2916    V /= np.sin(A/2.)
2917    return A,V
2918   
2919def makeQuat(A,B,C):
2920    ''' Make quaternion from rotation of A vector to B vector about C axis
2921
2922        :param np.array A,B,C: Cartesian 3-vectors
2923        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2924    '''
2925
2926    V1 = np.cross(A,C)
2927    V2 = np.cross(B,C)
2928    if nl.norm(V1)*nl.norm(V2):
2929        V1 /= nl.norm(V1)
2930        V2 /= nl.norm(V2)
2931        V3 = np.cross(V1,V2)
2932    else:
2933        V3 = np.zeros(3)
2934    Q = np.array([0.,0.,0.,1.])
2935    D = 0.
2936    if nl.norm(V3):
2937        V3 /= nl.norm(V3)
2938        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2939        D = np.arccos(D1)/2.0
2940        V1 = C-V3
2941        V2 = C+V3
2942        DM = nl.norm(V1)
2943        DP = nl.norm(V2)
2944        S = np.sin(D)
2945        Q[0] = np.cos(D)
2946        Q[1:] = V3*S
2947        D *= 2.
2948        if DM > DP:
2949            D *= -1.
2950    return Q,D
2951   
2952if __name__ == "__main__":
2953    from numpy import cos
2954    # minimum expected at ~-0.195
2955    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
2956    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
2957    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
2958    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
2959
2960    # minimum expected at ~[-0.195, -0.1]
2961    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
2962    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='cauchy')
2963    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='fast')
2964    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='boltzmann')
Note: See TracBrowser for help on using the repository browser.