source: trunk/GSASIImath.py @ 2072

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

replace local codes with routines for ZgZag/Block? wave calcs.

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