source: trunk/GSASIImath.py @ 2075

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

explicit plot destroy when GSAS2 closes
work on scrolling on phase data tabs - scroll back to place last used
fix(?) Block/ZigZag? derivatives
eliminate Refine & max index buttons for mod vector for single crystal data

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