source: trunk/GSASIIsasd.py @ 4782

Last change on this file since 4782 was 4095, checked in by vondreele, 6 years ago

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