source: trunk/GSASIImath.py @ 984

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

update Win2.7 g77 binaries
cov matrix now not modified by lam
no label plotting at all - fails
set MC/SA page after mc/sa run
fix 64 bit bug in switching to map peaks page

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