source: trunk/GSASIImath.py @ 2319

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

remove stray print & add \n to Dist argument for better column heading
fix symm problem in CalcDist? for seq distances

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