source: trunk/GSASIImath.py @ 1626

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

Changed "Transform atom " to "Transform draw atom" to fix menu problem
implement a 4D version of findOffset - works OK
implement facility to allow moving map on 4th dimension

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