source: trunk/GSASIImath.py @ 978

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

make MC/SA fortran stuff all REAL*8
add glFreeType.py - doesn't work correctly yet
implement glFreeType font in OpenGl? labels
add new menu items for RBs
MCSA anneal in a test mode no random start for fixed T0

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