source: trunk/GSASIImath.py @ 1060

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

mods to MC/SA stuff - some speedup done

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