source: trunk/GSASIImath.py @ 2070

Last change on this file since 2070 was 2070, checked in by vondreele, 6 years ago

remove extraneous argument from G2mth.posFourier
work on ZigZag/Block? SS functions; function OK, derivatives not.

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