source: trunk/GSASIImath.py @ 1995

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

fix SS atom plotting problem
SS map scale issue

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