source: trunk/GSASIIsasd.py @ 2546

Last change on this file since 2546 was 2546, checked in by vondreele, 5 years ago

cleanup of all GSASII main routines - remove unused variables, dead code, etc.

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