source: trunk/GSASIImath.py @ 1952

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

Add modulation of thermal parameters to drawing

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