source: trunk/GSASIImath.py @ 1055

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

working MC/SA

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