source: trunk/GSASIImath.py @ 1643

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

fix input of TOF hkl data
fix bond angle calculations
fix 2D & 3D hkl plotting

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