source: trunk/GSASIImath.py @ 2005

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

add print every 100 refl. in structurefactorderv routines
fix atom loc problem

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