source: trunk/GSASIImath.py @ 2061

Last change on this file since 2061 was 2061, checked in by vondreele, 7 years ago

modify export app peak lists to include d-spacing as a column
fix 2x error in U2Uij & Uij2U in G2spc - now same as one in G2lat
fix math errors for Uij modulation derivatives - all ok now

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