source: trunk/GSASIImath.py @ 1622

Last change on this file since 1622 was 1622, checked in by vondreele, 7 years ago

add 4D charge flipping - works to give 3D structure; only enabled for 4D data
fix problems with hkl limits in Fourier calcs. - see adjHKLmax
fix reflection table/2D & 3D display issues
disable MC/SA for 4D data

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