source: trunk/GSASIImath.py @ 1635

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

modulation functions & constraints

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