source: trunk/GSASIImath.py @ 2089

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

fix mod vector derivs
work on ZigZag/Block? derivs - seem ok on value & powder data, not pos or single xtal.

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