source: trunk/GSASIImath.py @ 1065

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

add MC refinement to MC/SA - works

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 99.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-09-25 18:04:46 +0000 (Wed, 25 Sep 2013) $
5# $Author: vondreele $
6# $Revision: 1065 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1065 2013-09-25 18:04:46Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIImath: computation module*
12================================
13
14Routines for least-squares minimization and other stuff
15
16'''
17import sys
18import os
19import os.path as ospath
20import random as rn
21import numpy as np
22import numpy.linalg as nl
23import numpy.ma as ma
24import cPickle
25import time
26import math
27import copy
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 1065 $")
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 set_range(self,x0,frac):
2091        delrange = frac*(self.upper-self.lower)
2092        self.upper = x0+delrange
2093        self.lower = x0-delrange
2094
2095    def accept_test(self, dE):
2096        T = self.T
2097        self.tests += 1
2098        if dE < 0:
2099            self.accepted += 1
2100            return 1
2101        p = exp(-dE*1.0/self.boltzmann/T)
2102        if (p > random.uniform(0.0, 1.0)):
2103            self.accepted += 1
2104            return 1
2105        return 0
2106
2107    def update_guess(self, x0):
2108        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2109
2110    def update_temp(self, x0):
2111        pass
2112
2113
2114#  A schedule due to Lester Ingber modified to use bounds - OK
2115class fast_sa(base_schedule):
2116    def init(self, **options):
2117        self.__dict__.update(options)
2118        if self.m is None:
2119            self.m = 1.0
2120        if self.n is None:
2121            self.n = 1.0
2122        self.c = self.m * exp(-self.n * self.quench)
2123
2124    def update_guess(self, x0):
2125        x0 = asarray(x0)
2126        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2127        T = self.T
2128        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2129        xnew = xc*(self.upper - self.lower)+self.lower
2130        return xnew
2131#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2132#        xc = y*(self.upper - self.lower)
2133#        xnew = x0 + xc
2134#        return xnew
2135
2136    def update_temp(self):
2137        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2138        self.k += 1
2139        return
2140
2141class cauchy_sa(base_schedule):     #modified to use bounds - not good
2142    def update_guess(self, x0):
2143        x0 = asarray(x0)
2144        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2145        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2146        xnew = xc*(self.upper - self.lower)+self.lower
2147        return xnew
2148#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2149#        xc = self.learn_rate * self.T * tan(numbers)
2150#        xnew = x0 + xc
2151#        return xnew
2152
2153    def update_temp(self):
2154        self.T = self.T0/(1+self.k)
2155        self.k += 1
2156        return
2157
2158class boltzmann_sa(base_schedule):
2159#    def update_guess(self, x0):
2160#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2161#        x0 = asarray(x0)
2162#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2163#
2164#        xnew = x0 + xc*std*self.learn_rate
2165#        return xnew
2166
2167    def update_temp(self):
2168        self.k += 1
2169        self.T = self.T0 / log(self.k+1.0)
2170        return
2171
2172class log_sa(base_schedule):        #OK
2173
2174    def init(self,**options):
2175        self.__dict__.update(options)
2176       
2177    def update_guess(self,x0):     #same as default
2178        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2179       
2180    def update_temp(self):
2181        self.k += 1
2182        self.T = self.T0*self.slope**self.k
2183       
2184class _state(object):
2185    def __init__(self):
2186        self.x = None
2187        self.cost = None
2188
2189# TODO:
2190#     allow for general annealing temperature profile
2191#     in that case use update given by alpha and omega and
2192#     variation of all previous updates and temperature?
2193
2194# Simulated annealing
2195
2196def anneal(func, x0, args=(), schedule='fast', full_output=0,
2197           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2198           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2199           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2200           ranRange=0.10,autoRan=False,dlg=None):
2201    """Minimize a function using simulated annealing.
2202
2203    Schedule is a schedule class implementing the annealing schedule.
2204    Available ones are 'fast', 'cauchy', 'boltzmann'
2205
2206    :param callable func: f(x, \*args)
2207        Function to be optimized.
2208    :param ndarray x0:
2209        Initial guess.
2210    :param tuple args:
2211        Extra parameters to `func`.
2212    :param base_schedule schedule:
2213        Annealing schedule to use (a class).
2214    :param bool full_output:
2215        Whether to return optional outputs.
2216    :param float T0:
2217        Initial Temperature (estimated as 1.2 times the largest
2218        cost-function deviation over random points in the range).
2219    :param float Tf:
2220        Final goal temperature.
2221    :param int maxeval:
2222        Maximum function evaluations.
2223    :param int maxaccept:
2224        Maximum changes to accept.
2225    :param int maxiter:
2226        Maximum cooling iterations.
2227    :param float learn_rate:
2228        Scale constant for adjusting guesses.
2229    :param float boltzmann:
2230        Boltzmann constant in acceptance test
2231        (increase for less stringent test at each temperature).
2232    :param float feps:
2233        Stopping relative error tolerance for the function value in
2234        last four coolings.
2235    :param float quench,m,n:
2236        Parameters to alter fast_sa schedule.
2237    :param float/ndarray lower,upper:
2238        Lower and upper bounds on `x`.
2239    :param int dwell:
2240        The number of times to search the space at each temperature.
2241    :param float slope:
2242        Parameter for log schedule
2243    :param bool ranStart=False:
2244        True for set 10% of ranges about x
2245
2246    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2247
2248     * xmin (ndarray): Point giving smallest value found.
2249     * Jmin (float): Minimum value of function found.
2250     * T (float): Final temperature.
2251     * feval (int): Number of function evaluations.
2252     * iters (int): Number of cooling iterations.
2253     * accept (int): Number of tests accepted.
2254     * retval (int): Flag indicating stopping condition:
2255
2256              *  0: Points no longer changing
2257              *  1: Cooled to final temperature
2258              *  2: Maximum function evaluations
2259              *  3: Maximum cooling iterations reached
2260              *  4: Maximum accepted query locations reached
2261              *  5: Final point not the minimum amongst encountered points
2262
2263    *Notes*:
2264    Simulated annealing is a random algorithm which uses no derivative
2265    information from the function being optimized. In practice it has
2266    been more useful in discrete optimization than continuous
2267    optimization, as there are usually better algorithms for continuous
2268    optimization problems.
2269
2270    Some experimentation by trying the difference temperature
2271    schedules and altering their parameters is likely required to
2272    obtain good performance.
2273
2274    The randomness in the algorithm comes from random sampling in numpy.
2275    To obtain the same results you can call numpy.random.seed with the
2276    same seed immediately before calling scipy.optimize.anneal.
2277
2278    We give a brief description of how the three temperature schedules
2279    generate new points and vary their temperature. Temperatures are
2280    only updated with iterations in the outer loop. The inner loop is
2281    over xrange(dwell), and new points are generated for
2282    every iteration in the inner loop. (Though whether the proposed
2283    new points are accepted is probabilistic.)
2284
2285    For readability, let d denote the dimension of the inputs to func.
2286    Also, let x_old denote the previous state, and k denote the
2287    iteration number of the outer loop. All other variables not
2288    defined below are input variables to scipy.optimize.anneal itself.
2289
2290    In the 'fast' schedule the updates are ::
2291
2292        u ~ Uniform(0, 1, size=d)
2293        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2294        xc = y * (upper - lower)
2295        x_new = x_old + xc
2296
2297        c = n * exp(-n * quench)
2298        T_new = T0 * exp(-c * k**quench)
2299
2300
2301    In the 'cauchy' schedule the updates are ::
2302
2303        u ~ Uniform(-pi/2, pi/2, size=d)
2304        xc = learn_rate * T * tan(u)
2305        x_new = x_old + xc
2306
2307        T_new = T0 / (1+k)
2308
2309    In the 'boltzmann' schedule the updates are ::
2310
2311        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2312        y ~ Normal(0, std, size=d)
2313        x_new = x_old + learn_rate * y
2314
2315        T_new = T0 / log(1+k)
2316
2317    """
2318    x0 = asarray(x0)
2319    lower = asarray(lower)
2320    upper = asarray(upper)
2321
2322    schedule = eval(schedule+'_sa()')
2323    #   initialize the schedule
2324    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2325                  learn_rate=learn_rate, lower=lower, upper=upper,
2326                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2327
2328    current_state, last_state, best_state = _state(), _state(), _state()
2329    if ranStart:
2330        schedule.set_range(x0,ranRange)
2331    if T0 is None:
2332        x0 = schedule.getstart_temp(best_state)
2333    else:
2334        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
2335        best_state.x = None
2336        best_state.cost = numpy.Inf
2337
2338    last_state.x = asarray(x0).copy()
2339    fval = func(x0,*args)
2340    schedule.feval += 1
2341    last_state.cost = fval
2342    if last_state.cost < best_state.cost:
2343        best_state.cost = fval
2344        best_state.x = asarray(x0).copy()
2345    schedule.T = schedule.T0
2346    fqueue = [100, 300, 500, 700]
2347    iters = 1
2348    keepGoing = True
2349    bestn = 0
2350    while keepGoing:
2351        retval = 0
2352        for n in xrange(dwell):
2353            current_state.x = schedule.update_guess(last_state.x)
2354            current_state.cost = func(current_state.x,*args)
2355            schedule.feval += 1
2356
2357            dE = current_state.cost - last_state.cost
2358            if schedule.accept_test(dE):
2359                last_state.x = current_state.x.copy()
2360                last_state.cost = current_state.cost
2361                if last_state.cost < best_state.cost:
2362                    best_state.x = last_state.x.copy()
2363                    best_state.cost = last_state.cost
2364                    bestn = n
2365                    if best_state.cost < 1.0 and autoRan:
2366                        schedule.set_range(x0,best_state.cost/2.)                       
2367        if dlg:
2368            GoOn = dlg.Update(min(100.,best_state.cost*100),
2369                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2370                    'Best trial:',bestn,  \
2371                    'MC/SA Residual:',best_state.cost*100,'%', \
2372                    ))[0]
2373            if not GoOn:
2374                best_state.x = last_state.x.copy()
2375                best_state.cost = last_state.cost
2376                retval = 5
2377        schedule.update_temp()
2378        iters += 1
2379        # Stopping conditions
2380        # 0) last saved values of f from each cooling step
2381        #     are all very similar (effectively cooled)
2382        # 1) Tf is set and we are below it
2383        # 2) maxeval is set and we are past it
2384        # 3) maxiter is set and we are past it
2385        # 4) maxaccept is set and we are past it
2386        # 5) user canceled run via progress bar
2387
2388        fqueue.append(squeeze(last_state.cost))
2389        fqueue.pop(0)
2390        af = asarray(fqueue)*1.0
2391        if retval == 5:
2392            print ' User terminated run; incomplete MC/SA'
2393            keepGoing = False
2394            break
2395        if all(abs((af-af[0])/af[0]) < feps):
2396            retval = 0
2397            if abs(af[-1]-best_state.cost) > feps*10:
2398                retval = 5
2399#                print "Warning: Cooled to %f at %s but this is not" \
2400#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2401#                      + " the smallest point found."
2402            break
2403        if (Tf is not None) and (schedule.T < Tf):
2404            retval = 1
2405            break
2406        if (maxeval is not None) and (schedule.feval > maxeval):
2407            retval = 2
2408            break
2409        if (iters > maxiter):
2410            print "Warning: Maximum number of iterations exceeded."
2411            retval = 3
2412            break
2413        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2414            retval = 4
2415            break
2416
2417    if full_output:
2418        return best_state.x, best_state.cost, schedule.T, \
2419               schedule.feval, iters, schedule.accepted, retval
2420    else:
2421        return best_state.x, retval
2422
2423def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q):
2424    outlist = []
2425    for n in range(iCyc):
2426        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2427        outlist.append(result[0])
2428        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2429    out_q.put(outlist)
2430
2431def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2432    import multiprocessing as mp
2433   
2434    nprocs = mp.cpu_count()
2435    out_q = mp.Queue()
2436    procs = []
2437    iCyc = np.zeros(nprocs)
2438    for i in range(nCyc):
2439        iCyc[i%nprocs] += 1
2440    for i in range(nprocs):
2441        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q))
2442        procs.append(p)
2443        p.start()
2444    resultlist = []
2445    for i in range(nprocs):
2446        resultlist += out_q.get()
2447    for p in procs:
2448        p.join()
2449    return resultlist
2450
2451def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2452    '''default doc string
2453   
2454    :param type name: description
2455   
2456    :returns: type name: description
2457   
2458    '''
2459   
2460    twopi = 2.0*np.pi
2461    global tsum
2462    tsum = 0.
2463   
2464    def getMDparms(item,pfx,parmDict,varyList):
2465        parmDict[pfx+'MDaxis'] = item['axis']
2466        parmDict[pfx+'MDval'] = item['Coef'][0]
2467        if item['Coef'][1]:
2468            varyList += [pfx+'MDval',]
2469            limits = item['Coef'][2]
2470            lower.append(limits[0])
2471            upper.append(limits[1])
2472                       
2473    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2474        parmDict[pfx+'Atype'] = item['atType']
2475        aTypes |= set([item['atType'],]) 
2476        pstr = ['Ax','Ay','Az']
2477        XYZ = [0,0,0]
2478        for i in range(3):
2479            name = pfx+pstr[i]
2480            parmDict[name] = item['Pos'][0][i]
2481            XYZ[i] = parmDict[name]
2482            if item['Pos'][1][i]:
2483                varyList += [name,]
2484                limits = item['Pos'][2][i]
2485                lower.append(limits[0])
2486                upper.append(limits[1])
2487        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2488           
2489    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2490        parmDict[mfx+'MolCent'] = item['MolCent']
2491        parmDict[mfx+'RBId'] = item['RBId']
2492        pstr = ['Px','Py','Pz']
2493        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
2494        for i in range(3):
2495            name = mfx+pstr[i]
2496            parmDict[name] = item['Pos'][0][i]
2497            if item['Pos'][1][i]:
2498                varyList += [name,]
2499                limits = item['Pos'][2][i]
2500                lower.append(limits[0])
2501                upper.append(limits[1])
2502        AV = item['Ori'][0]
2503        A = AV[0]
2504        V = AV[1:]
2505        for i in range(4):
2506            name = mfx+ostr[i]
2507            if i:
2508                parmDict[name] = V[i-1]
2509            else:
2510                parmDict[name] = A
2511            if item['Ovar'] == 'AV':
2512                varyList += [name,]
2513                limits = item['Ori'][2][i]
2514                lower.append(limits[0])
2515                upper.append(limits[1])
2516            elif item['Ovar'] == 'A' and not i:
2517                varyList += [name,]
2518                limits = item['Ori'][2][i]
2519                lower.append(limits[0])
2520                upper.append(limits[1])
2521        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2522            for i in range(len(item['Tor'][0])):
2523                name = mfx+'Tor'+str(i)
2524                parmDict[name] = item['Tor'][0][i]
2525                if item['Tor'][1][i]:
2526                    varyList += [name,]
2527                    limits = item['Tor'][2][i]
2528                    lower.append(limits[0])
2529                    upper.append(limits[1])
2530        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2531        aTypes |= set(atypes)
2532        atNo += len(atypes)
2533        return atNo
2534               
2535    def GetAtomM(Xdata,SGData):
2536        Mdata = []
2537        for xyz in Xdata:
2538            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2539        return np.array(Mdata)
2540       
2541    def GetAtomT(RBdata,parmDict):
2542        'Needs a doc string'
2543        atNo = parmDict['atNo']
2544        nfixAt = parmDict['nfixAt']
2545        Tdata = atNo*[' ',]
2546        for iatm in range(nfixAt):
2547            parm = ':'+str(iatm)+':Atype'
2548            if parm in parmDict:
2549                Tdata[iatm] = aTypes.index(parmDict[parm])
2550        iatm = nfixAt
2551        for iObj in range(parmDict['nObj']):
2552            pfx = str(iObj)+':'
2553            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2554                if parmDict[pfx+'Type'] == 'Vector':
2555                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2556                    nAtm = len(RBRes['rbVect'][0])
2557                else:       #Residue
2558                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2559                    nAtm = len(RBRes['rbXYZ'])
2560                for i in range(nAtm):
2561                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
2562                    iatm += 1
2563            elif parmDict[pfx+'Type'] == 'Atom':
2564                atNo = parmDict[pfx+'atNo']
2565                parm = pfx+'Atype'              #remove extra ':'
2566                if parm in parmDict:
2567                    Tdata[atNo] = aTypes.index(parmDict[parm])
2568                iatm += 1
2569            else:
2570                continue        #skips March Dollase
2571        return Tdata
2572       
2573    def GetAtomX(RBdata,parmDict):
2574        'Needs a doc string'
2575        Bmat = parmDict['Bmat']
2576        atNo = parmDict['atNo']
2577        nfixAt = parmDict['nfixAt']
2578        Xdata = np.zeros((3,atNo))
2579        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2580        for iatm in range(nfixAt):
2581            for key in keys:
2582                parm = ':'+str(iatm)+key
2583                if parm in parmDict:
2584                    keys[key][iatm] = parmDict[parm]
2585        iatm = nfixAt
2586        for iObj in range(parmDict['nObj']):
2587            pfx = str(iObj)+':'
2588            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2589                if parmDict[pfx+'Type'] == 'Vector':
2590                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2591                    vecs = RBRes['rbVect']
2592                    mags = RBRes['VectMag']
2593                    Cart = np.zeros_like(vecs[0])
2594                    for vec,mag in zip(vecs,mags):
2595                        Cart += vec*mag
2596                elif parmDict[pfx+'Type'] == 'Residue':
2597                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2598                    Cart = np.array(RBRes['rbXYZ'])
2599                    for itor,seq in enumerate(RBRes['rbSeq']):
2600                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
2601                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
2602                if parmDict[pfx+'MolCent'][1]:
2603                    Cart -= parmDict[pfx+'MolCent'][0]
2604                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2605                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2606                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
2607                iatm += len(Cart)
2608            elif parmDict[pfx+'Type'] == 'Atom':
2609                atNo = parmDict[pfx+'atNo']
2610                for key in keys:
2611                    parm = pfx+key[1:]              #remove extra ':'
2612                    if parm in parmDict:
2613                        keys[key][atNo] = parmDict[parm]
2614                iatm += 1
2615            else:
2616                continue        #skips March Dollase
2617        return Xdata.T
2618       
2619    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
2620        allX = np.inner(Xdata,SGM)+SGT
2621        allT = np.repeat(Tdata,allX.shape[1])
2622        allM = np.repeat(Mdata,allX.shape[1])
2623        allX = np.reshape(allX,(-1,3))
2624        return allT,allM,allX
2625
2626    def getAllX(Xdata,SGM,SGT):
2627        allX = np.inner(Xdata,SGM)+SGT
2628        allX = np.reshape(allX,(-1,3))
2629        return allX
2630       
2631    def normQuaternions(RBdata,parmDict,varyList,values):
2632        for iObj in range(parmDict['nObj']):
2633            pfx = str(iObj)+':'
2634            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2635                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2636                A,V = Q2AVdeg(Qori)
2637                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
2638                    if i:
2639                        parmDict[pfx+name] = V[i-1]
2640                    else:
2641                        parmDict[pfx+name] = A
2642       
2643    def mcsaCalc(values,refList,rcov,ifInv,allFF,RBdata,varyList,parmDict):
2644        ''' Compute structure factors for all h,k,l for phase
2645        input:
2646            refList: [ref] where each ref = h,k,l,m,d,...
2647            rcov:   array[nref,nref] covariance terms between Fo^2 values
2648            ifInv:  bool True if centrosymmetric
2649            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
2650            RBdata: [dict] rigid body dictionary
2651            varyList: [list] names of varied parameters in MC/SA (not used here)           
2652            ParmDict: [dict] problem parameters
2653        puts result F^2 in each ref[5] in refList
2654        returns:
2655            delt-F*rcov*delt-F/sum(Fo^2)^2
2656        '''       
2657        global tsum
2658        t0 = time.time()
2659        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
2660        Xdata = GetAtomX(RBdata,parmDict)                   #get new atom coords from RB
2661        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
2662        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
2663        MDaxis = parmDict['0:MDaxis']
2664        Gmat = parmDict['Gmat']
2665        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
2666        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
2667        if ifInv:
2668            refList[5] = Aterm
2669        else:
2670            Bterm = refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
2671            refList[5] = Aterm+Bterm
2672        sumFcsq = np.sum(refList[5])
2673        scale = parmDict['sumFosq']/sumFcsq
2674        refList[5] *= scale
2675        refList[6] = refList[4]-refList[5]
2676        M = np.inner(refList[6],np.inner(rcov,refList[6]))
2677        tsum += (time.time()-t0)
2678        return M/np.sum(refList[4]**2)
2679
2680    sq8ln2 = np.sqrt(8*np.log(2))
2681    sq2pi = np.sqrt(2*np.pi)
2682    sq4pi = np.sqrt(4*np.pi)
2683    generalData = data['General']
2684    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2685    Gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])[0]
2686    SGData = generalData['SGData']
2687    SGM = [SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))]
2688    SGT = [SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))]
2689    fixAtoms = data['Atoms']                       #if any
2690    cx,ct,cs = generalData['AtomPtrs'][:3]
2691    aTypes = set([])
2692    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2693    varyList = []
2694    atNo = 0
2695    for atm in fixAtoms:
2696        pfx = ':'+str(atNo)+':'
2697        parmDict[pfx+'Atype'] = atm[ct]
2698        aTypes |= set([atm[ct],])
2699        pstr = ['Ax','Ay','Az']
2700        parmDict[pfx+'Amul'] = atm[cs+1]
2701        for i in range(3):
2702            name = pfx+pstr[i]
2703            parmDict[name] = atm[cx+i]
2704        atNo += 1
2705    parmDict['nfixAt'] = len(fixAtoms)       
2706    MCSA = generalData['MCSA controls']
2707    reflName = MCSA['Data source']
2708    phaseName = generalData['Name']
2709    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2710    upper = []
2711    lower = []
2712    for i,item in enumerate(MCSAObjs):
2713        mfx = str(i)+':'
2714        parmDict[mfx+'Type'] = item['Type']
2715        if item['Type'] == 'MD':
2716            getMDparms(item,mfx,parmDict,varyList)
2717        elif item['Type'] == 'Atom':
2718            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2719            parmDict[mfx+'atNo'] = atNo
2720            atNo += 1
2721        elif item['Type'] in ['Residue','Vector']:
2722            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2723    parmDict['atNo'] = atNo                 #total no. of atoms
2724    parmDict['nObj'] = len(MCSAObjs)
2725    aTypes = list(aTypes)
2726    Tdata = GetAtomT(RBdata,parmDict)
2727    Xdata = GetAtomX(RBdata,parmDict)
2728    Mdata = GetAtomM(Xdata,SGData)
2729    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
2730    FFtables = G2el.GetFFtable(aTypes)
2731    refs = []
2732    allFF = []
2733    sumFosq = 0
2734    if 'PWDR' in reflName:
2735        for ref in reflData:
2736            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2737            if d >= MCSA['dmin']:
2738                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2739                SQ = 0.25/d**2
2740                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2741                refs.append([h,k,l,m,f*m,pos,sig])
2742                sumFosq += f*m
2743        nRef = len(refs)
2744        rcov = np.zeros((nRef,nRef))
2745        for iref,refI in enumerate(refs):
2746            rcov[iref][iref] = 1./(sq4pi*refI[6])
2747            for jref,refJ in enumerate(refs[:iref]):
2748                t1 = refI[6]**2+refJ[6]**2
2749                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2750                if t2 > 10.:
2751                    rcov[iref][jref] = 0.
2752                else:
2753                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2754        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2755        Rdiag = np.sqrt(np.diag(rcov))
2756        Rnorm = np.outer(Rdiag,Rdiag)
2757        rcov /= Rnorm
2758    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
2759        vList = covData['varyList']
2760        for iref,refI in enumerate(reflData):
2761            h,k,l,m,d,v,f,s = refI
2762            if d >= MCSA['dmin'] and v:       #skip unrefined ones
2763                SQ = 0.25/d**2
2764                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2765                refs.append([h,k,l,m,f*m,iref,0.])
2766                sumFosq += f*m
2767        nRef = len(refs)
2768        pfx = str(data['pId'])+'::PWLref:'
2769        if covData['freshCOV'] and generalData['doPawley'] and MCSA.get('newDmin',True):
2770            covMatrix = covData['covMatrix']
2771            rcov = np.zeros((nRef,nRef))       
2772            for iref,refI in enumerate(refs):
2773                I = refI[5]
2774                nameI = pfx+str(I)
2775                if nameI in vList:
2776                    Iindx = vList.index(nameI)
2777                    rcov[iref][iref] = covMatrix[Iindx][Iindx]
2778                    for jref,refJ in enumerate(refs[:iref]):
2779                        J = refJ[5]
2780                        nameJ = pfx+str(J)
2781                        try:
2782                            rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2783                        except ValueError:
2784                            rcov[iref][jref] = rcov[iref][jref-1]
2785                else:
2786                    rcov[iref] = rcov[iref-1]
2787                    rcov[iref][iref] = rcov[iref-1][iref-1]
2788            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2789            Rdiag = np.sqrt(np.diag(rcov))
2790            Rnorm = np.outer(Rdiag,Rdiag)
2791            rcov /= Rnorm
2792            MCSA['rcov'] = rcov
2793            covData['freshCOV'] = False
2794            MCSA['newDmin'] = False
2795        else:
2796            rcov = MCSA['rcov']
2797    elif 'HKLF' in reflName:
2798        for ref in reflData:
2799            [h,k,l,m,d],f = ref[:5],ref[6]
2800            if d >= MCSA['dmin']:
2801                SQ = 0.25/d**2
2802                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2803                refs.append([h,k,l,m,f*m,0.,0.])
2804                sumFosq += f*m
2805        nRef = len(refs)
2806        rcov = np.identity(len(refs))
2807    allFF = np.array(allFF).T
2808    refs = np.array(refs).T
2809    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
2810    print ' Number of parameters varied: %d'%(len(varyList))
2811    parmDict['sumFosq'] = sumFosq
2812    x0 = [parmDict[val] for val in varyList]
2813    ifInv = SGData['SGInv']
2814    results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,allFF,RBdata,varyList,parmDict),
2815        schedule=MCSA['Algorithm'], full_output=True,
2816        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2817        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
2818        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2819        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
2820        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
2821    M = mcsaCalc(results[0],refs,rcov,ifInv,allFF,RBdata,varyList,parmDict)
2822#    for ref in refs.T:
2823#        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])
2824#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
2825    Result = [False,False,results[1],results[2],]+list(results[0])
2826    Result.append(varyList)
2827    return Result,tsum
2828
2829       
2830################################################################################
2831##### Quaternion stuff
2832################################################################################
2833
2834def prodQQ(QA,QB):
2835    ''' Grassman quaternion product
2836        QA,QB quaternions; q=r+ai+bj+ck
2837    '''
2838    D = np.zeros(4)
2839    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2840    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2841    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2842    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2843   
2844#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
2845#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
2846   
2847    return D
2848   
2849def normQ(QA):
2850    ''' get length of quaternion & normalize it
2851        q=r+ai+bj+ck
2852    '''
2853    n = np.sqrt(np.sum(np.array(QA)**2))
2854    return QA/n
2855   
2856def invQ(Q):
2857    '''
2858        get inverse of quaternion
2859        q=r+ai+bj+ck; q* = r-ai-bj-ck
2860    '''
2861    return Q*np.array([1,-1,-1,-1])
2862   
2863def prodQVQ(Q,V):
2864    """
2865    compute the quaternion vector rotation qvq-1 = v'
2866    q=r+ai+bj+ck
2867    """
2868    T2 = Q[0]*Q[1]
2869    T3 = Q[0]*Q[2]
2870    T4 = Q[0]*Q[3]
2871    T5 = -Q[1]*Q[1]
2872    T6 = Q[1]*Q[2]
2873    T7 = Q[1]*Q[3]
2874    T8 = -Q[2]*Q[2]
2875    T9 = Q[2]*Q[3]
2876    T10 = -Q[3]*Q[3]
2877    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
2878    VP = 2.*np.inner(V,M)
2879    return VP+V
2880   
2881def Q2Mat(Q):
2882    ''' make rotation matrix from quaternion
2883        q=r+ai+bj+ck
2884    '''
2885    QN = normQ(Q)
2886    aa = QN[0]**2
2887    ab = QN[0]*QN[1]
2888    ac = QN[0]*QN[2]
2889    ad = QN[0]*QN[3]
2890    bb = QN[1]**2
2891    bc = QN[1]*QN[2]
2892    bd = QN[1]*QN[3]
2893    cc = QN[2]**2
2894    cd = QN[2]*QN[3]
2895    dd = QN[3]**2
2896    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2897        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2898        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2899    return np.array(M)
2900   
2901def AV2Q(A,V):
2902    ''' convert angle (radians) & vector to quaternion
2903        q=r+ai+bj+ck
2904    '''
2905    Q = np.zeros(4)
2906    d = nl.norm(np.array(V))
2907    if d:
2908        V /= d
2909        if not A:       #==0.
2910            A = 2.*np.pi
2911        p = A/2.
2912        Q[0] = np.cos(p)
2913        Q[1:4] = V*np.sin(p)
2914    else:
2915        Q[3] = 1.
2916    return Q
2917   
2918def AVdeg2Q(A,V):
2919    ''' convert angle (degrees) & vector to quaternion
2920        q=r+ai+bj+ck
2921    '''
2922    Q = np.zeros(4)
2923    d = nl.norm(np.array(V))
2924    if not A:       #== 0.!
2925        A = 360.
2926    if d:
2927        V /= d
2928        p = A/2.
2929        Q[0] = cosd(p)
2930        Q[1:4] = V*sind(p)
2931    else:
2932        Q[3] = 1.
2933    return Q
2934   
2935def Q2AVdeg(Q):
2936    ''' convert quaternion to angle (degrees 0-360) & normalized vector
2937        q=r+ai+bj+ck
2938    '''
2939    A = 2.*acosd(Q[0])
2940    V = np.array(Q[1:])
2941    V /= sind(A/2.)
2942    return A,V
2943   
2944def Q2AV(Q):
2945    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
2946        q=r+ai+bj+ck
2947    '''
2948    A = 2.*np.arccos(Q[0])
2949    V = np.array(Q[1:])
2950    V /= np.sin(A/2.)
2951    return A,V
2952   
2953def randomQ(r0,r1,r2,r3):
2954    ''' create random quaternion from 4 random numbers in range (-1,1)
2955    '''
2956    sum = 0
2957    Q = np.array(4)
2958    Q[0] = r0
2959    sum += Q[0]**2
2960    Q[1] = np.sqrt(1.-sum)*r1
2961    sum += Q[1]**2
2962    Q[2] = np.sqrt(1.-sum)*r2
2963    sum += Q[2]**2
2964    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
2965    return Q
2966   
2967def randomAVdeg(r0,r1,r2,r3):
2968    ''' create random angle (deg),vector from 4 random number in range (-1,1)
2969    '''
2970    return Q2AVdeg(randomQ(r0,r1,r2,r3))
2971   
2972def makeQuat(A,B,C):
2973    ''' Make quaternion from rotation of A vector to B vector about C axis
2974
2975        :param np.array A,B,C: Cartesian 3-vectors
2976        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2977    '''
2978
2979    V1 = np.cross(A,C)
2980    V2 = np.cross(B,C)
2981    if nl.norm(V1)*nl.norm(V2):
2982        V1 /= nl.norm(V1)
2983        V2 /= nl.norm(V2)
2984        V3 = np.cross(V1,V2)
2985    else:
2986        V3 = np.zeros(3)
2987    Q = np.array([0.,0.,0.,1.])
2988    D = 0.
2989    if nl.norm(V3):
2990        V3 /= nl.norm(V3)
2991        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2992        D = np.arccos(D1)/2.0
2993        V1 = C-V3
2994        V2 = C+V3
2995        DM = nl.norm(V1)
2996        DP = nl.norm(V2)
2997        S = np.sin(D)
2998        Q[0] = np.cos(D)
2999        Q[1:] = V3*S
3000        D *= 2.
3001        if DM > DP:
3002            D *= -1.
3003    return Q,D
3004   
3005if __name__ == "__main__":
3006    from numpy import cos
3007    # minimum expected at ~-0.195
3008    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3009    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3010    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3011    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3012
3013    # minimum expected at ~[-0.195, -0.1]
3014    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
3015    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')
3016    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')
3017    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.