source: trunk/GSASIImath.py @ 2073

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

numeric derivatives for ZigZag/Block? waves

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