source: trunk/GSASIImath.py @ 981

Last change on this file since 981 was 981, checked in by toby, 8 years ago

introduce regress option; fix esd printing; more docs; new Mac app with drag & drop for open; control reset of ref list on load

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