source: trunk/GSASIImath.py @ 2040

Last change on this file since 2040 was 2040, checked in by vondreele, 7 years ago

make sure XSC is reflection import default Type
convert Biso to Uiso for J2K atoms import
get plots to show hklm for SS 2D structure factor plots on status bar & tool tip
make wave coeff 5 places not 4
now get correct Fcalc2 for SS reflections; proper derivatives for x & Uij
(not for x waves or U waves)
correct error in single crystal Fcalc
2 & it's derivatives

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