source: trunk/GSASIImath.py @ 1646

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

allow search for negative density peaks in charge flip and Fourier maps
plot negative density in neutron maps & negative peaks in red

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