source: trunk/GSASIIsasd.py @ 1309

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

Implement structure factors for non-dilute small angle systems
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 56.9 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: 2014-04-30 19:05:57 +0000 (Wed, 30 Apr 2014) $
10# $Author: vondreele $
11# $Revision: 1309 $
12# $URL: trunk/GSASIIsasd.py $
13# $Id: GSASIIsasd.py 1309 2014-04-30 19:05:57Z 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: 1309 $")
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),212
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,tau = args
522#//      Input (fitting) variables are:
523#       //radius = w[0]
524#       //volume fraction = w[1]
525#       //epsilon (perurbation param) = w[2]
526#       //tau (stickiness) = w[3]
527#       Variable rad,phi,eps,tau,eta
528#       Variable sig,aa,etam1,qa,qb,qc,radic
529#       Variable lam,lam2,test,mu,alpha,BetaVar
530#       Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq
531#       rad = w[0]
532#       phi = w[1]
533#       eps = w[2]
534#       tau = w[3]
535#       
536#       eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps)
537#       
538#       sig = 2.0 * rad
539#       aa = sig/(1.0 - eps)
540#       etam1 = 1.0 - eta
541#//C
542#//C  SOLVE QUADRATIC FOR LAMBDA
543#//C
544#       qa = eta/12.0
545#       qb = -1.0*(tau + eta/(etam1))
546#       qc = (1.0 + eta/2.0)/(etam1*etam1)
547#       radic = qb*qb - 4.0*qa*qc
548#       if(radic<0)
549#               if(x>0.01 && x<0.015)
550#                       Print "Lambda unphysical - both roots imaginary"
551#               endif
552#               return(-1)
553#       endif
554#//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
555#       lam = (-1.0*qb-sqrt(radic))/(2.0*qa)
556#       lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa)
557#       if(lam2<lam)
558#               lam = lam2
559#       endif
560#       test = 1.0 + 2.0*eta
561#       mu = lam*eta*etam1
562#       if(mu>test)
563#               if(x>0.01 && x<0.015)
564#                Print "Lambda unphysical mu>test"
565#               endif
566#               return(-1)
567#       endif
568#       alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1)
569#       BetaVar = (mu - 3.0*eta)/(2.0*etam1*etam1)
570#       
571#//C
572#//C   CALCULATE THE STRUCTURE FACTOR
573#//C
574#
575#       qv = x
576#       kk = qv*aa
577#       k2 = kk*kk
578#       k3 = kk*k2
579#       ds = sin(kk)
580#       dc = cos(kk)
581#       aq1 = ((ds - kk*dc)*alpha)/k3
582#       aq2 = (BetaVar*(1.0-dc))/k2
583#       aq3 = (lam*ds)/(12.0*kk)
584#       aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
585#//
586#       bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
587#       bq2 = BetaVar*(1.0/kk - ds/k2)
588#       bq3 = (lam/12.0)*((1.0 - dc)/kk)
589#       bq = 12.0*eta*(bq1+bq2-bq3)
590#//
591#       sq = 1.0/(aq*aq +bq*bq)
592#
593#       Return (sq)
594    pass
595   
596#def HayterPenfoldSF(Q,args): #big & ugly function - do later (if ever)
597#    pass
598   
599def InterPrecipitateSF(Q,args):
600    ''' Computes structure factor for precipitates in a matrix
601    Refs.: E-Wen Huang, Peter K. Liaw, Lionel Porcar, Yun Liu, Yee-Lang Liu,
602    Ji-Jung Kai, and Wei-Ren Chen,APPLIED PHYSICS LETTERS 93, 161904 (2008)
603    R. Giordano, A. Grasso, and J. Teixeira, Phys. Rev. A 43, 6894    1991   
604    param float Q: Q value array (A-1)
605    param array args: [float R, float VolFr]: "radius" & volume fraction
606    returns numpy array S(Q)
607    '''
608    R,VolFr = args
609    QV2 = Q**2*VolFr**2
610    top = 1 - np.exp(-QV2/4)*np.cos(2.*Q*R)
611    bot = 1-2*np.exp(-QV2/4)*np.cos(2.*Q*R) + np.exp(-QV2/2)
612    return 2*(top/bot) - 1
613         
614################################################################################
615##### SB-MaxEnt
616################################################################################
617
618def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
619    '''Calculates the response matrix :math:`G(Q,r)`
620   
621    :param float q: :math:`Q`
622    :param float r: :math:`r`
623    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
624    :param function FFfxn: form factor function FF(q,r,args)
625    :param function Volfxn: volume function Vol(r,args)
626    :returns float: G(Q,r)
627    '''
628    FF = FFfxn(q,r,args)
629    Vol = Volfxn(r,args)
630    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
631   
632'''
633sbmaxent
634
635Entropy maximization routine as described in the article
636J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
637("MNRAS": "Monthly Notices of the Royal Astronomical Society")
638
639:license: Copyright (c) 2013, UChicago Argonne, LLC
640:license: This file is distributed subject to a Software License Agreement found
641     in the file LICENSE that is included with this distribution.
642
643References:
644
6451. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
6462. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
647   Neutron Scattering Data Analysis, Rutherford
648   Appleton Laboratory, UK, 1986; ed. MW Johnson,
649   IOP Conference Series 81 (1986) 81 - 86, Institute
650   of Physics, Bristol, UK.
6513. ID Culverwell and GP Clarke; Ibid. 87 - 96.
6524. JA Potton, GK Daniell, & BD Rainford,
653   J APPL CRYST 21 (1988) 663 - 668.
6545. JA Potton, GJ Daniell, & BD Rainford,
655   J APPL CRYST 21 (1988) 891 - 897.
656
657'''
658
659class MaxEntException(Exception): 
660    '''Any exception from this module'''
661    pass
662
663def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
664    '''
665    do the complete Maximum Entropy algorithm of Skilling and Bryan
666   
667    :param float datum[]:
668    :param float sigma[]:
669    :param float[][] G: transformation matrix
670    :param float base[]:
671    :param int IterMax:
672    :param obj image_to_data: opus function (defaults to opus)
673    :param obj data_to_image: tropus function (defaults to tropus)
674   
675    :returns float[]: :math:`f(r) dr`
676    '''
677       
678    TEST_LIMIT        = 0.05                    # for convergence
679    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
680    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
681    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
682    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
683   
684    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
685    MOVE_PASSES       = 0.001                   # convergence test in routine: move
686
687    def tropus (data, G):
688        '''
689        tropus: transform data-space -> solution-space:  [G] * data
690       
691        default definition, caller can use this definition or provide an alternative
692       
693        :param float[M] data: observations, ndarray of shape (M)
694        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
695        :returns float[N]: calculated image, ndarray of shape (N)
696        '''
697        return G.dot(data)
698
699    def opus (image, G):
700        '''
701        opus: transform solution-space -> data-space:  [G]^tr * image
702       
703        default definition, caller can use this definition or provide an alternative
704       
705        :param float[N] image: solution, ndarray of shape (N)
706        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
707        :returns float[M]: calculated data, ndarray of shape (M)
708        '''
709        return np.dot(G.T,image)    #G.transpose().dot(image)
710
711    def Dist(s2, beta):
712        '''measure the distance of this possible solution'''
713        w = 0
714        n = beta.shape[0]
715        for k in range(n):
716            z = -sum(s2[k] * beta)
717            w += beta[k] * z
718        return w
719   
720    def ChiNow(ax, c1, c2, s1, s2):
721        '''
722        ChiNow
723       
724        :returns tuple: (ChiNow computation of ``w``, beta)
725        '''
726       
727        bx = 1 - ax
728        a =   bx * c2  -  ax * s2
729        b = -(bx * c1  -  ax * s1)
730   
731        beta = ChoSol(a, b)
732        w = 1.0
733        for k in range(SEARCH_DIRECTIONS):
734            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
735        return w, beta
736   
737    def ChoSol(a, b):
738        '''
739        ChoSol: ? chop the solution vectors ?
740       
741        :returns: new vector beta
742        '''
743        n = b.shape[0]
744        fl = np.zeros((n,n))
745        bl = np.zeros_like(b)
746       
747        #print_arr("ChoSol: a", a)
748        #print_vec("ChoSol: b", b)
749   
750        if (a[0][0] <= 0):
751            msg = "ChoSol: a[0][0] = " 
752            msg += str(a[0][0])
753            msg += '  Value must be positive'
754            raise MaxEntException(msg)
755   
756        # first, compute fl from a
757        # note fl is a lower triangular matrix
758        fl[0][0] = math.sqrt (a[0][0])
759        for i in (1, 2):
760            fl[i][0] = a[i][0] / fl[0][0]
761            for j in range(1, i+1):
762                z = 0.0
763                for k in range(j):
764                    z += fl[i][k] * fl[j][k]
765                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
766                z = a[i][j] - z
767                if j == i:
768                    y = math.sqrt(z)
769                else:
770                    y = z / fl[j][j]
771                fl[i][j] = y
772        #print_arr("ChoSol: fl", fl)
773   
774        # next, compute bl from fl and b
775        bl[0] = b[0] / fl[0][0]
776        for i in (1, 2):
777            z = 0.0
778            for k in range(i):
779                z += fl[i][k] * bl[k]
780                #print "\t", i, k, z
781            bl[i] = (b[i] - z) / fl[i][i]
782        #print_vec("ChoSol: bl", bl)
783   
784        # last, compute beta from bl and fl
785        beta = np.empty((n))
786        beta[-1] = bl[-1] / fl[-1][-1]
787        for i in (1, 0):
788            z = 0.0
789            for k in range(i+1, n):
790                z += fl[k][i] * beta[k]
791                #print "\t\t", i, k, 'z=', z
792            beta[i] = (bl[i] - z) / fl[i][i]
793        #print_vec("ChoSol: beta", beta)
794   
795        return beta
796
797    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
798        '''
799        move beta one step closer towards the solution
800        '''
801        a_lower, a_upper = 0., 1.          # bracket  "a"
802        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
803        #print "MaxEntMove: cmin = %g" % cmin
804        if cmin*chisq > chizer:
805            ctarg = (1.0 + cmin)/2
806        else:
807            ctarg = chizer/chisq
808        f_lower = cmin - ctarg
809        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
810        f_upper = c_upper - ctarg
811   
812        fx = 2*MOVE_PASSES      # just to start off
813        loop = 1
814        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
815            a_new = (a_lower + a_upper) * 0.5           # search by bisection
816            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
817            fx = c_new - ctarg
818            # tighten the search range for the next pass
819            if f_lower*fx > 0:
820                a_lower, f_lower = a_new, fx
821            if f_upper*fx > 0:
822                a_upper, f_upper = a_new, fx
823            loop += 1
824   
825        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
826            msg = "MaxEntMove: Loop counter = " 
827            msg += str(MAX_MOVE_LOOPS)
828            msg += '  No convergence in alpha chop'
829            raise MaxEntException(msg)
830   
831        w = Dist (s2, beta);
832        m = SEARCH_DIRECTIONS
833        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
834            for k in range(m):
835                beta[k] *= math.sqrt (fSum/(blank*w))
836        chtarg = ctarg * chisq
837        return w, chtarg, loop, a_new, fx, beta
838       
839#MaxEnt_SB starts here   
840
841    if image_to_data == None:
842        image_to_data = opus
843    if data_to_image == None:
844        data_to_image = tropus
845    n   = len(base)
846    npt = len(datum)
847
848    # Note that the order of subscripts for
849    # "xi" and "eta" has been reversed from
850    # the convention used in the FORTRAN version
851    # to enable parts of them to be passed as
852    # as vectors to "image_to_data" and "data_to_image".
853    xi      = np.zeros((SEARCH_DIRECTIONS, n))
854    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
855    beta    = np.zeros((SEARCH_DIRECTIONS))
856    # s1      = np.zeros((SEARCH_DIRECTIONS))
857    # c1      = np.zeros((SEARCH_DIRECTIONS))
858    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
859    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
860
861    # TODO: replace blank (scalar) with base (vector)
862    blank = sum(base) / len(base)   # use the average value of base
863
864    chizer, chtarg = npt*1.0, npt*1.0
865    f = base * 1.0                              # starting distribution is base
866
867    fSum  = sum(f)                              # find the sum of the f-vector
868    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
869    chisq = sum(z*z)                            # Chi^2, SB eq. 4
870
871    for iter in range(IterMax):
872        ox = -2 * z / sigma                        # gradient of Chi^2
873
874        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
875        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
876        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
877        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
878        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
879
880        a = 1.0
881        b = 1.0 / cnorm
882        if iter == 0:
883            test = 0.0     # mismatch between entropy and ChiSquared gradients
884        else:
885            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
886            a = 0.5 / (snorm * test)
887            b *= 0.5 / test
888        xi[0] = f * cgrad / cnorm
889        xi[1] = f * (a * sgrad - b * cgrad)
890
891        eta[0] = image_to_data (xi[0], G);          # image --> data
892        eta[1] = image_to_data (xi[1], G);          # image --> data
893        ox = eta[1] / (sigma * sigma)
894        xi[2] = data_to_image (ox, G);              # data --> image
895        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
896        xi[2] = f * xi[2] * a
897        eta[2] = image_to_data (xi[2], G)           # image --> data
898       
899#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
900#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
901
902        # prepare the search directions for the conjugate gradient technique
903        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
904        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
905#         print_vec("MaxEnt: c1", c1)
906#         print_vec("MaxEnt: s1", s1)
907
908        for k in range(SEARCH_DIRECTIONS):
909            for l in range(k+1):
910                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
911                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
912#         print_arr("MaxEnt: c2", c2)
913#         print_arr("MaxEnt: s2", s2)
914
915        # reflect across the body diagonal
916        for k, l in ((0,1), (0,2), (1,2)):
917            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
918            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
919 
920        beta[0] = -0.5 * c1[0] / c2[0][0]
921        beta[1] = 0.0
922        beta[2] = 0.0
923        if (iter > 0):
924            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
925
926        f_old = f.copy()    # preserve the last image
927        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
928       
929        # As mentioned at the top of p.119,
930        # need to protect against stray negative values.
931        # In this case, set them to RESET_STRAYS * base[i]
932        #f = f.clip(RESET_STRAYS * blank, f.max())
933        for i in range(n):
934            if f[i] <= 0.0:
935                f[i] = RESET_STRAYS * base[i]
936        df = f - f_old
937        fSum = sum(f)
938        fChange = sum(df)
939
940        # calculate the normalized entropy
941        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
942        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
943        chisq = sum(z*z)                            # report this ChiSq
944
945        if report:
946            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
947            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
948            if iter > 0:
949                value = 100*( math.sqrt(chisq/chtarg)-1)
950            else:
951                value = 0
952      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
953            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
954
955        # See if we have finished our task.
956        # do the hardest test first
957        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
958            print ' Convergence achieved.'
959            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
960    print ' No convergence! Try increasing Error multiplier.'
961    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
962
963   
964###############################################################################
965#### IPG/TNNLS Routines
966###############################################################################
967
968def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
969    ''' An implementation of the Interior-Point Gradient method of
970    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
971    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
972    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
973    Problem addressed: Total Non-Negative Least Squares (TNNLS)
974    :param float datum[]:
975    :param float sigma[]:
976    :param float[][] G: transformation matrix
977    :param int IterMax:
978    :param float Qvec: data positions for Power = 0-4
979    :param float approach: 0.8 default fitting parameter
980    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
981   
982    '''
983    if Power < 0: 
984        GmatE = G/sigma[:np.newaxis]
985        IntE = datum/sigma
986        pwr = 0
987        QvecP = np.ones_like(datum)
988    else:
989        GmatE = G[:]
990        IntE = datum[:]
991        pwr = Power
992        QvecP = Qvec**pwr
993    Amat = GmatE*QvecP[:np.newaxis]
994    AAmat = np.inner(Amat,Amat)
995    Bvec = IntE*QvecP
996    Xw = np.ones_like(Bins)*1.e-32
997    calc = np.dot(G.T,Xw)
998    nIter = 0
999    err = 10.
1000    while (nIter<IterMax) and (err > 1.):
1001        #Step 1 in M&Z paper:
1002        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
1003        Dk = Xw/np.inner(AAmat,Xw)
1004        Pk = -Dk*Qk
1005        #Step 2 in M&Z paper:
1006        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
1007        alpWv = -Xw/Pk
1008        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
1009        #Step 3 in M&Z paper:
1010        shift = AkSt*Pk
1011        Xw += shift
1012        #done IPG; now results
1013        nIter += 1
1014        calc = np.dot(G.T,Xw)
1015        chisq = np.sum(((datum-calc)/sigma)**2)
1016        err = chisq/len(datum)
1017        if report:
1018            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
1019    return chisq,Xw,calc
1020
1021###############################################################################
1022#### SASD Utilities
1023###############################################################################
1024
1025def SetScale(Data,refData):
1026    Profile,Limits,Sample = Data
1027    refProfile,refLimits,refSample = refData
1028    x,y = Profile[:2]
1029    rx,ry = refProfile[:2]
1030    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
1031    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
1032    iBeg = np.searchsorted(x,Beg)
1033    iFin = np.searchsorted(x,Fin)
1034    sum = np.sum(y[iBeg:iFin])
1035    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
1036    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
1037   
1038def Bestimate(G,Rg,P):
1039    return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
1040   
1041###############################################################################
1042#### Size distribution
1043###############################################################################
1044
1045def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
1046    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1047        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1048        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1049        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1050        'Cylinder diam':[CylinderDFF,CylinderDVol]}
1051    Shape = data['Size']['Shape'][0]
1052    Parms = data['Size']['Shape'][1:]
1053    if data['Size']['logBins']:
1054        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
1055            data['Size']['Nbins']+1,True)/2.        #make radii
1056    else:
1057        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
1058            data['Size']['Nbins']+1,True)/2.        #make radii
1059    Dbins = np.diff(Bins)
1060    Bins = Bins[:-1]+Dbins/2.
1061    Contrast = Sample['Contrast'][1]
1062    Scale = Sample['Scale'][0]
1063    Sky = 10**data['Size']['MaxEnt']['Sky']
1064    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
1065    Back = data['Back']
1066    Q,Io,wt,Ic,Ib = Profile[:5]
1067    Qmin = Limits[1][0]
1068    Qmax = Limits[1][1]
1069    wtFactor = ProfDict['wtFactor']
1070    Ibeg = np.searchsorted(Q,Qmin)
1071    Ifin = np.searchsorted(Q,Qmax)
1072    BinMag = np.zeros_like(Bins)
1073    Ic[:] = 0.
1074    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
1075    if 'MaxEnt' == data['Size']['Method']:
1076        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
1077            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
1078            data['Size']['MaxEnt']['Niter'],report=True)
1079    elif 'IPG' == data['Size']['Method']:
1080        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
1081            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
1082            Power=data['Size']['IPG']['Power'],report=True)
1083    Ib[:] = Back[0]
1084    Ic[Ibeg:Ifin] += Back[0]
1085    print ' Final chi^2: %.3f'%(chisq)
1086    Vols = shapes[Shape][1](Bins,Parms)
1087    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
1088       
1089################################################################################
1090#### Modelling
1091################################################################################
1092
1093def ModelFit(Profile,ProfDict,Limits,Substances,Sample,Model):
1094    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1095        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1096        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1097        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1098        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
1099           
1100    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
1101            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
1102               
1103    parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
1104        'PkInt','PkPos','PkSig','PkGam',]
1105       
1106    FFparmOrder = ['Aspect ratio','Length','Diameter','Thickness']
1107   
1108    SFparmOrder = ['Dist','VolFr','epis','Sticky','Depth','Width']
1109
1110    def GetModelParms():
1111        parmDict = {'Scale':Sample['Scale'][0]}
1112        varyList = []
1113        values = []
1114        levelTypes = []
1115        Back = Model['Back']
1116        if Back[1]:
1117            varyList += ['Back',]
1118            values.append(Back[0])
1119        parmDict['Back'] = Back[0]
1120        partData = Model['Particle']
1121        parmDict['Matrix density'] = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
1122        for i,level in enumerate(partData['Levels']):
1123            cid = str(i)+':'
1124            controls = level['Controls']
1125            Type = controls['DistType']
1126            levelTypes.append(Type)
1127            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1128                if Type not in ['Monodosperse',]:
1129                    parmDict[cid+'NumPoints'] = controls['NumPoints']
1130                    parmDict[cid+'Cutoff'] = controls['Cutoff']
1131                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
1132                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
1133                parmDict[cid+'StrFact'] = sfxns[controls['StrFact']]
1134                parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
1135                for item in FFparmOrder:
1136                    if item in controls['FFargs']:
1137                        parmDict[cid+item] = controls['FFargs'][item][0]
1138                        if controls['FFargs'][item][1]:
1139                            varyList.append(cid+item)
1140                            values.append(controls['FFargs'][item][0])
1141                for item in SFparmOrder:
1142                    if item in controls.get('SFargs',{}):
1143                        parmDict[cid+item] = controls['SFargs'][item][0]
1144                        if controls['SFargs'][item][1]:
1145                            varyList.append(cid+item)
1146                            values.append(controls['SFargs'][item][0])
1147            distDict = controls['DistType']
1148            for item in parmOrder:
1149                if item in level[distDict]:
1150                    parmDict[cid+item] = level[distDict][item][0]
1151                    if level[distDict][item][1]:
1152                        values.append(level[distDict][item][0])
1153                        varyList.append(cid+item)
1154        return levelTypes,parmDict,varyList,values
1155       
1156    def SetModelParms():
1157        print ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
1158        if 'Back' in varyList:
1159            Model['Back'][0] = parmDict['Back']
1160            print '  %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
1161        partData = Model['Particle']
1162        for i,level in enumerate(partData['Levels']):
1163            controls = level['Controls']
1164            Type = controls['DistType']
1165            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1166                print ' Component %d: Type: %s: Structure Factor: %s'%(i,Type,controls['StrFact'])               
1167            else:
1168                print ' Component %d: Type: %s: '%(i,Type,)
1169            cid = str(i)+':'
1170            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
1171                for item in FFparmOrder:
1172                    if cid+item in varyList:
1173                        controls['FFargs'][item][0] = parmDict[cid+item]
1174                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1175                for item in SFparmOrder:
1176                    if cid+item in varyList:
1177                        controls['SFargs'][item][0] = parmDict[cid+item]
1178                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1179            distDict = controls['DistType']
1180            for item in level[distDict]:
1181                if cid+item in varyList:
1182                    level[distDict][item][0] = parmDict[cid+item]
1183                    print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
1184                   
1185    def calcSASD(values,Q,Io,wt,Ifb,levelTypes,parmDict,varyList):
1186        parmDict.update(zip(varyList,values))
1187        M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)+Ifb-parmDict['Scale']*Io)
1188        return M
1189       
1190    def getSASD(Q,levelTypes,parmDict):
1191        Ic = np.zeros_like(Q)
1192        rhoMat = parmDict['Matrix density']
1193        for i,Type in enumerate(levelTypes):
1194            cid = str(i)+':'
1195            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1196                FFfxn = parmDict[cid+'FormFact']
1197                Volfxn = parmDict[cid+'FFVolume']
1198                SFfxn = parmDict[cid+'StrFact']
1199                FFargs = []
1200                SFargs = []
1201                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter']:
1202                    if item in parmDict: 
1203                        FFargs.append(parmDict[item])
1204                for item in [cid+'Dist',cid+'VolFr',cid+'epis',cid+'Sticky',cid+'Depth',cid+'Width']:
1205                    if item in parmDict: 
1206                        SFargs.append(parmDict[item])
1207                distDict = {}
1208                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
1209                    if item in parmDict:
1210                        distDict[item.split(':')[1]] = parmDict[item]
1211                rho = parmDict[cid+'XAnom density']
1212                contrast = rho**2-rhoMat**2
1213                rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
1214                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
1215                dist *= parmDict[cid+'Volume']
1216                Ic += np.dot(Gmat,dist)*SFfxn(Q,args=SFargs)
1217            elif 'Unified' in Type:
1218                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
1219                    parmDict[cid+'P'],parmDict[cid+'Cutoff']
1220                Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
1221                Guin = G*np.exp(-(Q*Rg)**2/3)
1222                Porod = (B/Qstar**P)*np.exp(-(Q*Rgco)**2/3)
1223                Ic += Guin+Porod
1224            elif 'Porod' in Type:
1225                B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
1226                Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
1227                Ic += Porod
1228            elif 'Mono' in Type:
1229                FFfxn = parmDict[cid+'FormFact']
1230                Volfxn = parmDict[cid+'FFVolume']
1231                SFfxn = parmDict[cid+'StrFact']
1232                FFargs = []
1233                SFargs = []
1234                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
1235                    if item in parmDict: 
1236                        FFargs.append(parmDict[item])
1237                rho = parmDict[cid+'XAnom density']
1238                contrast = rho**2-rhoMat**2
1239                R = parmDict[cid+'Radius']
1240                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
1241                Ic += Gmat[0]*parmDict[cid+'Volume']*SFfxn(Q,args=SFargs)
1242            elif 'Bragg' in distFxn:
1243                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
1244                    parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
1245        Ic += parmDict['Back']  #/parmDict['Scale']
1246        return Ic
1247       
1248    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1249    Qmin = Limits[1][0]
1250    Qmax = Limits[1][1]
1251    wtFactor = ProfDict['wtFactor']
1252    Ibeg = np.searchsorted(Q,Qmin)
1253    Ifin = np.searchsorted(Q,Qmax)
1254    Ic[:] = 0
1255    levelTypes,parmDict,varyList,values = GetModelParms()
1256    if varyList:
1257        result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
1258            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList))
1259        ncyc = int(result[2]['nfev']/2)
1260        parmDict.update(zip(varyList,result[0]))
1261        chisq = np.sum(result[2]['fvec']**2)
1262        ncalc = result[2]['nfev']
1263        covM = result[1]
1264    else:   #nothing varied
1265        M = calcSASD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList)
1266        chisq = np.sum(M**2)
1267        ncalc = 0
1268        covM = []
1269        sig = []
1270        sigDict = {}
1271        result = []
1272    Rvals = {}
1273    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
1274    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
1275    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
1276    try:
1277        if len(covM):
1278            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
1279            sigDict = dict(zip(varyList,sig))
1280        print ' Results of small angle data modelling fit:'
1281        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
1282        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1283        SetModelParms()
1284        covMatrix = covM*Rvals['GOF']
1285        return True,result,varyList,sig,Rvals,covMatrix
1286    except ValueError:
1287        return False,0,0,0,0,0
1288   
1289def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
1290   
1291    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1292        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1293        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1294        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1295        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
1296    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
1297            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
1298
1299#    pdb.set_trace()
1300    partData = sasdData['Particle']
1301    rhoMat = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
1302    matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
1303    Scale = Sample['Scale']     #[value,flag]
1304    Back = sasdData['Back']
1305    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1306    Qmin = Limits[1][0]
1307    Qmax = Limits[1][1]
1308    wtFactor = ProfDict['wtFactor']
1309    Ibeg = np.searchsorted(Q,Qmin)
1310    Ifin = np.searchsorted(Q,Qmax)
1311    Ib[:] = Back[0]
1312    Ic[:] = 0
1313    Rbins = []
1314    Dist = []
1315    for level in partData['Levels']:
1316        controls = level['Controls']
1317        distFxn = controls['DistType']
1318        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1319            parmDict = level[controls['DistType']]
1320            FFfxn = shapes[controls['FormFact']][0]
1321            Volfxn = shapes[controls['FormFact']][1]
1322            SFfxn = sfxns[controls['StrFact']]
1323            FFargs = []
1324            SFargs = []
1325            for item in ['Dist','VolFr','epis','Sticky','Depth','Width',]:
1326                if item in controls.get('SFargs',{}):
1327                    SFargs.append(controls['SFargs'][item][0])
1328            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1329                if item in controls['FFargs']: 
1330                    FFargs.append(controls['FFargs'][item][0])
1331            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
1332            contrast = rho**2-rhoMat**2
1333            distDict = {}
1334            for item in parmDict:
1335                distDict[item] = parmDict[item][0]
1336            rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
1337            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
1338            dist *= level[distFxn]['Volume'][0]
1339            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)*SFfxn(Q[Ibeg:Ifin],args=SFargs)
1340            Rbins.append(rBins)
1341            Dist.append(dist/(4.*dBins))
1342        elif 'Unified' in distFxn:
1343            parmDict = level[controls['DistType']]
1344            Rg,G,B,P,Rgco = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],    \
1345                parmDict['P'][0],parmDict['Cutoff'][0]
1346            Qstar = Q[Ibeg:Ifin]/(scsp.erf(Q[Ibeg:Ifin]*Rg/np.sqrt(6)))**3
1347            Guin = G*np.exp(-(Q[Ibeg:Ifin]*Rg)**2/3)
1348            Porod = (B/Qstar**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1349            Ic[Ibeg:Ifin] += Guin+Porod
1350            Rbins.append([])
1351            Dist.append([])
1352        elif 'Porod' in distFxn:
1353            parmDict = level[controls['DistType']]
1354            B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
1355            Porod = (B/Q[Ibeg:Ifin]**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1356            Ic[Ibeg:Ifin] += Porod
1357            Rbins.append([])
1358            Dist.append([])
1359        elif 'Mono' in distFxn:
1360            parmDict = level[controls['DistType']]
1361            R = level[controls['DistType']]['Radius'][0]
1362            FFfxn = shapes[controls['FormFact']][0]
1363            Volfxn = shapes[controls['FormFact']][1]
1364            SFfxn = sfxns[controls['StrFact']]
1365            FFargs = []
1366            SFargs = [parmDict['Volume'][0],R]
1367            for item in ['epis','Sticky','Depth','Width',]:
1368                if item in controls.get('SFargs',{}):
1369                    SFargs.append(controls['SFargs'][item][0])
1370            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1371                if item in controls['FFargs']: 
1372                    FFargs.append(controls['FFargs'][item][0])
1373            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
1374            contrast = rho**2-rhoMat**2
1375            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
1376            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]*SFfxn(Q[Ibeg:Ifin],args=SFargs)
1377            Rbins.append([])
1378            Dist.append([])
1379        elif 'Bragg' in distFxn:
1380            parmDict = level[controls['DistType']]
1381            Ic[Ibeg:Ifin] += parmDict['PkInt'][0]*G2pwd.getPsVoigt(parmDict['PkPos'][0],
1382                parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
1383            Rbins.append([])
1384            Dist.append([])
1385    Ic[Ibeg:Ifin] += Back[0]
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        M = len(buf)
1461        buf = zip(*buf)         # transpose rows and columns
1462        q  = np.array(buf[0], dtype=np.float64)
1463        I  = np.array(buf[1], dtype=np.float64)
1464        dI = np.array(buf[2], dtype=np.float64)
1465        return q, I, dI
1466    print "MaxEnt_SB: "
1467    test_data_file = os.path.join( 'testinp', 'test.sas')
1468    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
1469    bkg   = 0.1     #   I = I - bkg
1470    dMin, dMax, nRadii = 25, 9000, 40
1471    defaultDistLevel = 1.0e-6
1472    IterMax = 40
1473    errFac = 1.05
1474   
1475    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
1476    dr   = r * (r[1]/r[0] - 1)          # step size
1477    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
1478    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
1479   
1480    qVec, I, dI = readTextData(test_data_file)
1481    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
1482   
1483    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
1484    if f_dr is None:
1485        print "no solution"
1486        return
1487   
1488    print "solution reached"
1489    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
1490        print '%10.4f %10.4f %12.4g'%(a,b,c)
1491
1492def tests():
1493    test_MaxEnt_SB(report=True)
1494
1495if __name__ == '__main__':
1496    tests()
1497
Note: See TracBrowser for help on using the repository browser.