source: trunk/GSASIIsasd.py @ 1743

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

fix to data selection issue
formatting error for histogram scale
change ':' to ';' for SASD parameter names
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 57.4 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII small angle calculation module*
5=======================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2015-03-19 16:44:06 +0000 (Thu, 19 Mar 2015) $
10# $Author: vondreele $
11# $Revision: 1743 $
12# $URL: trunk/GSASIIsasd.py $
13# $Id: GSASIIsasd.py 1743 2015-03-19 16:44:06Z vondreele $
14########### SVN repository information ###################
15import os
16import sys
17import math
18import time
19
20import numpy as np
21import scipy as sp
22import numpy.linalg as nl
23from numpy.fft import ifft, fft, fftshift
24import scipy.special as scsp
25import scipy.interpolate as si
26import scipy.stats as st
27import scipy.optimize as so
28#import pdb
29
30import GSASIIpath
31GSASIIpath.SetVersionNumber("$Revision: 1743 $")
32import GSASIIlattice as G2lat
33import GSASIIspc as G2spc
34import GSASIIElem as G2elem
35import GSASIIgrid as G2gd
36import GSASIIIO as G2IO
37import GSASIImath as G2mth
38import GSASIIpwd as G2pwd
39
40# trig functions in degrees
41sind = lambda x: math.sin(x*math.pi/180.)
42asind = lambda x: 180.*math.asin(x)/math.pi
43tand = lambda x: math.tan(x*math.pi/180.)
44atand = lambda x: 180.*math.atan(x)/math.pi
45atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
46cosd = lambda x: math.cos(x*math.pi/180.)
47acosd = lambda x: 180.*math.acos(x)/math.pi
48rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
49#numpy versions
50npsind = lambda x: np.sin(x*np.pi/180.)
51npasind = lambda x: 180.*np.arcsin(x)/math.pi
52npcosd = lambda x: np.cos(x*math.pi/180.)
53npacosd = lambda x: 180.*np.arccos(x)/math.pi
54nptand = lambda x: np.tan(x*math.pi/180.)
55npatand = lambda x: 180.*np.arctan(x)/np.pi
56npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
57npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
58npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
59   
60###############################################################################
61#### Particle form factors
62###############################################################################
63
64def SphereFF(Q,R,args=()):
65    ''' Compute hard sphere form factor - can use numpy arrays
66    param float Q: Q value array (usually in A-1)
67    param float R: sphere radius (Usually in A - must match Q-1 units)
68    param array args: ignored
69    returns float: form factors as array as needed
70    '''
71    QR = Q[:,np.newaxis]*R
72    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
73   
74def SpheroidFF(Q,R,args):
75    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
76    - can use numpy arrays for R & AR; will return corresponding numpy array
77    param float Q : Q value array (usually in A-1)
78    param float R: radius along 2 axes of spheroid
79    param array args: [float AR]: aspect ratio so 3rd axis = R*AR
80    returns float: form factors as array as needed
81    '''
82    NP = 50
83    AR = args[0]
84    if 0.99 < AR < 1.01:
85        return SphereFF(Q,R,0)
86    else:
87        cth = np.linspace(0,1.,NP)
88        Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
89        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
90           
91def CylinderFF(Q,R,args):
92    ''' Compute form factor for cylinders - can use numpy arrays
93    param float Q: Q value array (A-1)
94    param float R: cylinder radius (A)
95    param array args: [float L]: cylinder length (A)
96    returns float: form factor
97    '''
98    L = args[0]
99    NP = 200
100    alp = np.linspace(0,np.pi/2.,NP)
101    QL = Q[:,np.newaxis]*L
102    LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
103    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
104    QR = Q[:,np.newaxis]*R
105    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
106    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
107    LSBess = LBess*SBess
108    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
109   
110def CylinderDFF(Q,L,args):
111    ''' Compute form factor for cylinders - can use numpy arrays
112    param float Q: Q value array (A-1)
113    param float L: cylinder half length (A)
114    param array args: [float R]: cylinder radius (A)
115    returns float: form factor
116    '''
117    R = args[0]
118    return CylinderFF(Q,R,args=[2.*L])   
119   
120def CylinderARFF(Q,R,args): 
121    ''' Compute form factor for cylinders - can use numpy arrays
122    param float Q: Q value array (A-1)
123    param float R: cylinder radius (A)
124    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
125    returns float: form factor
126    '''
127    AR = args[0]
128    return CylinderFF(Q,R,args=[2.*R*AR])   
129   
130def UniSphereFF(Q,R,args=0):
131    ''' Compute form factor for unified sphere - can use numpy arrays
132    param float Q: Q value array (A-1)
133    param float R: cylinder radius (A)
134    param array args: ignored
135    returns float: form factor
136    '''
137    Rg = np.sqrt(3./5.)*R
138    B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
139    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
140    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
141   
142def UniRodFF(Q,R,args):
143    ''' Compute form factor for unified rod - can use numpy arrays
144    param float Q: Q value array (A-1)
145    param float R: cylinder radius (A)
146    param array args: [float R]: cylinder radius (A)
147    returns float: form factor
148    '''
149    L = args[0]
150    Rg2 = np.sqrt(R**2/2+L**2/12)
151    B2 = np.pi/L
152    Rg1 = np.sqrt(3.)*R/2.
153    G1 = (2./3.)*R/L
154    B1 = 4.*(L+R)/(R**3*L**2)
155    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
156    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
157    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
158    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
159    return np.sqrt(FF)
160   
161def UniRodARFF(Q,R,args):
162    ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
163    param float Q: Q value array (A-1)
164    param float R: cylinder radius (A)
165    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
166    returns float: form factor
167    '''
168    AR = args[0]
169    return UniRodFF(Q,R,args=[2.*AR*R,])
170   
171def UniDiskFF(Q,R,args):
172    ''' Compute form factor for unified disk - can use numpy arrays
173    param float Q: Q value array (A-1)
174    param float R: cylinder radius (A)
175    param array args: [float T]: disk thickness (A)
176    returns float: form factor
177    '''
178    T = args[0]
179    Rg2 = np.sqrt(R**2/2.+T**2/12.)
180    B2 = 2./R**2
181    Rg1 = np.sqrt(3.)*T/2.
182    RgC2 = 1.1*Rg1
183    G1 = (2./3.)*(T/R)**2
184    B1 = 4.*(T+R)/(R**3*T**2)
185    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
186    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
187    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
188    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
189    return np.sqrt(FF)
190   
191def UniTubeFF(Q,R,args):
192    ''' Compute form factor for unified tube - can use numpy arrays
193    assumes that core of tube is same as the matrix/solvent so contrast
194    is from tube wall vs matrix
195    param float Q: Q value array (A-1)
196    param float R: cylinder radius (A)
197    param array args: [float L,T]: tube length & wall thickness(A)
198    returns float: form factor
199    '''
200    L,T = args[:2]
201    Ri = R-T
202    DR2 = R**2-Ri**2
203    Vt = np.pi*DR2*L
204    Rg3 = np.sqrt(DR2/2.+L**2/12.)
205    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
206    B2 = np.pi**2*T/Vt
207    B3 = np.pi/L
208    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
209    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
210    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
211    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
212    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
213    FF += B1/QstV**4
214    return np.sqrt(FF)
215
216###############################################################################
217#### Particle volumes
218###############################################################################
219
220def SphereVol(R,args=()):
221    ''' Compute volume of sphere
222    - numpy array friendly
223    param float R: sphere radius
224    param array args: ignored
225    returns float: volume
226    '''
227    return (4./3.)*np.pi*R**3
228
229def SpheroidVol(R,args):
230    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
231    - numpy array friendly
232    param float R: radius along 2 axes of spheroid
233    param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
234    returns float: volume
235    '''
236    AR = args[0]
237    return AR*SphereVol(R)
238   
239def CylinderVol(R,args):
240    ''' Compute cylinder volume for radius & length
241    - numpy array friendly
242    param float R: diameter (A)
243    param array args: [float L]: length (A)
244    returns float:volume (A^3)
245    '''
246    L = args[0]
247    return np.pi*L*R**2
248   
249def CylinderDVol(L,args):
250    ''' Compute cylinder volume for length & diameter
251    - numpy array friendly
252    param float: L half length (A)
253    param array args: [float D]: diameter (A)
254    returns float:volume (A^3)
255    '''
256    D = args[0]
257    return CylinderVol(D/2.,[2.*L,])
258   
259def CylinderARVol(R,args):
260    ''' Compute cylinder volume for radius & aspect ratio = L/D
261    - numpy array friendly
262    param float: R radius (A)
263    param array args: [float AR]: =L/D=L/2R aspect ratio
264    returns float:volume
265    '''
266    AR = args[0]
267    return CylinderVol(R,[2.*R*AR,])
268   
269def UniSphereVol(R,args=()):
270    ''' Compute volume of sphere
271    - numpy array friendly
272    param float R: sphere radius
273    param array args: ignored
274    returns float: volume
275    '''
276    return SphereVol(R)
277   
278def UniRodVol(R,args):
279    ''' Compute cylinder volume for radius & length
280    - numpy array friendly
281    param float R: diameter (A)
282    param array args: [float L]: length (A)
283    returns float:volume (A^3)
284    '''
285    L = args[0]
286    return CylinderVol(R,[L,])
287   
288def UniRodARVol(R,args):
289    ''' Compute rod volume for radius & aspect ratio
290    - numpy array friendly
291    param float R: diameter (A)
292    param array args: [float AR]: =L/D=L/2R aspect ratio
293    returns float:volume (A^3)
294    '''
295    AR = args[0]
296    return CylinderARVol(R,[AR,])
297   
298def UniDiskVol(R,args):
299    ''' Compute disk volume for radius & thickness
300    - numpy array friendly
301    param float R: diameter (A)
302    param array args: [float T]: thickness
303    returns float:volume (A^3)
304    '''
305    T = args[0]
306    return CylinderVol(R,[T,])
307   
308def UniTubeVol(R,args):
309    ''' Compute tube volume for radius, length & wall thickness
310    - numpy array friendly
311    param float R: diameter (A)
312    param array args: [float L,T]: tube length & wall thickness(A)
313    returns float: volume (A^3) of tube wall
314    '''
315    L,T = args[:2]
316    return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
317   
318################################################################################
319#### Distribution functions & their cumulative fxns
320################################################################################
321
322def LogNormalDist(x,pos,args):
323    ''' Standard LogNormal distribution - numpy friendly on x axis
324    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
325    param float x: independent axis (can be numpy array)
326    param float pos: location of distribution
327    param float scale: width of distribution (m)
328    param float shape: shape - (sigma of log(LogNormal))
329    returns float: LogNormal distribution
330    '''
331    scale,shape = args
332    return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
333   
334def GaussDist(x,pos,args):
335    ''' Standard Normal distribution - numpy friendly on x axis
336    param float x: independent axis (can be numpy array)
337    param float pos: location of distribution
338    param float scale: width of distribution (sigma)
339    param float shape: not used
340    returns float: Normal distribution
341    '''
342    scale = args[0]
343    return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
344   
345def LSWDist(x,pos,args=[]):
346    ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
347    ref:
348    param float x: independent axis (can be numpy array)
349    param float pos: location of distribution
350    param float scale: not used
351    param float shape: not used
352    returns float: LSW distribution   
353    '''
354    redX = x/pos
355    result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
356   
357    return np.nan_to_num(result/pos)
358   
359def SchulzZimmDist(x,pos,args):
360    ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
361    ref: http://goldbook.iupac.org/S05502.html
362    param float x: independent axis (can be numpy array)
363    param float pos: location of distribution
364    param float scale: width of distribution (sigma)
365    param float shape: not used
366    returns float: Schulz-Zimm distribution
367    '''
368    scale = args[0]
369    b = (2.*pos/scale)**2
370    a = b/pos
371    if b<70:    #why bother?
372        return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
373    else:
374        return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
375           
376def LogNormalCume(x,pos,args):
377    ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
378    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
379    param float x: independent axis (can be numpy array)
380    param float pos: location of distribution
381    param float scale: width of distribution (sigma)
382    param float shape: shape parameter
383    returns float: LogNormal cumulative distribution
384    '''
385    scale,shape = args
386    lredX = np.log((x-pos)/scale)
387    return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
388   
389def GaussCume(x,pos,args):
390    ''' Standard Normal cumulative distribution - numpy friendly on x axis
391    param float x: independent axis (can be numpy array)
392    param float pos: location of distribution
393    param float scale: width of distribution (sigma)
394    param float shape: not used
395    returns float: Normal cumulative distribution
396    '''
397    scale = args[0]
398    redX = (x-pos)/scale
399    return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
400   
401def LSWCume(x,pos,args=[]):
402    ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
403    param float x: independent axis (can be numpy array)
404    param float pos: location of distribution
405    param float scale: not used
406    param float shape: not used
407    returns float: LSW cumulative distribution
408    '''
409    nP = 500
410    xMin,xMax = [0.,2.*pos]
411    X = np.linspace(xMin,xMax,nP,True)
412    fxn = LSWDist(X,pos)
413    mat = np.outer(np.ones(nP),fxn)
414    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
415    return np.interp(x,X,cume,0,1)
416   
417def SchulzZimmCume(x,pos,args):
418    ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
419    param float x: independent axis (can be numpy array)
420    param float pos: location of distribution
421    param float scale: width of distribution (sigma)
422    param float shape: not used
423    returns float: Normal distribution
424    '''
425    scale = args[0]
426    nP = 500
427    xMin = np.max([0.,pos-4.*scale])
428    xMax = np.min([pos+4.*scale,1.e5])
429    X = np.linspace(xMin,xMax,nP,True)
430    fxn = LSWDist(X,pos)
431    mat = np.outer(np.ones(nP),fxn)
432    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
433    return np.interp(x,X,cume,0,1)
434   
435    return []
436   
437################################################################################
438#### Structure factors for condensed systems
439################################################################################
440
441def DiluteSF(Q,args=[]):
442    ''' Default: no structure factor correction for dilute system
443    '''
444    return np.ones_like(Q)  #or return 1.0
445
446def HardSpheresSF(Q,args):
447    ''' Computes structure factor for not dilute monodisperse hard spheres
448    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
449    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
450   
451    param float Q: Q value array (A-1)
452    param array args: [float R, float VolFrac]: interparticle distance & volume fraction
453    returns numpy array S(Q)
454    '''
455   
456    R,VolFr = args
457    denom = (1.-VolFr)**4
458    num = (1.+2.*VolFr)**2
459    alp = num/denom
460    bet = -6.*VolFr*(1.+VolFr/2.)**2/denom
461    gamm = 0.5*VolFr*num/denom       
462    A = 2.*Q*R
463    A2 = A**2
464    A3 = A**3
465    A4 = A**4
466    Rca = np.cos(A)
467    Rsa = np.sin(A)
468    calp = alp*(Rsa/A2-Rca/A)
469    cbet = bet*(2.*Rsa/A2-(A2-2.)*Rca/A3-2./A3)
470    cgam = gamm*(-Rca/A+(4./A)*((3.*A2-6.)*Rca/A4+(A2-6.)*Rsa/A3+6./A4))
471    pfac = -24.*VolFr/A
472    C = pfac*(calp+cbet+cgam)
473    return 1./(1.-C)
474       
475def SquareWellSF(Q,args):
476    '''Computes structure factor for not dilute monodisperse hard sphere with a
477    square well potential interaction.
478    Refs.: SHARMA,SHARMA, PHYSICA 89A,(1977),213-
479   
480    :param float Q: Q value array (A-1)
481    :param array args: [float R, float VolFrac, float depth, float width]:
482        interparticle distance, volume fraction (<0.08), well depth (e/kT<1.5kT),
483        well width
484    :returns: numpy array S(Q)
485      well depth > 0 attractive & values outside above limits nonphysical cf.
486      Monte Carlo simulations
487    '''
488    R,VolFr,Depth,Width = args
489    eta = VolFr
490    eta2 = eta*eta
491    eta3 = eta*eta2
492    eta4 = eta*eta3       
493    etam1 = 1. - eta
494    etam14 = etam1**4
495    alp = (  ( (1. + 2.*eta)**2 ) + eta3*( eta-4.0 )  )/etam14
496    bet = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14
497    gam = 0.5*eta*( (1. + 2.*eta)**2 + eta3*(eta-4.) )/etam14
498
499    SK = 2.*Q*R
500    SK2 = SK*SK
501    SK3 = SK*SK2
502    SK4 = SK3*SK
503    T1 = alp * SK3 * ( np.sin(SK) - SK * np.cos(SK) )
504    T2 = bet * SK2 * ( 2.*SK*np.sin(SK) - (SK2-2.)*np.cos(SK) - 2.0 )
505    T3 =   ( 4.0*SK3 - 24.*SK ) * np.sin(SK) 
506    T3 = T3 - ( SK4 - 12.0*SK2 + 24.0 )*np.cos(SK) + 24.0   
507    T3 = gam*T3
508    T4 = -Depth*SK3*(np.sin(Width*SK) - Width*SK*np.cos(Width*SK)+ SK*np.cos(SK) - np.sin(SK) )
509    CK =  -24.0*eta*( T1 + T2 + T3 + T4 )/SK3/SK3
510    return 1./(1.-CK)
511   
512def StickyHardSpheresSF(Q,args):
513    ''' Computes structure factor for not dilute monodisperse hard spheres
514    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
515    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
516   
517    param float Q: Q value array (A-1)
518    param array args: [float R, float VolFrac]: sphere radius & volume fraction
519    returns numpy array S(Q)
520    '''
521    R,VolFr,epis,sticky = args
522    eta = VolFr/(1.0-epis)/(1.0-epis)/(1.0-epis)       
523    sig = 2.0 * R
524    aa = sig/(1.0 - epis)
525    etam1 = 1.0 - eta
526#  SOLVE QUADRATIC FOR LAMBDA
527    qa = eta/12.0
528    qb = -1.0*(sticky + eta/(etam1))
529    qc = (1.0 + eta/2.0)/(etam1*etam1)
530    radic = qb*qb - 4.0*qa*qc
531#   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
532    lam1 = (-1.0*qb-np.sqrt(radic))/(2.0*qa)
533    lam2 = (-1.0*qb+np.sqrt(radic))/(2.0*qa)
534    lam = min(lam1,lam2)
535    test = 1.0 + 2.0*eta
536    mu = lam*eta*etam1
537    alp = (1.0 + 2.0*eta - mu)/(etam1*etam1)
538    bet = (mu - 3.0*eta)/(2.0*etam1*etam1)
539#   CALCULATE THE STRUCTURE FACTOR<P></P>
540    kk = Q*aa
541    k2 = kk*kk
542    k3 = kk*k2
543    ds = np.sin(kk)
544    dc = np.cos(kk)
545    aq1 = ((ds - kk*dc)*alp)/k3
546    aq2 = (bet*(1.0-dc))/k2
547    aq3 = (lam*ds)/(12.0*kk)
548    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
549
550    bq1 = alp*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
551    bq2 = bet*(1.0/kk - ds/k2)
552    bq3 = (lam/12.0)*((1.0 - dc)/kk)
553    bq = 12.0*eta*(bq1+bq2-bq3)
554    sq = 1.0/(aq*aq +bq*bq)
555
556    return sq
557
558#def HayterPenfoldSF(Q,args): #big & ugly function - do later (if ever)
559#    pass
560   
561def InterPrecipitateSF(Q,args):
562    ''' Computes structure factor for precipitates in a matrix
563    Refs.: E-Wen Huang, Peter K. Liaw, Lionel Porcar, Yun Liu, Yee-Lang Liu,
564    Ji-Jung Kai, and Wei-Ren Chen,APPLIED PHYSICS LETTERS 93, 161904 (2008)
565    R. Giordano, A. Grasso, and J. Teixeira, Phys. Rev. A 43, 6894    1991   
566    param float Q: Q value array (A-1)
567    param array args: [float R, float VolFr]: "radius" & volume fraction
568    returns numpy array S(Q)
569    '''
570    R,VolFr = args
571    QV2 = Q**2*VolFr**2
572    top = 1 - np.exp(-QV2/4)*np.cos(2.*Q*R)
573    bot = 1-2*np.exp(-QV2/4)*np.cos(2.*Q*R) + np.exp(-QV2/2)
574    return 2*(top/bot) - 1
575         
576################################################################################
577##### SB-MaxEnt
578################################################################################
579
580def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
581    '''Calculates the response matrix :math:`G(Q,r)`
582   
583    :param float q: :math:`Q`
584    :param float r: :math:`r`
585    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
586    :param function FFfxn: form factor function FF(q,r,args)
587    :param function Volfxn: volume function Vol(r,args)
588    :returns float: G(Q,r)
589    '''
590    FF = FFfxn(q,r,args)
591    Vol = Volfxn(r,args)
592    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
593   
594'''
595sbmaxent
596
597Entropy maximization routine as described in the article
598J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
599("MNRAS": "Monthly Notices of the Royal Astronomical Society")
600
601:license: Copyright (c) 2013, UChicago Argonne, LLC
602:license: This file is distributed subject to a Software License Agreement found
603     in the file LICENSE that is included with this distribution.
604
605References:
606
6071. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
6082. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
609   Neutron Scattering Data Analysis, Rutherford
610   Appleton Laboratory, UK, 1986; ed. MW Johnson,
611   IOP Conference Series 81 (1986) 81 - 86, Institute
612   of Physics, Bristol, UK.
6133. ID Culverwell and GP Clarke; Ibid. 87 - 96.
6144. JA Potton, GK Daniell, & BD Rainford,
615   J APPL CRYST 21 (1988) 663 - 668.
6165. JA Potton, GJ Daniell, & BD Rainford,
617   J APPL CRYST 21 (1988) 891 - 897.
618
619'''
620
621class MaxEntException(Exception): 
622    '''Any exception from this module'''
623    pass
624
625def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
626    '''
627    do the complete Maximum Entropy algorithm of Skilling and Bryan
628   
629    :param float datum[]:
630    :param float sigma[]:
631    :param float[][] G: transformation matrix
632    :param float base[]:
633    :param int IterMax:
634    :param obj image_to_data: opus function (defaults to opus)
635    :param obj data_to_image: tropus function (defaults to tropus)
636   
637    :returns float[]: :math:`f(r) dr`
638    '''
639       
640    TEST_LIMIT        = 0.05                    # for convergence
641    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
642    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
643    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
644    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
645   
646    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
647    MOVE_PASSES       = 0.001                   # convergence test in routine: move
648
649    def tropus (data, G):
650        '''
651        tropus: transform data-space -> solution-space:  [G] * data
652       
653        default definition, caller can use this definition or provide an alternative
654       
655        :param float[M] data: observations, ndarray of shape (M)
656        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
657        :returns float[N]: calculated image, ndarray of shape (N)
658        '''
659        return G.dot(data)
660
661    def opus (image, G):
662        '''
663        opus: transform solution-space -> data-space:  [G]^tr * image
664       
665        default definition, caller can use this definition or provide an alternative
666       
667        :param float[N] image: solution, ndarray of shape (N)
668        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
669        :returns float[M]: calculated data, ndarray of shape (M)
670        '''
671        return np.dot(G.T,image)    #G.transpose().dot(image)
672
673    def Dist(s2, beta):
674        '''measure the distance of this possible solution'''
675        w = 0
676        n = beta.shape[0]
677        for k in range(n):
678            z = -sum(s2[k] * beta)
679            w += beta[k] * z
680        return w
681   
682    def ChiNow(ax, c1, c2, s1, s2):
683        '''
684        ChiNow
685       
686        :returns tuple: (ChiNow computation of ``w``, beta)
687        '''
688       
689        bx = 1 - ax
690        a =   bx * c2  -  ax * s2
691        b = -(bx * c1  -  ax * s1)
692   
693        beta = ChoSol(a, b)
694        w = 1.0
695        for k in range(SEARCH_DIRECTIONS):
696            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
697        return w, beta
698   
699    def ChoSol(a, b):
700        '''
701        ChoSol: ? chop the solution vectors ?
702       
703        :returns: new vector beta
704        '''
705        n = b.shape[0]
706        fl = np.zeros((n,n))
707        bl = np.zeros_like(b)
708       
709        #print_arr("ChoSol: a", a)
710        #print_vec("ChoSol: b", b)
711   
712        if (a[0][0] <= 0):
713            msg = "ChoSol: a[0][0] = " 
714            msg += str(a[0][0])
715            msg += '  Value must be positive'
716            raise MaxEntException(msg)
717   
718        # first, compute fl from a
719        # note fl is a lower triangular matrix
720        fl[0][0] = math.sqrt (a[0][0])
721        for i in (1, 2):
722            fl[i][0] = a[i][0] / fl[0][0]
723            for j in range(1, i+1):
724                z = 0.0
725                for k in range(j):
726                    z += fl[i][k] * fl[j][k]
727                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
728                z = a[i][j] - z
729                if j == i:
730                    y = math.sqrt(z)
731                else:
732                    y = z / fl[j][j]
733                fl[i][j] = y
734        #print_arr("ChoSol: fl", fl)
735   
736        # next, compute bl from fl and b
737        bl[0] = b[0] / fl[0][0]
738        for i in (1, 2):
739            z = 0.0
740            for k in range(i):
741                z += fl[i][k] * bl[k]
742                #print "\t", i, k, z
743            bl[i] = (b[i] - z) / fl[i][i]
744        #print_vec("ChoSol: bl", bl)
745   
746        # last, compute beta from bl and fl
747        beta = np.empty((n))
748        beta[-1] = bl[-1] / fl[-1][-1]
749        for i in (1, 0):
750            z = 0.0
751            for k in range(i+1, n):
752                z += fl[k][i] * beta[k]
753                #print "\t\t", i, k, 'z=', z
754            beta[i] = (bl[i] - z) / fl[i][i]
755        #print_vec("ChoSol: beta", beta)
756   
757        return beta
758
759    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
760        '''
761        move beta one step closer towards the solution
762        '''
763        a_lower, a_upper = 0., 1.          # bracket  "a"
764        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
765        #print "MaxEntMove: cmin = %g" % cmin
766        if cmin*chisq > chizer:
767            ctarg = (1.0 + cmin)/2
768        else:
769            ctarg = chizer/chisq
770        f_lower = cmin - ctarg
771        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
772        f_upper = c_upper - ctarg
773   
774        fx = 2*MOVE_PASSES      # just to start off
775        loop = 1
776        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
777            a_new = (a_lower + a_upper) * 0.5           # search by bisection
778            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
779            fx = c_new - ctarg
780            # tighten the search range for the next pass
781            if f_lower*fx > 0:
782                a_lower, f_lower = a_new, fx
783            if f_upper*fx > 0:
784                a_upper, f_upper = a_new, fx
785            loop += 1
786   
787        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
788            msg = "MaxEntMove: Loop counter = " 
789            msg += str(MAX_MOVE_LOOPS)
790            msg += '  No convergence in alpha chop'
791            raise MaxEntException(msg)
792   
793        w = Dist (s2, beta);
794        m = SEARCH_DIRECTIONS
795        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
796            for k in range(m):
797                beta[k] *= math.sqrt (fSum/(blank*w))
798        chtarg = ctarg * chisq
799        return w, chtarg, loop, a_new, fx, beta
800       
801#MaxEnt_SB starts here   
802
803    if image_to_data == None:
804        image_to_data = opus
805    if data_to_image == None:
806        data_to_image = tropus
807    n   = len(base)
808    npt = len(datum)
809
810    # Note that the order of subscripts for
811    # "xi" and "eta" has been reversed from
812    # the convention used in the FORTRAN version
813    # to enable parts of them to be passed as
814    # as vectors to "image_to_data" and "data_to_image".
815    xi      = np.zeros((SEARCH_DIRECTIONS, n))
816    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
817    beta    = np.zeros((SEARCH_DIRECTIONS))
818    # s1      = np.zeros((SEARCH_DIRECTIONS))
819    # c1      = np.zeros((SEARCH_DIRECTIONS))
820    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
821    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
822
823    # TODO: replace blank (scalar) with base (vector)
824    blank = sum(base) / len(base)   # use the average value of base
825
826    chizer, chtarg = npt*1.0, npt*1.0
827    f = base * 1.0                              # starting distribution is base
828
829    fSum  = sum(f)                              # find the sum of the f-vector
830    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
831    chisq = sum(z*z)                            # Chi^2, SB eq. 4
832
833    for iter in range(IterMax):
834        ox = -2 * z / sigma                        # gradient of Chi^2
835
836        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
837        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
838        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
839        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
840        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
841
842        a = 1.0
843        b = 1.0 / cnorm
844        if iter == 0:
845            test = 0.0     # mismatch between entropy and ChiSquared gradients
846        else:
847            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
848            a = 0.5 / (snorm * test)
849            b *= 0.5 / test
850        xi[0] = f * cgrad / cnorm
851        xi[1] = f * (a * sgrad - b * cgrad)
852
853        eta[0] = image_to_data (xi[0], G);          # image --> data
854        eta[1] = image_to_data (xi[1], G);          # image --> data
855        ox = eta[1] / (sigma * sigma)
856        xi[2] = data_to_image (ox, G);              # data --> image
857        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
858        xi[2] = f * xi[2] * a
859        eta[2] = image_to_data (xi[2], G)           # image --> data
860       
861#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
862#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
863
864        # prepare the search directions for the conjugate gradient technique
865        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
866        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
867#         print_vec("MaxEnt: c1", c1)
868#         print_vec("MaxEnt: s1", s1)
869
870        for k in range(SEARCH_DIRECTIONS):
871            for l in range(k+1):
872                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
873                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
874#         print_arr("MaxEnt: c2", c2)
875#         print_arr("MaxEnt: s2", s2)
876
877        # reflect across the body diagonal
878        for k, l in ((0,1), (0,2), (1,2)):
879            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
880            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
881 
882        beta[0] = -0.5 * c1[0] / c2[0][0]
883        beta[1] = 0.0
884        beta[2] = 0.0
885        if (iter > 0):
886            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
887
888        f_old = f.copy()    # preserve the last image
889        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
890       
891        # As mentioned at the top of p.119,
892        # need to protect against stray negative values.
893        # In this case, set them to RESET_STRAYS * base[i]
894        #f = f.clip(RESET_STRAYS * blank, f.max())
895        for i in range(n):
896            if f[i] <= 0.0:
897                f[i] = RESET_STRAYS * base[i]
898        df = f - f_old
899        fSum = sum(f)
900        fChange = sum(df)
901
902        # calculate the normalized entropy
903        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
904        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
905        chisq = sum(z*z)                            # report this ChiSq
906
907        if report:
908            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
909            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
910            if iter > 0:
911                value = 100*( math.sqrt(chisq/chtarg)-1)
912            else:
913                value = 0
914      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
915            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
916
917        # See if we have finished our task.
918        # do the hardest test first
919        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
920            print ' Convergence achieved.'
921            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
922    print ' No convergence! Try increasing Error multiplier.'
923    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
924
925   
926###############################################################################
927#### IPG/TNNLS Routines
928###############################################################################
929
930def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
931    ''' An implementation of the Interior-Point Gradient method of
932    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
933    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
934    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
935    Problem addressed: Total Non-Negative Least Squares (TNNLS)
936    :param float datum[]:
937    :param float sigma[]:
938    :param float[][] G: transformation matrix
939    :param int IterMax:
940    :param float Qvec: data positions for Power = 0-4
941    :param float approach: 0.8 default fitting parameter
942    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
943   
944    '''
945    if Power < 0: 
946        GmatE = G/sigma[:np.newaxis]
947        IntE = datum/sigma
948        pwr = 0
949        QvecP = np.ones_like(datum)
950    else:
951        GmatE = G[:]
952        IntE = datum[:]
953        pwr = Power
954        QvecP = Qvec**pwr
955    Amat = GmatE*QvecP[:np.newaxis]
956    AAmat = np.inner(Amat,Amat)
957    Bvec = IntE*QvecP
958    Xw = np.ones_like(Bins)*1.e-32
959    calc = np.dot(G.T,Xw)
960    nIter = 0
961    err = 10.
962    while (nIter<IterMax) and (err > 1.):
963        #Step 1 in M&Z paper:
964        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
965        Dk = Xw/np.inner(AAmat,Xw)
966        Pk = -Dk*Qk
967        #Step 2 in M&Z paper:
968        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
969        alpWv = -Xw/Pk
970        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
971        #Step 3 in M&Z paper:
972        shift = AkSt*Pk
973        Xw += shift
974        #done IPG; now results
975        nIter += 1
976        calc = np.dot(G.T,Xw)
977        chisq = np.sum(((datum-calc)/sigma)**2)
978        err = chisq/len(datum)
979        if report:
980            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
981    return chisq,Xw,calc
982
983###############################################################################
984#### SASD Utilities
985###############################################################################
986
987def SetScale(Data,refData):
988    Profile,Limits,Sample = Data
989    refProfile,refLimits,refSample = refData
990    x,y = Profile[:2]
991    rx,ry = refProfile[:2]
992    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
993    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
994    iBeg = np.searchsorted(x,Beg)
995    iFin = np.searchsorted(x,Fin)+1        #include last point
996    sum = np.sum(y[iBeg:iFin])
997    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
998    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
999   
1000def Bestimate(G,Rg,P):
1001    return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
1002   
1003def SmearData(Ic,Q,slitLen,Back):
1004    Np = Q.shape[0]
1005    Qtemp = np.concatenate([Q,Q[-1]+20*Q])
1006    Ictemp = np.concatenate([Ic,Ic[-1]*(1-(Qtemp[Np:]-Qtemp[Np])/(20*Qtemp[Np-1]))])
1007    Icsm = np.zeros_like(Q)
1008    Qsm = 2*slitLen*(np.interp(np.arange(2*Np)/2.,np.arange(Np),Q)-Q[0])/(Q[-1]-Q[0])
1009    Sp = np.searchsorted(Qsm,slitLen)
1010    DQsm = np.diff(Qsm)[:Sp]
1011    Ism = np.interp(np.sqrt(Q[:,np.newaxis]**2+Qsm**2),Qtemp,Ictemp)
1012    Icsm = np.sum((Ism[:,:Sp]*DQsm),axis=1)
1013    Icsm /= slitLen
1014    return Icsm
1015   
1016###############################################################################
1017#### Size distribution
1018###############################################################################
1019
1020def SizeDistribution(Profile,ProfDict,Limits,Sample,data):
1021    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1022        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1023        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1024        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1025        'Cylinder diam':[CylinderDFF,CylinderDVol]}
1026    Shape = data['Size']['Shape'][0]
1027    SlitLen = Sample.get('SlitLen',0.0)
1028    Parms = data['Size']['Shape'][1:]
1029    if data['Size']['logBins']:
1030        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
1031            data['Size']['Nbins']+1,True)/2.        #make radii
1032    else:
1033        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
1034            data['Size']['Nbins']+1,True)/2.        #make radii
1035    Dbins = np.diff(Bins)
1036    Bins = Bins[:-1]+Dbins/2.
1037    Contrast = Sample['Contrast'][1]
1038    Scale = Sample['Scale'][0]
1039    Sky = 10**data['Size']['MaxEnt']['Sky']
1040    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
1041    Back = data['Back']
1042    Q,Io,wt,Ic,Ib = Profile[:5]
1043    Qmin = Limits[1][0]
1044    Qmax = Limits[1][1]
1045    wtFactor = ProfDict['wtFactor']
1046    Ibeg = np.searchsorted(Q,Qmin)
1047    Ifin = np.searchsorted(Q,Qmax)+1        #include last point
1048    BinMag = np.zeros_like(Bins)
1049    Ic[:] = 0.
1050    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
1051    if 'MaxEnt' == data['Size']['Method']:
1052        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
1053            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
1054            data['Size']['MaxEnt']['Niter'],report=True)
1055    elif 'IPG' == data['Size']['Method']:
1056        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
1057            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
1058            Power=data['Size']['IPG']['Power'],report=True)
1059    Ib[:] = Back[0]
1060    Ic[Ibeg:Ifin] += Back[0]
1061    print ' Final chi^2: %.3f'%(chisq)
1062    Vols = shapes[Shape][1](Bins,Parms)
1063    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
1064       
1065################################################################################
1066#### Modelling
1067################################################################################
1068
1069def ModelFit(Profile,ProfDict,Limits,Sample,Model):
1070    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1071        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1072        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1073        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1074        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
1075           
1076    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
1077            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
1078               
1079    parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
1080        'PkInt','PkPos','PkSig','PkGam',]
1081       
1082    FFparmOrder = ['Aspect ratio','Length','Diameter','Thickness']
1083   
1084    SFparmOrder = ['Dist','VolFr','epis','Sticky','Depth','Width']
1085
1086    def GetModelParms():
1087        parmDict = {'Scale':Sample['Scale'][0],'SlitLen':Sample.get('SlitLen',0.0),}
1088        varyList = []
1089        values = []
1090        levelTypes = []
1091        Back = Model['Back']
1092        if Back[1]:
1093            varyList += ['Back',]
1094            values.append(Back[0])
1095        parmDict['Back'] = Back[0]
1096        partData = Model['Particle']
1097        for i,level in enumerate(partData['Levels']):
1098            cid = str(i)+';'
1099            controls = level['Controls']
1100            Type = controls['DistType']
1101            levelTypes.append(Type)
1102            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1103                if Type not in ['Monodosperse',]:
1104                    parmDict[cid+'NumPoints'] = controls['NumPoints']
1105                    parmDict[cid+'Cutoff'] = controls['Cutoff']
1106                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
1107                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
1108                parmDict[cid+'StrFact'] = sfxns[controls['StrFact']]
1109                parmDict[cid+'Contrast'] = controls['Contrast']
1110                for item in FFparmOrder:
1111                    if item in controls['FFargs']:
1112                        parmDict[cid+item] = controls['FFargs'][item][0]
1113                        if controls['FFargs'][item][1]:
1114                            varyList.append(cid+item)
1115                            values.append(controls['FFargs'][item][0])
1116                for item in SFparmOrder:
1117                    if item in controls.get('SFargs',{}):
1118                        parmDict[cid+item] = controls['SFargs'][item][0]
1119                        if controls['SFargs'][item][1]:
1120                            varyList.append(cid+item)
1121                            values.append(controls['SFargs'][item][0])
1122            distDict = controls['DistType']
1123            for item in parmOrder:
1124                if item in level[distDict]:
1125                    parmDict[cid+item] = level[distDict][item][0]
1126                    if level[distDict][item][1]:
1127                        values.append(level[distDict][item][0])
1128                        varyList.append(cid+item)
1129        return levelTypes,parmDict,varyList,values
1130       
1131    def SetModelParms():
1132        print ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1133        if 'Back' in varyList:
1134            Model['Back'][0] = parmDict['Back']
1135            print '  %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
1136        partData = Model['Particle']
1137        for i,level in enumerate(partData['Levels']):
1138            controls = level['Controls']
1139            Type = controls['DistType']
1140            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1141                print ' Component %d: Type: %s: Structure Factor: %s Contrast: %12.3f'  \
1142                    %(i,Type,controls['StrFact'],controls['Contrast'])               
1143            else:
1144                print ' Component %d: Type: %s: '%(i,Type,)
1145            cid = str(i)+';'
1146            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1147                for item in FFparmOrder:
1148                    if cid+item in varyList:
1149                        controls['FFargs'][item][0] = parmDict[cid+item]
1150                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1151                for item in SFparmOrder:
1152                    if cid+item in varyList:
1153                        controls['SFargs'][item][0] = parmDict[cid+item]
1154                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1155            distDict = controls['DistType']
1156            for item in level[distDict]:
1157                if cid+item in varyList:
1158                    level[distDict][item][0] = parmDict[cid+item]
1159                    print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1160                   
1161    def calcSASD(values,Q,Io,wt,Ifb,levelTypes,parmDict,varyList):
1162        parmDict.update(zip(varyList,values))
1163        M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)+Ifb-parmDict['Scale']*Io)
1164        return M
1165       
1166    def getSASD(Q,levelTypes,parmDict):
1167        Ic = np.zeros_like(Q)
1168        for i,Type in enumerate(levelTypes):
1169            cid = str(i)+';'
1170            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1171                FFfxn = parmDict[cid+'FormFact']
1172                Volfxn = parmDict[cid+'FFVolume']
1173                SFfxn = parmDict[cid+'StrFact']
1174                FFargs = []
1175                SFargs = []
1176                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter']:
1177                    if item in parmDict: 
1178                        FFargs.append(parmDict[item])
1179                for item in [cid+'Dist',cid+'VolFr',cid+'epis',cid+'Sticky',cid+'Depth',cid+'Width']:
1180                    if item in parmDict: 
1181                        SFargs.append(parmDict[item])
1182                distDict = {}
1183                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
1184                    if item in parmDict:
1185                        distDict[item.split(';')[1]] = parmDict[item]
1186                contrast = parmDict[cid+'Contrast']
1187                rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
1188                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
1189                dist *= parmDict[cid+'Volume']
1190                Ic += np.dot(Gmat,dist)*SFfxn(Q,args=SFargs)
1191            elif 'Unified' in Type:
1192                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
1193                    parmDict[cid+'P'],parmDict[cid+'Cutoff']
1194                Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
1195                Guin = G*np.exp(-(Q*Rg)**2/3)
1196                Porod = (B/Qstar**P)*np.exp(-(Q*Rgco)**2/3)
1197                Ic += Guin+Porod
1198            elif 'Porod' in Type:
1199                B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
1200                Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
1201                Ic += Porod
1202            elif 'Mono' in Type:
1203                FFfxn = parmDict[cid+'FormFact']
1204                Volfxn = parmDict[cid+'FFVolume']
1205                SFfxn = parmDict[cid+'StrFact']
1206                FFargs = []
1207                SFargs = []
1208                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
1209                    if item in parmDict: 
1210                        FFargs.append(parmDict[item])
1211                for item in [cid+'Dist',cid+'VolFr',cid+'epis',cid+'Sticky',cid+'Depth',cid+'Width']:
1212                    if item in parmDict: 
1213                        SFargs.append(parmDict[item])
1214                contrast = parmDict[cid+'Contrast']
1215                R = parmDict[cid+'Radius']
1216                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
1217                Ic += Gmat[0]*parmDict[cid+'Volume']*SFfxn(Q,args=SFargs)
1218            elif 'Bragg' in Type:
1219                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
1220                    parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
1221        Ic += parmDict['Back']  #/parmDict['Scale']
1222        slitLen = Sample['SlitLen']
1223        if slitLen:
1224            Ic = SmearData(Ic,Q,slitLen,parmDict['Back'])
1225        return Ic
1226       
1227    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1228    Qmin = Limits[1][0]
1229    Qmax = Limits[1][1]
1230    wtFactor = ProfDict['wtFactor']
1231    Ibeg = np.searchsorted(Q,Qmin)
1232    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
1233    Ic[:] = 0
1234    levelTypes,parmDict,varyList,values = GetModelParms()
1235    if varyList:
1236        result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
1237            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList))
1238        ncyc = int(result[2]['nfev']/2)
1239        parmDict.update(zip(varyList,result[0]))
1240        chisq = np.sum(result[2]['fvec']**2)
1241        ncalc = result[2]['nfev']
1242        covM = result[1]
1243    else:   #nothing varied
1244        M = calcSASD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList)
1245        chisq = np.sum(M**2)
1246        ncalc = 0
1247        covM = []
1248        sig = []
1249        sigDict = {}
1250        result = []
1251    Rvals = {}
1252    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
1253    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
1254    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
1255    Msg = 'Failed to converge'
1256    try:
1257        Nans = np.isnan(result[0])
1258        if np.any(Nans):
1259            name = varyList[Nans.nonzero(True)[0]]
1260            Msg = 'Nan result for '+name+'!'
1261            raise ValueError
1262        Negs = np.less_equal(result[0],0.)
1263        if np.any(Negs):
1264            name = varyList[Negs.nonzero(True)[0]]
1265            Msg = 'negative coefficient for '+name+'!'
1266            raise ValueError
1267        if len(covM):
1268            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
1269            sigDict = dict(zip(varyList,sig))
1270        print ' Results of small angle data modelling fit:'
1271        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
1272        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1273        SetModelParms()
1274        covMatrix = covM*Rvals['GOF']
1275        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
1276    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
1277        return False,0,0,0,0,0,0,Msg
1278   
1279def ModelFxn(Profile,ProfDict,Limits,Sample,sasdData):
1280   
1281    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1282        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1283        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1284        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1285        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
1286    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
1287            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
1288
1289#    pdb.set_trace()
1290    partData = sasdData['Particle']
1291    matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
1292    Scale = Sample['Scale']     #[value,flag]
1293    SlitLen = Sample.get('SlitLen',0.0)
1294    Back = sasdData['Back']
1295    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1296    Qmin = Limits[1][0]
1297    Qmax = Limits[1][1]
1298    wtFactor = ProfDict['wtFactor']
1299    Ibeg = np.searchsorted(Q,Qmin)
1300    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
1301    Ib[:] = Back[0]
1302    Ic[:] = 0
1303    Rbins = []
1304    Dist = []
1305    for level in partData['Levels']:
1306        controls = level['Controls']
1307        distFxn = controls['DistType']
1308        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1309            parmDict = level[controls['DistType']]
1310            FFfxn = shapes[controls['FormFact']][0]
1311            Volfxn = shapes[controls['FormFact']][1]
1312            SFfxn = sfxns[controls['StrFact']]
1313            FFargs = []
1314            SFargs = []
1315            for item in ['Dist','VolFr','epis','Sticky','Depth','Width',]:
1316                if item in controls.get('SFargs',{}):
1317                    SFargs.append(controls['SFargs'][item][0])
1318            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1319                if item in controls['FFargs']: 
1320                    FFargs.append(controls['FFargs'][item][0])
1321            contrast = controls['Contrast']
1322            distDict = {}
1323            for item in parmDict:
1324                distDict[item] = parmDict[item][0]
1325            rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
1326            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
1327            dist *= level[distFxn]['Volume'][0]
1328            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)*SFfxn(Q[Ibeg:Ifin],args=SFargs)
1329            Rbins.append(rBins)
1330            Dist.append(dist/(4.*dBins))
1331        elif 'Unified' in distFxn:
1332            parmDict = level[controls['DistType']]
1333            Rg,G,B,P,Rgco = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],    \
1334                parmDict['P'][0],parmDict['Cutoff'][0]
1335            Qstar = Q[Ibeg:Ifin]/(scsp.erf(Q[Ibeg:Ifin]*Rg/np.sqrt(6)))**3
1336            Guin = G*np.exp(-(Q[Ibeg:Ifin]*Rg)**2/3)
1337            Porod = (B/Qstar**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1338            Ic[Ibeg:Ifin] += Guin+Porod
1339            Rbins.append([])
1340            Dist.append([])
1341        elif 'Porod' in distFxn:
1342            parmDict = level[controls['DistType']]
1343            B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
1344            Porod = (B/Q[Ibeg:Ifin]**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1345            Ic[Ibeg:Ifin] += Porod
1346            Rbins.append([])
1347            Dist.append([])
1348        elif 'Mono' in distFxn:
1349            parmDict = level[controls['DistType']]
1350            R = level[controls['DistType']]['Radius'][0]
1351            FFfxn = shapes[controls['FormFact']][0]
1352            Volfxn = shapes[controls['FormFact']][1]
1353            SFfxn = sfxns[controls['StrFact']]
1354            FFargs = []
1355            SFargs = []
1356            for item in ['Dist','VolFr','epis','Sticky','Depth','Width',]:
1357                if item in controls.get('SFargs',{}):
1358                    SFargs.append(controls['SFargs'][item][0])
1359            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1360                if item in controls['FFargs']: 
1361                    FFargs.append(controls['FFargs'][item][0])
1362            contrast = controls['Contrast']
1363            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
1364            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]*SFfxn(Q[Ibeg:Ifin],args=SFargs)
1365            Rbins.append([])
1366            Dist.append([])
1367        elif 'Bragg' in distFxn:
1368            parmDict = level[controls['DistType']]
1369            Ic[Ibeg:Ifin] += parmDict['PkInt'][0]*G2pwd.getPsVoigt(parmDict['PkPos'][0],
1370                parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
1371            Rbins.append([])
1372            Dist.append([])
1373    Ic[Ibeg:Ifin] += Back[0]
1374    slitLen = Sample.get('SlitLen',0.)
1375    if slitLen:
1376        Ic[Ibeg:Ifin] = SmearData(Ic,Q,slitLen,Back[0])[Ibeg:Ifin]
1377    sasdData['Size Calc'] = [Rbins,Dist]
1378   
1379def MakeDiamDist(DistName,nPoints,cutoff,distDict):
1380   
1381    if 'LogNormal' in DistName:
1382        distFxn = 'LogNormalDist'
1383        cumeFxn = 'LogNormalCume'
1384        pos = distDict['MinSize']
1385        args = [distDict['Mean'],distDict['StdDev']]
1386        step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))
1387        mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)
1388        minX = 1. #pos
1389        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
1390    elif 'Gauss' in DistName:
1391        distFxn = 'GaussDist'
1392        cumeFxn = 'GaussCume'
1393        pos = distDict['Mean']
1394        args = [distDict['StdDev']]
1395        step = 0.02*distDict['StdDev']
1396        mode = distDict['Mean']
1397        minX = np.max([mode-4.*distDict['StdDev'],1.])
1398        maxX = np.min([mode+4.*distDict['StdDev'],1.e5])
1399    elif 'LSW' in DistName:
1400        distFxn = 'LSWDist'
1401        cumeFxn = 'LSWCume'
1402        pos = distDict['Mean']
1403        args = []
1404        minX,maxX = [0.,2.*pos]
1405    elif 'Schulz' in DistName:
1406        distFxn = 'SchulzZimmDist'
1407        cumeFxn = 'SchulzZimmCume'
1408        pos = distDict['Mean']
1409        args = [distDict['StdDev']]
1410        minX = np.max([1.,pos-4.*distDict['StdDev']])
1411        maxX = np.min([pos+4.*distDict['StdDev'],1.e5])
1412    nP = 500
1413    Diam = np.logspace(0.,5.,nP,True)
1414    TCW = eval(cumeFxn+'(Diam,pos,args)')
1415    CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints+1,True)
1416    rBins = np.interp(CumeTgts,TCW,Diam,0,0)
1417    dBins = np.diff(rBins)
1418    rBins = rBins[:-1]+dBins/2.
1419    return rBins,dBins,eval(distFxn+'(rBins,pos,args)')
1420
1421################################################################################
1422#### MaxEnt testing stuff
1423################################################################################
1424
1425def print_vec(text, a):
1426    '''print the contents of a vector to the console'''
1427    n = a.shape[0]
1428    print "%s[ = (" % text,
1429    for i in range(n):
1430        s = " %g, " % a[i]
1431        print s,
1432    print ")"
1433
1434def print_arr(text, a):
1435    '''print the contents of an array to the console'''
1436    n, m = a.shape
1437    print "%s[][] = (" % text
1438    for i in range(n):
1439        print " (",
1440        for j in range(m):
1441            print " %g, " % a[i][j],
1442        print "),"
1443    print ")"
1444
1445def test_MaxEnt_SB(report=True):
1446    def readTextData(filename):
1447        '''return q, I, dI from a 3-column text file'''
1448        if not os.path.exists(filename):
1449            raise Exception("file not found: " + filename)
1450        buf = [line.split() for line in open(filename, 'r').readlines()]
1451        M = len(buf)
1452        buf = zip(*buf)         # transpose rows and columns
1453        q  = np.array(buf[0], dtype=np.float64)
1454        I  = np.array(buf[1], dtype=np.float64)
1455        dI = np.array(buf[2], dtype=np.float64)
1456        return q, I, dI
1457    print "MaxEnt_SB: "
1458    test_data_file = os.path.join( 'testinp', 'test.sas')
1459    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
1460    bkg   = 0.1     #   I = I - bkg
1461    dMin, dMax, nRadii = 25, 9000, 40
1462    defaultDistLevel = 1.0e-6
1463    IterMax = 40
1464    errFac = 1.05
1465   
1466    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
1467    dr   = r * (r[1]/r[0] - 1)          # step size
1468    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
1469    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
1470   
1471    qVec, I, dI = readTextData(test_data_file)
1472    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
1473   
1474    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
1475    if f_dr is None:
1476        print "no solution"
1477        return
1478   
1479    print "solution reached"
1480    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
1481        print '%10.4f %10.4f %12.4g'%(a,b,c)
1482
1483def tests():
1484    test_MaxEnt_SB(report=True)
1485
1486if __name__ == '__main__':
1487    tests()
1488
Note: See TracBrowser for help on using the repository browser.