source: trunk/GSASIImath.py @ 989

Last change on this file since 989 was 989, checked in by toby, 9 years ago

add constraint derivs for single xtal; allow cellFill to work without esds; update docs; fix controls copy bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 96.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-07-14 03:21:12 +0000 (Sun, 14 Jul 2013) $
5# $Author: toby $
6# $Revision: 989 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 989 2013-07-14 03:21:12Z toby $
9########### SVN repository information ###################
10'''
11*GSASIImath: computation module*
12================================
13
14Routines for least-squares minimization and other stuff
15
16'''
17import sys
18import os
19import os.path as ospath
20import random as rn
21import numpy as np
22import numpy.linalg as nl
23import numpy.ma as ma
24import cPickle
25import time
26import math
27import copy
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 989 $")
30import GSASIIElem as G2el
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import numpy.fft as fft
34import pypowder as pyd
35
36sind = lambda x: np.sin(x*np.pi/180.)
37cosd = lambda x: np.cos(x*np.pi/180.)
38tand = lambda x: np.tan(x*np.pi/180.)
39asind = lambda x: 180.*np.arcsin(x)/np.pi
40acosd = lambda x: 180.*np.arccos(x)/np.pi
41atand = lambda x: 180.*np.arctan(x)/np.pi
42atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
43
44def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,Print=False):
45   
46    """
47    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
48    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`
49   
50    ::
51   
52                    Nobs
53        x = arg min(sum(func(y)**2,axis=0))
54                    y=0
55
56    :param function func: callable method or function
57        should take at least one (possibly length N vector) argument and
58        returns M floating point numbers.
59    :param np.ndarray x0: The starting estimate for the minimization of length N
60    :param function Hess: callable method or function
61        A required function or method to compute the weighted vector and Hessian for func.
62        It must be a symmetric NxN array
63    :param tuple args: Any extra arguments to func are placed in this tuple.
64    :param float ftol: Relative error desired in the sum of squares.
65    :param float xtol: Relative error desired in the approximate solution.
66    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
67        until other limits are met (ftol, xtol)
68    :param bool Print: True for printing results (residuals & times) by cycle
69
70    :returns: (x,cov_x,infodict) where
71
72      * x : ndarray
73        The solution (or the result of the last iteration for an unsuccessful
74        call).
75      * cov_x : ndarray
76        Uses the fjac and ipvt optional outputs to construct an
77        estimate of the jacobian around the solution.  ``None`` if a
78        singular matrix encountered (indicates very flat curvature in
79        some direction).  This matrix must be multiplied by the
80        residual standard deviation to get the covariance of the
81        parameter estimates -- see curve_fit.
82      * infodict : dict
83        a dictionary of optional outputs with the keys:
84
85         * 'fvec' : the function evaluated at the output
86         * 'num cyc':
87         * 'nfev':
88         * 'lamMax':
89         * 'psing':
90           
91    """
92               
93    x0 = np.array(x0, ndmin=1)      #might be redundant?
94    n = len(x0)
95    if type(args) != type(()):
96        args = (args,)
97       
98    icycle = 0
99    One = np.ones((n,n))
100    lam = 0.001
101    lamMax = lam
102    nfev = 0
103    if Print:
104        print ' Hessian refinement on %d variables:'%(n)
105    while icycle < maxcyc:
106        time0 = time.time()
107        M = func(x0,*args)
108        nfev += 1
109        chisq0 = np.sum(M**2)
110        Yvec,Amat = Hess(x0,*args)
111        Adiag = np.sqrt(np.diag(Amat))
112        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
113        if np.any(psing):                #hard singularity in matrix
114            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
115        Anorm = np.outer(Adiag,Adiag)
116        Yvec /= Adiag
117        Amat /= Anorm
118        while True:
119            Lam = np.eye(Amat.shape[0])*lam
120            Amatlam = Amat*(One+Lam)
121            try:
122                Xvec = nl.solve(Amatlam,Yvec)
123            except nl.LinAlgError:
124                print 'ouch #1'
125                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
126                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
127            Xvec /= Adiag
128            M2 = func(x0+Xvec,*args)
129            nfev += 1
130            chisq1 = np.sum(M2**2)
131            if chisq1 > chisq0:
132                lam *= 10.
133            else:
134                x0 += Xvec
135                lam /= 10.
136                break
137            if lam > 10.e5:
138                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
139                break
140        lamMax = max(lamMax,lam)
141        if Print:
142            print ' Cycle: %d, Time: %.2fs, Chi**2: %.3g, Lambda: %.3g'%(icycle,time.time()-time0,chisq1,lam)
143        if (chisq0-chisq1)/chisq0 < ftol:
144            break
145        icycle += 1
146    M = func(x0,*args)
147    nfev += 1
148    Yvec,Amat = Hess(x0,*args)
149    lam = 1.e-8          #unmodified matrix - comment out for matrix modification
150    Lam = np.eye(Amat.shape[0])*lam
151    Amatlam = Amat*(One+Lam)
152    try:
153        Bmat = nl.inv(Amatlam)
154        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
155    except nl.LinAlgError:
156        print 'ouch #2 linear algebra error in LS'
157        psing = []
158        if maxcyc:
159            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
160        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
161   
162def getVCov(varyNames,varyList,covMatrix):
163    '''obtain variance-covariance terms for a set of variables. NB: the varyList
164    and covMatrix were saved by the last least squares refinement so they must match.
165   
166    :param list varyNames: variable names to find v-cov matric for
167    :param list varyList: full list of all variables in v-cov matrix
168    :param nparray covMatrix: full variance-covariance matrix from the last
169     least squares refinement
170   
171    :returns: nparray vcov: variance-covariance matrix for the variables given
172     in varyNames
173   
174    '''
175    vcov = np.zeros((len(varyNames),len(varyNames)))
176    for i1,name1 in enumerate(varyNames):
177        for i2,name2 in enumerate(varyNames):
178            try:
179                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
180            except ValueError:
181                vcov[i1][i2] = 0.0
182    return vcov
183
184def FindAtomIndexByIDs(atomData,IDs,Draw=True):
185    '''finds the set of atom array indices for a list of atom IDs. Will search
186    either the Atom table or the drawAtom table.
187   
188    :param list atomData: Atom or drawAtom table containting coordinates, etc.
189    :param list IDs: atom IDs to be found
190    :param bool Draw: True if drawAtom table to be searched; False if Atom table
191      is searched
192   
193    :returns: list indx: atom (or drawAtom) indices
194   
195    '''
196    indx = []
197    for i,atom in enumerate(atomData):
198        if Draw and atom[-3] in IDs:
199            indx.append(i)
200        elif atom[-1] in IDs:
201            indx.append(i)
202    return indx
203
204def FillAtomLookUp(atomData):
205    '''create a dictionary of atom indexes with atom IDs as keys
206   
207    :param list atomData: Atom table to be used
208   
209    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
210   
211    '''
212    atomLookUp = {}
213    for iatm,atom in enumerate(atomData):
214        atomLookUp[atom[-1]] = iatm
215    return atomLookUp
216
217def GetAtomsById(atomData,atomLookUp,IdList):
218    '''gets a list of atoms from Atom table that match a set of atom IDs
219   
220    :param list atomData: Atom table to be used
221    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
222    :param list IdList: atom IDs to be found
223   
224    :returns: list atoms: list of atoms found
225   
226    '''
227    atoms = []
228    for id in IdList:
229        atoms.append(atomData[atomLookUp[id]])
230    return atoms
231   
232def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
233    '''gets atom parameters for atoms using atom IDs
234   
235    :param list atomData: Atom table to be used
236    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
237    :param list IdList: atom IDs to be found
238    :param int itemLoc: pointer to desired 1st item in an atom table entry
239    :param int numItems: number of items to be retrieved
240   
241    :returns: type name: description
242   
243    '''
244    Items = []
245    if not isinstance(IdList,list):
246        IdList = [IdList,]
247    for id in IdList:
248        if numItems == 1:
249            Items.append(atomData[atomLookUp[id]][itemLoc])
250        else:
251            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
252    return Items
253   
254def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
255    '''default doc string
256   
257    :param type name: description
258   
259    :returns: type name: description
260   
261    '''
262    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
263    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
264    XYZ = []
265    for ind in indx:
266        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
267        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
268        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
269    return XYZ
270
271def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
272    '''default doc string
273   
274    :param type name: description
275   
276    :returns: type name: description
277   
278    '''
279    for atom in atomData:
280        XYZ = np.inner(Amat,atom[cx:cx+3])
281        if atom[cia] == 'A':
282            UIJ = atom[cia+2:cia+8]
283               
284def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
285    '''default doc string
286   
287    :param type name: description
288   
289    :returns: type name: description
290   
291    '''
292    TLStype,TLS = rbObj['ThermalMotion'][:2]
293    Tmat = np.zeros((3,3))
294    Lmat = np.zeros((3,3))
295    Smat = np.zeros((3,3))
296    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
297        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
298    if 'T' in TLStype:
299        Tmat = G2lat.U6toUij(TLS[:6])
300    if 'L' in TLStype:
301        Lmat = G2lat.U6toUij(TLS[6:12])
302    if 'S' in TLStype:
303        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
304    XYZ = np.inner(Amat,xyz)
305    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
306    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
307    beta = np.inner(np.inner(g,Umat),g)
308    return G2lat.UijtoU6(beta)*gvec   
309       
310def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
311    '''default doc string
312   
313    :param type name: description
314   
315    :returns: type name: description
316   
317    '''
318    cx,ct,cs,cia = atPtrs
319    TLStype,TLS = rbObj['ThermalMotion'][:2]
320    Tmat = np.zeros((3,3))
321    Lmat = np.zeros((3,3))
322    Smat = np.zeros((3,3))
323    G,g = G2lat.A2Gmat(Amat)
324    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
325    if 'T' in TLStype:
326        Tmat = G2lat.U6toUij(TLS[:6])
327    if 'L' in TLStype:
328        Lmat = G2lat.U6toUij(TLS[6:12])
329    if 'S' in TLStype:
330        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
331    for atom in atomData:
332        XYZ = np.inner(Amat,atom[cx:cx+3])
333        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
334        if 'U' in TSLtype:
335            atom[cia+1] = TLS[0]
336            atom[cia] = 'I'
337        else:
338            atom[cia] = 'A'
339            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
340            beta = np.inner(np.inner(g,Umat),g)
341            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
342
343def GetXYZDist(xyz,XYZ,Amat):
344    '''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
345        and are in crystal coordinates; Amat is crystal to Cart matrix
346   
347    :param type name: description
348   
349    :returns: type name: description
350   
351    '''
352    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
353
354def getAtomXYZ(atoms,cx):
355    '''default doc string
356   
357    :param type name: description
358   
359    :returns: type name: description
360   
361    '''
362    XYZ = []
363    for atom in atoms:
364        XYZ.append(atom[cx:cx+3])
365    return np.array(XYZ)
366
367def RotateRBXYZ(Bmat,Cart,oriQ):
368    '''rotate & transform cartesian coordinates to crystallographic ones
369    no translation applied. To be used for numerical derivatives
370   
371    :param type name: description
372   
373    :returns: type name: description
374   
375    '''
376    ''' returns crystal coordinates for atoms described by RBObj
377    '''
378    XYZ = np.zeros_like(Cart)
379    for i,xyz in enumerate(Cart):
380        X = prodQVQ(oriQ,xyz)
381        XYZ[i] = np.inner(Bmat,X)
382    return XYZ
383
384def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
385    '''default doc string
386   
387    :param type name: description
388   
389    :returns: type name: description
390   
391    '''
392    ''' returns crystal coordinates for atoms described by RBObj
393    '''
394    RBRes = RBData[RBType][RBObj['RBId']]
395    if RBType == 'Vector':
396        vecs = RBRes['rbVect']
397        mags = RBRes['VectMag']
398        Cart = np.zeros_like(vecs[0])
399        for vec,mag in zip(vecs,mags):
400            Cart += vec*mag
401    elif RBType == 'Residue':
402        Cart = np.array(RBRes['rbXYZ'])
403        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
404            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
405            for ride in seq[3]:
406                Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
407    XYZ = np.zeros_like(Cart)
408    for i,xyz in enumerate(Cart):
409        X = prodQVQ(RBObj['Orient'][0],xyz)
410        XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
411    return XYZ,Cart
412
413def UpdateMCSAxyz(Bmat,MCSA):
414    '''default doc string
415   
416    :param type name: description
417   
418    :returns: type name: description
419   
420    '''
421    xyz = []
422    atTypes = []
423    iatm = 0
424    for model in MCSA['Models'][1:]:        #skip the MD model
425        if model['Type'] == 'Atom':
426            xyz.append(model['Pos'][0])
427            atTypes.append(model['atType'])
428            iatm += 1
429        else:
430            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
431            Pos = np.array(model['Pos'][0])
432            Qori = np.array(model['Ori'][0])
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        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1302    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1303        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1304    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1305        extra = -math.log10(abs(value))
1306        if extra > 0: extra += 1
1307        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1308    if esd > 0:
1309        out += ("({:d})").format(intesd)  # add the esd
1310    elif nTZ and '.' in out:
1311        out = out.rstrip('0')  # strip zeros to right of decimal
1312        out = out.rstrip('.')  # and decimal place when not needed
1313    if valoff != 0:
1314        out += ("e{:d}").format(valoff) # add an exponent, when needed
1315    return out
1316
1317################################################################################
1318##### Fourier & charge flip stuff
1319################################################################################
1320
1321def adjHKLmax(SGData,Hmax):
1322    '''default doc string
1323   
1324    :param type name: description
1325   
1326    :returns: type name: description
1327   
1328    '''
1329    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1330        Hmax[0] = ((Hmax[0]+3)/6)*6
1331        Hmax[1] = ((Hmax[1]+3)/6)*6
1332        Hmax[2] = ((Hmax[2]+1)/4)*4
1333    else:
1334        Hmax[0] = ((Hmax[0]+2)/4)*4
1335        Hmax[1] = ((Hmax[1]+2)/4)*4
1336        Hmax[2] = ((Hmax[2]+1)/4)*4
1337
1338def OmitMap(data,reflData):
1339    '''default doc string
1340   
1341    :param type name: description
1342   
1343    :returns: type name: description
1344   
1345    '''
1346    generalData = data['General']
1347    if not generalData['Map']['MapType']:
1348        print '**** ERROR - Fourier map not defined'
1349        return
1350    mapData = generalData['Map']
1351    dmin = mapData['Resolution']
1352    SGData = generalData['SGData']
1353    cell = generalData['Cell'][1:8]       
1354    A = G2lat.cell2A(cell[:6])
1355    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1356    adjHKLmax(SGData,Hmax)
1357    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1358    time0 = time.time()
1359    for ref in reflData:
1360        if ref[4] >= dmin:
1361            Fosq,Fcsq,ph = ref[8:11]
1362            for i,hkl in enumerate(ref[11]):
1363                hkl = np.asarray(hkl,dtype='i')
1364                dp = 360.*ref[12][i]
1365                a = cosd(ph+dp)
1366                b = sind(ph+dp)
1367                phasep = complex(a,b)
1368                phasem = complex(a,-b)
1369                F = np.sqrt(Fosq)
1370                h,k,l = hkl+Hmax
1371                Fhkl[h,k,l] = F*phasep
1372                h,k,l = -hkl+Hmax
1373                Fhkl[h,k,l] = F*phasem
1374    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1375    print 'NB: this is just an Fobs map for now - under development'
1376    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1377    print rho.shape
1378    mapData['rho'] = np.real(rho)
1379    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1380    return mapData
1381
1382def FourierMap(data,reflData):
1383    '''default doc string
1384   
1385    :param type name: description
1386   
1387    :returns: type name: description
1388   
1389    '''
1390    generalData = data['General']
1391    if not generalData['Map']['MapType']:
1392        print '**** ERROR - Fourier map not defined'
1393        return
1394    mapData = generalData['Map']
1395    dmin = mapData['Resolution']
1396    SGData = generalData['SGData']
1397    cell = generalData['Cell'][1:8]       
1398    A = G2lat.cell2A(cell[:6])
1399    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1400    adjHKLmax(SGData,Hmax)
1401    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1402#    Fhkl[0,0,0] = generalData['F000X']
1403    time0 = time.time()
1404    for ref in reflData:
1405        if ref[4] >= dmin:
1406            Fosq,Fcsq,ph = ref[8:11]
1407            for i,hkl in enumerate(ref[11]):
1408                hkl = np.asarray(hkl,dtype='i')
1409                dp = 360.*ref[12][i]
1410                a = cosd(ph+dp)
1411                b = sind(ph+dp)
1412                phasep = complex(a,b)
1413                phasem = complex(a,-b)
1414                if 'Fobs' in mapData['MapType']:
1415                    F = np.sqrt(Fosq)
1416                    h,k,l = hkl+Hmax
1417                    Fhkl[h,k,l] = F*phasep
1418                    h,k,l = -hkl+Hmax
1419                    Fhkl[h,k,l] = F*phasem
1420                elif 'Fcalc' in mapData['MapType']:
1421                    F = np.sqrt(Fcsq)
1422                    h,k,l = hkl+Hmax
1423                    Fhkl[h,k,l] = F*phasep
1424                    h,k,l = -hkl+Hmax
1425                    Fhkl[h,k,l] = F*phasem
1426                elif 'delt-F' in mapData['MapType']:
1427                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
1428                    h,k,l = hkl+Hmax
1429                    Fhkl[h,k,l] = dF*phasep
1430                    h,k,l = -hkl+Hmax
1431                    Fhkl[h,k,l] = dF*phasem
1432                elif '2*Fo-Fc' in mapData['MapType']:
1433                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1434                    h,k,l = hkl+Hmax
1435                    Fhkl[h,k,l] = F*phasep
1436                    h,k,l = -hkl+Hmax
1437                    Fhkl[h,k,l] = F*phasem
1438                elif 'Patterson' in mapData['MapType']:
1439                    h,k,l = hkl+Hmax
1440                    Fhkl[h,k,l] = complex(Fosq,0.)
1441                    h,k,l = -hkl+Hmax
1442                    Fhkl[h,k,l] = complex(Fosq,0.)
1443    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1444    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1445    mapData['rho'] = np.real(rho)
1446    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1447    return mapData
1448   
1449# map printing for testing purposes
1450def printRho(SGLaue,rho,rhoMax):                         
1451    '''default doc string
1452   
1453    :param type name: description
1454   
1455    :returns: type name: description
1456   
1457    '''
1458    dim = len(rho.shape)
1459    if dim == 2:
1460        ix,jy = rho.shape
1461        for j in range(jy):
1462            line = ''
1463            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1464                line += (jy-j)*'  '
1465            for i in range(ix):
1466                r = int(100*rho[i,j]/rhoMax)
1467                line += '%4d'%(r)
1468            print line+'\n'
1469    else:
1470        ix,jy,kz = rho.shape
1471        for k in range(kz):
1472            print 'k = ',k
1473            for j in range(jy):
1474                line = ''
1475                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1476                    line += (jy-j)*'  '
1477                for i in range(ix):
1478                    r = int(100*rho[i,j,k]/rhoMax)
1479                    line += '%4d'%(r)
1480                print line+'\n'
1481## keep this
1482               
1483def findOffset(SGData,A,Fhkl):   
1484    '''default doc string
1485   
1486    :param type name: description
1487   
1488    :returns: type name: description
1489   
1490    '''
1491    if SGData['SpGrp'] == 'P 1':
1492        return [0,0,0]   
1493    hklShape = Fhkl.shape
1494    hklHalf = np.array(hklShape)/2
1495    sortHKL = np.argsort(Fhkl.flatten())
1496    Fdict = {}
1497    for hkl in sortHKL:
1498        HKL = np.unravel_index(hkl,hklShape)
1499        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1500        if F == 0.:
1501            break
1502        Fdict['%.6f'%(np.absolute(F))] = hkl
1503    Flist = np.flipud(np.sort(Fdict.keys()))
1504    F = str(1.e6)
1505    i = 0
1506    DH = []
1507    Dphi = []
1508    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1509    for F in Flist:
1510        hkl = np.unravel_index(Fdict[F],hklShape)
1511        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1512        Uniq = np.array(Uniq,dtype='i')
1513        Phi = np.array(Phi)
1514        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1515        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1516        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1517        ang0 = np.angle(Fh0,deg=True)/360.
1518        for H,phi in zip(Uniq,Phi)[1:]:
1519            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1520            dH = H-hkl
1521            dang = ang-ang0
1522            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
1523                continue
1524            DH.append(dH)
1525            Dphi.append((dang+.5) % 1.0)
1526        if i > 20 or len(DH) > 30:
1527            break
1528        i += 1
1529    DH = np.array(DH)
1530    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1531    Dphi = np.array(Dphi)
1532    steps = np.array(hklShape)
1533    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1534    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1535    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1536    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1537    hist,bins = np.histogram(Mmap,bins=1000)
1538#    for i,item in enumerate(hist[:10]):
1539#        print item,bins[i]
1540    chisq = np.min(Mmap)
1541    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1542    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1543#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1544    return DX
1545   
1546def ChargeFlip(data,reflData,pgbar):
1547    '''default doc string
1548   
1549    :param type name: description
1550   
1551    :returns: type name: description
1552   
1553    '''
1554    generalData = data['General']
1555    mapData = generalData['Map']
1556    flipData = generalData['Flip']
1557    FFtable = {}
1558    if 'None' not in flipData['Norm element']:
1559        normElem = flipData['Norm element'].upper()
1560        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1561        for ff in FFs:
1562            if ff['Symbol'] == normElem:
1563                FFtable.update(ff)
1564    dmin = flipData['Resolution']
1565    SGData = generalData['SGData']
1566    cell = generalData['Cell'][1:8]       
1567    A = G2lat.cell2A(cell[:6])
1568    Vol = cell[6]
1569    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1570    adjHKLmax(SGData,Hmax)
1571    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1572    time0 = time.time()
1573    for ref in reflData:
1574        dsp = ref[4]
1575        if dsp >= dmin:
1576            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1577            if FFtable:
1578                SQ = 0.25/dsp**2
1579                ff *= G2el.ScatFac(FFtable,SQ)[0]
1580            if ref[8] > 0.:         #use only +ve Fobs**2
1581                E = np.sqrt(ref[8])/ff
1582            else:
1583                E = 0.
1584            ph = ref[10]
1585            ph = rn.uniform(0.,360.)
1586            for i,hkl in enumerate(ref[11]):
1587                hkl = np.asarray(hkl,dtype='i')
1588                dp = 360.*ref[12][i]
1589                a = cosd(ph+dp)
1590                b = sind(ph+dp)
1591                phasep = complex(a,b)
1592                phasem = complex(a,-b)
1593                h,k,l = hkl+Hmax
1594                Ehkl[h,k,l] = E*phasep
1595                h,k,l = -hkl+Hmax       #Friedel pair refl.
1596                Ehkl[h,k,l] = E*phasem
1597#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1598    CEhkl = copy.copy(Ehkl)
1599    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1600    Emask = ma.getmask(MEhkl)
1601    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1602    Ncyc = 0
1603    old = np.seterr(all='raise')
1604    while True:       
1605        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1606        CEsig = np.std(CErho)
1607        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1608        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
1609        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1610        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1611        phase = CFhkl/np.absolute(CFhkl)
1612        CEhkl = np.absolute(Ehkl)*phase
1613        Ncyc += 1
1614        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1615        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1616        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1617        if Rcf < 5.:
1618            break
1619        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1620        if not GoOn or Ncyc > 10000:
1621            break
1622    np.seterr(**old)
1623    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1624    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1625    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1626    roll = findOffset(SGData,A,CEhkl)
1627       
1628    mapData['Rcf'] = Rcf
1629    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1630    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1631    mapData['rollMap'] = [0,0,0]
1632    return mapData
1633   
1634def SearchMap(data):
1635    '''Does a search of a density map for peaks meeting the criterion of peak
1636    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
1637    mapData is data['General']['mapData']; the map is also in mapData.
1638
1639    :param data: the phase data structure
1640
1641    :returns: (peaks,mags,dzeros) where
1642
1643        * peaks : ndarray
1644          x,y,z positions of the peaks found in the map
1645        * mags : ndarray
1646          the magnitudes of the peaks
1647        * dzeros : ndarray
1648          the distance of the peaks from  the unit cell origin
1649
1650    '''       
1651    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1652   
1653    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
1654   
1655    def noDuplicate(xyz,peaks,Amat):
1656        XYZ = np.inner(Amat,xyz)
1657        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1658            print ' Peak',xyz,' <0.5A from another peak'
1659            return False
1660        return True
1661                           
1662    def fixSpecialPos(xyz,SGData,Amat):
1663        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1664        X = []
1665        xyzs = [equiv[0] for equiv in equivs]
1666        for x in xyzs:
1667            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1668                X.append(x)
1669        if len(X) > 1:
1670            return np.average(X,axis=0)
1671        else:
1672            return xyz
1673       
1674    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1675        Mag,x0,y0,z0,sig = parms
1676        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
1677#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
1678        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
1679       
1680    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1681        Mag,x0,y0,z0,sig = parms
1682        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1683        return M
1684       
1685    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1686        Mag,x0,y0,z0,sig = parms
1687        dMdv = np.zeros(([5,]+list(rX.shape)))
1688        delt = .01
1689        for i in range(5):
1690            parms[i] -= delt
1691            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1692            parms[i] += 2.*delt
1693            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1694            parms[i] -= delt
1695            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1696        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1697        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1698        dMdv = np.reshape(dMdv,(5,rX.size))
1699        Hess = np.inner(dMdv,dMdv)
1700       
1701        return Vec,Hess
1702       
1703    generalData = data['General']
1704    phaseName = generalData['Name']
1705    SGData = generalData['SGData']
1706    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1707    drawingData = data['Drawing']
1708    peaks = []
1709    mags = []
1710    dzeros = []
1711    try:
1712        mapData = generalData['Map']
1713        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1714        rho = copy.copy(mapData['rho'])     #don't mess up original
1715        mapHalf = np.array(rho.shape)/2
1716        res = mapData['Resolution']
1717        incre = np.array(rho.shape,dtype=np.float)
1718        step = max(1.0,1./res)+1
1719        steps = np.array(3*[step,])
1720    except KeyError:
1721        print '**** ERROR - Fourier map not defined'
1722        return peaks,mags
1723    rhoMask = ma.array(rho,mask=(rho<contLevel))
1724    indices = (-1,0,1)
1725    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1726    for roll in rolls:
1727        if np.any(roll):
1728            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1729    indx = np.transpose(rhoMask.nonzero())
1730    peaks = indx/incre
1731    mags = rhoMask[rhoMask.nonzero()]
1732    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1733        rho = rollMap(rho,ind)
1734        rMM = mapHalf-steps
1735        rMP = mapHalf+steps+1
1736        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1737        peakInt = np.sum(rhoPeak)*res**3
1738        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1739        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
1740        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1741            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1742        x1 = result[0]
1743        if not np.any(x1 < 0):
1744            mag = x1[0]
1745            peak = (np.array(x1[1:4])-ind)/incre
1746        peak = fixSpecialPos(peak,SGData,Amat)
1747        rho = rollMap(rho,-ind)       
1748    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1749    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1750   
1751def sortArray(data,pos,reverse=False):
1752    '''data is a list of items
1753    sort by pos in list; reverse if True
1754    '''
1755    T = []
1756    for i,M in enumerate(data):
1757        T.append((M[pos],i))
1758    D = dict(zip(T,data))
1759    T.sort()
1760    if reverse:
1761        T.reverse()
1762    X = []
1763    for key in T:
1764        X.append(D[key])
1765    return X
1766
1767def PeaksEquiv(data,Ind):
1768    '''Find the equivalent map peaks for those selected. Works on the
1769    contents of data['Map Peaks'].
1770
1771    :param data: the phase data structure
1772    :param list Ind: list of selected peak indices
1773    :returns: augmented list of peaks including those related by symmetry to the
1774      ones in Ind
1775
1776    '''       
1777    def Duplicate(xyz,peaks,Amat):
1778        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1779            return True
1780        return False
1781                           
1782    generalData = data['General']
1783    cell = generalData['Cell'][1:7]
1784    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1785    A = G2lat.cell2A(cell)
1786    SGData = generalData['SGData']
1787    mapPeaks = data['Map Peaks']
1788    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1789    Indx = {}
1790    for ind in Ind:
1791        xyz = np.array(mapPeaks[ind][1:4])
1792        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1793#        for x in xyzs: print x
1794        for jnd,xyz in enumerate(XYZ):       
1795            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1796    Ind = []
1797    for ind in Indx:
1798        if Indx[ind]:
1799            Ind.append(ind)
1800    return Ind
1801               
1802def PeaksUnique(data,Ind):
1803    '''Finds the symmetry unique set of peaks from those selected. Works on the
1804    contents of data['Map Peaks'].
1805
1806    :param data: the phase data structure
1807    :param list Ind: list of selected peak indices
1808    :returns: the list of symmetry unique peaks from among those given in Ind
1809
1810    '''       
1811#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1812
1813    def noDuplicate(xyz,peaks,Amat):
1814        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1815            return False
1816        return True
1817                           
1818    generalData = data['General']
1819    cell = generalData['Cell'][1:7]
1820    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1821    A = G2lat.cell2A(cell)
1822    SGData = generalData['SGData']
1823    mapPeaks = data['Map Peaks']
1824    Indx = {}
1825    XYZ = {}
1826    for ind in Ind:
1827        XYZ[ind] = np.array(mapPeaks[ind][1:4])
1828        Indx[ind] = True
1829    for ind in Ind:
1830        if Indx[ind]:
1831            xyz = XYZ[ind]
1832            for jnd in Ind:
1833                if ind != jnd and Indx[jnd]:                       
1834                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1835                    xyzs = np.array([equiv[0] for equiv in Equiv])
1836                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1837    Ind = []
1838    for ind in Indx:
1839        if Indx[ind]:
1840            Ind.append(ind)
1841    return Ind
1842   
1843################################################################################
1844##### single peak fitting profile fxn stuff
1845################################################################################
1846
1847def getCWsig(ins,pos):
1848    '''default doc string
1849   
1850    :param type name: description
1851   
1852    :returns: type name: description
1853   
1854    '''
1855    tp = tand(pos/2.0)
1856    return ins['U']*tp**2+ins['V']*tp+ins['W']
1857   
1858def getCWsigDeriv(pos):
1859    '''default doc string
1860   
1861    :param type name: description
1862   
1863    :returns: type name: description
1864   
1865    '''
1866    tp = tand(pos/2.0)
1867    return tp**2,tp,1.0
1868   
1869def getCWgam(ins,pos):
1870    '''default doc string
1871   
1872    :param type name: description
1873   
1874    :returns: type name: description
1875   
1876    '''
1877    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
1878   
1879def getCWgamDeriv(pos):
1880    '''default doc string
1881   
1882    :param type name: description
1883   
1884    :returns: type name: description
1885   
1886    '''
1887    return 1./cosd(pos/2.0),tand(pos/2.0)
1888   
1889def getTOFsig(ins,dsp):
1890    '''default doc string
1891   
1892    :param type name: description
1893   
1894    :returns: type name: description
1895   
1896    '''
1897    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-q']*dsp
1898   
1899def getTOFsigDeriv(dsp):
1900    '''default doc string
1901   
1902    :param type name: description
1903   
1904    :returns: type name: description
1905   
1906    '''
1907    return 1.0,dsp**2,dsp
1908   
1909def getTOFgamma(ins,dsp):
1910    '''default doc string
1911   
1912    :param type name: description
1913   
1914    :returns: type name: description
1915   
1916    '''
1917    return ins['X']*dsp+ins['Y']*dsp**2
1918   
1919def getTOFgammaDeriv(dsp):
1920    '''default doc string
1921   
1922    :param type name: description
1923   
1924    :returns: type name: description
1925   
1926    '''
1927    return dsp,dsp**2
1928   
1929def getTOFbeta(ins,dsp):
1930    '''default doc string
1931   
1932    :param type name: description
1933   
1934    :returns: type name: description
1935   
1936    '''
1937    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp
1938   
1939def getTOFbetaDeriv(dsp):
1940    '''default doc string
1941   
1942    :param type name: description
1943   
1944    :returns: type name: description
1945   
1946    '''
1947    return 1.0,1./dsp**4,1./dsp
1948   
1949def getTOFalpha(ins,dsp):
1950    '''default doc string
1951   
1952    :param type name: description
1953   
1954    :returns: type name: description
1955   
1956    '''
1957    return ins['alpha']/dsp
1958   
1959def getTOFalphaDeriv(dsp):
1960    '''default doc string
1961   
1962    :param type name: description
1963   
1964    :returns: type name: description
1965   
1966    '''
1967    return 1./dsp
1968   
1969def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1970    '''default doc string
1971   
1972    :param type name: description
1973   
1974    :returns: type name: description
1975   
1976    '''
1977    ind = 0
1978    if useFit:
1979        ind = 1
1980    ins = {}
1981    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
1982        for x in ['U','V','W','X','Y']:
1983            ins[x] = Parms[x][ind]
1984        if ifQ:                              #qplot - convert back to 2-theta
1985            pos = 2.0*asind(pos*wave/(4*math.pi))
1986        sig = getCWsig(ins,pos)
1987        gam = getCWgam(ins,pos)           
1988        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
1989    else:
1990        if ifQ:
1991            dsp = 2.*np.pi/pos
1992            pos = Parms['difC']*dsp
1993        else:
1994            dsp = pos/Parms['difC'][1]
1995        if 'Pdabc' in Parms2:
1996            for x in ['sig-0','sig-1','sig-q','X','Y']:
1997                ins[x] = Parms[x][ind]
1998            Pdabc = Parms2['Pdabc'].T
1999            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2000            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2001        else:
2002            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y']:
2003                ins[x] = Parms[x][ind]
2004            alp = getTOFalpha(ins,dsp)
2005            bet = getTOFbeta(ins,dsp)
2006        sig = getTOFsig(ins,dsp)
2007        gam = getTOFgamma(ins,dsp)
2008        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2009    return XY
2010   
2011################################################################################
2012##### MC/SA stuff
2013################################################################################
2014
2015#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2016# Original Author: Travis Oliphant 2002
2017# Bug-fixes in 2006 by Tim Leslie
2018
2019
2020import numpy
2021from numpy import asarray, tan, exp, ones, squeeze, sign, \
2022     all, log, sqrt, pi, shape, array, minimum, where
2023from numpy import random
2024
2025#__all__ = ['anneal']
2026
2027_double_min = numpy.finfo(float).min
2028_double_max = numpy.finfo(float).max
2029class base_schedule(object):
2030    def __init__(self):
2031        self.dwell = 20
2032        self.learn_rate = 0.5
2033        self.lower = -10
2034        self.upper = 10
2035        self.Ninit = 50
2036        self.accepted = 0
2037        self.tests = 0
2038        self.feval = 0
2039        self.k = 0
2040        self.T = None
2041
2042    def init(self, **options):
2043        self.__dict__.update(options)
2044        self.lower = asarray(self.lower)
2045        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2046        self.upper = asarray(self.upper)
2047        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2048        self.k = 0
2049        self.accepted = 0
2050        self.feval = 0
2051        self.tests = 0
2052
2053    def getstart_temp(self, best_state):
2054        """ Find a matching starting temperature and starting parameters vector
2055        i.e. find x0 such that func(x0) = T0.
2056
2057        Parameters
2058        ----------
2059        best_state : _state
2060            A _state object to store the function value and x0 found.
2061
2062        returns
2063        -------
2064        x0 : array
2065            The starting parameters vector.
2066        """
2067
2068        assert(not self.dims is None)
2069        lrange = self.lower
2070        urange = self.upper
2071        fmax = _double_min
2072        fmin = _double_max
2073        for _ in range(self.Ninit):
2074            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2075            fval = self.func(x0, *self.args)
2076            self.feval += 1
2077            if fval > fmax:
2078                fmax = fval
2079            if fval < fmin:
2080                fmin = fval
2081                best_state.cost = fval
2082                best_state.x = array(x0)
2083
2084        self.T0 = (fmax-fmin)*1.5
2085        return best_state.x
2086
2087    def accept_test(self, dE):
2088        T = self.T
2089        self.tests += 1
2090        if dE < 0:
2091            self.accepted += 1
2092            return 1
2093        p = exp(-dE*1.0/self.boltzmann/T)
2094        if (p > random.uniform(0.0, 1.0)):
2095            self.accepted += 1
2096            return 1
2097        return 0
2098
2099    def update_guess(self, x0):
2100        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2101
2102    def update_temp(self, x0):
2103        pass
2104
2105
2106#  A schedule due to Lester Ingber modified to use bounds - OK
2107class fast_sa(base_schedule):
2108    def init(self, **options):
2109        self.__dict__.update(options)
2110        if self.m is None:
2111            self.m = 1.0
2112        if self.n is None:
2113            self.n = 1.0
2114        self.c = self.m * exp(-self.n * self.quench)
2115
2116    def update_guess(self, x0):
2117        x0 = asarray(x0)
2118        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2119        T = self.T
2120        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2121        xnew = xc*(self.upper - self.lower)+self.lower
2122        return xnew
2123#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2124#        xc = y*(self.upper - self.lower)
2125#        xnew = x0 + xc
2126#        return xnew
2127
2128    def update_temp(self):
2129        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2130        self.k += 1
2131        return
2132
2133class cauchy_sa(base_schedule):     #modified to use bounds - not good
2134    def update_guess(self, x0):
2135        x0 = asarray(x0)
2136        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2137        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2138        xnew = xc*(self.upper - self.lower)+self.lower
2139        return xnew
2140#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2141#        xc = self.learn_rate * self.T * tan(numbers)
2142#        xnew = x0 + xc
2143#        return xnew
2144
2145    def update_temp(self):
2146        self.T = self.T0/(1+self.k)
2147        self.k += 1
2148        return
2149
2150class boltzmann_sa(base_schedule):
2151#    def update_guess(self, x0):
2152#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2153#        x0 = asarray(x0)
2154#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2155#
2156#        xnew = x0 + xc*std*self.learn_rate
2157#        return xnew
2158
2159    def update_temp(self):
2160        self.k += 1
2161        self.T = self.T0 / log(self.k+1.0)
2162        return
2163
2164class log_sa(base_schedule):        #OK
2165
2166    def init(self,**options):
2167        self.__dict__.update(options)
2168       
2169#    def update_guess(self,x0):
2170#        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2171       
2172    def update_temp(self):
2173        self.k += 1
2174        self.T = self.T0*self.slope**self.k
2175       
2176class Tremayne_sa(base_schedule):   #needs fixing for two steps
2177
2178    def init(self,**options):
2179        self.__dict__.update(options)
2180
2181#    def update_guess(self,x0):
2182#        x0 = np.asarray(x0)
2183#        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
2184#        xnew = x0+u
2185#        return xnew       
2186   
2187    def update_temp(self):
2188        self.k += 1
2189        self.T = self.T0*self.slope**k
2190   
2191class _state(object):
2192    def __init__(self):
2193        self.x = None
2194        self.cost = None
2195
2196# TODO:
2197#     allow for general annealing temperature profile
2198#     in that case use update given by alpha and omega and
2199#     variation of all previous updates and temperature?
2200
2201# Simulated annealing
2202
2203def anneal(func, x0, args=(), schedule='fast', full_output=0,
2204           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2205           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2206           lower=-100, upper=100, dwell=50, slope=0.9,dlg=None):
2207    """Minimize a function using simulated annealing.
2208
2209    Schedule is a schedule class implementing the annealing schedule.
2210    Available ones are 'fast', 'cauchy', 'boltzmann'
2211
2212    :param callable func: f(x, \*args)
2213        Function to be optimized.
2214    :param ndarray x0:
2215        Initial guess.
2216    :param tuple args:
2217        Extra parameters to `func`.
2218    :param base_schedule schedule:
2219        Annealing schedule to use (a class).
2220    :param bool full_output:
2221        Whether to return optional outputs.
2222    :param float T0:
2223        Initial Temperature (estimated as 1.2 times the largest
2224        cost-function deviation over random points in the range).
2225    :param float Tf:
2226        Final goal temperature.
2227    :param int maxeval:
2228        Maximum function evaluations.
2229    :param int maxaccept:
2230        Maximum changes to accept.
2231    :param int maxiter:
2232        Maximum cooling iterations.
2233    :param float learn_rate:
2234        Scale constant for adjusting guesses.
2235    :param float boltzmann:
2236        Boltzmann constant in acceptance test
2237        (increase for less stringent test at each temperature).
2238    :param float feps:
2239        Stopping relative error tolerance for the function value in
2240        last four coolings.
2241    :param float quench,m,n:
2242        Parameters to alter fast_sa schedule.
2243    :param float/ndarray lower,upper:
2244        Lower and upper bounds on `x`.
2245    :param int dwell:
2246        The number of times to search the space at each temperature.
2247    :param float slope:
2248        Parameter for log schedule
2249
2250    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2251
2252     * xmin (ndarray): Point giving smallest value found.
2253     * Jmin (float): Minimum value of function found.
2254     * T (float): Final temperature.
2255     * feval (int): Number of function evaluations.
2256     * iters (int): Number of cooling iterations.
2257     * accept (int): Number of tests accepted.
2258     * retval (int): Flag indicating stopping condition:
2259
2260              *  0: Points no longer changing
2261              *  1: Cooled to final temperature
2262              *  2: Maximum function evaluations
2263              *  3: Maximum cooling iterations reached
2264              *  4: Maximum accepted query locations reached
2265              *  5: Final point not the minimum amongst encountered points
2266
2267    *Notes*:
2268    Simulated annealing is a random algorithm which uses no derivative
2269    information from the function being optimized. In practice it has
2270    been more useful in discrete optimization than continuous
2271    optimization, as there are usually better algorithms for continuous
2272    optimization problems.
2273
2274    Some experimentation by trying the difference temperature
2275    schedules and altering their parameters is likely required to
2276    obtain good performance.
2277
2278    The randomness in the algorithm comes from random sampling in numpy.
2279    To obtain the same results you can call numpy.random.seed with the
2280    same seed immediately before calling scipy.optimize.anneal.
2281
2282    We give a brief description of how the three temperature schedules
2283    generate new points and vary their temperature. Temperatures are
2284    only updated with iterations in the outer loop. The inner loop is
2285    over xrange(dwell), and new points are generated for
2286    every iteration in the inner loop. (Though whether the proposed
2287    new points are accepted is probabilistic.)
2288
2289    For readability, let d denote the dimension of the inputs to func.
2290    Also, let x_old denote the previous state, and k denote the
2291    iteration number of the outer loop. All other variables not
2292    defined below are input variables to scipy.optimize.anneal itself.
2293
2294    In the 'fast' schedule the updates are ::
2295
2296        u ~ Uniform(0, 1, size=d)
2297        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2298        xc = y * (upper - lower)
2299        x_new = x_old + xc
2300
2301        c = n * exp(-n * quench)
2302        T_new = T0 * exp(-c * k**quench)
2303
2304
2305    In the 'cauchy' schedule the updates are ::
2306
2307        u ~ Uniform(-pi/2, pi/2, size=d)
2308        xc = learn_rate * T * tan(u)
2309        x_new = x_old + xc
2310
2311        T_new = T0 / (1+k)
2312
2313    In the 'boltzmann' schedule the updates are ::
2314
2315        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2316        y ~ Normal(0, std, size=d)
2317        x_new = x_old + learn_rate * y
2318
2319        T_new = T0 / log(1+k)
2320
2321    """
2322    x0 = asarray(x0)
2323    lower = asarray(lower)
2324    upper = asarray(upper)
2325
2326    schedule = eval(schedule+'_sa()')
2327    #   initialize the schedule
2328    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2329                  learn_rate=learn_rate, lower=lower, upper=upper,
2330                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2331
2332    current_state, last_state, best_state = _state(), _state(), _state()
2333    if T0 is None:
2334        x0 = schedule.getstart_temp(best_state)
2335    else:
2336        x0 = random.uniform(size=len(x0))*(upper-lower) + lower #comment to avoid random start
2337        best_state.x = None
2338        best_state.cost = numpy.Inf
2339
2340    last_state.x = asarray(x0).copy()
2341    fval = func(x0,*args)
2342    schedule.feval += 1
2343    last_state.cost = fval
2344    if last_state.cost < best_state.cost:
2345        best_state.cost = fval
2346        best_state.x = asarray(x0).copy()
2347    schedule.T = schedule.T0
2348    fqueue = [100, 300, 500, 700]
2349    iters = 0
2350    keepGoing = True
2351    while keepGoing:
2352        retval = 0
2353        for n in xrange(dwell):
2354            current_state.x = schedule.update_guess(last_state.x)
2355            current_state.cost = func(current_state.x,*args)
2356            schedule.feval += 1
2357
2358            dE = current_state.cost - last_state.cost
2359            if schedule.accept_test(dE):
2360                last_state.x = current_state.x.copy()
2361                last_state.cost = current_state.cost
2362                if last_state.cost < best_state.cost:
2363                    best_state.x = last_state.x.copy()
2364                    best_state.cost = last_state.cost
2365        if dlg:
2366            GoOn = dlg.Update(best_state.cost*100,
2367                newmsg='%s%8.5f\n%s%8.4f%s'%('Temperature =',schedule.T,'MC/SA Residual =',best_state.cost*100,'%'))[0]
2368            if not GoOn:
2369                best_state.x = last_state.x.copy()
2370                best_state.cost = last_state.cost
2371                retval = 5
2372        schedule.update_temp()
2373        iters += 1
2374        # Stopping conditions
2375        # 0) last saved values of f from each cooling step
2376        #     are all very similar (effectively cooled)
2377        # 1) Tf is set and we are below it
2378        # 2) maxeval is set and we are past it
2379        # 3) maxiter is set and we are past it
2380        # 4) maxaccept is set and we are past it
2381        # 5) user canceled run via progress bar
2382
2383        fqueue.append(squeeze(last_state.cost))
2384        fqueue.pop(0)
2385        af = asarray(fqueue)*1.0
2386        if retval == 5:
2387            print ' User terminated run; incomplete MC/SA'
2388            keepGoing = False
2389            break
2390        if all(abs((af-af[0])/af[0]) < feps):
2391            retval = 0
2392            if abs(af[-1]-best_state.cost) > feps*10:
2393                retval = 5
2394#                print "Warning: Cooled to %f at %s but this is not" \
2395#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2396#                      + " the smallest point found."
2397            break
2398        if (Tf is not None) and (schedule.T < Tf):
2399            retval = 1
2400            break
2401        if (maxeval is not None) and (schedule.feval > maxeval):
2402            retval = 2
2403            break
2404        if (iters > maxiter):
2405            print "Warning: Maximum number of iterations exceeded."
2406            retval = 3
2407            break
2408        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2409            retval = 4
2410            break
2411
2412    if full_output:
2413        return best_state.x, best_state.cost, schedule.T, \
2414               schedule.feval, iters, schedule.accepted, retval
2415    else:
2416        return best_state.x, retval
2417
2418def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q):
2419    outlist = []
2420    for n in range(iCyc):
2421        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2422        outlist.append(result[0])
2423        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2424    out_q.put(outlist)
2425
2426def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2427    import multiprocessing as mp
2428   
2429    nprocs = mp.cpu_count()
2430    out_q = mp.Queue()
2431    procs = []
2432    iCyc = np.zeros(nprocs)
2433    for i in range(nCyc):
2434        iCyc[i%nprocs] += 1
2435    for i in range(nprocs):
2436        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q))
2437        procs.append(p)
2438        p.start()
2439    resultlist = []
2440    for i in range(nprocs):
2441        resultlist += out_q.get()
2442    for p in procs:
2443        p.join()
2444    return resultlist
2445
2446def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2447    '''default doc string
2448   
2449    :param type name: description
2450   
2451    :returns: type name: description
2452   
2453    '''
2454    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
2455   
2456    twopi = 2.0*np.pi
2457    global tsum
2458    tsum = 0.
2459   
2460    def getMDparms(item,pfx,parmDict,varyList):
2461        parmDict[pfx+'MDaxis'] = item['axis']
2462        parmDict[pfx+'MDval'] = item['Coef'][0]
2463        if item['Coef'][1]:
2464            varyList += [pfx+'MDval',]
2465            limits = item['Coef'][2]
2466            lower.append(limits[0])
2467            upper.append(limits[1])
2468                       
2469    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2470        parmDict[pfx+'Atype'] = item['atType']
2471        aTypes |= set([item['atType'],]) 
2472        pstr = ['Ax','Ay','Az']
2473        XYZ = [0,0,0]
2474        for i in range(3):
2475            name = pfx+pstr[i]
2476            parmDict[name] = item['Pos'][0][i]
2477            XYZ[i] = parmDict[name]
2478            if item['Pos'][1][i]:
2479                varyList += [name,]
2480                limits = item['Pos'][2][i]
2481                lower.append(limits[0])
2482                upper.append(limits[1])
2483        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2484           
2485    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2486        parmDict[mfx+'MolCent'] = item['MolCent']
2487        parmDict[mfx+'RBId'] = item['RBId']
2488        pstr = ['Px','Py','Pz']
2489        ostr = ['Qa','Qi','Qj','Qk']
2490        for i in range(3):
2491            name = mfx+pstr[i]
2492            parmDict[name] = item['Pos'][0][i]
2493            if item['Pos'][1][i]:
2494                varyList += [name,]
2495                limits = item['Pos'][2][i]
2496                lower.append(limits[0])
2497                upper.append(limits[1])
2498        for i in range(4):
2499            name = mfx+ostr[i]
2500            parmDict[name] = item['Ori'][0][i]
2501            if item['Ovar'] == 'AV' and i:
2502                varyList += [name,]
2503                limits = item['Ori'][2][i]
2504                lower.append(limits[0])
2505                upper.append(limits[1])
2506            elif item['Ovar'] == 'A' and not i:
2507                varyList += [name,]
2508                limits = item['Ori'][2][i]
2509                lower.append(limits[0])
2510                upper.append(limits[1])
2511        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2512            for i in range(len(item['Tor'][0])):
2513                name = mfx+'Tor'+str(i)
2514                parmDict[name] = item['Tor'][0][i]
2515                if item['Tor'][1][i]:
2516                    varyList += [name,]
2517                    limits = item['Tor'][2][i]
2518                    lower.append(limits[0])
2519                    upper.append(limits[1])
2520        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2521        aTypes |= set(atypes)
2522        atNo += len(atypes)
2523        return atNo
2524       
2525    def GetAtomM(Xdata,SGData):
2526        Mdata = []
2527        for xyz in Xdata:
2528            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2529        return np.array(Mdata)
2530       
2531    def GetAtomTX(RBdata,parmDict):
2532        'Needs a doc string'
2533        Bmat = parmDict['Bmat']
2534        atNo = parmDict['atNo']
2535        nfixAt = parmDict['nfixAt']
2536        Tdata = atNo*[' ',]
2537        Xdata = np.zeros((3,atNo))
2538        keys = {':Atype':Tdata,':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2539        for iatm in range(nfixAt):
2540            for key in keys:
2541                parm = ':'+str(iatm)+key
2542                if parm in parmDict:
2543                    if key == ':Atype':
2544                        keys[key][iatm] = aTypes.index(parmDict[parm])
2545                    else:
2546                        keys[key][iatm] = parmDict[parm]
2547        iatm = nfixAt
2548        for iObj in range(parmDict['nObj']):
2549            pfx = str(iObj)+':'
2550            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2551                if parmDict[pfx+'Type'] == 'Vector':
2552                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2553                    vecs = RBRes['rbVect']
2554                    mags = RBRes['VectMag']
2555                    Cart = np.zeros_like(vecs[0])
2556                    for vec,mag in zip(vecs,mags):
2557                        Cart += vec*mag
2558                elif parmDict[pfx+'Type'] == 'Residue':
2559                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2560                    Cart = np.array(RBRes['rbXYZ'])
2561                    for itor,seq in enumerate(RBRes['rbSeq']):
2562                        tName = pfx+'Tor'+str(itor)
2563                        QuatA = AVdeg2Q(parmDict[tName],Cart[seq[0]]-Cart[seq[1]])
2564                        for ride in seq[3]:
2565                            Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
2566                if parmDict[pfx+'MolCent'][1]:
2567                    Cart -= parmDict[pfx+'MolCent'][0]
2568                Qori = normQ(np.array([parmDict[pfx+'Qa'],parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']]))
2569                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2570                for i,x in enumerate(Cart):
2571                    X = np.inner(Bmat,prodQVQ(Qori,x))+Pos
2572                    for j in range(3):
2573                        Xdata[j][iatm] = X[j]
2574                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
2575                    iatm += 1
2576            elif parmDict[pfx+'Type'] == 'Atom':
2577                atNo = parmDict[pfx+'atNo']
2578                for key in keys:
2579                    parm = pfx+key[1:]              #remove extra ':'
2580                    if parm in parmDict:
2581                        if key == ':Atype':
2582                            keys[key][atNo] = aTypes.index(parmDict[parm])
2583                        else:
2584                            keys[key][atNo] = parmDict[parm]
2585                iatm += 1
2586            else:
2587                continue        #skips March Dollase
2588        return Tdata,Xdata.T
2589   
2590    def calcMDcorr(MDval,MDaxis,Uniq,G):
2591        ''' Calls fortran routine'''
2592        MDcorr = pyd.pymdcalc(MDval,MDaxis,len(Uniq),Uniq.flatten(),G)
2593        return MDcorr
2594       
2595    def mcsaMDSFcalc(ifInv,Tdata,Mdata,Xdata,MDval,MDaxis,G,mul,FFs,Uniq,Phi):
2596        ''' Calls fortran routine'''
2597        Icalc = pyd.pymcsamdsfcalc(ifInv,len(Tdata),Tdata,Mdata,Xdata.flatten(),
2598            MDval,MDaxis,G,mul,len(FFs),FFs,len(Uniq),Uniq.flatten(),Phi)
2599        return Icalc
2600
2601    def mcsaSFcalc(ifInv,Tdata,Mdata,Xdata,mul,FFs,Uniq,Phi):
2602        ''' Calls fortran routine'''
2603        Icalc = pyd.pymcsasfcalc(ifInv,len(Tdata),Tdata,Mdata,Xdata.flatten(),
2604            mul,len(FFs),FFs,len(Uniq),Uniq.flatten(),Phi)
2605        return Icalc
2606
2607    def mcsaCalc(values,refList,rcov,ifInv,RBdata,varyList,parmDict):
2608        ''' Compute structure factors for all h,k,l for phase
2609        input:
2610            refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2611            ParmDict:
2612        puts result F^2 in each ref[8] in refList
2613        '''       
2614        global tsum
2615        parmDict.update(dict(zip(varyList,values)))
2616        Tdata,Xdata = GetAtomTX(RBdata,parmDict)
2617        Mdata = parmDict['Mdata']
2618        MDval = parmDict['0:MDval']
2619        MDaxis = parmDict['0:MDaxis']
2620        Gmat = parmDict['Gmat']
2621        Srefs = np.zeros(len(refList))
2622        sumFcsq = 0.
2623        for refl in refList:
2624            t0 = time.time()
2625            refl[5] = mcsaMDSFcalc(ifInv,Tdata,Mdata,Xdata,MDval,MDaxis,Gmat,
2626                refl[3],refl[7],refl[8],refl[9])
2627#            refl[5] = mcsaSFcalc(ifInv,Tdata,Mdata,Xdata,refl[3],refl[7],refl[8],refl[9])
2628#            refl[5] *= calcMDcorr(MDval,MDaxis,refl[8],Gmat)
2629            tsum += (time.time()-t0)
2630            sumFcsq += refl[5]
2631        scale = (parmDict['sumFosq']/sumFcsq)
2632        for iref,refl in enumerate(refList):
2633            refl[5] *= scale
2634            Srefs[iref] = refl[4]-refl[5]
2635        M = np.inner(Srefs,np.inner(rcov,Srefs))
2636        return M/parmDict['sumFosq']**2
2637   
2638    sq8ln2 = np.sqrt(8*np.log(2))
2639    sq2pi = np.sqrt(2*np.pi)
2640    sq4pi = np.sqrt(4*np.pi)
2641    generalData = data['General']
2642    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2643    Gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])[0]
2644    SGData = generalData['SGData']
2645    fixAtoms = data['Atoms']                       #if any
2646    cx,ct,cs = generalData['AtomPtrs'][:3]
2647    aTypes = set([])
2648    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2649    varyList = []
2650    atNo = 0
2651    for atm in fixAtoms:
2652        pfx = ':'+str(atNo)+':'
2653        parmDict[pfx+'Atype'] = atm[ct]
2654        aTypes |= set([atm[ct],])
2655        pstr = ['Ax','Ay','Az']
2656        parmDict[pfx+'Amul'] = atm[cs+1]
2657        for i in range(3):
2658            name = pfx+pstr[i]
2659            parmDict[name] = atm[cx+i]
2660        atNo += 1
2661    parmDict['nfixAt'] = len(fixAtoms)       
2662    MCSA = generalData['MCSA controls']
2663    reflName = MCSA['Data source']
2664    phaseName = generalData['Name']
2665    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2666    upper = []
2667    lower = []
2668    for i,item in enumerate(MCSAObjs):
2669        mfx = str(i)+':'
2670        parmDict[mfx+'Type'] = item['Type']
2671        if item['Type'] == 'MD':
2672            getMDparms(item,mfx,parmDict,varyList)
2673        elif item['Type'] == 'Atom':
2674            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2675            parmDict[mfx+'atNo'] = atNo
2676            atNo += 1
2677        elif item['Type'] in ['Residue','Vector']:
2678            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2679    parmDict['atNo'] = atNo                 #total no. of atoms
2680    parmDict['nObj'] = len(MCSAObjs)
2681    aTypes = list(aTypes)
2682    Tdata,Xdata = GetAtomTX(RBdata,parmDict)
2683    parmDict['Mdata'] = GetAtomM(Xdata,SGData)
2684    FFtables = G2el.GetFFtable(aTypes)
2685    refs = []
2686    sumFosq = 0
2687    if 'PWDR' in reflName:
2688        for ref in reflData:
2689            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2690            if d >= MCSA['dmin']:
2691                sig = gamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2692                SQ = 0.25/d**2
2693                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2694                FFs = G2el.getFFvalues(FFtables,SQ,True)
2695                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
2696                sumFosq += f*m
2697        nRef = len(refs)
2698        rcov = np.zeros((nRef,nRef))
2699        for iref,refI in enumerate(refs):
2700            rcov[iref][iref] = 1./(sq4pi*refI[6])
2701            for jref,refJ in enumerate(refs[:iref]):
2702                t1 = refI[6]**2+refJ[6]**2
2703                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2704                if t2 > 10.:
2705                    rcov[iref][jref] = 0.
2706                else:
2707                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2708        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2709        Rdiag = np.sqrt(np.diag(rcov))
2710        Rnorm = np.outer(Rdiag,Rdiag)
2711        rcov /= Rnorm
2712    elif 'Pawley' in reflName:
2713        covMatrix = covData['covMatrix']
2714        vList = covData['varyList']
2715        for iref,refI in enumerate(reflData):
2716            h,k,l,m,d,v,f,s = refI
2717            if d >= MCSA['dmin'] and v:       #skip unrefined ones
2718                SQ = 0.25/d**2
2719                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2720                FFs = G2el.getFFvalues(FFtables,SQ,True)
2721                refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
2722                sumFosq += f*m
2723        nRef = len(refs)
2724        pfx = str(data['pId'])+'::PWLref:'
2725        rcov = np.zeros((nRef,nRef))       
2726        for iref,refI in enumerate(refs):
2727            I = refI[5]
2728            nameI = pfx+str(I)
2729            if nameI in vList:
2730                Iindx = vList.index(nameI)
2731                rcov[iref][iref] = covMatrix[Iindx][Iindx]
2732                for jref,refJ in enumerate(refs[:iref]):
2733                    J = refJ[5]
2734                    nameJ = pfx+str(J)
2735                    try:
2736                        rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2737                    except ValueError:
2738                        rcov[iref][jref] = rcov[iref][jref-1]
2739            else:
2740                rcov[iref] = rcov[iref-1]
2741                rcov[iref][iref] = rcov[iref-1][iref-1]
2742        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2743        Rdiag = np.sqrt(np.diag(rcov))
2744        Rnorm = np.outer(Rdiag,Rdiag)
2745        rcov /= Rnorm
2746    elif 'HKLF' in reflName:
2747        for ref in reflData:
2748            [h,k,l,m,d],f = ref[:5],ref[6]
2749            if d >= MCSA['dmin']:
2750                SQ = 0.25/d**2
2751                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2752                FFs = G2el.getFFvalues(FFtables,SQ,True)
2753                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
2754                sumFosq += f*m
2755        nRef = len(refs)
2756        rcov = np.identity(len(refs))
2757    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
2758    print ' Number of parameters varied: %d'%(len(varyList))
2759    parmDict['sumFosq'] = sumFosq
2760    x0 = [parmDict[val] for val in varyList]
2761    ifInv = SGData['SGInv']
2762    results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,RBdata,varyList,parmDict),
2763        schedule=MCSA['Algorithm'], full_output=True,
2764        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2765        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
2766        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2767        lower=lower, upper=upper, slope=MCSA['log slope'],dlg=pgbar)
2768    Result = [False,False,results[1],results[2],]+list(results[0])
2769    Result.append(varyList)
2770    return Result,tsum
2771
2772       
2773################################################################################
2774##### Quaternion stuff
2775################################################################################
2776
2777def prodQQ(QA,QB):
2778    ''' Grassman quaternion product
2779        QA,QB quaternions; q=r+ai+bj+ck
2780    '''
2781    D = np.zeros(4)
2782    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2783    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2784    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2785    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2786   
2787#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
2788#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
2789   
2790    return D
2791   
2792def normQ(QA):
2793    ''' get length of quaternion & normalize it
2794        q=r+ai+bj+ck
2795    '''
2796    n = np.sqrt(np.sum(np.array(QA)**2))
2797    return QA/n
2798   
2799def invQ(Q):
2800    '''
2801        get inverse of quaternion
2802        q=r+ai+bj+ck; q* = r-ai-bj-ck
2803    '''
2804    return Q*np.array([1,-1,-1,-1])
2805   
2806def prodQVQ(Q,V):
2807    """
2808    compute the quaternion vector rotation qvq-1 = v'
2809    q=r+ai+bj+ck
2810    """
2811    VP = np.zeros(3)
2812    T2 = Q[0]*Q[1]
2813    T3 = Q[0]*Q[2]
2814    T4 = Q[0]*Q[3]
2815    T5 = -Q[1]*Q[1]
2816    T6 = Q[1]*Q[2]
2817    T7 = Q[1]*Q[3]
2818    T8 = -Q[2]*Q[2]
2819    T9 = Q[2]*Q[3]
2820    T10 = -Q[3]*Q[3]
2821    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
2822    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
2823    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
2824    return VP   
2825   
2826def Q2Mat(Q):
2827    ''' make rotation matrix from quaternion
2828        q=r+ai+bj+ck
2829    '''
2830    QN = normQ(Q)
2831    aa = QN[0]**2
2832    ab = QN[0]*QN[1]
2833    ac = QN[0]*QN[2]
2834    ad = QN[0]*QN[3]
2835    bb = QN[1]**2
2836    bc = QN[1]*QN[2]
2837    bd = QN[1]*QN[3]
2838    cc = QN[2]**2
2839    cd = QN[2]*QN[3]
2840    dd = QN[3]**2
2841    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2842        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2843        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2844    return np.array(M)
2845   
2846def AV2Q(A,V):
2847    ''' convert angle (radians) & vector to quaternion
2848        q=r+ai+bj+ck
2849    '''
2850    Q = np.zeros(4)
2851    d = np.sqrt(np.sum(np.array(V)**2))
2852    if d:
2853        V /= d
2854    else:
2855        V = np.array([0.,0.,1.])
2856    p = A/2.
2857    Q[0] = np.cos(p)
2858    Q[1:4] = V*np.sin(p)
2859    return Q
2860   
2861def AVdeg2Q(A,V):
2862    ''' convert angle (degrees) & vector to quaternion
2863        q=r+ai+bj+ck
2864    '''
2865    Q = np.zeros(4)
2866    d = np.sqrt(np.sum(np.array(V)**2))
2867    if d:
2868        V /= d
2869    else:
2870        V = np.array([0.,0.,1.])
2871    p = A/2.
2872    Q[0] = cosd(p)
2873    Q[1:4] = V*sind(p)
2874    return Q
2875   
2876def Q2AVdeg(Q):
2877    ''' convert quaternion to angle (degrees 0-360) & normalized vector
2878        q=r+ai+bj+ck
2879    '''
2880    A = 2.*acosd(Q[0])
2881    V = np.array(Q[1:])
2882    d = np.sqrt(np.sum(V**2))
2883    if d:
2884        V /= d
2885    else:
2886        A = 0.
2887        V = np.array([0.,0.,0.])
2888    return A,V
2889   
2890def Q2AV(Q):
2891    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
2892        q=r+ai+bj+ck
2893    '''
2894    A = 2.*np.arccos(Q[0])
2895    V = np.array(Q[1:])
2896    d = np.sqrt(np.sum(V**2))
2897    if d:
2898        V /= d
2899    else:
2900        A = 0.
2901        V = np.array([0.,0.,0.])
2902    return A,V
2903   
2904def makeQuat(A,B,C):
2905    ''' Make quaternion from rotation of A vector to B vector about C axis
2906
2907        :param np.array A,B,C: Cartesian 3-vectors
2908        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2909    '''
2910
2911    V1 = np.cross(A,C)
2912    V2 = np.cross(B,C)
2913    if nl.norm(V1)*nl.norm(V2):
2914        V1 /= nl.norm(V1)
2915        V2 /= nl.norm(V2)
2916        V3 = np.cross(V1,V2)
2917    else:
2918        V3 = np.zeros(3)
2919    Q = np.array([0.,0.,0.,1.])
2920    D = 0.
2921    if nl.norm(V3):
2922        V3 /= nl.norm(V3)
2923        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2924        D = np.arccos(D1)/2.0
2925        V1 = C-V3
2926        V2 = C+V3
2927        DM = nl.norm(V1)
2928        DP = nl.norm(V2)
2929        S = np.sin(D)
2930        Q[0] = np.cos(D)
2931        Q[1:] = V3*S
2932        D *= 2.
2933        if DM > DP:
2934            D *= -1.
2935    return Q,D
2936   
2937if __name__ == "__main__":
2938    from numpy import cos
2939    # minimum expected at ~-0.195
2940    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
2941    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
2942    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
2943    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
2944
2945    # minimum expected at ~[-0.195, -0.1]
2946    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
2947    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')
2948    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')
2949    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.