source: trunk/GSASIImath.py @ 1639

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

remove a couple of error messages from G2math.
implement use of multiple reflection sets for Fourier maps & charge flipping; both 3 & 4D cases - needed for neutron TOF single crystal data.

  • 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-09 21:11:16 +0000 (Mon, 09 Feb 2015) $
5# $Author: vondreele $
6# $Revision: 1639 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1639 2015-02-09 21:11:16Z 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: 1639 $")
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
1093        TxT = G2spc.MoveToUnitCell(TxT)+U
1094        return np.inner(Amat,(TxT-Ox))
1095       
1096    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
1097        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
1098        VecA /= np.sqrt(np.sum(VecA**2))
1099        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
1100        VecB /= np.sqrt(np.sum(VecB**2))
1101        edge = VecB-VecA
1102        edge = np.sum(edge**2)
1103        angle = (2.-edge)/2.
1104        angle = max(angle,-1.)
1105        return acosd(angle)
1106       
1107    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
1108    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
1109    invA = invB = 1
1110    invA = TopA/abs(TopA)
1111    invB = TopB/abs(TopB)
1112    centA = abs(TopA)/100
1113    centB = abs(TopB)/100
1114    opA = abs(TopA)%100-1
1115    opB = abs(TopB)%100-1
1116    MA,TA = SGData['SGOps'][opA]
1117    MB,TB = SGData['SGOps'][opB]
1118    CA = SGData['SGCen'][centA]
1119    CB = SGData['SGCen'][centB]
1120    if 'covMatrix' in covData:
1121        covMatrix = covData['covMatrix']
1122        varyList = covData['varyList']
1123        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
1124        dx = .00001
1125        dadx = np.zeros(9)
1126        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1127        for i in [0,1,2]:
1128            OxA[i] -= dx
1129            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1130            OxA[i] += 2*dx
1131            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1132            OxA[i] -= dx
1133           
1134            TxA[i] -= dx
1135            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1136            TxA[i] += 2*dx
1137            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1138            TxA[i] -= dx
1139           
1140            TxB[i] -= dx
1141            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1142            TxB[i] += 2*dx
1143            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1144            TxB[i] -= dx
1145           
1146        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1147        if sigAng < 0.01:
1148            sigAng = 0.0
1149        return Ang,sigAng
1150    else:
1151        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
1152
1153def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
1154    '''default doc string
1155   
1156    :param type name: description
1157   
1158    :returns: type name: description
1159   
1160    '''
1161    def calcDist(Atoms,SyOps,Amat):
1162        XYZ = []
1163        for i,atom in enumerate(Atoms):
1164            Inv,M,T,C,U = SyOps[i]
1165            XYZ.append(np.array(atom[1:4]))
1166            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1167            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1168        V1 = XYZ[1]-XYZ[0]
1169        return np.sqrt(np.sum(V1**2))
1170       
1171    Inv = []
1172    SyOps = []
1173    names = []
1174    for i,atom in enumerate(Oatoms):
1175        names += atom[-1]
1176        Op,unit = Atoms[i][-1]
1177        inv = Op/abs(Op)
1178        m,t = SGData['SGOps'][abs(Op)%100-1]
1179        c = SGData['SGCen'][abs(Op)/100]
1180        SyOps.append([inv,m,t,c,unit])
1181    Dist = calcDist(Oatoms,SyOps,Amat)
1182   
1183    sig = -0.001
1184    if 'covMatrix' in covData:
1185        parmNames = []
1186        dx = .00001
1187        dadx = np.zeros(6)
1188        for i in range(6):
1189            ia = i/3
1190            ix = i%3
1191            Oatoms[ia][ix+1] += dx
1192            a0 = calcDist(Oatoms,SyOps,Amat)
1193            Oatoms[ia][ix+1] -= 2*dx
1194            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1195        covMatrix = covData['covMatrix']
1196        varyList = covData['varyList']
1197        DistVcov = getVCov(names,varyList,covMatrix)
1198        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1199        if sig < 0.001:
1200            sig = -0.001
1201   
1202    return Dist,sig
1203
1204def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
1205    '''default doc string
1206   
1207    :param type name: description
1208   
1209    :returns: type name: description
1210   
1211    '''
1212
1213    def calcAngle(Atoms,SyOps,Amat):
1214        XYZ = []
1215        for i,atom in enumerate(Atoms):
1216            Inv,M,T,C,U = SyOps[i]
1217            XYZ.append(np.array(atom[1:4]))
1218            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1219            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1220        V1 = XYZ[1]-XYZ[0]
1221        V1 /= np.sqrt(np.sum(V1**2))
1222        V2 = XYZ[1]-XYZ[2]
1223        V2 /= np.sqrt(np.sum(V2**2))
1224        V3 = V2-V1
1225        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1226        return acosd(cang)
1227
1228    Inv = []
1229    SyOps = []
1230    names = []
1231    for i,atom in enumerate(Oatoms):
1232        names += atom[-1]
1233        Op,unit = Atoms[i][-1]
1234        inv = Op/abs(Op)
1235        m,t = SGData['SGOps'][abs(Op)%100-1]
1236        c = SGData['SGCen'][abs(Op)/100]
1237        SyOps.append([inv,m,t,c,unit])
1238    Angle = calcAngle(Oatoms,SyOps,Amat)
1239   
1240    sig = -0.01
1241    if 'covMatrix' in covData:
1242        parmNames = []
1243        dx = .00001
1244        dadx = np.zeros(9)
1245        for i in range(9):
1246            ia = i/3
1247            ix = i%3
1248            Oatoms[ia][ix+1] += dx
1249            a0 = calcAngle(Oatoms,SyOps,Amat)
1250            Oatoms[ia][ix+1] -= 2*dx
1251            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1252        covMatrix = covData['covMatrix']
1253        varyList = covData['varyList']
1254        AngVcov = getVCov(names,varyList,covMatrix)
1255        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1256        if sig < 0.01:
1257            sig = -0.01
1258   
1259    return Angle,sig
1260
1261def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1262    '''default doc string
1263   
1264    :param type name: description
1265   
1266    :returns: type name: description
1267   
1268    '''
1269
1270    def calcTorsion(Atoms,SyOps,Amat):
1271       
1272        XYZ = []
1273        for i,atom in enumerate(Atoms):
1274            Inv,M,T,C,U = SyOps[i]
1275            XYZ.append(np.array(atom[1:4]))
1276            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1277            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1278        V1 = XYZ[1]-XYZ[0]
1279        V2 = XYZ[2]-XYZ[1]
1280        V3 = XYZ[3]-XYZ[2]
1281        V1 /= np.sqrt(np.sum(V1**2))
1282        V2 /= np.sqrt(np.sum(V2**2))
1283        V3 /= np.sqrt(np.sum(V3**2))
1284        M = np.array([V1,V2,V3])
1285        D = nl.det(M)
1286        Ang = 1.0
1287        P12 = np.dot(V1,V2)
1288        P13 = np.dot(V1,V3)
1289        P23 = np.dot(V2,V3)
1290        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1291        return Tors
1292           
1293    Inv = []
1294    SyOps = []
1295    names = []
1296    for i,atom in enumerate(Oatoms):
1297        names += atom[-1]
1298        Op,unit = Atoms[i][-1]
1299        inv = Op/abs(Op)
1300        m,t = SGData['SGOps'][abs(Op)%100-1]
1301        c = SGData['SGCen'][abs(Op)/100]
1302        SyOps.append([inv,m,t,c,unit])
1303    Tors = calcTorsion(Oatoms,SyOps,Amat)
1304   
1305    sig = -0.01
1306    if 'covMatrix' in covData:
1307        parmNames = []
1308        dx = .00001
1309        dadx = np.zeros(12)
1310        for i in range(12):
1311            ia = i/3
1312            ix = i%3
1313            Oatoms[ia][ix+1] -= dx
1314            a0 = calcTorsion(Oatoms,SyOps,Amat)
1315            Oatoms[ia][ix+1] += 2*dx
1316            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1317            Oatoms[ia][ix+1] -= dx           
1318        covMatrix = covData['covMatrix']
1319        varyList = covData['varyList']
1320        TorVcov = getVCov(names,varyList,covMatrix)
1321        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1322        if sig < 0.01:
1323            sig = -0.01
1324   
1325    return Tors,sig
1326       
1327def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1328    '''default doc string
1329   
1330    :param type name: description
1331   
1332    :returns: type name: description
1333   
1334    '''
1335
1336    def calcDist(Atoms,SyOps,Amat):
1337        XYZ = []
1338        for i,atom in enumerate(Atoms):
1339            Inv,M,T,C,U = SyOps[i]
1340            XYZ.append(np.array(atom[1:4]))
1341            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1342            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1343        V1 = XYZ[1]-XYZ[0]
1344        return np.sqrt(np.sum(V1**2))
1345       
1346    def calcAngle(Atoms,SyOps,Amat):
1347        XYZ = []
1348        for i,atom in enumerate(Atoms):
1349            Inv,M,T,C,U = SyOps[i]
1350            XYZ.append(np.array(atom[1:4]))
1351            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1352            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1353        V1 = XYZ[1]-XYZ[0]
1354        V1 /= np.sqrt(np.sum(V1**2))
1355        V2 = XYZ[1]-XYZ[2]
1356        V2 /= np.sqrt(np.sum(V2**2))
1357        V3 = V2-V1
1358        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1359        return acosd(cang)
1360
1361    def calcTorsion(Atoms,SyOps,Amat):
1362       
1363        XYZ = []
1364        for i,atom in enumerate(Atoms):
1365            Inv,M,T,C,U = SyOps[i]
1366            XYZ.append(np.array(atom[1:4]))
1367            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1368            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1369        V1 = XYZ[1]-XYZ[0]
1370        V2 = XYZ[2]-XYZ[1]
1371        V3 = XYZ[3]-XYZ[2]
1372        V1 /= np.sqrt(np.sum(V1**2))
1373        V2 /= np.sqrt(np.sum(V2**2))
1374        V3 /= np.sqrt(np.sum(V3**2))
1375        M = np.array([V1,V2,V3])
1376        D = nl.det(M)
1377        Ang = 1.0
1378        P12 = np.dot(V1,V2)
1379        P13 = np.dot(V1,V3)
1380        P23 = np.dot(V2,V3)
1381        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1382        return Tors
1383           
1384    Inv = []
1385    SyOps = []
1386    names = []
1387    for i,atom in enumerate(Oatoms):
1388        names += atom[-1]
1389        Op,unit = Atoms[i][-1]
1390        inv = Op/abs(Op)
1391        m,t = SGData['SGOps'][abs(Op)%100-1]
1392        c = SGData['SGCen'][abs(Op)/100]
1393        SyOps.append([inv,m,t,c,unit])
1394    M = len(Oatoms)
1395    if M == 2:
1396        Val = calcDist(Oatoms,SyOps,Amat)
1397    elif M == 3:
1398        Val = calcAngle(Oatoms,SyOps,Amat)
1399    else:
1400        Val = calcTorsion(Oatoms,SyOps,Amat)
1401   
1402    sigVals = [-0.001,-0.01,-0.01]
1403    sig = sigVals[M-3]
1404    if 'covMatrix' in covData:
1405        parmNames = []
1406        dx = .00001
1407        N = M*3
1408        dadx = np.zeros(N)
1409        for i in range(N):
1410            ia = i/3
1411            ix = i%3
1412            Oatoms[ia][ix+1] += dx
1413            if M == 2:
1414                a0 = calcDist(Oatoms,SyOps,Amat)
1415            elif M == 3:
1416                a0 = calcAngle(Oatoms,SyOps,Amat)
1417            else:
1418                a0 = calcTorsion(Oatoms,SyOps,Amat)
1419            Oatoms[ia][ix+1] -= 2*dx
1420            if M == 2:
1421                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1422            elif M == 3:
1423                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1424            else:
1425                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1426        covMatrix = covData['covMatrix']
1427        varyList = covData['varyList']
1428        Vcov = getVCov(names,varyList,covMatrix)
1429        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1430        if sig < sigVals[M-3]:
1431            sig = sigVals[M-3]
1432   
1433    return Val,sig
1434       
1435def ValEsd(value,esd=0,nTZ=False):
1436    '''Format a floating point number with a given level of precision or
1437    with in crystallographic format with a "esd", as value(esd). If esd is
1438    negative the number is formatted with the level of significant figures
1439    appropriate if abs(esd) were the esd, but the esd is not included.
1440    if the esd is zero, approximately 6 significant figures are printed.
1441    nTZ=True causes "extra" zeros to be removed after the decimal place.
1442    for example:
1443
1444      * "1.235(3)" for value=1.2346 & esd=0.003
1445      * "1.235(3)e4" for value=12346. & esd=30
1446      * "1.235(3)e6" for value=0.12346e7 & esd=3000
1447      * "1.235" for value=1.2346 & esd=-0.003
1448      * "1.240" for value=1.2395 & esd=-0.003
1449      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
1450      * "1.23460" for value=1.2346 & esd=0.0
1451
1452    :param float value: number to be formatted
1453    :param float esd: uncertainty or if esd < 0, specifies level of
1454      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
1455    :param bool nTZ: True to remove trailing zeros (default is False)
1456    :returns: value(esd) or value as a string
1457
1458    '''
1459    # Note: this routine is Python 3 compatible -- I think
1460    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
1461    if math.isnan(value): # invalid value, bail out
1462        return '?'
1463    if math.isnan(esd): # invalid esd, treat as zero
1464        esd = 0
1465        esdoff = 5
1466    elif esd != 0:
1467        # transform the esd to a one or two digit integer
1468        l = math.log10(abs(esd)) % 1.
1469        if l < math.log10(cutoff): l+= 1.       
1470        intesd = int(round(10**l)) # esd as integer
1471        # determine the number of digits offset for the esd
1472        esdoff = int(round(math.log10(intesd*1./abs(esd))))
1473    else:
1474        esdoff = 5
1475    valoff = 0
1476    if abs(value) < abs(esdoff): # value is effectively zero
1477        pass
1478    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
1479        # where the digit offset is to the left of the decimal place or where too many
1480        # digits are needed
1481        if abs(value) > 1:
1482            valoff = int(math.log10(abs(value)))
1483        elif abs(value) > 0:
1484            valoff = int(math.log10(abs(value))-0.9999999)
1485        else:
1486            valoff = 0
1487    if esd != 0:
1488        if valoff+esdoff < 0:
1489            valoff = esdoff = 0
1490        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1491    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1492        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1493    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1494        if abs(value) > 0:           
1495            extra = -math.log10(abs(value))
1496        else:
1497            extra = 0
1498        if extra > 0: extra += 1
1499        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1500    if esd > 0:
1501        out += ("({:d})").format(intesd)  # add the esd
1502    elif nTZ and '.' in out:
1503        out = out.rstrip('0')  # strip zeros to right of decimal
1504        out = out.rstrip('.')  # and decimal place when not needed
1505    if valoff != 0:
1506        out += ("e{:d}").format(valoff) # add an exponent, when needed
1507    return out
1508
1509################################################################################
1510##### Fourier & charge flip stuff
1511################################################################################
1512
1513def adjHKLmax(SGData,Hmax):
1514    '''default doc string
1515   
1516    :param type name: description
1517   
1518    :returns: type name: description
1519   
1520    '''
1521    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1522        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
1523        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
1524        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1525    else:
1526        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
1527        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
1528        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1529
1530def OmitMap(data,reflDict,pgbar=None):
1531    '''default doc string
1532   
1533    :param type name: description
1534   
1535    :returns: type name: description
1536   
1537    '''
1538    generalData = data['General']
1539    if not generalData['Map']['MapType']:
1540        print '**** ERROR - Fourier map not defined'
1541        return
1542    mapData = generalData['Map']
1543    dmin = mapData['Resolution']
1544    SGData = generalData['SGData']
1545    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1546    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1547    cell = generalData['Cell'][1:8]       
1548    A = G2lat.cell2A(cell[:6])
1549    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1550    adjHKLmax(SGData,Hmax)
1551    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1552    time0 = time.time()
1553    for iref,ref in enumerate(reflDict['RefList']):
1554        if ref[4] >= dmin:
1555            Fosq,Fcsq,ph = ref[8:11]
1556            Uniq = np.inner(ref[:3],SGMT)
1557            Phi = np.inner(ref[:3],SGT)
1558            for i,hkl in enumerate(Uniq):        #uses uniq
1559                hkl = np.asarray(hkl,dtype='i')
1560                dp = 360.*Phi[i]                #and phi
1561                a = cosd(ph+dp)
1562                b = sind(ph+dp)
1563                phasep = complex(a,b)
1564                phasem = complex(a,-b)
1565                Fo = np.sqrt(Fosq)
1566                if '2Fo-Fc' in mapData['MapType']:
1567                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1568                else:
1569                    F = np.sqrt(Fosq)
1570                h,k,l = hkl+Hmax
1571                Fhkl[h,k,l] = F*phasep
1572                h,k,l = -hkl+Hmax
1573                Fhkl[h,k,l] = F*phasem
1574    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1575    M = np.mgrid[0:4,0:4,0:4]
1576    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
1577    iBeg = blkIds*rho0.shape/4
1578    iFin = (blkIds+1)*rho0.shape/4
1579    rho_omit = np.zeros_like(rho0)
1580    nBlk = 0
1581    for iB,iF in zip(iBeg,iFin):
1582        rho1 = np.copy(rho0)
1583        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
1584        Fnew = fft.ifftshift(fft.ifftn(rho1))
1585        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
1586        phase = Fnew/np.absolute(Fnew)
1587        OFhkl = np.absolute(Fhkl)*phase
1588        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
1589        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]])
1590        nBlk += 1
1591        pgbar.Update(nBlk)
1592    mapData['rho'] = np.real(rho_omit)/cell[6]
1593    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1594    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1595    return mapData
1596   
1597def Fourier4DMap(data,reflDict):
1598    '''default doc string
1599   
1600    :param type name: description
1601   
1602    :returns: type name: description
1603   
1604    '''
1605    generalData = data['General']
1606    mapData = generalData['4DmapData']
1607    dmin = mapData['Resolution']
1608    SGData = generalData['SGData']
1609    SSGData = generalData['SSGData']
1610    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1611    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1612    cell = generalData['Cell'][1:8]       
1613    A = G2lat.cell2A(cell[:6])
1614    maxM = generalData['SuperVec'][2]
1615    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
1616    adjHKLmax(SGData,Hmax)
1617    Hmax = np.asarray(Hmax,dtype='i')+1
1618    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1619    time0 = time.time()
1620    for iref,ref in enumerate(reflDict['RefList']):
1621        if ref[5] > dmin:
1622            Fosq,Fcsq,ph = ref[9:12]
1623            Uniq = np.inner(ref[:4],SSGMT)
1624            Phi = np.inner(ref[:4],SSGT)
1625            for i,hkl in enumerate(Uniq):        #uses uniq
1626                hkl = np.asarray(hkl,dtype='i')
1627                dp = 360.*Phi[i]                #and phi
1628                a = cosd(ph+dp)
1629                b = sind(ph+dp)
1630                phasep = complex(a,b)
1631                phasem = complex(a,-b)
1632                if 'Fobs' in mapData['MapType']:
1633                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1634                    h,k,l,m = hkl+Hmax
1635                    Fhkl[h,k,l,m] = F*phasep
1636                    h,k,l,m = -hkl+Hmax
1637                    Fhkl[h,k,l,m] = F*phasem
1638                elif 'delt-F' in mapData['MapType']:
1639                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1640                    h,k,l,m = hkl+Hmax
1641                    Fhkl[h,k,l,m] = dF*phasep
1642                    h,k,l,m = -hkl+Hmax
1643                    Fhkl[h,k,l,m] = dF*phasem
1644    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1645    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1646    mapData['rho'] = np.real(rho)
1647    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1648    return mapData
1649
1650def FourierMap(data,reflDict):
1651    '''default doc string
1652   
1653    :param type name: description
1654   
1655    :returns: type name: description
1656   
1657    '''
1658    generalData = data['General']
1659    mapData = generalData['Map']
1660    dmin = mapData['Resolution']
1661    SGData = generalData['SGData']
1662    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1663    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1664    cell = generalData['Cell'][1:8]       
1665    A = G2lat.cell2A(cell[:6])
1666    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1667    adjHKLmax(SGData,Hmax)
1668    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1669#    Fhkl[0,0,0] = generalData['F000X']
1670    time0 = time.time()
1671    for iref,ref in enumerate(reflDict['RefList']):
1672        if ref[4] > dmin:
1673            Fosq,Fcsq,ph = ref[8:11]
1674            Uniq = np.inner(ref[:3],SGMT)
1675            Phi = np.inner(ref[:3],SGT)
1676            for i,hkl in enumerate(Uniq):        #uses uniq
1677                hkl = np.asarray(hkl,dtype='i')
1678                dp = 360.*Phi[i]                #and phi
1679                a = cosd(ph+dp)
1680                b = sind(ph+dp)
1681                phasep = complex(a,b)
1682                phasem = complex(a,-b)
1683                if 'Fobs' in mapData['MapType']:
1684                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1685                    h,k,l = hkl+Hmax
1686                    Fhkl[h,k,l] = F*phasep
1687                    h,k,l = -hkl+Hmax
1688                    Fhkl[h,k,l] = F*phasem
1689                elif 'Fcalc' in mapData['MapType']:
1690                    F = np.sqrt(Fcsq)
1691                    h,k,l = hkl+Hmax
1692                    Fhkl[h,k,l] = F*phasep
1693                    h,k,l = -hkl+Hmax
1694                    Fhkl[h,k,l] = F*phasem
1695                elif 'delt-F' in mapData['MapType']:
1696                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1697                    h,k,l = hkl+Hmax
1698                    Fhkl[h,k,l] = dF*phasep
1699                    h,k,l = -hkl+Hmax
1700                    Fhkl[h,k,l] = dF*phasem
1701                elif '2*Fo-Fc' in mapData['MapType']:
1702                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1703                    h,k,l = hkl+Hmax
1704                    Fhkl[h,k,l] = F*phasep
1705                    h,k,l = -hkl+Hmax
1706                    Fhkl[h,k,l] = F*phasem
1707                elif 'Patterson' in mapData['MapType']:
1708                    h,k,l = hkl+Hmax
1709                    Fhkl[h,k,l] = complex(Fosq,0.)
1710                    h,k,l = -hkl+Hmax
1711                    Fhkl[h,k,l] = complex(Fosq,0.)
1712    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1713    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1714    mapData['rho'] = np.real(rho)
1715    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1716    return mapData
1717   
1718# map printing for testing purposes
1719def printRho(SGLaue,rho,rhoMax):                         
1720    '''default doc string
1721   
1722    :param type name: description
1723   
1724    :returns: type name: description
1725   
1726    '''
1727    dim = len(rho.shape)
1728    if dim == 2:
1729        ix,jy = rho.shape
1730        for j in range(jy):
1731            line = ''
1732            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1733                line += (jy-j)*'  '
1734            for i in range(ix):
1735                r = int(100*rho[i,j]/rhoMax)
1736                line += '%4d'%(r)
1737            print line+'\n'
1738    else:
1739        ix,jy,kz = rho.shape
1740        for k in range(kz):
1741            print 'k = ',k
1742            for j in range(jy):
1743                line = ''
1744                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1745                    line += (jy-j)*'  '
1746                for i in range(ix):
1747                    r = int(100*rho[i,j,k]/rhoMax)
1748                    line += '%4d'%(r)
1749                print line+'\n'
1750## keep this
1751               
1752def findOffset(SGData,A,Fhkl):   
1753    '''default doc string
1754   
1755    :param type name: description
1756   
1757    :returns: type name: description
1758   
1759    '''
1760    if SGData['SpGrp'] == 'P 1':
1761        return [0,0,0]   
1762    hklShape = Fhkl.shape
1763    hklHalf = np.array(hklShape)/2
1764    sortHKL = np.argsort(Fhkl.flatten())
1765    Fdict = {}
1766    for hkl in sortHKL:
1767        HKL = np.unravel_index(hkl,hklShape)
1768        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1769        if F == 0.:
1770            break
1771        Fdict['%.6f'%(np.absolute(F))] = hkl
1772    Flist = np.flipud(np.sort(Fdict.keys()))
1773    F = str(1.e6)
1774    i = 0
1775    DH = []
1776    Dphi = []
1777    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1778    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1779    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1780    for F in Flist:
1781        hkl = np.unravel_index(Fdict[F],hklShape)
1782        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
1783            continue
1784        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1785        Uniq = np.array(Uniq,dtype='i')
1786        Phi = np.array(Phi)
1787        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1788        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1789        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1790        ang0 = np.angle(Fh0,deg=True)/360.
1791        for H,phi in zip(Uniq,Phi)[1:]:
1792            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1793            dH = H-hkl
1794            dang = ang-ang0
1795            DH.append(dH)
1796            Dphi.append((dang+.5) % 1.0)
1797        if i > 20 or len(DH) > 30:
1798            break
1799        i += 1
1800    DH = np.array(DH)
1801    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1802    Dphi = np.array(Dphi)
1803    steps = np.array(hklShape)
1804    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1805    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1806    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1807    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1808    hist,bins = np.histogram(Mmap,bins=1000)
1809#    for i,item in enumerate(hist[:10]):
1810#        print item,bins[i]
1811    chisq = np.min(Mmap)
1812    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1813    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1814#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1815    return DX
1816   
1817def ChargeFlip(data,reflDict,pgbar):
1818    '''default doc string
1819   
1820    :param type name: description
1821   
1822    :returns: type name: description
1823   
1824    '''
1825    generalData = data['General']
1826    mapData = generalData['Map']
1827    flipData = generalData['Flip']
1828    FFtable = {}
1829    if 'None' not in flipData['Norm element']:
1830        normElem = flipData['Norm element'].upper()
1831        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1832        for ff in FFs:
1833            if ff['Symbol'] == normElem:
1834                FFtable.update(ff)
1835    dmin = flipData['Resolution']
1836    SGData = generalData['SGData']
1837    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1838    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1839    cell = generalData['Cell'][1:8]       
1840    A = G2lat.cell2A(cell[:6])
1841    Vol = cell[6]
1842    im = 0
1843    if generalData['Type'] in ['modulated','magnetic',]:
1844        im = 1
1845    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1846    adjHKLmax(SGData,Hmax)
1847    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1848    time0 = time.time()
1849    for iref,ref in enumerate(reflDict['RefList']):
1850        dsp = ref[4+im]
1851        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
1852            continue
1853        if dsp > dmin:
1854            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1855            if FFtable:
1856                SQ = 0.25/dsp**2
1857                ff *= G2el.ScatFac(FFtable,SQ)[0]
1858            if ref[8+im] > 0.:         #use only +ve Fobs**2
1859                E = np.sqrt(ref[8+im])/ff
1860            else:
1861                E = 0.
1862            ph = ref[10]
1863            ph = rn.uniform(0.,360.)
1864            Uniq = np.inner(ref[:3],SGMT)
1865            Phi = np.inner(ref[:3],SGT)
1866            for i,hkl in enumerate(Uniq):        #uses uniq
1867                hkl = np.asarray(hkl,dtype='i')
1868                dp = 360.*Phi[i]                #and phi
1869                a = cosd(ph+dp)
1870                b = sind(ph+dp)
1871                phasep = complex(a,b)
1872                phasem = complex(a,-b)
1873                h,k,l = hkl+Hmax
1874                Ehkl[h,k,l] = E*phasep
1875                h,k,l = -hkl+Hmax       #Friedel pair refl.
1876                Ehkl[h,k,l] = E*phasem
1877#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1878    CEhkl = copy.copy(Ehkl)
1879    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1880    Emask = ma.getmask(MEhkl)
1881    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1882    Ncyc = 0
1883    old = np.seterr(all='raise')
1884    while True:       
1885        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1886        CEsig = np.std(CErho)
1887        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1888        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
1889        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1890        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1891        phase = CFhkl/np.absolute(CFhkl)
1892        CEhkl = np.absolute(Ehkl)*phase
1893        Ncyc += 1
1894        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1895        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1896        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1897        if Rcf < 5.:
1898            break
1899        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1900        if not GoOn or Ncyc > 10000:
1901            break
1902    np.seterr(**old)
1903    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1904    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1905    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1906    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
1907       
1908    mapData['Rcf'] = Rcf
1909    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1910    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1911    return mapData
1912   
1913def findSSOffset(SGData,SSGData,A,Fhklm):   
1914    '''default doc string
1915   
1916    :param type name: description
1917   
1918    :returns: type name: description
1919   
1920    '''
1921    if SGData['SpGrp'] == 'P 1':
1922        return [0,0,0,0]   
1923    hklmShape = Fhklm.shape
1924    hklmHalf = np.array(hklmShape)/2
1925    sortHKLM = np.argsort(Fhklm.flatten())
1926    Fdict = {}
1927    for hklm in sortHKLM:
1928        HKLM = np.unravel_index(hklm,hklmShape)
1929        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
1930        if F == 0.:
1931            break
1932        Fdict['%.6f'%(np.absolute(F))] = hklm
1933    Flist = np.flipud(np.sort(Fdict.keys()))
1934    F = str(1.e6)
1935    i = 0
1936    DH = []
1937    Dphi = []
1938    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1939    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1940    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1941    for F in Flist:
1942        hklm = np.unravel_index(Fdict[F],hklmShape)
1943        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
1944            continue
1945        Uniq = np.inner(hklm-hklmHalf,SSGMT)
1946        Phi = np.inner(hklm-hklmHalf,SSGT)
1947        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
1948        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1949        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
1950        ang0 = np.angle(Fh0,deg=True)/360.
1951        for H,phi in zip(Uniq,Phi)[1:]:
1952            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
1953            dH = H-hklm
1954            dang = ang-ang0
1955            DH.append(dH)
1956            Dphi.append((dang+.5) % 1.0)
1957        if i > 20 or len(DH) > 30:
1958            break
1959        i += 1
1960    DH = np.array(DH)
1961    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1962    Dphi = np.array(Dphi)
1963    steps = np.array(hklmShape)
1964    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]]
1965    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
1966    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
1967    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1968    hist,bins = np.histogram(Mmap,bins=1000)
1969#    for i,item in enumerate(hist[:10]):
1970#        print item,bins[i]
1971    chisq = np.min(Mmap)
1972    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1973    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
1974#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1975    return DX
1976   
1977def SSChargeFlip(data,reflDict,pgbar):
1978    '''default doc string
1979   
1980    :param type name: description
1981   
1982    :returns: type name: description
1983   
1984    '''
1985    generalData = data['General']
1986    mapData = generalData['Map']
1987    map4DData = {}
1988    flipData = generalData['Flip']
1989    FFtable = {}
1990    if 'None' not in flipData['Norm element']:
1991        normElem = flipData['Norm element'].upper()
1992        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1993        for ff in FFs:
1994            if ff['Symbol'] == normElem:
1995                FFtable.update(ff)
1996    dmin = flipData['Resolution']
1997    SGData = generalData['SGData']
1998    SSGData = generalData['SSGData']
1999    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2000    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2001    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2002    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2003    cell = generalData['Cell'][1:8]       
2004    A = G2lat.cell2A(cell[:6])
2005    Vol = cell[6]
2006    maxM = generalData['SuperVec'][2]
2007    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
2008    adjHKLmax(SGData,Hmax)
2009    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2010    time0 = time.time()
2011    for iref,ref in enumerate(reflDict['RefList']):
2012        dsp = ref[5]
2013        if dsp > dmin:
2014            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2015            if FFtable:
2016                SQ = 0.25/dsp**2
2017                ff *= G2el.ScatFac(FFtable,SQ)[0]
2018            if ref[9] > 0.:         #use only +ve Fobs**2
2019                E = np.sqrt(ref[9])/ff
2020            else:
2021                E = 0.
2022            ph = ref[11]
2023            ph = rn.uniform(0.,360.)
2024            Uniq = np.inner(ref[:4],SSGMT)
2025            Phi = np.inner(ref[:4],SSGT)
2026            for i,hklm in enumerate(Uniq):        #uses uniq
2027                hklm = np.asarray(hklm,dtype='i')
2028                dp = 360.*Phi[i]                #and phi
2029                a = cosd(ph+dp)
2030                b = sind(ph+dp)
2031                phasep = complex(a,b)
2032                phasem = complex(a,-b)
2033                h,k,l,m = hklm+Hmax
2034                Ehkl[h,k,l,m] = E*phasep
2035                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
2036                Ehkl[h,k,l,m] = E*phasem
2037#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2038    CEhkl = copy.copy(Ehkl)
2039    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2040    Emask = ma.getmask(MEhkl)
2041    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2042    Ncyc = 0
2043    old = np.seterr(all='raise')
2044    while True:       
2045        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2046        CEsig = np.std(CErho)
2047        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2048        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2049        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2050        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2051        phase = CFhkl/np.absolute(CFhkl)
2052        CEhkl = np.absolute(Ehkl)*phase
2053        Ncyc += 1
2054        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2055        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2056        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2057        if Rcf < 5.:
2058            break
2059        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2060        if not GoOn or Ncyc > 10000:
2061            break
2062    np.seterr(**old)
2063    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2064    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))
2065    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2066    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2067    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2068
2069    mapData['Rcf'] = Rcf
2070    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2071    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2072
2073    map4DData['Rcf'] = Rcf
2074    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))
2075    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
2076    return mapData,map4DData
2077   
2078def SearchMap(generalData,drawingData):
2079    '''Does a search of a density map for peaks meeting the criterion of peak
2080    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
2081    mapData is data['General']['mapData']; the map is also in mapData.
2082
2083    :param data: the phase data structure
2084
2085    :returns: (peaks,mags,dzeros) where
2086
2087        * peaks : ndarray
2088          x,y,z positions of the peaks found in the map
2089        * mags : ndarray
2090          the magnitudes of the peaks
2091        * dzeros : ndarray
2092          the distance of the peaks from  the unit cell origin
2093
2094    '''       
2095    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2096   
2097    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
2098   
2099#    def noDuplicate(xyz,peaks,Amat):
2100#        XYZ = np.inner(Amat,xyz)
2101#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2102#            print ' Peak',xyz,' <0.5A from another peak'
2103#            return False
2104#        return True
2105#                           
2106    def fixSpecialPos(xyz,SGData,Amat):
2107        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
2108        X = []
2109        xyzs = [equiv[0] for equiv in equivs]
2110        for x in xyzs:
2111            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
2112                X.append(x)
2113        if len(X) > 1:
2114            return np.average(X,axis=0)
2115        else:
2116            return xyz
2117       
2118    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
2119        Mag,x0,y0,z0,sig = parms
2120        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
2121#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
2122        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
2123       
2124    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
2125        Mag,x0,y0,z0,sig = parms
2126        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2127        return M
2128       
2129    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
2130        Mag,x0,y0,z0,sig = parms
2131        dMdv = np.zeros(([5,]+list(rX.shape)))
2132        delt = .01
2133        for i in range(5):
2134            parms[i] -= delt
2135            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2136            parms[i] += 2.*delt
2137            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2138            parms[i] -= delt
2139            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
2140        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2141        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
2142        dMdv = np.reshape(dMdv,(5,rX.size))
2143        Hess = np.inner(dMdv,dMdv)
2144       
2145        return Vec,Hess
2146       
2147    phaseName = generalData['Name']
2148    SGData = generalData['SGData']
2149    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2150    peaks = []
2151    mags = []
2152    dzeros = []
2153    try:
2154        mapData = generalData['Map']
2155        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
2156        rho = copy.copy(mapData['rho'])     #don't mess up original
2157        mapHalf = np.array(rho.shape)/2
2158        res = mapData['Resolution']
2159        incre = np.array(rho.shape,dtype=np.float)
2160        step = max(1.0,1./res)+1
2161        steps = np.array(3*[step,])
2162    except KeyError:
2163        print '**** ERROR - Fourier map not defined'
2164        return peaks,mags
2165    rhoMask = ma.array(rho,mask=(rho<contLevel))
2166    indices = (-1,0,1)
2167    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2168    for roll in rolls:
2169        if np.any(roll):
2170            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
2171    indx = np.transpose(rhoMask.nonzero())
2172    peaks = indx/incre
2173    mags = rhoMask[rhoMask.nonzero()]
2174    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
2175        rho = rollMap(rho,ind)
2176        rMM = mapHalf-steps
2177        rMP = mapHalf+steps+1
2178        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2179        peakInt = np.sum(rhoPeak)*res**3
2180        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2181        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
2182        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
2183            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
2184        x1 = result[0]
2185        if not np.any(x1 < 0):
2186            mag = x1[0]
2187            peak = (np.array(x1[1:4])-ind)/incre
2188        peak = fixSpecialPos(peak,SGData,Amat)
2189        rho = rollMap(rho,-ind)       
2190    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
2191    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
2192   
2193def sortArray(data,pos,reverse=False):
2194    '''data is a list of items
2195    sort by pos in list; reverse if True
2196    '''
2197    T = []
2198    for i,M in enumerate(data):
2199        T.append((M[pos],i))
2200    D = dict(zip(T,data))
2201    T.sort()
2202    if reverse:
2203        T.reverse()
2204    X = []
2205    for key in T:
2206        X.append(D[key])
2207    return X
2208
2209def PeaksEquiv(data,Ind):
2210    '''Find the equivalent map peaks for those selected. Works on the
2211    contents of data['Map Peaks'].
2212
2213    :param data: the phase data structure
2214    :param list Ind: list of selected peak indices
2215    :returns: augmented list of peaks including those related by symmetry to the
2216      ones in Ind
2217
2218    '''       
2219    def Duplicate(xyz,peaks,Amat):
2220        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2221            return True
2222        return False
2223                           
2224    generalData = data['General']
2225    cell = generalData['Cell'][1:7]
2226    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2227    A = G2lat.cell2A(cell)
2228    SGData = generalData['SGData']
2229    mapPeaks = data['Map Peaks']
2230    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
2231    Indx = {}
2232    for ind in Ind:
2233        xyz = np.array(mapPeaks[ind][1:4])
2234        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
2235        for jnd,xyz in enumerate(XYZ):       
2236            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
2237    Ind = []
2238    for ind in Indx:
2239        if Indx[ind]:
2240            Ind.append(ind)
2241    return Ind
2242               
2243def PeaksUnique(data,Ind):
2244    '''Finds the symmetry unique set of peaks from those selected. Works on the
2245    contents of data['Map Peaks'].
2246
2247    :param data: the phase data structure
2248    :param list Ind: list of selected peak indices
2249    :returns: the list of symmetry unique peaks from among those given in Ind
2250
2251    '''       
2252#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
2253
2254    def noDuplicate(xyz,peaks,Amat):
2255        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=1.0) for peak in peaks]:
2256            return False
2257        return True
2258                           
2259    generalData = data['General']
2260    cell = generalData['Cell'][1:7]
2261    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2262    A = G2lat.cell2A(cell)
2263    SGData = generalData['SGData']
2264    mapPeaks = data['Map Peaks']
2265    Indx = {}
2266    XYZ = {}
2267    for ind in Ind:
2268        XYZ[ind] = np.array(mapPeaks[ind][1:4])
2269        Indx[ind] = True
2270    for ind in Ind:
2271        if Indx[ind]:
2272            xyz = XYZ[ind]
2273            for jnd in Ind:
2274                if ind != jnd and Indx[jnd]:                       
2275                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
2276                    xyzs = np.array([equiv[0] for equiv in Equiv])
2277                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
2278    Ind = []
2279    for ind in Indx:
2280        if Indx[ind]:
2281            Ind.append(ind)
2282    return Ind
2283   
2284################################################################################
2285##### single peak fitting profile fxn stuff
2286################################################################################
2287
2288def getCWsig(ins,pos):
2289    '''get CW peak profile sigma
2290   
2291    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
2292      as values only
2293    :param float pos: 2-theta of peak
2294    :returns: float getCWsig: peak sigma
2295   
2296    '''
2297    tp = tand(pos/2.0)
2298    return ins['U']*tp**2+ins['V']*tp+ins['W']
2299   
2300def getCWsigDeriv(pos):
2301    '''get derivatives of CW peak profile sigma wrt U,V, & W
2302   
2303    :param float pos: 2-theta of peak
2304   
2305    :returns: list getCWsigDeriv: d(sig)/dU, d(sig)/dV & d(sig)/dW
2306   
2307    '''
2308    tp = tand(pos/2.0)
2309    return tp**2,tp,1.0
2310   
2311def getCWgam(ins,pos):
2312    '''get CW peak profile gamma
2313   
2314    :param dict ins: instrument parameters with at least 'X' & 'Y'
2315      as values only
2316    :param float pos: 2-theta of peak
2317    :returns: float getCWgam: peak gamma
2318   
2319    '''
2320    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2321   
2322def getCWgamDeriv(pos):
2323    '''get derivatives of CW peak profile gamma wrt X & Y
2324   
2325    :param float pos: 2-theta of peak
2326   
2327    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
2328   
2329    '''
2330    return 1./cosd(pos/2.0),tand(pos/2.0)
2331   
2332def getTOFsig(ins,dsp):
2333    '''get TOF peak profile sigma
2334   
2335    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
2336      as values only
2337    :param float dsp: d-spacing of peak
2338   
2339    :returns: float getTOFsig: peak sigma
2340   
2341    '''
2342    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
2343   
2344def getTOFsigDeriv(dsp):
2345    '''get derivatives of TOF peak profile gamma wrt sig-0, sig-1, & sig-q
2346   
2347    :param float dsp: d-spacing of peak
2348   
2349    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
2350   
2351    '''
2352    return 1.0,dsp**2,dsp**4,1./dsp**2
2353   
2354def getTOFgamma(ins,dsp):
2355    '''get TOF peak profile gamma
2356   
2357    :param dict ins: instrument parameters with at least 'X' & 'Y'
2358      as values only
2359    :param float dsp: d-spacing of peak
2360   
2361    :returns: float getTOFgamma: peak gamma
2362   
2363    '''
2364    return ins['X']*dsp+ins['Y']*dsp**2
2365   
2366def getTOFgammaDeriv(dsp):
2367    '''get derivatives of TOF peak profile gamma wrt X & Y
2368   
2369    :param float dsp: d-spacing of peak
2370   
2371    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
2372   
2373    '''
2374    return dsp,dsp**2
2375   
2376def getTOFbeta(ins,dsp):
2377    '''get TOF peak profile beta
2378   
2379    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
2380      as values only
2381    :param float dsp: d-spacing of peak
2382   
2383    :returns: float getTOFbeta: peak beat
2384   
2385    '''
2386    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
2387   
2388def getTOFbetaDeriv(dsp):
2389    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
2390   
2391    :param float dsp: d-spacing of peak
2392   
2393    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
2394   
2395    '''
2396    return 1.0,1./dsp**4,1./dsp**2
2397   
2398def getTOFalpha(ins,dsp):
2399    '''get TOF peak profile alpha
2400   
2401    :param dict ins: instrument parameters with at least 'alpha'
2402      as values only
2403    :param float dsp: d-spacing of peak
2404   
2405    :returns: flaot getTOFalpha: peak alpha
2406   
2407    '''
2408    return ins['alpha']/dsp
2409   
2410def getTOFalphaDeriv(dsp):
2411    '''get derivatives of TOF peak profile beta wrt alpha
2412   
2413    :param float dsp: d-spacing of peak
2414   
2415    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
2416   
2417    '''
2418    return 1./dsp
2419   
2420def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2421    '''set starting peak parameters for single peak fits from plot selection or auto selection
2422   
2423    :param dict Parms: instrument parameters dictionary
2424    :param dict Parms2: table lookup for TOF profile coefficients
2425    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
2426    :param float mag: peak top magnitude from pick
2427    :param bool ifQ: True if pos in Q
2428    :param bool useFit: True if use fitted CW Parms values (not defaults)
2429   
2430    :returns: list XY: peak list entry:
2431        for CW: [pos,0,mag,1,sig,0,gam,0]
2432        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2433        NB: mag refinement set by default, all others off
2434   
2435    '''
2436    ind = 0
2437    if useFit:
2438        ind = 1
2439    ins = {}
2440    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
2441        for x in ['U','V','W','X','Y']:
2442            ins[x] = Parms[x][ind]
2443        if ifQ:                              #qplot - convert back to 2-theta
2444            pos = 2.0*asind(pos*wave/(4*math.pi))
2445        sig = getCWsig(ins,pos)
2446        gam = getCWgam(ins,pos)           
2447        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
2448    else:
2449        if ifQ:
2450            dsp = 2.*np.pi/pos
2451            pos = Parms['difC']*dsp
2452        else:
2453            dsp = pos/Parms['difC'][1]
2454        if 'Pdabc' in Parms2:
2455            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
2456                ins[x] = Parms[x][ind]
2457            Pdabc = Parms2['Pdabc'].T
2458            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2459            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2460        else:
2461            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
2462                ins[x] = Parms[x][ind]
2463            alp = getTOFalpha(ins,dsp)
2464            bet = getTOFbeta(ins,dsp)
2465        sig = getTOFsig(ins,dsp)
2466        gam = getTOFgamma(ins,dsp)
2467        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2468    return XY
2469   
2470################################################################################
2471##### MC/SA stuff
2472################################################################################
2473
2474#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2475# Original Author: Travis Oliphant 2002
2476# Bug-fixes in 2006 by Tim Leslie
2477
2478
2479import numpy
2480from numpy import asarray, tan, exp, ones, squeeze, sign, \
2481     all, log, sqrt, pi, shape, array, minimum, where
2482from numpy import random
2483
2484#__all__ = ['anneal']
2485
2486_double_min = numpy.finfo(float).min
2487_double_max = numpy.finfo(float).max
2488class base_schedule(object):
2489    def __init__(self):
2490        self.dwell = 20
2491        self.learn_rate = 0.5
2492        self.lower = -10
2493        self.upper = 10
2494        self.Ninit = 50
2495        self.accepted = 0
2496        self.tests = 0
2497        self.feval = 0
2498        self.k = 0
2499        self.T = None
2500
2501    def init(self, **options):
2502        self.__dict__.update(options)
2503        self.lower = asarray(self.lower)
2504        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2505        self.upper = asarray(self.upper)
2506        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2507        self.k = 0
2508        self.accepted = 0
2509        self.feval = 0
2510        self.tests = 0
2511
2512    def getstart_temp(self, best_state):
2513        """ Find a matching starting temperature and starting parameters vector
2514        i.e. find x0 such that func(x0) = T0.
2515
2516        Parameters
2517        ----------
2518        best_state : _state
2519            A _state object to store the function value and x0 found.
2520
2521        returns
2522        -------
2523        x0 : array
2524            The starting parameters vector.
2525        """
2526
2527        assert(not self.dims is None)
2528        lrange = self.lower
2529        urange = self.upper
2530        fmax = _double_min
2531        fmin = _double_max
2532        for _ in range(self.Ninit):
2533            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2534            fval = self.func(x0, *self.args)
2535            self.feval += 1
2536            if fval > fmax:
2537                fmax = fval
2538            if fval < fmin:
2539                fmin = fval
2540                best_state.cost = fval
2541                best_state.x = array(x0)
2542
2543        self.T0 = (fmax-fmin)*1.5
2544        return best_state.x
2545       
2546    def set_range(self,x0,frac):
2547        delrange = frac*(self.upper-self.lower)
2548        self.upper = x0+delrange
2549        self.lower = x0-delrange
2550
2551    def accept_test(self, dE):
2552        T = self.T
2553        self.tests += 1
2554        if dE < 0:
2555            self.accepted += 1
2556            return 1
2557        p = exp(-dE*1.0/self.boltzmann/T)
2558        if (p > random.uniform(0.0, 1.0)):
2559            self.accepted += 1
2560            return 1
2561        return 0
2562
2563    def update_guess(self, x0):
2564        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2565
2566    def update_temp(self, x0):
2567        pass
2568
2569
2570#  A schedule due to Lester Ingber modified to use bounds - OK
2571class fast_sa(base_schedule):
2572    def init(self, **options):
2573        self.__dict__.update(options)
2574        if self.m is None:
2575            self.m = 1.0
2576        if self.n is None:
2577            self.n = 1.0
2578        self.c = self.m * exp(-self.n * self.quench)
2579
2580    def update_guess(self, x0):
2581        x0 = asarray(x0)
2582        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2583        T = self.T
2584        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2585        xnew = xc*(self.upper - self.lower)+self.lower
2586        return xnew
2587#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2588#        xc = y*(self.upper - self.lower)
2589#        xnew = x0 + xc
2590#        return xnew
2591
2592    def update_temp(self):
2593        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2594        self.k += 1
2595        return
2596
2597class cauchy_sa(base_schedule):     #modified to use bounds - not good
2598    def update_guess(self, x0):
2599        x0 = asarray(x0)
2600        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2601        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2602        xnew = xc*(self.upper - self.lower)+self.lower
2603        return xnew
2604#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2605#        xc = self.learn_rate * self.T * tan(numbers)
2606#        xnew = x0 + xc
2607#        return xnew
2608
2609    def update_temp(self):
2610        self.T = self.T0/(1+self.k)
2611        self.k += 1
2612        return
2613
2614class boltzmann_sa(base_schedule):
2615#    def update_guess(self, x0):
2616#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2617#        x0 = asarray(x0)
2618#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2619#
2620#        xnew = x0 + xc*std*self.learn_rate
2621#        return xnew
2622
2623    def update_temp(self):
2624        self.k += 1
2625        self.T = self.T0 / log(self.k+1.0)
2626        return
2627
2628class log_sa(base_schedule):        #OK
2629
2630    def init(self,**options):
2631        self.__dict__.update(options)
2632       
2633    def update_guess(self,x0):     #same as default
2634        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2635       
2636    def update_temp(self):
2637        self.k += 1
2638        self.T = self.T0*self.slope**self.k
2639       
2640class _state(object):
2641    def __init__(self):
2642        self.x = None
2643        self.cost = None
2644
2645# TODO:
2646#     allow for general annealing temperature profile
2647#     in that case use update given by alpha and omega and
2648#     variation of all previous updates and temperature?
2649
2650# Simulated annealing
2651
2652def anneal(func, x0, args=(), schedule='fast', full_output=0,
2653           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2654           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2655           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2656           ranRange=0.10,autoRan=False,dlg=None):
2657    """Minimize a function using simulated annealing.
2658
2659    Schedule is a schedule class implementing the annealing schedule.
2660    Available ones are 'fast', 'cauchy', 'boltzmann'
2661
2662    :param callable func: f(x, \*args)
2663        Function to be optimized.
2664    :param ndarray x0:
2665        Initial guess.
2666    :param tuple args:
2667        Extra parameters to `func`.
2668    :param base_schedule schedule:
2669        Annealing schedule to use (a class).
2670    :param bool full_output:
2671        Whether to return optional outputs.
2672    :param float T0:
2673        Initial Temperature (estimated as 1.2 times the largest
2674        cost-function deviation over random points in the range).
2675    :param float Tf:
2676        Final goal temperature.
2677    :param int maxeval:
2678        Maximum function evaluations.
2679    :param int maxaccept:
2680        Maximum changes to accept.
2681    :param int maxiter:
2682        Maximum cooling iterations.
2683    :param float learn_rate:
2684        Scale constant for adjusting guesses.
2685    :param float boltzmann:
2686        Boltzmann constant in acceptance test
2687        (increase for less stringent test at each temperature).
2688    :param float feps:
2689        Stopping relative error tolerance for the function value in
2690        last four coolings.
2691    :param float quench,m,n:
2692        Parameters to alter fast_sa schedule.
2693    :param float/ndarray lower,upper:
2694        Lower and upper bounds on `x`.
2695    :param int dwell:
2696        The number of times to search the space at each temperature.
2697    :param float slope:
2698        Parameter for log schedule
2699    :param bool ranStart=False:
2700        True for set 10% of ranges about x
2701
2702    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2703
2704     * xmin (ndarray): Point giving smallest value found.
2705     * Jmin (float): Minimum value of function found.
2706     * T (float): Final temperature.
2707     * feval (int): Number of function evaluations.
2708     * iters (int): Number of cooling iterations.
2709     * accept (int): Number of tests accepted.
2710     * retval (int): Flag indicating stopping condition:
2711
2712              *  0: Points no longer changing
2713              *  1: Cooled to final temperature
2714              *  2: Maximum function evaluations
2715              *  3: Maximum cooling iterations reached
2716              *  4: Maximum accepted query locations reached
2717              *  5: Final point not the minimum amongst encountered points
2718
2719    *Notes*:
2720    Simulated annealing is a random algorithm which uses no derivative
2721    information from the function being optimized. In practice it has
2722    been more useful in discrete optimization than continuous
2723    optimization, as there are usually better algorithms for continuous
2724    optimization problems.
2725
2726    Some experimentation by trying the difference temperature
2727    schedules and altering their parameters is likely required to
2728    obtain good performance.
2729
2730    The randomness in the algorithm comes from random sampling in numpy.
2731    To obtain the same results you can call numpy.random.seed with the
2732    same seed immediately before calling scipy.optimize.anneal.
2733
2734    We give a brief description of how the three temperature schedules
2735    generate new points and vary their temperature. Temperatures are
2736    only updated with iterations in the outer loop. The inner loop is
2737    over xrange(dwell), and new points are generated for
2738    every iteration in the inner loop. (Though whether the proposed
2739    new points are accepted is probabilistic.)
2740
2741    For readability, let d denote the dimension of the inputs to func.
2742    Also, let x_old denote the previous state, and k denote the
2743    iteration number of the outer loop. All other variables not
2744    defined below are input variables to scipy.optimize.anneal itself.
2745
2746    In the 'fast' schedule the updates are ::
2747
2748        u ~ Uniform(0, 1, size=d)
2749        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2750        xc = y * (upper - lower)
2751        x_new = x_old + xc
2752
2753        c = n * exp(-n * quench)
2754        T_new = T0 * exp(-c * k**quench)
2755
2756
2757    In the 'cauchy' schedule the updates are ::
2758
2759        u ~ Uniform(-pi/2, pi/2, size=d)
2760        xc = learn_rate * T * tan(u)
2761        x_new = x_old + xc
2762
2763        T_new = T0 / (1+k)
2764
2765    In the 'boltzmann' schedule the updates are ::
2766
2767        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2768        y ~ Normal(0, std, size=d)
2769        x_new = x_old + learn_rate * y
2770
2771        T_new = T0 / log(1+k)
2772
2773    """
2774    x0 = asarray(x0)
2775    lower = asarray(lower)
2776    upper = asarray(upper)
2777
2778    schedule = eval(schedule+'_sa()')
2779    #   initialize the schedule
2780    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2781                  learn_rate=learn_rate, lower=lower, upper=upper,
2782                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2783
2784    current_state, last_state, best_state = _state(), _state(), _state()
2785    if ranStart:
2786        schedule.set_range(x0,ranRange)
2787    if T0 is None:
2788        x0 = schedule.getstart_temp(best_state)
2789    else:
2790        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
2791        best_state.x = None
2792        best_state.cost = numpy.Inf
2793
2794    last_state.x = asarray(x0).copy()
2795    fval = func(x0,*args)
2796    schedule.feval += 1
2797    last_state.cost = fval
2798    if last_state.cost < best_state.cost:
2799        best_state.cost = fval
2800        best_state.x = asarray(x0).copy()
2801    schedule.T = schedule.T0
2802    fqueue = [100, 300, 500, 700]
2803    iters = 1
2804    keepGoing = True
2805    bestn = 0
2806    while keepGoing:
2807        retval = 0
2808        for n in xrange(dwell):
2809            current_state.x = schedule.update_guess(last_state.x)
2810            current_state.cost = func(current_state.x,*args)
2811            schedule.feval += 1
2812
2813            dE = current_state.cost - last_state.cost
2814            if schedule.accept_test(dE):
2815                last_state.x = current_state.x.copy()
2816                last_state.cost = current_state.cost
2817                if last_state.cost < best_state.cost:
2818                    best_state.x = last_state.x.copy()
2819                    best_state.cost = last_state.cost
2820                    bestn = n
2821                    if best_state.cost < 1.0 and autoRan:
2822                        schedule.set_range(x0,best_state.cost/2.)                       
2823        if dlg:
2824            GoOn = dlg.Update(min(100.,best_state.cost*100),
2825                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2826                    'Best trial:',bestn,  \
2827                    'MC/SA Residual:',best_state.cost*100,'%', \
2828                    ))[0]
2829            if not GoOn:
2830                best_state.x = last_state.x.copy()
2831                best_state.cost = last_state.cost
2832                retval = 5
2833        schedule.update_temp()
2834        iters += 1
2835        # Stopping conditions
2836        # 0) last saved values of f from each cooling step
2837        #     are all very similar (effectively cooled)
2838        # 1) Tf is set and we are below it
2839        # 2) maxeval is set and we are past it
2840        # 3) maxiter is set and we are past it
2841        # 4) maxaccept is set and we are past it
2842        # 5) user canceled run via progress bar
2843
2844        fqueue.append(squeeze(last_state.cost))
2845        fqueue.pop(0)
2846        af = asarray(fqueue)*1.0
2847        if retval == 5:
2848            print ' User terminated run; incomplete MC/SA'
2849            keepGoing = False
2850            break
2851        if all(abs((af-af[0])/af[0]) < feps):
2852            retval = 0
2853            if abs(af[-1]-best_state.cost) > feps*10:
2854                retval = 5
2855#                print "Warning: Cooled to %f at %s but this is not" \
2856#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2857#                      + " the smallest point found."
2858            break
2859        if (Tf is not None) and (schedule.T < Tf):
2860            retval = 1
2861            break
2862        if (maxeval is not None) and (schedule.feval > maxeval):
2863            retval = 2
2864            break
2865        if (iters > maxiter):
2866            print "Warning: Maximum number of iterations exceeded."
2867            retval = 3
2868            break
2869        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2870            retval = 4
2871            break
2872
2873    if full_output:
2874        return best_state.x, best_state.cost, schedule.T, \
2875               schedule.feval, iters, schedule.accepted, retval
2876    else:
2877        return best_state.x, retval
2878
2879def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
2880    outlist = []
2881    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
2882    for n in range(iCyc):
2883        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2884        outlist.append(result[0])
2885        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2886    out_q.put(outlist)
2887
2888def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2889    import multiprocessing as mp
2890   
2891    nprocs = mp.cpu_count()
2892    out_q = mp.Queue()
2893    procs = []
2894    iCyc = np.zeros(nprocs)
2895    for i in range(nCyc):
2896        iCyc[i%nprocs] += 1
2897    for i in range(nprocs):
2898        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
2899        procs.append(p)
2900        p.start()
2901    resultlist = []
2902    for i in range(nprocs):
2903        resultlist += out_q.get()
2904    for p in procs:
2905        p.join()
2906    return resultlist
2907
2908def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2909    '''default doc string
2910   
2911    :param type name: description
2912   
2913    :returns: type name: description
2914   
2915    '''
2916   
2917    global tsum
2918    tsum = 0.
2919   
2920    def getMDparms(item,pfx,parmDict,varyList):
2921        parmDict[pfx+'MDaxis'] = item['axis']
2922        parmDict[pfx+'MDval'] = item['Coef'][0]
2923        if item['Coef'][1]:
2924            varyList += [pfx+'MDval',]
2925            limits = item['Coef'][2]
2926            lower.append(limits[0])
2927            upper.append(limits[1])
2928                       
2929    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2930        parmDict[pfx+'Atype'] = item['atType']
2931        aTypes |= set([item['atType'],]) 
2932        pstr = ['Ax','Ay','Az']
2933        XYZ = [0,0,0]
2934        for i in range(3):
2935            name = pfx+pstr[i]
2936            parmDict[name] = item['Pos'][0][i]
2937            XYZ[i] = parmDict[name]
2938            if item['Pos'][1][i]:
2939                varyList += [name,]
2940                limits = item['Pos'][2][i]
2941                lower.append(limits[0])
2942                upper.append(limits[1])
2943        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2944           
2945    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2946        parmDict[mfx+'MolCent'] = item['MolCent']
2947        parmDict[mfx+'RBId'] = item['RBId']
2948        pstr = ['Px','Py','Pz']
2949        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
2950        for i in range(3):
2951            name = mfx+pstr[i]
2952            parmDict[name] = item['Pos'][0][i]
2953            if item['Pos'][1][i]:
2954                varyList += [name,]
2955                limits = item['Pos'][2][i]
2956                lower.append(limits[0])
2957                upper.append(limits[1])
2958        AV = item['Ori'][0]
2959        A = AV[0]
2960        V = AV[1:]
2961        for i in range(4):
2962            name = mfx+ostr[i]
2963            if i:
2964                parmDict[name] = V[i-1]
2965            else:
2966                parmDict[name] = A
2967            if item['Ovar'] == 'AV':
2968                varyList += [name,]
2969                limits = item['Ori'][2][i]
2970                lower.append(limits[0])
2971                upper.append(limits[1])
2972            elif item['Ovar'] == 'A' and not i:
2973                varyList += [name,]
2974                limits = item['Ori'][2][i]
2975                lower.append(limits[0])
2976                upper.append(limits[1])
2977        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2978            for i in range(len(item['Tor'][0])):
2979                name = mfx+'Tor'+str(i)
2980                parmDict[name] = item['Tor'][0][i]
2981                if item['Tor'][1][i]:
2982                    varyList += [name,]
2983                    limits = item['Tor'][2][i]
2984                    lower.append(limits[0])
2985                    upper.append(limits[1])
2986        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2987        aTypes |= set(atypes)
2988        atNo += len(atypes)
2989        return atNo
2990               
2991    def GetAtomM(Xdata,SGData):
2992        Mdata = []
2993        for xyz in Xdata:
2994            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2995        return np.array(Mdata)
2996       
2997    def GetAtomT(RBdata,parmDict):
2998        'Needs a doc string'
2999        atNo = parmDict['atNo']
3000        nfixAt = parmDict['nfixAt']
3001        Tdata = atNo*[' ',]
3002        for iatm in range(nfixAt):
3003            parm = ':'+str(iatm)+':Atype'
3004            if parm in parmDict:
3005                Tdata[iatm] = aTypes.index(parmDict[parm])
3006        iatm = nfixAt
3007        for iObj in range(parmDict['nObj']):
3008            pfx = str(iObj)+':'
3009            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3010                if parmDict[pfx+'Type'] == 'Vector':
3011                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3012                    nAtm = len(RBRes['rbVect'][0])
3013                else:       #Residue
3014                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3015                    nAtm = len(RBRes['rbXYZ'])
3016                for i in range(nAtm):
3017                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
3018                    iatm += 1
3019            elif parmDict[pfx+'Type'] == 'Atom':
3020                atNo = parmDict[pfx+'atNo']
3021                parm = pfx+'Atype'              #remove extra ':'
3022                if parm in parmDict:
3023                    Tdata[atNo] = aTypes.index(parmDict[parm])
3024                iatm += 1
3025            else:
3026                continue        #skips March Dollase
3027        return Tdata
3028       
3029    def GetAtomX(RBdata,parmDict):
3030        'Needs a doc string'
3031        Bmat = parmDict['Bmat']
3032        atNo = parmDict['atNo']
3033        nfixAt = parmDict['nfixAt']
3034        Xdata = np.zeros((3,atNo))
3035        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
3036        for iatm in range(nfixAt):
3037            for key in keys:
3038                parm = ':'+str(iatm)+key
3039                if parm in parmDict:
3040                    keys[key][iatm] = parmDict[parm]
3041        iatm = nfixAt
3042        for iObj in range(parmDict['nObj']):
3043            pfx = str(iObj)+':'
3044            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3045                if parmDict[pfx+'Type'] == 'Vector':
3046                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3047                    vecs = RBRes['rbVect']
3048                    mags = RBRes['VectMag']
3049                    Cart = np.zeros_like(vecs[0])
3050                    for vec,mag in zip(vecs,mags):
3051                        Cart += vec*mag
3052                elif parmDict[pfx+'Type'] == 'Residue':
3053                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3054                    Cart = np.array(RBRes['rbXYZ'])
3055                    for itor,seq in enumerate(RBRes['rbSeq']):
3056                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
3057                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
3058                if parmDict[pfx+'MolCent'][1]:
3059                    Cart -= parmDict[pfx+'MolCent'][0]
3060                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3061                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
3062                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
3063                iatm += len(Cart)
3064            elif parmDict[pfx+'Type'] == 'Atom':
3065                atNo = parmDict[pfx+'atNo']
3066                for key in keys:
3067                    parm = pfx+key[1:]              #remove extra ':'
3068                    if parm in parmDict:
3069                        keys[key][atNo] = parmDict[parm]
3070                iatm += 1
3071            else:
3072                continue        #skips March Dollase
3073        return Xdata.T
3074       
3075    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
3076        allX = np.inner(Xdata,SGM)+SGT
3077        allT = np.repeat(Tdata,allX.shape[1])
3078        allM = np.repeat(Mdata,allX.shape[1])
3079        allX = np.reshape(allX,(-1,3))
3080        return allT,allM,allX
3081
3082    def getAllX(Xdata,SGM,SGT):
3083        allX = np.inner(Xdata,SGM)+SGT
3084        allX = np.reshape(allX,(-1,3))
3085        return allX
3086       
3087    def normQuaternions(RBdata,parmDict,varyList,values):
3088        for iObj in range(parmDict['nObj']):
3089            pfx = str(iObj)+':'
3090            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3091                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3092                A,V = Q2AVdeg(Qori)
3093                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
3094                    if i:
3095                        parmDict[pfx+name] = V[i-1]
3096                    else:
3097                        parmDict[pfx+name] = A
3098       
3099    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
3100        ''' Compute structure factors for all h,k,l for phase
3101        input:
3102            refList: [ref] where each ref = h,k,l,m,d,...
3103            rcov:   array[nref,nref] covariance terms between Fo^2 values
3104            ifInv:  bool True if centrosymmetric
3105            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
3106            RBdata: [dict] rigid body dictionary
3107            varyList: [list] names of varied parameters in MC/SA (not used here)           
3108            ParmDict: [dict] problem parameters
3109        puts result F^2 in each ref[5] in refList
3110        returns:
3111            delt-F*rcov*delt-F/sum(Fo^2)^2
3112        '''       
3113        global tsum
3114        t0 = time.time()
3115        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
3116        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
3117        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
3118        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
3119        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
3120        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
3121        refList[5] = Aterm
3122        if not ifInv:
3123            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
3124        if len(cosTable):        #apply MD correction
3125            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
3126        sumFcsq = np.sum(refList[5])
3127        scale = parmDict['sumFosq']/sumFcsq
3128        refList[5] *= scale
3129        refList[6] = refList[4]-refList[5]
3130        M = np.inner(refList[6],np.inner(rcov,refList[6]))
3131        tsum += (time.time()-t0)
3132        return M/np.sum(refList[4]**2)
3133
3134    sq8ln2 = np.sqrt(8*np.log(2))
3135    sq2pi = np.sqrt(2*np.pi)
3136    sq4pi = np.sqrt(4*np.pi)
3137    generalData = data['General']
3138    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3139    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
3140    SGData = generalData['SGData']
3141    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
3142    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
3143    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
3144    fixAtoms = data['Atoms']                       #if any
3145    cx,ct,cs = generalData['AtomPtrs'][:3]
3146    aTypes = set([])
3147    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
3148    varyList = []
3149    atNo = 0
3150    for atm in fixAtoms:
3151        pfx = ':'+str(atNo)+':'
3152        parmDict[pfx+'Atype'] = atm[ct]
3153        aTypes |= set([atm[ct],])
3154        pstr = ['Ax','Ay','Az']
3155        parmDict[pfx+'Amul'] = atm[cs+1]
3156        for i in range(3):
3157            name = pfx+pstr[i]
3158            parmDict[name] = atm[cx+i]
3159        atNo += 1
3160    parmDict['nfixAt'] = len(fixAtoms)       
3161    MCSA = generalData['MCSA controls']
3162    reflName = MCSA['Data source']
3163    phaseName = generalData['Name']
3164    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
3165    upper = []
3166    lower = []
3167    MDvec = np.zeros(3)
3168    for i,item in enumerate(MCSAObjs):
3169        mfx = str(i)+':'
3170        parmDict[mfx+'Type'] = item['Type']
3171        if item['Type'] == 'MD':
3172            getMDparms(item,mfx,parmDict,varyList)
3173            MDvec = np.array(item['axis'])
3174        elif item['Type'] == 'Atom':
3175            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
3176            parmDict[mfx+'atNo'] = atNo
3177            atNo += 1
3178        elif item['Type'] in ['Residue','Vector']:
3179            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
3180    parmDict['atNo'] = atNo                 #total no. of atoms
3181    parmDict['nObj'] = len(MCSAObjs)
3182    aTypes = list(aTypes)
3183    Tdata = GetAtomT(RBdata,parmDict)
3184    Xdata = GetAtomX(RBdata,parmDict)
3185    Mdata = GetAtomM(Xdata,SGData)
3186    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
3187    FFtables = G2el.GetFFtable(aTypes)
3188    refs = []
3189    allFF = []
3190    cosTable = []
3191    sumFosq = 0
3192    if 'PWDR' in reflName:
3193        for ref in reflData:
3194            h,k,l,m,d,pos,sig,gam,f = ref[:9]
3195            if d >= MCSA['dmin']:
3196                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
3197                SQ = 0.25/d**2
3198                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3199                refs.append([h,k,l,m,f*m,pos,sig])
3200                sumFosq += f*m
3201                Heqv = np.inner(np.array([h,k,l]),SGMT)
3202                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3203        nRef = len(refs)
3204        cosTable = np.array(cosTable)**2
3205        rcov = np.zeros((nRef,nRef))
3206        for iref,refI in enumerate(refs):
3207            rcov[iref][iref] = 1./(sq4pi*refI[6])
3208            for jref,refJ in enumerate(refs[:iref]):
3209                t1 = refI[6]**2+refJ[6]**2
3210                t2 = (refJ[5]-refI[5])**2/(2.*t1)
3211                if t2 > 10.:
3212                    rcov[iref][jref] = 0.
3213                else:
3214                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
3215        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3216        Rdiag = np.sqrt(np.diag(rcov))
3217        Rnorm = np.outer(Rdiag,Rdiag)
3218        rcov /= Rnorm
3219    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
3220        vNames = []
3221        pfx = str(data['pId'])+'::PWLref:'
3222        for iref,refI in enumerate(reflData):           #Pawley reflection set
3223            h,k,l,m,d,v,f,s = refI
3224            if d >= MCSA['dmin'] and v:       #skip unrefined ones
3225                vNames.append(pfx+str(iref))
3226                SQ = 0.25/d**2
3227                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3228                refs.append([h,k,l,m,f*m,iref,0.])
3229                sumFosq += f*m
3230                Heqv = np.inner(np.array([h,k,l]),SGMT)
3231                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3232        cosTable = np.array(cosTable)**2
3233        nRef = len(refs)
3234#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
3235        if covData['freshCOV'] or  MCSA['newDmin']:
3236            vList = covData['varyList']
3237            covMatrix = covData['covMatrix']
3238            rcov = getVCov(vNames,vList,covMatrix)
3239            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3240            Rdiag = np.sqrt(np.diag(rcov))
3241            Rdiag = np.where(Rdiag,Rdiag,1.0)
3242            Rnorm = np.outer(Rdiag,Rdiag)
3243            rcov /= Rnorm
3244            MCSA['rcov'] = rcov
3245            covData['freshCOV'] = False
3246            MCSA['newDmin'] = False
3247        else:
3248            rcov = MCSA['rcov']
3249    elif 'HKLF' in reflName:
3250        for ref in reflData:
3251            [h,k,l,m,d],f = ref[:5],ref[6]
3252            if d >= MCSA['dmin']:
3253                SQ = 0.25/d**2
3254                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3255                refs.append([h,k,l,m,f*m,0.,0.])
3256                sumFosq += f*m
3257        nRef = len(refs)
3258        rcov = np.identity(len(refs))
3259    allFF = np.array(allFF).T
3260    refs = np.array(refs).T
3261    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
3262    print ' Number of parameters varied: %d'%(len(varyList))
3263    parmDict['sumFosq'] = sumFosq
3264    x0 = [parmDict[val] for val in varyList]
3265    ifInv = SGData['SGInv']
3266    # consider replacing anneal with scipy.optimize.basinhopping
3267    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
3268        schedule=MCSA['Algorithm'], full_output=True,
3269        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
3270        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
3271        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
3272        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
3273        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
3274    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
3275#    for ref in refs.T:
3276#        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])
3277#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
3278    Result = [False,False,results[1],results[2],]+list(results[0])
3279    Result.append(varyList)
3280    return Result,tsum
3281
3282       
3283################################################################################
3284##### Quaternion stuff
3285################################################################################
3286
3287def prodQQ(QA,QB):
3288    ''' Grassman quaternion product
3289        QA,QB quaternions; q=r+ai+bj+ck
3290    '''
3291    D = np.zeros(4)
3292    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
3293    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
3294    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
3295    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
3296   
3297#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
3298#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
3299   
3300    return D
3301   
3302def normQ(QA):
3303    ''' get length of quaternion & normalize it
3304        q=r+ai+bj+ck
3305    '''
3306    n = np.sqrt(np.sum(np.array(QA)**2))
3307    return QA/n
3308   
3309def invQ(Q):
3310    '''
3311        get inverse of quaternion
3312        q=r+ai+bj+ck; q* = r-ai-bj-ck
3313    '''
3314    return Q*np.array([1,-1,-1,-1])
3315   
3316def prodQVQ(Q,V):
3317    """
3318    compute the quaternion vector rotation qvq-1 = v'
3319    q=r+ai+bj+ck
3320    """
3321    T2 = Q[0]*Q[1]
3322    T3 = Q[0]*Q[2]
3323    T4 = Q[0]*Q[3]
3324    T5 = -Q[1]*Q[1]
3325    T6 = Q[1]*Q[2]
3326    T7 = Q[1]*Q[3]
3327    T8 = -Q[2]*Q[2]
3328    T9 = Q[2]*Q[3]
3329    T10 = -Q[3]*Q[3]
3330    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3331    VP = 2.*np.inner(V,M)
3332    return VP+V
3333   
3334def Q2Mat(Q):
3335    ''' make rotation matrix from quaternion
3336        q=r+ai+bj+ck
3337    '''
3338    QN = normQ(Q)
3339    aa = QN[0]**2
3340    ab = QN[0]*QN[1]
3341    ac = QN[0]*QN[2]
3342    ad = QN[0]*QN[3]
3343    bb = QN[1]**2
3344    bc = QN[1]*QN[2]
3345    bd = QN[1]*QN[3]
3346    cc = QN[2]**2
3347    cd = QN[2]*QN[3]
3348    dd = QN[3]**2
3349    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
3350        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
3351        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
3352    return np.array(M)
3353   
3354def AV2Q(A,V):
3355    ''' convert angle (radians) & vector to quaternion
3356        q=r+ai+bj+ck
3357    '''
3358    Q = np.zeros(4)
3359    d = nl.norm(np.array(V))
3360    if d:
3361        V /= d
3362        if not A:       #==0.
3363            A = 2.*np.pi
3364        p = A/2.
3365        Q[0] = np.cos(p)
3366        Q[1:4] = V*np.sin(p)
3367    else:
3368        Q[3] = 1.
3369    return Q
3370   
3371def AVdeg2Q(A,V):
3372    ''' convert angle (degrees) & vector to quaternion
3373        q=r+ai+bj+ck
3374    '''
3375    Q = np.zeros(4)
3376    d = nl.norm(np.array(V))
3377    if not A:       #== 0.!
3378        A = 360.
3379    if d:
3380        V /= d
3381        p = A/2.
3382        Q[0] = cosd(p)
3383        Q[1:4] = V*sind(p)
3384    else:
3385        Q[3] = 1.
3386    return Q
3387   
3388def Q2AVdeg(Q):
3389    ''' convert quaternion to angle (degrees 0-360) & normalized vector
3390        q=r+ai+bj+ck
3391    '''
3392    A = 2.*acosd(Q[0])
3393    V = np.array(Q[1:])
3394    V /= sind(A/2.)
3395    return A,V
3396   
3397def Q2AV(Q):
3398    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
3399        q=r+ai+bj+ck
3400    '''
3401    A = 2.*np.arccos(Q[0])
3402    V = np.array(Q[1:])
3403    V /= np.sin(A/2.)
3404    return A,V
3405   
3406def randomQ(r0,r1,r2,r3):
3407    ''' create random quaternion from 4 random numbers in range (-1,1)
3408    '''
3409    sum = 0
3410    Q = np.array(4)
3411    Q[0] = r0
3412    sum += Q[0]**2
3413    Q[1] = np.sqrt(1.-sum)*r1
3414    sum += Q[1]**2
3415    Q[2] = np.sqrt(1.-sum)*r2
3416    sum += Q[2]**2
3417    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3418    return Q
3419   
3420def randomAVdeg(r0,r1,r2,r3):
3421    ''' create random angle (deg),vector from 4 random number in range (-1,1)
3422    '''
3423    return Q2AVdeg(randomQ(r0,r1,r2,r3))
3424   
3425def makeQuat(A,B,C):
3426    ''' Make quaternion from rotation of A vector to B vector about C axis
3427
3428        :param np.array A,B,C: Cartesian 3-vectors
3429        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3430    '''
3431
3432    V1 = np.cross(A,C)
3433    V2 = np.cross(B,C)
3434    if nl.norm(V1)*nl.norm(V2):
3435        V1 /= nl.norm(V1)
3436        V2 /= nl.norm(V2)
3437        V3 = np.cross(V1,V2)
3438    else:
3439        V3 = np.zeros(3)
3440    Q = np.array([0.,0.,0.,1.])
3441    D = 0.
3442    if nl.norm(V3):
3443        V3 /= nl.norm(V3)
3444        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3445        D = np.arccos(D1)/2.0
3446        V1 = C-V3
3447        V2 = C+V3
3448        DM = nl.norm(V1)
3449        DP = nl.norm(V2)
3450        S = np.sin(D)
3451        Q[0] = np.cos(D)
3452        Q[1:] = V3*S
3453        D *= 2.
3454        if DM > DP:
3455            D *= -1.
3456    return Q,D
3457   
3458def annealtests():
3459    from numpy import cos
3460    # minimum expected at ~-0.195
3461    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3462    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3463    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3464    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3465
3466    # minimum expected at ~[-0.195, -0.1]
3467    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
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='cauchy')
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='fast')
3470    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')
3471
3472
3473if __name__ == '__main__':
3474    annealtests()
Note: See TracBrowser for help on using the repository browser.