source: trunk/GSASIImath.py @ 1066

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

more MC/SA fixes

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