source: trunk/GSASIIsasd.py @ 2802

Last change on this file since 2802 was 2802, checked in by toby, 6 years ago

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