source: trunk/GSASIImath.py @ 1051

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

remove all mcsa fortran code - not needed
replace 32-bit Windows libraries
complete rework of MC/SA structure factor calcs - no loops!
all numpy/python significant speedup >6x over fortran & >100x over
reflection looped numpy/python

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