source: trunk/GSASIImath.py @ 1935

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

apply Uiso (equiv) for H-atom generation/update
revise H-atom addition to allow consecutive addition
deleted H atoms treated correctly

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