source: trunk/GSASIImath.py @ 1645

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

define a mapDataType? for the data type used to make the Fourier/charge flip map & then show negative density for neutron maps.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 118.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2015-02-10 22:16:02 +0000 (Tue, 10 Feb 2015) $
5# $Author: vondreele $
6# $Revision: 1645 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1645 2015-02-10 22:16:02Z 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: 1645 $")
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):
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 data: the phase data structure
2088
2089    :returns: (peaks,mags,dzeros) where
2090
2091        * peaks : ndarray
2092          x,y,z positions of the peaks found in the map
2093        * mags : ndarray
2094          the magnitudes of the peaks
2095        * dzeros : ndarray
2096          the distance of the peaks from  the unit cell origin
2097
2098    '''       
2099    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2100   
2101    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
2102   
2103#    def noDuplicate(xyz,peaks,Amat):
2104#        XYZ = np.inner(Amat,xyz)
2105#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2106#            print ' Peak',xyz,' <0.5A from another peak'
2107#            return False
2108#        return True
2109#                           
2110    def fixSpecialPos(xyz,SGData,Amat):
2111        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
2112        X = []
2113        xyzs = [equiv[0] for equiv in equivs]
2114        for x in xyzs:
2115            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
2116                X.append(x)
2117        if len(X) > 1:
2118            return np.average(X,axis=0)
2119        else:
2120            return xyz
2121       
2122    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
2123        Mag,x0,y0,z0,sig = parms
2124        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
2125#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
2126        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
2127       
2128    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
2129        Mag,x0,y0,z0,sig = parms
2130        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2131        return M
2132       
2133    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
2134        Mag,x0,y0,z0,sig = parms
2135        dMdv = np.zeros(([5,]+list(rX.shape)))
2136        delt = .01
2137        for i in range(5):
2138            parms[i] -= delt
2139            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2140            parms[i] += 2.*delt
2141            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2142            parms[i] -= delt
2143            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
2144        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2145        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
2146        dMdv = np.reshape(dMdv,(5,rX.size))
2147        Hess = np.inner(dMdv,dMdv)
2148       
2149        return Vec,Hess
2150       
2151    phaseName = generalData['Name']
2152    SGData = generalData['SGData']
2153    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2154    peaks = []
2155    mags = []
2156    dzeros = []
2157    try:
2158        mapData = generalData['Map']
2159        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
2160        rho = copy.copy(mapData['rho'])     #don't mess up original
2161        mapHalf = np.array(rho.shape)/2
2162        res = mapData['Resolution']
2163        incre = np.array(rho.shape,dtype=np.float)
2164        step = max(1.0,1./res)+1
2165        steps = np.array(3*[step,])
2166    except KeyError:
2167        print '**** ERROR - Fourier map not defined'
2168        return peaks,mags
2169    rhoMask = ma.array(rho,mask=(rho<contLevel))
2170    indices = (-1,0,1)
2171    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2172    for roll in rolls:
2173        if np.any(roll):
2174            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
2175    indx = np.transpose(rhoMask.nonzero())
2176    peaks = indx/incre
2177    mags = rhoMask[rhoMask.nonzero()]
2178    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
2179        rho = rollMap(rho,ind)
2180        rMM = mapHalf-steps
2181        rMP = mapHalf+steps+1
2182        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2183        peakInt = np.sum(rhoPeak)*res**3
2184        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2185        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
2186        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
2187            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
2188        x1 = result[0]
2189        if not np.any(x1 < 0):
2190            mag = x1[0]
2191            peak = (np.array(x1[1:4])-ind)/incre
2192        peak = fixSpecialPos(peak,SGData,Amat)
2193        rho = rollMap(rho,-ind)       
2194    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
2195    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
2196   
2197def sortArray(data,pos,reverse=False):
2198    '''data is a list of items
2199    sort by pos in list; reverse if True
2200    '''
2201    T = []
2202    for i,M in enumerate(data):
2203        T.append((M[pos],i))
2204    D = dict(zip(T,data))
2205    T.sort()
2206    if reverse:
2207        T.reverse()
2208    X = []
2209    for key in T:
2210        X.append(D[key])
2211    return X
2212
2213def PeaksEquiv(data,Ind):
2214    '''Find the equivalent map peaks for those selected. Works on the
2215    contents of data['Map Peaks'].
2216
2217    :param data: the phase data structure
2218    :param list Ind: list of selected peak indices
2219    :returns: augmented list of peaks including those related by symmetry to the
2220      ones in Ind
2221
2222    '''       
2223    def Duplicate(xyz,peaks,Amat):
2224        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2225            return True
2226        return False
2227                           
2228    generalData = data['General']
2229    cell = generalData['Cell'][1:7]
2230    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2231    A = G2lat.cell2A(cell)
2232    SGData = generalData['SGData']
2233    mapPeaks = data['Map Peaks']
2234    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
2235    Indx = {}
2236    for ind in Ind:
2237        xyz = np.array(mapPeaks[ind][1:4])
2238        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
2239        for jnd,xyz in enumerate(XYZ):       
2240            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
2241    Ind = []
2242    for ind in Indx:
2243        if Indx[ind]:
2244            Ind.append(ind)
2245    return Ind
2246               
2247def PeaksUnique(data,Ind):
2248    '''Finds the symmetry unique set of peaks from those selected. Works on the
2249    contents of data['Map Peaks'].
2250
2251    :param data: the phase data structure
2252    :param list Ind: list of selected peak indices
2253    :returns: the list of symmetry unique peaks from among those given in Ind
2254
2255    '''       
2256#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
2257
2258    def noDuplicate(xyz,peaks,Amat):
2259        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=1.0) for peak in peaks]:
2260            return False
2261        return True
2262                           
2263    generalData = data['General']
2264    cell = generalData['Cell'][1:7]
2265    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2266    A = G2lat.cell2A(cell)
2267    SGData = generalData['SGData']
2268    mapPeaks = data['Map Peaks']
2269    Indx = {}
2270    XYZ = {}
2271    for ind in Ind:
2272        XYZ[ind] = np.array(mapPeaks[ind][1:4])
2273        Indx[ind] = True
2274    for ind in Ind:
2275        if Indx[ind]:
2276            xyz = XYZ[ind]
2277            for jnd in Ind:
2278                if ind != jnd and Indx[jnd]:                       
2279                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
2280                    xyzs = np.array([equiv[0] for equiv in Equiv])
2281                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
2282    Ind = []
2283    for ind in Indx:
2284        if Indx[ind]:
2285            Ind.append(ind)
2286    return Ind
2287   
2288################################################################################
2289##### single peak fitting profile fxn stuff
2290################################################################################
2291
2292def getCWsig(ins,pos):
2293    '''get CW peak profile sigma
2294   
2295    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
2296      as values only
2297    :param float pos: 2-theta of peak
2298    :returns: float getCWsig: peak sigma
2299   
2300    '''
2301    tp = tand(pos/2.0)
2302    return ins['U']*tp**2+ins['V']*tp+ins['W']
2303   
2304def getCWsigDeriv(pos):
2305    '''get derivatives of CW peak profile sigma wrt U,V, & W
2306   
2307    :param float pos: 2-theta of peak
2308   
2309    :returns: list getCWsigDeriv: d(sig)/dU, d(sig)/dV & d(sig)/dW
2310   
2311    '''
2312    tp = tand(pos/2.0)
2313    return tp**2,tp,1.0
2314   
2315def getCWgam(ins,pos):
2316    '''get CW peak profile gamma
2317   
2318    :param dict ins: instrument parameters with at least 'X' & 'Y'
2319      as values only
2320    :param float pos: 2-theta of peak
2321    :returns: float getCWgam: peak gamma
2322   
2323    '''
2324    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2325   
2326def getCWgamDeriv(pos):
2327    '''get derivatives of CW peak profile gamma wrt X & Y
2328   
2329    :param float pos: 2-theta of peak
2330   
2331    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
2332   
2333    '''
2334    return 1./cosd(pos/2.0),tand(pos/2.0)
2335   
2336def getTOFsig(ins,dsp):
2337    '''get TOF peak profile sigma
2338   
2339    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
2340      as values only
2341    :param float dsp: d-spacing of peak
2342   
2343    :returns: float getTOFsig: peak sigma
2344   
2345    '''
2346    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
2347   
2348def getTOFsigDeriv(dsp):
2349    '''get derivatives of TOF peak profile gamma wrt sig-0, sig-1, & sig-q
2350   
2351    :param float dsp: d-spacing of peak
2352   
2353    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
2354   
2355    '''
2356    return 1.0,dsp**2,dsp**4,1./dsp**2
2357   
2358def getTOFgamma(ins,dsp):
2359    '''get TOF peak profile gamma
2360   
2361    :param dict ins: instrument parameters with at least 'X' & 'Y'
2362      as values only
2363    :param float dsp: d-spacing of peak
2364   
2365    :returns: float getTOFgamma: peak gamma
2366   
2367    '''
2368    return ins['X']*dsp+ins['Y']*dsp**2
2369   
2370def getTOFgammaDeriv(dsp):
2371    '''get derivatives of TOF peak profile gamma wrt X & Y
2372   
2373    :param float dsp: d-spacing of peak
2374   
2375    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
2376   
2377    '''
2378    return dsp,dsp**2
2379   
2380def getTOFbeta(ins,dsp):
2381    '''get TOF peak profile beta
2382   
2383    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
2384      as values only
2385    :param float dsp: d-spacing of peak
2386   
2387    :returns: float getTOFbeta: peak beat
2388   
2389    '''
2390    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
2391   
2392def getTOFbetaDeriv(dsp):
2393    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
2394   
2395    :param float dsp: d-spacing of peak
2396   
2397    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
2398   
2399    '''
2400    return 1.0,1./dsp**4,1./dsp**2
2401   
2402def getTOFalpha(ins,dsp):
2403    '''get TOF peak profile alpha
2404   
2405    :param dict ins: instrument parameters with at least 'alpha'
2406      as values only
2407    :param float dsp: d-spacing of peak
2408   
2409    :returns: flaot getTOFalpha: peak alpha
2410   
2411    '''
2412    return ins['alpha']/dsp
2413   
2414def getTOFalphaDeriv(dsp):
2415    '''get derivatives of TOF peak profile beta wrt alpha
2416   
2417    :param float dsp: d-spacing of peak
2418   
2419    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
2420   
2421    '''
2422    return 1./dsp
2423   
2424def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2425    '''set starting peak parameters for single peak fits from plot selection or auto selection
2426   
2427    :param dict Parms: instrument parameters dictionary
2428    :param dict Parms2: table lookup for TOF profile coefficients
2429    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
2430    :param float mag: peak top magnitude from pick
2431    :param bool ifQ: True if pos in Q
2432    :param bool useFit: True if use fitted CW Parms values (not defaults)
2433   
2434    :returns: list XY: peak list entry:
2435        for CW: [pos,0,mag,1,sig,0,gam,0]
2436        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2437        NB: mag refinement set by default, all others off
2438   
2439    '''
2440    ind = 0
2441    if useFit:
2442        ind = 1
2443    ins = {}
2444    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
2445        for x in ['U','V','W','X','Y']:
2446            ins[x] = Parms[x][ind]
2447        if ifQ:                              #qplot - convert back to 2-theta
2448            pos = 2.0*asind(pos*wave/(4*math.pi))
2449        sig = getCWsig(ins,pos)
2450        gam = getCWgam(ins,pos)           
2451        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
2452    else:
2453        if ifQ:
2454            dsp = 2.*np.pi/pos
2455            pos = Parms['difC']*dsp
2456        else:
2457            dsp = pos/Parms['difC'][1]
2458        if 'Pdabc' in Parms2:
2459            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
2460                ins[x] = Parms[x][ind]
2461            Pdabc = Parms2['Pdabc'].T
2462            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2463            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2464        else:
2465            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
2466                ins[x] = Parms[x][ind]
2467            alp = getTOFalpha(ins,dsp)
2468            bet = getTOFbeta(ins,dsp)
2469        sig = getTOFsig(ins,dsp)
2470        gam = getTOFgamma(ins,dsp)
2471        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2472    return XY
2473   
2474################################################################################
2475##### MC/SA stuff
2476################################################################################
2477
2478#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2479# Original Author: Travis Oliphant 2002
2480# Bug-fixes in 2006 by Tim Leslie
2481
2482
2483import numpy
2484from numpy import asarray, tan, exp, ones, squeeze, sign, \
2485     all, log, sqrt, pi, shape, array, minimum, where
2486from numpy import random
2487
2488#__all__ = ['anneal']
2489
2490_double_min = numpy.finfo(float).min
2491_double_max = numpy.finfo(float).max
2492class base_schedule(object):
2493    def __init__(self):
2494        self.dwell = 20
2495        self.learn_rate = 0.5
2496        self.lower = -10
2497        self.upper = 10
2498        self.Ninit = 50
2499        self.accepted = 0
2500        self.tests = 0
2501        self.feval = 0
2502        self.k = 0
2503        self.T = None
2504
2505    def init(self, **options):
2506        self.__dict__.update(options)
2507        self.lower = asarray(self.lower)
2508        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2509        self.upper = asarray(self.upper)
2510        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2511        self.k = 0
2512        self.accepted = 0
2513        self.feval = 0
2514        self.tests = 0
2515
2516    def getstart_temp(self, best_state):
2517        """ Find a matching starting temperature and starting parameters vector
2518        i.e. find x0 such that func(x0) = T0.
2519
2520        Parameters
2521        ----------
2522        best_state : _state
2523            A _state object to store the function value and x0 found.
2524
2525        returns
2526        -------
2527        x0 : array
2528            The starting parameters vector.
2529        """
2530
2531        assert(not self.dims is None)
2532        lrange = self.lower
2533        urange = self.upper
2534        fmax = _double_min
2535        fmin = _double_max
2536        for _ in range(self.Ninit):
2537            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2538            fval = self.func(x0, *self.args)
2539            self.feval += 1
2540            if fval > fmax:
2541                fmax = fval
2542            if fval < fmin:
2543                fmin = fval
2544                best_state.cost = fval
2545                best_state.x = array(x0)
2546
2547        self.T0 = (fmax-fmin)*1.5
2548        return best_state.x
2549       
2550    def set_range(self,x0,frac):
2551        delrange = frac*(self.upper-self.lower)
2552        self.upper = x0+delrange
2553        self.lower = x0-delrange
2554
2555    def accept_test(self, dE):
2556        T = self.T
2557        self.tests += 1
2558        if dE < 0:
2559            self.accepted += 1
2560            return 1
2561        p = exp(-dE*1.0/self.boltzmann/T)
2562        if (p > random.uniform(0.0, 1.0)):
2563            self.accepted += 1
2564            return 1
2565        return 0
2566
2567    def update_guess(self, x0):
2568        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2569
2570    def update_temp(self, x0):
2571        pass
2572
2573
2574#  A schedule due to Lester Ingber modified to use bounds - OK
2575class fast_sa(base_schedule):
2576    def init(self, **options):
2577        self.__dict__.update(options)
2578        if self.m is None:
2579            self.m = 1.0
2580        if self.n is None:
2581            self.n = 1.0
2582        self.c = self.m * exp(-self.n * self.quench)
2583
2584    def update_guess(self, x0):
2585        x0 = asarray(x0)
2586        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2587        T = self.T
2588        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2589        xnew = xc*(self.upper - self.lower)+self.lower
2590        return xnew
2591#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2592#        xc = y*(self.upper - self.lower)
2593#        xnew = x0 + xc
2594#        return xnew
2595
2596    def update_temp(self):
2597        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2598        self.k += 1
2599        return
2600
2601class cauchy_sa(base_schedule):     #modified to use bounds - not good
2602    def update_guess(self, x0):
2603        x0 = asarray(x0)
2604        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2605        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2606        xnew = xc*(self.upper - self.lower)+self.lower
2607        return xnew
2608#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2609#        xc = self.learn_rate * self.T * tan(numbers)
2610#        xnew = x0 + xc
2611#        return xnew
2612
2613    def update_temp(self):
2614        self.T = self.T0/(1+self.k)
2615        self.k += 1
2616        return
2617
2618class boltzmann_sa(base_schedule):
2619#    def update_guess(self, x0):
2620#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2621#        x0 = asarray(x0)
2622#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2623#
2624#        xnew = x0 + xc*std*self.learn_rate
2625#        return xnew
2626
2627    def update_temp(self):
2628        self.k += 1
2629        self.T = self.T0 / log(self.k+1.0)
2630        return
2631
2632class log_sa(base_schedule):        #OK
2633
2634    def init(self,**options):
2635        self.__dict__.update(options)
2636       
2637    def update_guess(self,x0):     #same as default
2638        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2639       
2640    def update_temp(self):
2641        self.k += 1
2642        self.T = self.T0*self.slope**self.k
2643       
2644class _state(object):
2645    def __init__(self):
2646        self.x = None
2647        self.cost = None
2648
2649# TODO:
2650#     allow for general annealing temperature profile
2651#     in that case use update given by alpha and omega and
2652#     variation of all previous updates and temperature?
2653
2654# Simulated annealing
2655
2656def anneal(func, x0, args=(), schedule='fast', full_output=0,
2657           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2658           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2659           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2660           ranRange=0.10,autoRan=False,dlg=None):
2661    """Minimize a function using simulated annealing.
2662
2663    Schedule is a schedule class implementing the annealing schedule.
2664    Available ones are 'fast', 'cauchy', 'boltzmann'
2665
2666    :param callable func: f(x, \*args)
2667        Function to be optimized.
2668    :param ndarray x0:
2669        Initial guess.
2670    :param tuple args:
2671        Extra parameters to `func`.
2672    :param base_schedule schedule:
2673        Annealing schedule to use (a class).
2674    :param bool full_output:
2675        Whether to return optional outputs.
2676    :param float T0:
2677        Initial Temperature (estimated as 1.2 times the largest
2678        cost-function deviation over random points in the range).
2679    :param float Tf:
2680        Final goal temperature.
2681    :param int maxeval:
2682        Maximum function evaluations.
2683    :param int maxaccept:
2684        Maximum changes to accept.
2685    :param int maxiter:
2686        Maximum cooling iterations.
2687    :param float learn_rate:
2688        Scale constant for adjusting guesses.
2689    :param float boltzmann:
2690        Boltzmann constant in acceptance test
2691        (increase for less stringent test at each temperature).
2692    :param float feps:
2693        Stopping relative error tolerance for the function value in
2694        last four coolings.
2695    :param float quench,m,n:
2696        Parameters to alter fast_sa schedule.
2697    :param float/ndarray lower,upper:
2698        Lower and upper bounds on `x`.
2699    :param int dwell:
2700        The number of times to search the space at each temperature.
2701    :param float slope:
2702        Parameter for log schedule
2703    :param bool ranStart=False:
2704        True for set 10% of ranges about x
2705
2706    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2707
2708     * xmin (ndarray): Point giving smallest value found.
2709     * Jmin (float): Minimum value of function found.
2710     * T (float): Final temperature.
2711     * feval (int): Number of function evaluations.
2712     * iters (int): Number of cooling iterations.
2713     * accept (int): Number of tests accepted.
2714     * retval (int): Flag indicating stopping condition:
2715
2716              *  0: Points no longer changing
2717              *  1: Cooled to final temperature
2718              *  2: Maximum function evaluations
2719              *  3: Maximum cooling iterations reached
2720              *  4: Maximum accepted query locations reached
2721              *  5: Final point not the minimum amongst encountered points
2722
2723    *Notes*:
2724    Simulated annealing is a random algorithm which uses no derivative
2725    information from the function being optimized. In practice it has
2726    been more useful in discrete optimization than continuous
2727    optimization, as there are usually better algorithms for continuous
2728    optimization problems.
2729
2730    Some experimentation by trying the difference temperature
2731    schedules and altering their parameters is likely required to
2732    obtain good performance.
2733
2734    The randomness in the algorithm comes from random sampling in numpy.
2735    To obtain the same results you can call numpy.random.seed with the
2736    same seed immediately before calling scipy.optimize.anneal.
2737
2738    We give a brief description of how the three temperature schedules
2739    generate new points and vary their temperature. Temperatures are
2740    only updated with iterations in the outer loop. The inner loop is
2741    over xrange(dwell), and new points are generated for
2742    every iteration in the inner loop. (Though whether the proposed
2743    new points are accepted is probabilistic.)
2744
2745    For readability, let d denote the dimension of the inputs to func.
2746    Also, let x_old denote the previous state, and k denote the
2747    iteration number of the outer loop. All other variables not
2748    defined below are input variables to scipy.optimize.anneal itself.
2749
2750    In the 'fast' schedule the updates are ::
2751
2752        u ~ Uniform(0, 1, size=d)
2753        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2754        xc = y * (upper - lower)
2755        x_new = x_old + xc
2756
2757        c = n * exp(-n * quench)
2758        T_new = T0 * exp(-c * k**quench)
2759
2760
2761    In the 'cauchy' schedule the updates are ::
2762
2763        u ~ Uniform(-pi/2, pi/2, size=d)
2764        xc = learn_rate * T * tan(u)
2765        x_new = x_old + xc
2766
2767        T_new = T0 / (1+k)
2768
2769    In the 'boltzmann' schedule the updates are ::
2770
2771        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2772        y ~ Normal(0, std, size=d)
2773        x_new = x_old + learn_rate * y
2774
2775        T_new = T0 / log(1+k)
2776
2777    """
2778    x0 = asarray(x0)
2779    lower = asarray(lower)
2780    upper = asarray(upper)
2781
2782    schedule = eval(schedule+'_sa()')
2783    #   initialize the schedule
2784    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2785                  learn_rate=learn_rate, lower=lower, upper=upper,
2786                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2787
2788    current_state, last_state, best_state = _state(), _state(), _state()
2789    if ranStart:
2790        schedule.set_range(x0,ranRange)
2791    if T0 is None:
2792        x0 = schedule.getstart_temp(best_state)
2793    else:
2794        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
2795        best_state.x = None
2796        best_state.cost = numpy.Inf
2797
2798    last_state.x = asarray(x0).copy()
2799    fval = func(x0,*args)
2800    schedule.feval += 1
2801    last_state.cost = fval
2802    if last_state.cost < best_state.cost:
2803        best_state.cost = fval
2804        best_state.x = asarray(x0).copy()
2805    schedule.T = schedule.T0
2806    fqueue = [100, 300, 500, 700]
2807    iters = 1
2808    keepGoing = True
2809    bestn = 0
2810    while keepGoing:
2811        retval = 0
2812        for n in xrange(dwell):
2813            current_state.x = schedule.update_guess(last_state.x)
2814            current_state.cost = func(current_state.x,*args)
2815            schedule.feval += 1
2816
2817            dE = current_state.cost - last_state.cost
2818            if schedule.accept_test(dE):
2819                last_state.x = current_state.x.copy()
2820                last_state.cost = current_state.cost
2821                if last_state.cost < best_state.cost:
2822                    best_state.x = last_state.x.copy()
2823                    best_state.cost = last_state.cost
2824                    bestn = n
2825                    if best_state.cost < 1.0 and autoRan:
2826                        schedule.set_range(x0,best_state.cost/2.)                       
2827        if dlg:
2828            GoOn = dlg.Update(min(100.,best_state.cost*100),
2829                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2830                    'Best trial:',bestn,  \
2831                    'MC/SA Residual:',best_state.cost*100,'%', \
2832                    ))[0]
2833            if not GoOn:
2834                best_state.x = last_state.x.copy()
2835                best_state.cost = last_state.cost
2836                retval = 5
2837        schedule.update_temp()
2838        iters += 1
2839        # Stopping conditions
2840        # 0) last saved values of f from each cooling step
2841        #     are all very similar (effectively cooled)
2842        # 1) Tf is set and we are below it
2843        # 2) maxeval is set and we are past it
2844        # 3) maxiter is set and we are past it
2845        # 4) maxaccept is set and we are past it
2846        # 5) user canceled run via progress bar
2847
2848        fqueue.append(squeeze(last_state.cost))
2849        fqueue.pop(0)
2850        af = asarray(fqueue)*1.0
2851        if retval == 5:
2852            print ' User terminated run; incomplete MC/SA'
2853            keepGoing = False
2854            break
2855        if all(abs((af-af[0])/af[0]) < feps):
2856            retval = 0
2857            if abs(af[-1]-best_state.cost) > feps*10:
2858                retval = 5
2859#                print "Warning: Cooled to %f at %s but this is not" \
2860#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2861#                      + " the smallest point found."
2862            break
2863        if (Tf is not None) and (schedule.T < Tf):
2864            retval = 1
2865            break
2866        if (maxeval is not None) and (schedule.feval > maxeval):
2867            retval = 2
2868            break
2869        if (iters > maxiter):
2870            print "Warning: Maximum number of iterations exceeded."
2871            retval = 3
2872            break
2873        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2874            retval = 4
2875            break
2876
2877    if full_output:
2878        return best_state.x, best_state.cost, schedule.T, \
2879               schedule.feval, iters, schedule.accepted, retval
2880    else:
2881        return best_state.x, retval
2882
2883def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
2884    outlist = []
2885    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
2886    for n in range(iCyc):
2887        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2888        outlist.append(result[0])
2889        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2890    out_q.put(outlist)
2891
2892def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2893    import multiprocessing as mp
2894   
2895    nprocs = mp.cpu_count()
2896    out_q = mp.Queue()
2897    procs = []
2898    iCyc = np.zeros(nprocs)
2899    for i in range(nCyc):
2900        iCyc[i%nprocs] += 1
2901    for i in range(nprocs):
2902        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
2903        procs.append(p)
2904        p.start()
2905    resultlist = []
2906    for i in range(nprocs):
2907        resultlist += out_q.get()
2908    for p in procs:
2909        p.join()
2910    return resultlist
2911
2912def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2913    '''default doc string
2914   
2915    :param type name: description
2916   
2917    :returns: type name: description
2918   
2919    '''
2920   
2921    global tsum
2922    tsum = 0.
2923   
2924    def getMDparms(item,pfx,parmDict,varyList):
2925        parmDict[pfx+'MDaxis'] = item['axis']
2926        parmDict[pfx+'MDval'] = item['Coef'][0]
2927        if item['Coef'][1]:
2928            varyList += [pfx+'MDval',]
2929            limits = item['Coef'][2]
2930            lower.append(limits[0])
2931            upper.append(limits[1])
2932                       
2933    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2934        parmDict[pfx+'Atype'] = item['atType']
2935        aTypes |= set([item['atType'],]) 
2936        pstr = ['Ax','Ay','Az']
2937        XYZ = [0,0,0]
2938        for i in range(3):
2939            name = pfx+pstr[i]
2940            parmDict[name] = item['Pos'][0][i]
2941            XYZ[i] = parmDict[name]
2942            if item['Pos'][1][i]:
2943                varyList += [name,]
2944                limits = item['Pos'][2][i]
2945                lower.append(limits[0])
2946                upper.append(limits[1])
2947        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2948           
2949    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2950        parmDict[mfx+'MolCent'] = item['MolCent']
2951        parmDict[mfx+'RBId'] = item['RBId']
2952        pstr = ['Px','Py','Pz']
2953        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
2954        for i in range(3):
2955            name = mfx+pstr[i]
2956            parmDict[name] = item['Pos'][0][i]
2957            if item['Pos'][1][i]:
2958                varyList += [name,]
2959                limits = item['Pos'][2][i]
2960                lower.append(limits[0])
2961                upper.append(limits[1])
2962        AV = item['Ori'][0]
2963        A = AV[0]
2964        V = AV[1:]
2965        for i in range(4):
2966            name = mfx+ostr[i]
2967            if i:
2968                parmDict[name] = V[i-1]
2969            else:
2970                parmDict[name] = A
2971            if item['Ovar'] == 'AV':
2972                varyList += [name,]
2973                limits = item['Ori'][2][i]
2974                lower.append(limits[0])
2975                upper.append(limits[1])
2976            elif item['Ovar'] == 'A' and not i:
2977                varyList += [name,]
2978                limits = item['Ori'][2][i]
2979                lower.append(limits[0])
2980                upper.append(limits[1])
2981        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2982            for i in range(len(item['Tor'][0])):
2983                name = mfx+'Tor'+str(i)
2984                parmDict[name] = item['Tor'][0][i]
2985                if item['Tor'][1][i]:
2986                    varyList += [name,]
2987                    limits = item['Tor'][2][i]
2988                    lower.append(limits[0])
2989                    upper.append(limits[1])
2990        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2991        aTypes |= set(atypes)
2992        atNo += len(atypes)
2993        return atNo
2994               
2995    def GetAtomM(Xdata,SGData):
2996        Mdata = []
2997        for xyz in Xdata:
2998            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2999        return np.array(Mdata)
3000       
3001    def GetAtomT(RBdata,parmDict):
3002        'Needs a doc string'
3003        atNo = parmDict['atNo']
3004        nfixAt = parmDict['nfixAt']
3005        Tdata = atNo*[' ',]
3006        for iatm in range(nfixAt):
3007            parm = ':'+str(iatm)+':Atype'
3008            if parm in parmDict:
3009                Tdata[iatm] = aTypes.index(parmDict[parm])
3010        iatm = nfixAt
3011        for iObj in range(parmDict['nObj']):
3012            pfx = str(iObj)+':'
3013            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3014                if parmDict[pfx+'Type'] == 'Vector':
3015                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3016                    nAtm = len(RBRes['rbVect'][0])
3017                else:       #Residue
3018                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3019                    nAtm = len(RBRes['rbXYZ'])
3020                for i in range(nAtm):
3021                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
3022                    iatm += 1
3023            elif parmDict[pfx+'Type'] == 'Atom':
3024                atNo = parmDict[pfx+'atNo']
3025                parm = pfx+'Atype'              #remove extra ':'
3026                if parm in parmDict:
3027                    Tdata[atNo] = aTypes.index(parmDict[parm])
3028                iatm += 1
3029            else:
3030                continue        #skips March Dollase
3031        return Tdata
3032       
3033    def GetAtomX(RBdata,parmDict):
3034        'Needs a doc string'
3035        Bmat = parmDict['Bmat']
3036        atNo = parmDict['atNo']
3037        nfixAt = parmDict['nfixAt']
3038        Xdata = np.zeros((3,atNo))
3039        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
3040        for iatm in range(nfixAt):
3041            for key in keys:
3042                parm = ':'+str(iatm)+key
3043                if parm in parmDict:
3044                    keys[key][iatm] = parmDict[parm]
3045        iatm = nfixAt
3046        for iObj in range(parmDict['nObj']):
3047            pfx = str(iObj)+':'
3048            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3049                if parmDict[pfx+'Type'] == 'Vector':
3050                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3051                    vecs = RBRes['rbVect']
3052                    mags = RBRes['VectMag']
3053                    Cart = np.zeros_like(vecs[0])
3054                    for vec,mag in zip(vecs,mags):
3055                        Cart += vec*mag
3056                elif parmDict[pfx+'Type'] == 'Residue':
3057                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3058                    Cart = np.array(RBRes['rbXYZ'])
3059                    for itor,seq in enumerate(RBRes['rbSeq']):
3060                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
3061                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
3062                if parmDict[pfx+'MolCent'][1]:
3063                    Cart -= parmDict[pfx+'MolCent'][0]
3064                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3065                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
3066                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
3067                iatm += len(Cart)
3068            elif parmDict[pfx+'Type'] == 'Atom':
3069                atNo = parmDict[pfx+'atNo']
3070                for key in keys:
3071                    parm = pfx+key[1:]              #remove extra ':'
3072                    if parm in parmDict:
3073                        keys[key][atNo] = parmDict[parm]
3074                iatm += 1
3075            else:
3076                continue        #skips March Dollase
3077        return Xdata.T
3078       
3079    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
3080        allX = np.inner(Xdata,SGM)+SGT
3081        allT = np.repeat(Tdata,allX.shape[1])
3082        allM = np.repeat(Mdata,allX.shape[1])
3083        allX = np.reshape(allX,(-1,3))
3084        return allT,allM,allX
3085
3086    def getAllX(Xdata,SGM,SGT):
3087        allX = np.inner(Xdata,SGM)+SGT
3088        allX = np.reshape(allX,(-1,3))
3089        return allX
3090       
3091    def normQuaternions(RBdata,parmDict,varyList,values):
3092        for iObj in range(parmDict['nObj']):
3093            pfx = str(iObj)+':'
3094            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3095                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3096                A,V = Q2AVdeg(Qori)
3097                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
3098                    if i:
3099                        parmDict[pfx+name] = V[i-1]
3100                    else:
3101                        parmDict[pfx+name] = A
3102       
3103    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
3104        ''' Compute structure factors for all h,k,l for phase
3105        input:
3106            refList: [ref] where each ref = h,k,l,m,d,...
3107            rcov:   array[nref,nref] covariance terms between Fo^2 values
3108            ifInv:  bool True if centrosymmetric
3109            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
3110            RBdata: [dict] rigid body dictionary
3111            varyList: [list] names of varied parameters in MC/SA (not used here)           
3112            ParmDict: [dict] problem parameters
3113        puts result F^2 in each ref[5] in refList
3114        returns:
3115            delt-F*rcov*delt-F/sum(Fo^2)^2
3116        '''       
3117        global tsum
3118        t0 = time.time()
3119        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
3120        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
3121        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
3122        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
3123        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
3124        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
3125        refList[5] = Aterm
3126        if not ifInv:
3127            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
3128        if len(cosTable):        #apply MD correction
3129            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
3130        sumFcsq = np.sum(refList[5])
3131        scale = parmDict['sumFosq']/sumFcsq
3132        refList[5] *= scale
3133        refList[6] = refList[4]-refList[5]
3134        M = np.inner(refList[6],np.inner(rcov,refList[6]))
3135        tsum += (time.time()-t0)
3136        return M/np.sum(refList[4]**2)
3137
3138    sq8ln2 = np.sqrt(8*np.log(2))
3139    sq2pi = np.sqrt(2*np.pi)
3140    sq4pi = np.sqrt(4*np.pi)
3141    generalData = data['General']
3142    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3143    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
3144    SGData = generalData['SGData']
3145    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
3146    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
3147    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
3148    fixAtoms = data['Atoms']                       #if any
3149    cx,ct,cs = generalData['AtomPtrs'][:3]
3150    aTypes = set([])
3151    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
3152    varyList = []
3153    atNo = 0
3154    for atm in fixAtoms:
3155        pfx = ':'+str(atNo)+':'
3156        parmDict[pfx+'Atype'] = atm[ct]
3157        aTypes |= set([atm[ct],])
3158        pstr = ['Ax','Ay','Az']
3159        parmDict[pfx+'Amul'] = atm[cs+1]
3160        for i in range(3):
3161            name = pfx+pstr[i]
3162            parmDict[name] = atm[cx+i]
3163        atNo += 1
3164    parmDict['nfixAt'] = len(fixAtoms)       
3165    MCSA = generalData['MCSA controls']
3166    reflName = MCSA['Data source']
3167    phaseName = generalData['Name']
3168    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
3169    upper = []
3170    lower = []
3171    MDvec = np.zeros(3)
3172    for i,item in enumerate(MCSAObjs):
3173        mfx = str(i)+':'
3174        parmDict[mfx+'Type'] = item['Type']
3175        if item['Type'] == 'MD':
3176            getMDparms(item,mfx,parmDict,varyList)
3177            MDvec = np.array(item['axis'])
3178        elif item['Type'] == 'Atom':
3179            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
3180            parmDict[mfx+'atNo'] = atNo
3181            atNo += 1
3182        elif item['Type'] in ['Residue','Vector']:
3183            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
3184    parmDict['atNo'] = atNo                 #total no. of atoms
3185    parmDict['nObj'] = len(MCSAObjs)
3186    aTypes = list(aTypes)
3187    Tdata = GetAtomT(RBdata,parmDict)
3188    Xdata = GetAtomX(RBdata,parmDict)
3189    Mdata = GetAtomM(Xdata,SGData)
3190    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
3191    FFtables = G2el.GetFFtable(aTypes)
3192    refs = []
3193    allFF = []
3194    cosTable = []
3195    sumFosq = 0
3196    if 'PWDR' in reflName:
3197        for ref in reflData:
3198            h,k,l,m,d,pos,sig,gam,f = ref[:9]
3199            if d >= MCSA['dmin']:
3200                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
3201                SQ = 0.25/d**2
3202                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3203                refs.append([h,k,l,m,f*m,pos,sig])
3204                sumFosq += f*m
3205                Heqv = np.inner(np.array([h,k,l]),SGMT)
3206                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3207        nRef = len(refs)
3208        cosTable = np.array(cosTable)**2
3209        rcov = np.zeros((nRef,nRef))
3210        for iref,refI in enumerate(refs):
3211            rcov[iref][iref] = 1./(sq4pi*refI[6])
3212            for jref,refJ in enumerate(refs[:iref]):
3213                t1 = refI[6]**2+refJ[6]**2
3214                t2 = (refJ[5]-refI[5])**2/(2.*t1)
3215                if t2 > 10.:
3216                    rcov[iref][jref] = 0.
3217                else:
3218                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
3219        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3220        Rdiag = np.sqrt(np.diag(rcov))
3221        Rnorm = np.outer(Rdiag,Rdiag)
3222        rcov /= Rnorm
3223    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
3224        vNames = []
3225        pfx = str(data['pId'])+'::PWLref:'
3226        for iref,refI in enumerate(reflData):           #Pawley reflection set
3227            h,k,l,m,d,v,f,s = refI
3228            if d >= MCSA['dmin'] and v:       #skip unrefined ones
3229                vNames.append(pfx+str(iref))
3230                SQ = 0.25/d**2
3231                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3232                refs.append([h,k,l,m,f*m,iref,0.])
3233                sumFosq += f*m
3234                Heqv = np.inner(np.array([h,k,l]),SGMT)
3235                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3236        cosTable = np.array(cosTable)**2
3237        nRef = len(refs)
3238#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
3239        if covData['freshCOV'] or  MCSA['newDmin']:
3240            vList = covData['varyList']
3241            covMatrix = covData['covMatrix']
3242            rcov = getVCov(vNames,vList,covMatrix)
3243            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3244            Rdiag = np.sqrt(np.diag(rcov))
3245            Rdiag = np.where(Rdiag,Rdiag,1.0)
3246            Rnorm = np.outer(Rdiag,Rdiag)
3247            rcov /= Rnorm
3248            MCSA['rcov'] = rcov
3249            covData['freshCOV'] = False
3250            MCSA['newDmin'] = False
3251        else:
3252            rcov = MCSA['rcov']
3253    elif 'HKLF' in reflName:
3254        for ref in reflData:
3255            [h,k,l,m,d],f = ref[:5],ref[6]
3256            if d >= MCSA['dmin']:
3257                SQ = 0.25/d**2
3258                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3259                refs.append([h,k,l,m,f*m,0.,0.])
3260                sumFosq += f*m
3261        nRef = len(refs)
3262        rcov = np.identity(len(refs))
3263    allFF = np.array(allFF).T
3264    refs = np.array(refs).T
3265    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
3266    print ' Number of parameters varied: %d'%(len(varyList))
3267    parmDict['sumFosq'] = sumFosq
3268    x0 = [parmDict[val] for val in varyList]
3269    ifInv = SGData['SGInv']
3270    # consider replacing anneal with scipy.optimize.basinhopping
3271    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
3272        schedule=MCSA['Algorithm'], full_output=True,
3273        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
3274        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
3275        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
3276        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
3277        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
3278    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
3279#    for ref in refs.T:
3280#        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])
3281#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
3282    Result = [False,False,results[1],results[2],]+list(results[0])
3283    Result.append(varyList)
3284    return Result,tsum
3285
3286       
3287################################################################################
3288##### Quaternion stuff
3289################################################################################
3290
3291def prodQQ(QA,QB):
3292    ''' Grassman quaternion product
3293        QA,QB quaternions; q=r+ai+bj+ck
3294    '''
3295    D = np.zeros(4)
3296    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
3297    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
3298    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
3299    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
3300   
3301#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
3302#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
3303   
3304    return D
3305   
3306def normQ(QA):
3307    ''' get length of quaternion & normalize it
3308        q=r+ai+bj+ck
3309    '''
3310    n = np.sqrt(np.sum(np.array(QA)**2))
3311    return QA/n
3312   
3313def invQ(Q):
3314    '''
3315        get inverse of quaternion
3316        q=r+ai+bj+ck; q* = r-ai-bj-ck
3317    '''
3318    return Q*np.array([1,-1,-1,-1])
3319   
3320def prodQVQ(Q,V):
3321    """
3322    compute the quaternion vector rotation qvq-1 = v'
3323    q=r+ai+bj+ck
3324    """
3325    T2 = Q[0]*Q[1]
3326    T3 = Q[0]*Q[2]
3327    T4 = Q[0]*Q[3]
3328    T5 = -Q[1]*Q[1]
3329    T6 = Q[1]*Q[2]
3330    T7 = Q[1]*Q[3]
3331    T8 = -Q[2]*Q[2]
3332    T9 = Q[2]*Q[3]
3333    T10 = -Q[3]*Q[3]
3334    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3335    VP = 2.*np.inner(V,M)
3336    return VP+V
3337   
3338def Q2Mat(Q):
3339    ''' make rotation matrix from quaternion
3340        q=r+ai+bj+ck
3341    '''
3342    QN = normQ(Q)
3343    aa = QN[0]**2
3344    ab = QN[0]*QN[1]
3345    ac = QN[0]*QN[2]
3346    ad = QN[0]*QN[3]
3347    bb = QN[1]**2
3348    bc = QN[1]*QN[2]
3349    bd = QN[1]*QN[3]
3350    cc = QN[2]**2
3351    cd = QN[2]*QN[3]
3352    dd = QN[3]**2
3353    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
3354        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
3355        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
3356    return np.array(M)
3357   
3358def AV2Q(A,V):
3359    ''' convert angle (radians) & vector to quaternion
3360        q=r+ai+bj+ck
3361    '''
3362    Q = np.zeros(4)
3363    d = nl.norm(np.array(V))
3364    if d:
3365        V /= d
3366        if not A:       #==0.
3367            A = 2.*np.pi
3368        p = A/2.
3369        Q[0] = np.cos(p)
3370        Q[1:4] = V*np.sin(p)
3371    else:
3372        Q[3] = 1.
3373    return Q
3374   
3375def AVdeg2Q(A,V):
3376    ''' convert angle (degrees) & vector to quaternion
3377        q=r+ai+bj+ck
3378    '''
3379    Q = np.zeros(4)
3380    d = nl.norm(np.array(V))
3381    if not A:       #== 0.!
3382        A = 360.
3383    if d:
3384        V /= d
3385        p = A/2.
3386        Q[0] = cosd(p)
3387        Q[1:4] = V*sind(p)
3388    else:
3389        Q[3] = 1.
3390    return Q
3391   
3392def Q2AVdeg(Q):
3393    ''' convert quaternion to angle (degrees 0-360) & normalized vector
3394        q=r+ai+bj+ck
3395    '''
3396    A = 2.*acosd(Q[0])
3397    V = np.array(Q[1:])
3398    V /= sind(A/2.)
3399    return A,V
3400   
3401def Q2AV(Q):
3402    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
3403        q=r+ai+bj+ck
3404    '''
3405    A = 2.*np.arccos(Q[0])
3406    V = np.array(Q[1:])
3407    V /= np.sin(A/2.)
3408    return A,V
3409   
3410def randomQ(r0,r1,r2,r3):
3411    ''' create random quaternion from 4 random numbers in range (-1,1)
3412    '''
3413    sum = 0
3414    Q = np.array(4)
3415    Q[0] = r0
3416    sum += Q[0]**2
3417    Q[1] = np.sqrt(1.-sum)*r1
3418    sum += Q[1]**2
3419    Q[2] = np.sqrt(1.-sum)*r2
3420    sum += Q[2]**2
3421    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3422    return Q
3423   
3424def randomAVdeg(r0,r1,r2,r3):
3425    ''' create random angle (deg),vector from 4 random number in range (-1,1)
3426    '''
3427    return Q2AVdeg(randomQ(r0,r1,r2,r3))
3428   
3429def makeQuat(A,B,C):
3430    ''' Make quaternion from rotation of A vector to B vector about C axis
3431
3432        :param np.array A,B,C: Cartesian 3-vectors
3433        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3434    '''
3435
3436    V1 = np.cross(A,C)
3437    V2 = np.cross(B,C)
3438    if nl.norm(V1)*nl.norm(V2):
3439        V1 /= nl.norm(V1)
3440        V2 /= nl.norm(V2)
3441        V3 = np.cross(V1,V2)
3442    else:
3443        V3 = np.zeros(3)
3444    Q = np.array([0.,0.,0.,1.])
3445    D = 0.
3446    if nl.norm(V3):
3447        V3 /= nl.norm(V3)
3448        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3449        D = np.arccos(D1)/2.0
3450        V1 = C-V3
3451        V2 = C+V3
3452        DM = nl.norm(V1)
3453        DP = nl.norm(V2)
3454        S = np.sin(D)
3455        Q[0] = np.cos(D)
3456        Q[1:] = V3*S
3457        D *= 2.
3458        if DM > DP:
3459            D *= -1.
3460    return Q,D
3461   
3462def annealtests():
3463    from numpy import cos
3464    # minimum expected at ~-0.195
3465    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3466    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3467    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3468    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3469
3470    # minimum expected at ~[-0.195, -0.1]
3471    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
3472    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')
3473    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')
3474    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')
3475
3476
3477if __name__ == '__main__':
3478    annealtests()
Note: See TracBrowser for help on using the repository browser.