source: trunk/GSASIImath.py @ 1996

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

fix Fade problem
ModulationDeriv? work
twin sf error fixed

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