source: trunk/GSASIImath.py @ 1063

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

change sf calc timing for MC/SA

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