source: trunk/GSASIImath.py @ 1951

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

implement drawing of incommensurate structures - live variation wrt tau
required careful attention to siteSym text & some reorganization of modules
more work on special pos code

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