source: trunk/GSASIImath.py @ 970

Last change on this file since 970 was 970, checked in by vondreele, 9 years ago

more fixups of MC/SA

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