source: trunk/GSASIIsasd.py @ 1301

Last change on this file since 1301 was 1301, checked in by vondreele, 8 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 98.1 KB
Line 
1<<<<<<< .mine
2#/usr/bin/env python
3# -*- coding: utf-8 -*-
4'''
5*GSASII small angle calculation module*
6==================================
7
8'''
9########### SVN repository information ###################
10# $Date: 2014-04-25 17:51:35 +0000 (Fri, 25 Apr 2014) $
11# $Author: vondreele $
12# $Revision: 1301 $
13# $URL: trunk/GSASIIsasd.py $
14# $Id: GSASIIsasd.py 1301 2014-04-25 17:51:35Z vondreele $
15########### SVN repository information ###################
16import os
17import sys
18import math
19import time
20
21import numpy as np
22import scipy as sp
23import numpy.linalg as nl
24from numpy.fft import ifft, fft, fftshift
25import scipy.special as scsp
26import scipy.interpolate as si
27import scipy.stats as st
28import scipy.optimize as so
29#import pdb
30
31import GSASIIpath
32GSASIIpath.SetVersionNumber("$Revision: 1301 $")
33import GSASIIlattice as G2lat
34import GSASIIspc as G2spc
35import GSASIIElem as G2elem
36import GSASIIgrid as G2gd
37import GSASIIIO as G2IO
38import GSASIImath as G2mth
39import GSASIIpwd as G2pwd
40
41# trig functions in degrees
42sind = lambda x: math.sin(x*math.pi/180.)
43asind = lambda x: 180.*math.asin(x)/math.pi
44tand = lambda x: math.tan(x*math.pi/180.)
45atand = lambda x: 180.*math.atan(x)/math.pi
46atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
47cosd = lambda x: math.cos(x*math.pi/180.)
48acosd = lambda x: 180.*math.acos(x)/math.pi
49rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
50#numpy versions
51npsind = lambda x: np.sin(x*np.pi/180.)
52npasind = lambda x: 180.*np.arcsin(x)/math.pi
53npcosd = lambda x: np.cos(x*math.pi/180.)
54npacosd = lambda x: 180.*np.arccos(x)/math.pi
55nptand = lambda x: np.tan(x*math.pi/180.)
56npatand = lambda x: 180.*np.arctan(x)/np.pi
57npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
58npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
59npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
60   
61###############################################################################
62#### Particle form factors
63###############################################################################
64
65def SphereFF(Q,R,args=()):
66    ''' Compute hard sphere form factor - can use numpy arrays
67    param float Q: Q value array (usually in A-1)
68    param float R: sphere radius (Usually in A - must match Q-1 units)
69    param array args: ignored
70    returns float: form factors as array as needed
71    '''
72    QR = Q[:,np.newaxis]*R
73    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
74   
75def SpheroidFF(Q,R,args):
76    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
77    - can use numpy arrays for R & AR; will return corresponding numpy array
78    param float Q : Q value array (usually in A-1)
79    param float R: radius along 2 axes of spheroid
80    param array args: [float AR]: aspect ratio so 3rd axis = R*AR
81    returns float: form factors as array as needed
82    '''
83    NP = 50
84    AR = args[0]
85    if 0.99 < AR < 1.01:
86        return SphereFF(Q,R,0)
87    else:
88        cth = np.linspace(0,1.,NP)
89        Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
90        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
91           
92def CylinderFF(Q,R,args):
93    ''' Compute form factor for cylinders - can use numpy arrays
94    param float Q: Q value array (A-1)
95    param float R: cylinder radius (A)
96    param array args: [float L]: cylinder length (A)
97    returns float: form factor
98    '''
99    L = args[0]
100    NP = 200
101    alp = np.linspace(0,np.pi/2.,NP)
102    QL = Q[:,np.newaxis]*L
103    LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
104    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
105    QR = Q[:,np.newaxis]*R
106    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
107    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
108    LSBess = LBess*SBess
109    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
110   
111def CylinderDFF(Q,L,args):
112    ''' Compute form factor for cylinders - can use numpy arrays
113    param float Q: Q value array (A-1)
114    param float L: cylinder half length (A)
115    param array args: [float R]: cylinder radius (A)
116    returns float: form factor
117    '''
118    R = args[0]
119    return CylinderFF(Q,R,args=[2.*L])   
120   
121def CylinderARFF(Q,R,args): 
122    ''' Compute form factor for cylinders - can use numpy arrays
123    param float Q: Q value array (A-1)
124    param float R: cylinder radius (A)
125    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
126    returns float: form factor
127    '''
128    AR = args[0]
129    return CylinderFF(Q,R,args=[2.*R*AR])   
130   
131def UniSphereFF(Q,R,args=0):
132    ''' Compute form factor for unified sphere - can use numpy arrays
133    param float Q: Q value array (A-1)
134    param float R: cylinder radius (A)
135    param array args: ignored
136    returns float: form factor
137    '''
138    Rg = np.sqrt(3./5.)*R
139    B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
140    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
141    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
142   
143def UniRodFF(Q,R,args):
144    ''' Compute form factor for unified rod - can use numpy arrays
145    param float Q: Q value array (A-1)
146    param float R: cylinder radius (A)
147    param array args: [float R]: cylinder radius (A)
148    returns float: form factor
149    '''
150    L = args[0]
151    Rg2 = np.sqrt(R**2/2+L**2/12)
152    B2 = np.pi/L
153    Rg1 = np.sqrt(3.)*R/2.
154    G1 = (2./3.)*R/L
155    B1 = 4.*(L+R)/(R**3*L**2)
156    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
157    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
158    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
159    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
160    return np.sqrt(FF)
161   
162def UniRodARFF(Q,R,args):
163    ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
164    param float Q: Q value array (A-1)
165    param float R: cylinder radius (A)
166    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
167    returns float: form factor
168    '''
169    AR = args[0]
170    return UniRodFF(Q,R,args=[2.*AR*R,])
171   
172def UniDiskFF(Q,R,args):
173    ''' Compute form factor for unified disk - can use numpy arrays
174    param float Q: Q value array (A-1)
175    param float R: cylinder radius (A)
176    param array args: [float T]: disk thickness (A)
177    returns float: form factor
178    '''
179    T = args[0]
180    Rg2 = np.sqrt(R**2/2.+T**2/12.)
181    B2 = 2./R**2
182    Rg1 = np.sqrt(3.)*T/2.
183    RgC2 = 1.1*Rg1
184    G1 = (2./3.)*(T/R)**2
185    B1 = 4.*(T+R)/(R**3*T**2)
186    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
187    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
188    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
189    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
190    return np.sqrt(FF)
191   
192def UniTubeFF(Q,R,args):
193    ''' Compute form factor for unified tube - can use numpy arrays
194    assumes that core of tube is same as the matrix/solvent so contrast
195    is from tube wall vs matrix
196    param float Q: Q value array (A-1)
197    param float R: cylinder radius (A)
198    param array args: [float L,T]: tube length & wall thickness(A)
199    returns float: form factor
200    '''
201    L,T = args[:2]
202    Ri = R-T
203    DR2 = R**2-Ri**2
204    Vt = np.pi*DR2*L
205    Rg3 = np.sqrt(DR2/2.+L**2/12.)
206    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
207    B2 = np.pi**2*T/Vt
208    B3 = np.pi/L
209    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
210    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
211    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
212    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
213    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
214    FF += B1/QstV**4
215    return np.sqrt(FF)
216
217###############################################################################
218#### Particle volumes
219###############################################################################
220
221def SphereVol(R,args=()):
222    ''' Compute volume of sphere
223    - numpy array friendly
224    param float R: sphere radius
225    param array args: ignored
226    returns float: volume
227    '''
228    return (4./3.)*np.pi*R**3
229
230def SpheroidVol(R,args):
231    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
232    - numpy array friendly
233    param float R: radius along 2 axes of spheroid
234    param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
235    returns float: volume
236    '''
237    AR = args[0]
238    return AR*SphereVol(R)
239   
240def CylinderVol(R,args):
241    ''' Compute cylinder volume for radius & length
242    - numpy array friendly
243    param float R: diameter (A)
244    param array args: [float L]: length (A)
245    returns float:volume (A^3)
246    '''
247    L = args[0]
248    return np.pi*L*R**2
249   
250def CylinderDVol(L,args):
251    ''' Compute cylinder volume for length & diameter
252    - numpy array friendly
253    param float: L half length (A)
254    param array args: [float D]: diameter (A)
255    returns float:volume (A^3)
256    '''
257    D = args[0]
258    return CylinderVol(D/2.,[2.*L,])
259   
260def CylinderARVol(R,args):
261    ''' Compute cylinder volume for radius & aspect ratio = L/D
262    - numpy array friendly
263    param float: R radius (A)
264    param array args: [float AR]: =L/D=L/2R aspect ratio
265    returns float:volume
266    '''
267    AR = args[0]
268    return CylinderVol(R,[2.*R*AR,])
269   
270def UniSphereVol(R,args=()):
271    ''' Compute volume of sphere
272    - numpy array friendly
273    param float R: sphere radius
274    param array args: ignored
275    returns float: volume
276    '''
277    return SphereVol(R)
278   
279def UniRodVol(R,args):
280    ''' Compute cylinder volume for radius & length
281    - numpy array friendly
282    param float R: diameter (A)
283    param array args: [float L]: length (A)
284    returns float:volume (A^3)
285    '''
286    L = args[0]
287    return CylinderVol(R,[L,])
288   
289def UniRodARVol(R,args):
290    ''' Compute rod volume for radius & aspect ratio
291    - numpy array friendly
292    param float R: diameter (A)
293    param array args: [float AR]: =L/D=L/2R aspect ratio
294    returns float:volume (A^3)
295    '''
296    AR = args[0]
297    return CylinderARVol(R,[AR,])
298   
299def UniDiskVol(R,args):
300    ''' Compute disk volume for radius & thickness
301    - numpy array friendly
302    param float R: diameter (A)
303    param array args: [float T]: thickness
304    returns float:volume (A^3)
305    '''
306    T = args[0]
307    return CylinderVol(R,[T,])
308   
309def UniTubeVol(R,args):
310    ''' Compute tube volume for radius, length & wall thickness
311    - numpy array friendly
312    param float R: diameter (A)
313    param array args: [float L,T]: tube length & wall thickness(A)
314    returns float: volume (A^3) of tube wall
315    '''
316    L,T = args[:2]
317    return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
318   
319################################################################################
320#### Distribution functions & their cumulative fxns
321################################################################################
322
323def LogNormalDist(x,pos,args):
324    ''' Standard LogNormal distribution - numpy friendly on x axis
325    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
326    param float x: independent axis (can be numpy array)
327    param float pos: location of distribution
328    param float scale: width of distribution (m)
329    param float shape: shape - (sigma of log(LogNormal))
330    returns float: LogNormal distribution
331    '''
332    scale,shape = args
333    return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
334   
335def GaussDist(x,pos,args):
336    ''' Standard Normal distribution - numpy friendly on x axis
337    param float x: independent axis (can be numpy array)
338    param float pos: location of distribution
339    param float scale: width of distribution (sigma)
340    param float shape: not used
341    returns float: Normal distribution
342    '''
343    scale = args[0]
344    return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
345   
346def LSWDist(x,pos,args=[]):
347    ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
348    ref:
349    param float x: independent axis (can be numpy array)
350    param float pos: location of distribution
351    param float scale: not used
352    param float shape: not used
353    returns float: LSW distribution   
354    '''
355    redX = x/pos
356    result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
357   
358    return np.nan_to_num(result/pos)
359   
360def SchulzZimmDist(x,pos,args):
361    ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
362    ref: http://goldbook.iupac.org/S05502.html
363    param float x: independent axis (can be numpy array)
364    param float pos: location of distribution
365    param float scale: width of distribution (sigma)
366    param float shape: not used
367    returns float: Schulz-Zimm distribution
368    '''
369    scale = args[0]
370    b = (2.*pos/scale)**2
371    a = b/pos
372    if b<70:    #why bother?
373        return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
374    else:
375        return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
376           
377def LogNormalCume(x,pos,args):
378    ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
379    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
380    param float x: independent axis (can be numpy array)
381    param float pos: location of distribution
382    param float scale: width of distribution (sigma)
383    param float shape: shape parameter
384    returns float: LogNormal cumulative distribution
385    '''
386    scale,shape = args
387    lredX = np.log((x-pos)/scale)
388    return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
389   
390def GaussCume(x,pos,args):
391    ''' Standard Normal cumulative distribution - numpy friendly on x axis
392    param float x: independent axis (can be numpy array)
393    param float pos: location of distribution
394    param float scale: width of distribution (sigma)
395    param float shape: not used
396    returns float: Normal cumulative distribution
397    '''
398    scale = args[0]
399    redX = (x-pos)/scale
400    return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
401   
402def LSWCume(x,pos,args=[]):
403    ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
404    param float x: independent axis (can be numpy array)
405    param float pos: location of distribution
406    param float scale: not used
407    param float shape: not used
408    returns float: LSW cumulative distribution
409    '''
410    nP = 500
411    xMin,xMax = [0.,2.*pos]
412    X = np.linspace(xMin,xMax,nP,True)
413    fxn = LSWDist(X,pos)
414    mat = np.outer(np.ones(nP),fxn)
415    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
416    return np.interp(x,X,cume,0,1)
417   
418def SchulzZimmCume(x,pos,args):
419    ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
420    param float x: independent axis (can be numpy array)
421    param float pos: location of distribution
422    param float scale: width of distribution (sigma)
423    param float shape: not used
424    returns float: Normal distribution
425    '''
426    scale = args[0]
427    nP = 500
428    xMin = np.max([0.,pos-4.*scale])
429    xMax = np.min([pos+4.*scale,1.e5])
430    X = np.linspace(xMin,xMax,nP,True)
431    fxn = LSWDist(X,pos)
432    mat = np.outer(np.ones(nP),fxn)
433    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
434    return np.interp(x,X,cume,0,1)
435   
436    return []
437   
438         
439################################################################################
440##### SB-MaxEnt
441################################################################################
442
443def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
444    '''Calculates the response matrix :math:`G(Q,r)`
445   
446    :param float q: :math:`Q`
447    :param float r: :math:`r`
448    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
449    :param function FFfxn: form factor function FF(q,r,args)
450    :param function Volfxn: volume function Vol(r,args)
451    :returns float: G(Q,r)
452    '''
453    FF = FFfxn(q,r,args)
454    Vol = Volfxn(r,args)
455    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
456   
457'''
458sbmaxent
459
460Entropy maximization routine as described in the article
461J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
462("MNRAS": "Monthly Notices of the Royal Astronomical Society")
463
464:license: Copyright (c) 2013, UChicago Argonne, LLC
465:license: This file is distributed subject to a Software License Agreement found
466     in the file LICENSE that is included with this distribution.
467
468References:
469
4701. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
4712. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
472   Neutron Scattering Data Analysis, Rutherford
473   Appleton Laboratory, UK, 1986; ed. MW Johnson,
474   IOP Conference Series 81 (1986) 81 - 86, Institute
475   of Physics, Bristol, UK.
4763. ID Culverwell and GP Clarke; Ibid. 87 - 96.
4774. JA Potton, GK Daniell, & BD Rainford,
478   J APPL CRYST 21 (1988) 663 - 668.
4795. JA Potton, GJ Daniell, & BD Rainford,
480   J APPL CRYST 21 (1988) 891 - 897.
481
482'''
483
484class MaxEntException(Exception): 
485    '''Any exception from this module'''
486    pass
487
488def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
489    '''
490    do the complete Maximum Entropy algorithm of Skilling and Bryan
491   
492    :param float datum[]:
493    :param float sigma[]:
494    :param float[][] G: transformation matrix
495    :param float base[]:
496    :param int IterMax:
497    :param obj image_to_data: opus function (defaults to opus)
498    :param obj data_to_image: tropus function (defaults to tropus)
499   
500    :returns float[]: :math:`f(r) dr`
501    '''
502       
503    TEST_LIMIT        = 0.05                    # for convergence
504    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
505    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
506    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
507    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
508   
509    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
510    MOVE_PASSES       = 0.001                   # convergence test in routine: move
511
512    def tropus (data, G):
513        '''
514        tropus: transform data-space -> solution-space:  [G] * data
515       
516        default definition, caller can use this definition or provide an alternative
517       
518        :param float[M] data: observations, ndarray of shape (M)
519        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
520        :returns float[N]: calculated image, ndarray of shape (N)
521        '''
522        return G.dot(data)
523
524    def opus (image, G):
525        '''
526        opus: transform solution-space -> data-space:  [G]^tr * image
527       
528        default definition, caller can use this definition or provide an alternative
529       
530        :param float[N] image: solution, ndarray of shape (N)
531        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
532        :returns float[M]: calculated data, ndarray of shape (M)
533        '''
534        return np.dot(G.T,image)    #G.transpose().dot(image)
535
536    def Dist(s2, beta):
537        '''measure the distance of this possible solution'''
538        w = 0
539        n = beta.shape[0]
540        for k in range(n):
541            z = -sum(s2[k] * beta)
542            w += beta[k] * z
543        return w
544   
545    def ChiNow(ax, c1, c2, s1, s2):
546        '''
547        ChiNow
548       
549        :returns tuple: (ChiNow computation of ``w``, beta)
550        '''
551       
552        bx = 1 - ax
553        a =   bx * c2  -  ax * s2
554        b = -(bx * c1  -  ax * s1)
555   
556        beta = ChoSol(a, b)
557        w = 1.0
558        for k in range(SEARCH_DIRECTIONS):
559            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
560        return w, beta
561   
562    def ChoSol(a, b):
563        '''
564        ChoSol: ? chop the solution vectors ?
565       
566        :returns: new vector beta
567        '''
568        n = b.shape[0]
569        fl = np.zeros((n,n))
570        bl = np.zeros_like(b)
571       
572        #print_arr("ChoSol: a", a)
573        #print_vec("ChoSol: b", b)
574   
575        if (a[0][0] <= 0):
576            msg = "ChoSol: a[0][0] = " 
577            msg += str(a[0][0])
578            msg += '  Value must be positive'
579            raise MaxEntException(msg)
580   
581        # first, compute fl from a
582        # note fl is a lower triangular matrix
583        fl[0][0] = math.sqrt (a[0][0])
584        for i in (1, 2):
585            fl[i][0] = a[i][0] / fl[0][0]
586            for j in range(1, i+1):
587                z = 0.0
588                for k in range(j):
589                    z += fl[i][k] * fl[j][k]
590                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
591                z = a[i][j] - z
592                if j == i:
593                    y = math.sqrt(z)
594                else:
595                    y = z / fl[j][j]
596                fl[i][j] = y
597        #print_arr("ChoSol: fl", fl)
598   
599        # next, compute bl from fl and b
600        bl[0] = b[0] / fl[0][0]
601        for i in (1, 2):
602            z = 0.0
603            for k in range(i):
604                z += fl[i][k] * bl[k]
605                #print "\t", i, k, z
606            bl[i] = (b[i] - z) / fl[i][i]
607        #print_vec("ChoSol: bl", bl)
608   
609        # last, compute beta from bl and fl
610        beta = np.empty((n))
611        beta[-1] = bl[-1] / fl[-1][-1]
612        for i in (1, 0):
613            z = 0.0
614            for k in range(i+1, n):
615                z += fl[k][i] * beta[k]
616                #print "\t\t", i, k, 'z=', z
617            beta[i] = (bl[i] - z) / fl[i][i]
618        #print_vec("ChoSol: beta", beta)
619   
620        return beta
621
622    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
623        '''
624        move beta one step closer towards the solution
625        '''
626        a_lower, a_upper = 0., 1.          # bracket  "a"
627        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
628        #print "MaxEntMove: cmin = %g" % cmin
629        if cmin*chisq > chizer:
630            ctarg = (1.0 + cmin)/2
631        else:
632            ctarg = chizer/chisq
633        f_lower = cmin - ctarg
634        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
635        f_upper = c_upper - ctarg
636   
637        fx = 2*MOVE_PASSES      # just to start off
638        loop = 1
639        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
640            a_new = (a_lower + a_upper) * 0.5           # search by bisection
641            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
642            fx = c_new - ctarg
643            # tighten the search range for the next pass
644            if f_lower*fx > 0:
645                a_lower, f_lower = a_new, fx
646            if f_upper*fx > 0:
647                a_upper, f_upper = a_new, fx
648            loop += 1
649   
650        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
651            msg = "MaxEntMove: Loop counter = " 
652            msg += str(MAX_MOVE_LOOPS)
653            msg += '  No convergence in alpha chop'
654            raise MaxEntException(msg)
655   
656        w = Dist (s2, beta);
657        m = SEARCH_DIRECTIONS
658        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
659            for k in range(m):
660                beta[k] *= math.sqrt (fSum/(blank*w))
661        chtarg = ctarg * chisq
662        return w, chtarg, loop, a_new, fx, beta
663       
664#MaxEnt_SB starts here   
665
666    if image_to_data == None:
667        image_to_data = opus
668    if data_to_image == None:
669        data_to_image = tropus
670    n   = len(base)
671    npt = len(datum)
672
673    # Note that the order of subscripts for
674    # "xi" and "eta" has been reversed from
675    # the convention used in the FORTRAN version
676    # to enable parts of them to be passed as
677    # as vectors to "image_to_data" and "data_to_image".
678    xi      = np.zeros((SEARCH_DIRECTIONS, n))
679    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
680    beta    = np.zeros((SEARCH_DIRECTIONS))
681    # s1      = np.zeros((SEARCH_DIRECTIONS))
682    # c1      = np.zeros((SEARCH_DIRECTIONS))
683    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
684    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
685
686    # TODO: replace blank (scalar) with base (vector)
687    blank = sum(base) / len(base)   # use the average value of base
688
689    chizer, chtarg = npt*1.0, npt*1.0
690    f = base * 1.0                              # starting distribution is base
691
692    fSum  = sum(f)                              # find the sum of the f-vector
693    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
694    chisq = sum(z*z)                            # Chi^2, SB eq. 4
695
696    for iter in range(IterMax):
697        ox = -2 * z / sigma                        # gradient of Chi^2
698
699        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
700        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
701        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
702        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
703        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
704
705        a = 1.0
706        b = 1.0 / cnorm
707        if iter == 0:
708            test = 0.0     # mismatch between entropy and ChiSquared gradients
709        else:
710            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
711            a = 0.5 / (snorm * test)
712            b *= 0.5 / test
713        xi[0] = f * cgrad / cnorm
714        xi[1] = f * (a * sgrad - b * cgrad)
715
716        eta[0] = image_to_data (xi[0], G);          # image --> data
717        eta[1] = image_to_data (xi[1], G);          # image --> data
718        ox = eta[1] / (sigma * sigma)
719        xi[2] = data_to_image (ox, G);              # data --> image
720        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
721        xi[2] = f * xi[2] * a
722        eta[2] = image_to_data (xi[2], G)           # image --> data
723       
724#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
725#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
726
727        # prepare the search directions for the conjugate gradient technique
728        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
729        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
730#         print_vec("MaxEnt: c1", c1)
731#         print_vec("MaxEnt: s1", s1)
732
733        for k in range(SEARCH_DIRECTIONS):
734            for l in range(k+1):
735                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
736                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
737#         print_arr("MaxEnt: c2", c2)
738#         print_arr("MaxEnt: s2", s2)
739
740        # reflect across the body diagonal
741        for k, l in ((0,1), (0,2), (1,2)):
742            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
743            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
744 
745        beta[0] = -0.5 * c1[0] / c2[0][0]
746        beta[1] = 0.0
747        beta[2] = 0.0
748        if (iter > 0):
749            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
750
751        f_old = f.copy()    # preserve the last image
752        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
753       
754        # As mentioned at the top of p.119,
755        # need to protect against stray negative values.
756        # In this case, set them to RESET_STRAYS * base[i]
757        #f = f.clip(RESET_STRAYS * blank, f.max())
758        for i in range(n):
759            if f[i] <= 0.0:
760                f[i] = RESET_STRAYS * base[i]
761        df = f - f_old
762        fSum = sum(f)
763        fChange = sum(df)
764
765        # calculate the normalized entropy
766        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
767        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
768        chisq = sum(z*z)                            # report this ChiSq
769
770        if report:
771            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
772            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
773            if iter > 0:
774                value = 100*( math.sqrt(chisq/chtarg)-1)
775            else:
776                value = 0
777      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
778            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
779
780        # See if we have finished our task.
781        # do the hardest test first
782        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
783            print ' Convergence achieved.'
784            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
785    print ' No convergence! Try increasing Error multiplier.'
786    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
787
788   
789###############################################################################
790#### IPG/TNNLS Routines
791###############################################################################
792
793def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
794    ''' An implementation of the Interior-Point Gradient method of
795    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
796    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
797    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
798    Problem addressed: Total Non-Negative Least Squares (TNNLS)
799    :param float datum[]:
800    :param float sigma[]:
801    :param float[][] G: transformation matrix
802    :param int IterMax:
803    :param float Qvec: data positions for Power = 0-4
804    :param float approach: 0.8 default fitting parameter
805    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
806   
807    '''
808    if Power < 0: 
809        GmatE = G/sigma[:np.newaxis]
810        IntE = datum/sigma
811        pwr = 0
812        QvecP = np.ones_like(datum)
813    else:
814        GmatE = G[:]
815        IntE = datum[:]
816        pwr = Power
817        QvecP = Qvec**pwr
818    Amat = GmatE*QvecP[:np.newaxis]
819    AAmat = np.inner(Amat,Amat)
820    Bvec = IntE*QvecP
821    Xw = np.ones_like(Bins)*1.e-32
822    calc = np.dot(G.T,Xw)
823    nIter = 0
824    err = 10.
825    while (nIter<IterMax) and (err > 1.):
826        #Step 1 in M&Z paper:
827        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
828        Dk = Xw/np.inner(AAmat,Xw)
829        Pk = -Dk*Qk
830        #Step 2 in M&Z paper:
831        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
832        alpWv = -Xw/Pk
833        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
834        #Step 3 in M&Z paper:
835        shift = AkSt*Pk
836        Xw += shift
837        #done IPG; now results
838        nIter += 1
839        calc = np.dot(G.T,Xw)
840        chisq = np.sum(((datum-calc)/sigma)**2)
841        err = chisq/len(datum)
842        if report:
843            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
844    return chisq,Xw,calc
845
846###############################################################################
847#### SASD Utilities
848###############################################################################
849
850def SetScale(Data,refData):
851    Profile,Limits,Sample = Data
852    refProfile,refLimits,refSample = refData
853    x,y = Profile[:2]
854    rx,ry = refProfile[:2]
855    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
856    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
857    iBeg = np.searchsorted(x,Beg)
858    iFin = np.searchsorted(x,Fin)
859    sum = np.sum(y[iBeg:iFin])
860    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
861    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
862   
863def Bestimate(G,Rg,P):
864    return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
865   
866###############################################################################
867#### Size distribution
868###############################################################################
869
870def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
871    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
872        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
873        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
874        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
875        'Cylinder diam':[CylinderDFF,CylinderDVol]}
876    Shape = data['Size']['Shape'][0]
877    Parms = data['Size']['Shape'][1:]
878    if data['Size']['logBins']:
879        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
880            data['Size']['Nbins']+1,True)/2.        #make radii
881    else:
882        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
883            data['Size']['Nbins']+1,True)/2.        #make radii
884    Dbins = np.diff(Bins)
885    Bins = Bins[:-1]+Dbins/2.
886    Contrast = Sample['Contrast'][1]
887    Scale = Sample['Scale'][0]
888    Sky = 10**data['Size']['MaxEnt']['Sky']
889    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
890    Back = data['Back']
891    Q,Io,wt,Ic,Ib = Profile[:5]
892    Qmin = Limits[1][0]
893    Qmax = Limits[1][1]
894    wtFactor = ProfDict['wtFactor']
895    Ibeg = np.searchsorted(Q,Qmin)
896    Ifin = np.searchsorted(Q,Qmax)
897    BinMag = np.zeros_like(Bins)
898    Ic[:] = 0.
899    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
900    if 'MaxEnt' == data['Size']['Method']:
901        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
902            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
903            data['Size']['MaxEnt']['Niter'],report=True)
904    elif 'IPG' == data['Size']['Method']:
905        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
906            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
907            Power=data['Size']['IPG']['Power'],report=True)
908    Ib[:] = Back[0]
909    Ic[Ibeg:Ifin] += Back[0]
910    print ' Final chi^2: %.3f'%(chisq)
911    Vols = shapes[Shape][1](Bins,Parms)
912    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
913       
914################################################################################
915#### Modelling
916################################################################################
917
918def ModelFit(Profile,ProfDict,Limits,Substances,Sample,Model):
919    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
920        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
921        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
922        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
923        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
924
925    def GetModelParms():
926        parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
927            'PkInt','PkPos','PkSig','PkGam',]
928        parmDict = {'Scale':Sample['Scale'][0]}
929        varyList = []
930        values = []
931        levelTypes = []
932        Back = Model['Back']
933        if Back[1]:
934            varyList += ['Back',]
935            values.append(Back[0])
936        parmDict['Back'] = Back[0]
937        partData = Model['Particle']
938        parmDict['Matrix density'] = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
939        for i,level in enumerate(partData['Levels']):
940            cid = str(i)+':'
941            controls = level['Controls']
942            Type = controls['DistType']
943            levelTypes.append(Type)
944            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
945                if Type not in ['Monodosperse',]:
946                    parmDict[cid+'NumPoints'] = controls['NumPoints']
947                    parmDict[cid+'Cutoff'] = controls['Cutoff']
948                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
949                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
950                parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
951                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
952                    if item in controls['FFargs']:
953                        parmDict[cid+item] = controls['FFargs'][item][0]
954                        if controls['FFargs'][item][1]:
955                            varyList.append(cid+item)
956                            values.append(controls['FFargs'][item][0])
957            distDict = controls['DistType']
958            for item in parmOrder:
959                if item in level[distDict]:
960                    parmDict[cid+item] = level[distDict][item][0]
961                    if level[distDict][item][1]:
962                        values.append(level[distDict][item][0])
963                        varyList.append(cid+item)
964        return levelTypes,parmDict,varyList,values
965       
966    def SetModelParms():
967        print ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
968        if 'Back' in varyList:
969            Model['Back'][0] = parmDict['Back']
970            print %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
971        partData = Model['Particle']
972        for i,level in enumerate(partData['Levels']):
973            Type = level['Controls']['DistType']
974            print ' Component %d: Type: %s:'%(i,Type)
975            cid = str(i)+':'
976            controls = level['Controls']
977            Type = controls['DistType']
978            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
979                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
980                    if cid+item in varyList:
981                        controls['FFargs'][item][0] = parmDict[cid+item]
982                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
983            distDict = controls['DistType']
984            for item in level[distDict]:
985                if cid+item in varyList:
986                    level[distDict][item][0] = parmDict[cid+item]
987                    print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
988                   
989    def calcSASD(values,Q,Io,wt,Ifb,levelTypes,parmDict,varyList):
990        parmDict.update(zip(varyList,values))
991        M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)+Ifb-parmDict['Scale']*Io)
992        return M
993       
994    def getSASD(Q,levelTypes,parmDict):
995        Ic = np.zeros_like(Q)
996        rhoMat = parmDict['Matrix density']
997        for i,Type in enumerate(levelTypes):
998            cid = str(i)+':'
999            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1000                FFfxn = parmDict[cid+'FormFact']
1001                Volfxn = parmDict[cid+'FFVolume']
1002                FFargs = []
1003                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
1004                    if item in parmDict: 
1005                        FFargs.append(parmDict[item])
1006                distDict = {}
1007                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
1008                    if item in parmDict:
1009                        distDict[item.split(':')[1]] = parmDict[item]
1010                rho = parmDict[cid+'XAnom density']
1011                contrast = rho**2-rhoMat**2
1012                rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
1013                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
1014                dist *= parmDict[cid+'Volume']
1015                Ic += np.dot(Gmat,dist)
1016            elif 'Unified' in Type:
1017                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
1018                    parmDict[cid+'P'],parmDict[cid+'Cutoff']
1019                Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
1020                Guin = G*np.exp(-(Q*Rg)**2/3)
1021                Porod = (B/Qstar**P)*np.exp(-(Q*Rgco)**2/3)
1022                Ic += Guin+Porod
1023            elif 'Porod' in Type:
1024                B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
1025                Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
1026                Ic += Porod
1027            elif 'Mono' in Type:
1028                FFfxn = parmDict[cid+'FormFact']
1029                Volfxn = parmDict[cid+'FFVolume']
1030                FFargs = []
1031                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
1032                    if item in parmDict: 
1033                        FFargs.append(parmDict[item])
1034                rho = parmDict[cid+'XAnom density']
1035                contrast = rho**2-rhoMat**2
1036                R = parmDict[cid+'Radius']
1037                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
1038                Ic += Gmat[0]*parmDict[cid+'Volume']
1039            elif 'Bragg' in distFxn:
1040                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
1041                    parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
1042        Ic += parmDict['Back']  #/parmDict['Scale']
1043        return Ic
1044       
1045    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1046    Qmin = Limits[1][0]
1047    Qmax = Limits[1][1]
1048    wtFactor = ProfDict['wtFactor']
1049    Ibeg = np.searchsorted(Q,Qmin)
1050    Ifin = np.searchsorted(Q,Qmax)
1051    Ic[:] = 0
1052    levelTypes,parmDict,varyList,values = GetModelParms()
1053    if varyList:
1054        result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
1055            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList))
1056        ncyc = int(result[2]['nfev']/2)
1057        parmDict.update(zip(varyList,result[0]))
1058        chisq = np.sum(result[2]['fvec']**2)
1059        ncalc = result[2]['nfev']
1060        covM = result[1]
1061    else:   #nothing varied
1062        M = calcSASD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList)
1063        chisq = np.sum(M**2)
1064        ncalc = 0
1065        covM = []
1066        sig = []
1067        sigDict = {}
1068        result = []
1069    Rvals = {}
1070    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
1071    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
1072    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
1073    try:
1074        if len(covM):
1075            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
1076            sigDict = dict(zip(varyList,sig))
1077        print ' Results of small angle data modelling fit:'
1078        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
1079        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1080        SetModelParms()
1081        covMatrix = covM*Rvals['GOF']
1082        return True,result,varyList,sig,Rvals,covMatrix
1083    except ValueError:
1084        return False,0,0,0,0,0
1085   
1086def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
1087    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1088        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1089        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1090        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1091        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
1092#    pdb.set_trace()
1093    partData = sasdData['Particle']
1094    rhoMat = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
1095    matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
1096    Scale = Sample['Scale']     #[value,flag]
1097    Back = sasdData['Back']
1098    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1099    Qmin = Limits[1][0]
1100    Qmax = Limits[1][1]
1101    wtFactor = ProfDict['wtFactor']
1102    Ibeg = np.searchsorted(Q,Qmin)
1103    Ifin = np.searchsorted(Q,Qmax)
1104    Ib[:] = Back[0]
1105    Ic[:] = 0
1106    Rbins = []
1107    Dist = []
1108    for level in partData['Levels']:
1109        controls = level['Controls']
1110        distFxn = controls['DistType']
1111        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1112            FFfxn = shapes[controls['FormFact']][0]
1113            Volfxn = shapes[controls['FormFact']][1]
1114            FFargs = []
1115            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1116                if item in controls['FFargs']: 
1117                    FFargs.append(controls['FFargs'][item][0])
1118            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
1119            contrast = rho**2-rhoMat**2
1120            parmDict = level[controls['DistType']]
1121            distDict = {}
1122            for item in parmDict:
1123                distDict[item] = parmDict[item][0]
1124            rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
1125            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
1126            dist *= level[distFxn]['Volume'][0]
1127            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)
1128            Rbins.append(rBins)
1129            Dist.append(dist/(4.*dBins))
1130        elif 'Unified' in distFxn:
1131            parmDict = level[controls['DistType']]
1132            Rg,G,B,P,Rgco = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],    \
1133                parmDict['P'][0],parmDict['Cutoff'][0]
1134            Qstar = Q[Ibeg:Ifin]/(scsp.erf(Q[Ibeg:Ifin]*Rg/np.sqrt(6)))**3
1135            Guin = G*np.exp(-(Q[Ibeg:Ifin]*Rg)**2/3)
1136            Porod = (B/Qstar**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1137            Ic[Ibeg:Ifin] += Guin+Porod
1138            Rbins.append([])
1139            Dist.append([])
1140        elif 'Porod' in distFxn:
1141            parmDict = level[controls['DistType']]
1142            B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
1143            Porod = (B/Q[Ibeg:Ifin]**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1144            Ic[Ibeg:Ifin] += Porod
1145            Rbins.append([])
1146            Dist.append([])
1147        elif 'Mono' in distFxn:
1148            FFfxn = shapes[controls['FormFact']][0]
1149            Volfxn = shapes[controls['FormFact']][1]
1150            FFargs = []
1151            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1152                if item in controls['FFargs']: 
1153                    FFargs.append(controls['FFargs'][item][0])
1154            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
1155            contrast = rho**2-rhoMat**2
1156            R = level[controls['DistType']]['Radius'][0]
1157            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
1158            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]
1159            Rbins.append([])
1160            Dist.append([])
1161        elif 'Bragg' in distFxn:
1162            parmDict = level[controls['DistType']]
1163            Ic[Ibeg:Ifin] += parmDict['PkInt'][0]*G2pwd.getPsVoigt(parmDict['PkPos'][0],
1164                parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
1165            Rbins.append([])
1166            Dist.append([])
1167    Ic[Ibeg:Ifin] += Back[0]
1168    sasdData['Size Calc'] = [Rbins,Dist]
1169   
1170def MakeDiamDist(DistName,nPoints,cutoff,distDict):
1171   
1172    if 'LogNormal' in DistName:
1173        distFxn = 'LogNormalDist'
1174        cumeFxn = 'LogNormalCume'
1175        pos = distDict['MinSize']
1176        args = [distDict['Mean'],distDict['StdDev']]
1177        step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))
1178        mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)
1179        minX = 1. #pos
1180        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
1181    elif 'Gauss' in DistName:
1182        distFxn = 'GaussDist'
1183        cumeFxn = 'GaussCume'
1184        pos = distDict['Mean']
1185        args = [distDict['StdDev']]
1186        step = 0.02*distDict['StdDev']
1187        mode = distDict['Mean']
1188        minX = np.max([mode-4.*distDict['StdDev'],1.])
1189        maxX = np.min([mode+4.*distDict['StdDev'],1.e5])
1190    elif 'LSW' in DistName:
1191        distFxn = 'LSWDist'
1192        cumeFxn = 'LSWCume'
1193        pos = distDict['Mean']
1194        args = []
1195        minX,maxX = [0.,2.*pos]
1196    elif 'Schulz' in DistName:
1197        distFxn = 'SchulzZimmDist'
1198        cumeFxn = 'SchulzZimmCume'
1199        pos = distDict['Mean']
1200        args = [distDict['StdDev']]
1201        minX = np.max([1.,pos-4.*distDict['StdDev']])
1202        maxX = np.min([pos+4.*distDict['StdDev'],1.e5])
1203    nP = 500
1204    Diam = np.logspace(0.,5.,nP,True)
1205    TCW = eval(cumeFxn+'(Diam,pos,args)')
1206    CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints+1,True)
1207    rBins = np.interp(CumeTgts,TCW,Diam,0,0)
1208    dBins = np.diff(rBins)
1209    rBins = rBins[:-1]+dBins/2.
1210    return rBins,dBins,eval(distFxn+'(rBins,pos,args)')
1211
1212################################################################################
1213#### MaxEnt testing stuff
1214################################################################################
1215
1216def print_vec(text, a):
1217    '''print the contents of a vector to the console'''
1218    n = a.shape[0]
1219    print "%s[ = (" % text,
1220    for i in range(n):
1221        s = " %g, " % a[i]
1222        print s,
1223    print ")"
1224
1225def print_arr(text, a):
1226    '''print the contents of an array to the console'''
1227    n, m = a.shape
1228    print "%s[][] = (" % text
1229    for i in range(n):
1230        print " (",
1231        for j in range(m):
1232            print " %g, " % a[i][j],
1233        print "),"
1234    print ")"
1235
1236def test_MaxEnt_SB(report=True):
1237    def readTextData(filename):
1238        '''return q, I, dI from a 3-column text file'''
1239        if not os.path.exists(filename):
1240            raise Exception("file not found: " + filename)
1241        buf = [line.split() for line in open(filename, 'r').readlines()]
1242        M = len(buf)
1243        buf = zip(*buf)         # transpose rows and columns
1244        q  = np.array(buf[0], dtype=np.float64)
1245        I  = np.array(buf[1], dtype=np.float64)
1246        dI = np.array(buf[2], dtype=np.float64)
1247        return q, I, dI
1248    print "MaxEnt_SB: "
1249    test_data_file = os.path.join( 'testinp', 'test.sas')
1250    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
1251    bkg   = 0.1     #   I = I - bkg
1252    dMin, dMax, nRadii = 25, 9000, 40
1253    defaultDistLevel = 1.0e-6
1254    IterMax = 40
1255    errFac = 1.05
1256   
1257    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
1258    dr   = r * (r[1]/r[0] - 1)          # step size
1259    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
1260    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
1261   
1262    qVec, I, dI = readTextData(test_data_file)
1263    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
1264   
1265    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
1266    if f_dr is None:
1267        print "no solution"
1268        return
1269   
1270    print "solution reached"
1271    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
1272        print '%10.4f %10.4f %12.4g'%(a,b,c)
1273
1274def tests():
1275    test_MaxEnt_SB(report=True)
1276
1277if __name__ == '__main__':
1278    tests()
1279
1280=======
1281#/usr/bin/env python
1282# -*- coding: utf-8 -*-
1283'''
1284*GSASII small angle calculation module*
1285========================================
1286
1287'''
1288########### SVN repository information ###################
1289# $Date: 2014-04-25 17:51:35 +0000 (Fri, 25 Apr 2014) $
1290# $Author: vondreele $
1291# $Revision: 1301 $
1292# $URL: trunk/GSASIIsasd.py $
1293# $Id: GSASIIsasd.py 1301 2014-04-25 17:51:35Z vondreele $
1294########### SVN repository information ###################
1295import os
1296import sys
1297import math
1298import time
1299
1300import numpy as np
1301import scipy as sp
1302import numpy.linalg as nl
1303from numpy.fft import ifft, fft, fftshift
1304import scipy.special as scsp
1305import scipy.interpolate as si
1306import scipy.stats as st
1307import scipy.optimize as so
1308#import pdb
1309
1310import GSASIIpath
1311GSASIIpath.SetVersionNumber("$Revision: 1301 $")
1312import GSASIIlattice as G2lat
1313import GSASIIspc as G2spc
1314import GSASIIElem as G2elem
1315import GSASIIgrid as G2gd
1316import GSASIIIO as G2IO
1317import GSASIImath as G2mth
1318import GSASIIpwd as G2pwd
1319
1320# trig functions in degrees
1321sind = lambda x: math.sin(x*math.pi/180.)
1322asind = lambda x: 180.*math.asin(x)/math.pi
1323tand = lambda x: math.tan(x*math.pi/180.)
1324atand = lambda x: 180.*math.atan(x)/math.pi
1325atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
1326cosd = lambda x: math.cos(x*math.pi/180.)
1327acosd = lambda x: 180.*math.acos(x)/math.pi
1328rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
1329#numpy versions
1330npsind = lambda x: np.sin(x*np.pi/180.)
1331npasind = lambda x: 180.*np.arcsin(x)/math.pi
1332npcosd = lambda x: np.cos(x*math.pi/180.)
1333npacosd = lambda x: 180.*np.arccos(x)/math.pi
1334nptand = lambda x: np.tan(x*math.pi/180.)
1335npatand = lambda x: 180.*np.arctan(x)/np.pi
1336npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
1337npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
1338npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
1339   
1340###############################################################################
1341#### Particle form factors
1342###############################################################################
1343
1344def SphereFF(Q,R,args=()):
1345    ''' Compute hard sphere form factor - can use numpy arrays
1346    param float Q: Q value array (usually in A-1)
1347    param float R: sphere radius (Usually in A - must match Q-1 units)
1348    param array args: ignored
1349    returns float: form factors as array as needed
1350    '''
1351    QR = Q[:,np.newaxis]*R
1352    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
1353   
1354def SpheroidFF(Q,R,args):
1355    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
1356    - can use numpy arrays for R & AR; will return corresponding numpy array
1357    param float Q : Q value array (usually in A-1)
1358    param float R: radius along 2 axes of spheroid
1359    param array args: [float AR]: aspect ratio so 3rd axis = R*AR
1360    returns float: form factors as array as needed
1361    '''
1362    NP = 50
1363    AR = args[0]
1364    if 0.99 < AR < 1.01:
1365        return SphereFF(Q,R,0)
1366    else:
1367        cth = np.linspace(0,1.,NP)
1368        Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
1369        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
1370           
1371def CylinderFF(Q,R,args):
1372    ''' Compute form factor for cylinders - can use numpy arrays
1373    param float Q: Q value array (A-1)
1374    param float R: cylinder radius (A)
1375    param array args: [float L]: cylinder length (A)
1376    returns float: form factor
1377    '''
1378    L = args[0]
1379    NP = 200
1380    alp = np.linspace(0,np.pi/2.,NP)
1381    QL = Q[:,np.newaxis]*L
1382    LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
1383    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
1384    QR = Q[:,np.newaxis]*R
1385    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
1386    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
1387    LSBess = LBess*SBess
1388    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
1389   
1390def CylinderDFF(Q,L,args):
1391    ''' Compute form factor for cylinders - can use numpy arrays
1392    param float Q: Q value array (A-1)
1393    param float L: cylinder half length (A)
1394    param array args: [float R]: cylinder radius (A)
1395    returns float: form factor
1396    '''
1397    R = args[0]
1398    return CylinderFF(Q,R,args=[2.*L])   
1399   
1400def CylinderARFF(Q,R,args): 
1401    ''' Compute form factor for cylinders - can use numpy arrays
1402    param float Q: Q value array (A-1)
1403    param float R: cylinder radius (A)
1404    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
1405    returns float: form factor
1406    '''
1407    AR = args[0]
1408    return CylinderFF(Q,R,args=[2.*R*AR])   
1409   
1410def UniSphereFF(Q,R,args=0):
1411    ''' Compute form factor for unified sphere - can use numpy arrays
1412    param float Q: Q value array (A-1)
1413    param float R: cylinder radius (A)
1414    param array args: ignored
1415    returns float: form factor
1416    '''
1417    Rg = np.sqrt(3./5.)*R
1418    B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
1419    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
1420    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
1421   
1422def UniRodFF(Q,R,args):
1423    ''' Compute form factor for unified rod - can use numpy arrays
1424    param float Q: Q value array (A-1)
1425    param float R: cylinder radius (A)
1426    param array args: [float R]: cylinder radius (A)
1427    returns float: form factor
1428    '''
1429    L = args[0]
1430    Rg2 = np.sqrt(R**2/2+L**2/12)
1431    B2 = np.pi/L
1432    Rg1 = np.sqrt(3.)*R/2.
1433    G1 = (2./3.)*R/L
1434    B1 = 4.*(L+R)/(R**3*L**2)
1435    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
1436    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
1437    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
1438    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
1439    return np.sqrt(FF)
1440   
1441def UniRodARFF(Q,R,args):
1442    ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
1443    param float Q: Q value array (A-1)
1444    param float R: cylinder radius (A)
1445    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
1446    returns float: form factor
1447    '''
1448    AR = args[0]
1449    return UniRodFF(Q,R,args=[2.*AR*R,])
1450   
1451def UniDiskFF(Q,R,args):
1452    ''' Compute form factor for unified disk - can use numpy arrays
1453    param float Q: Q value array (A-1)
1454    param float R: cylinder radius (A)
1455    param array args: [float T]: disk thickness (A)
1456    returns float: form factor
1457    '''
1458    T = args[0]
1459    Rg2 = np.sqrt(R**2/2.+T**2/12.)
1460    B2 = 2./R**2
1461    Rg1 = np.sqrt(3.)*T/2.
1462    RgC2 = 1.1*Rg1
1463    G1 = (2./3.)*(T/R)**2
1464    B1 = 4.*(T+R)/(R**3*T**2)
1465    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
1466    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
1467    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
1468    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
1469    return np.sqrt(FF)
1470   
1471def UniTubeFF(Q,R,args):
1472    ''' Compute form factor for unified tube - can use numpy arrays
1473    assumes that core of tube is same as the matrix/solvent so contrast
1474    is from tube wall vs matrix
1475    param float Q: Q value array (A-1)
1476    param float R: cylinder radius (A)
1477    param array args: [float L,T]: tube length & wall thickness(A)
1478    returns float: form factor
1479    '''
1480    L,T = args[:2]
1481    Ri = R-T
1482    DR2 = R**2-Ri**2
1483    Vt = np.pi*DR2*L
1484    Rg3 = np.sqrt(DR2/2.+L**2/12.)
1485    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
1486    B2 = np.pi**2*T/Vt
1487    B3 = np.pi/L
1488    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
1489    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
1490    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
1491    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
1492    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
1493    FF += B1/QstV**4
1494    return np.sqrt(FF)
1495
1496###############################################################################
1497#### Particle volumes
1498###############################################################################
1499
1500def SphereVol(R,args=()):
1501    ''' Compute volume of sphere
1502    - numpy array friendly
1503    param float R: sphere radius
1504    param array args: ignored
1505    returns float: volume
1506    '''
1507    return (4./3.)*np.pi*R**3
1508
1509def SpheroidVol(R,args):
1510    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
1511    - numpy array friendly
1512    param float R: radius along 2 axes of spheroid
1513    param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
1514    returns float: volume
1515    '''
1516    AR = args[0]
1517    return AR*SphereVol(R)
1518   
1519def CylinderVol(R,args):
1520    ''' Compute cylinder volume for radius & length
1521    - numpy array friendly
1522    param float R: diameter (A)
1523    param array args: [float L]: length (A)
1524    returns float:volume (A^3)
1525    '''
1526    L = args[0]
1527    return np.pi*L*R**2
1528   
1529def CylinderDVol(L,args):
1530    ''' Compute cylinder volume for length & diameter
1531    - numpy array friendly
1532    param float: L half length (A)
1533    param array args: [float D]: diameter (A)
1534    returns float:volume (A^3)
1535    '''
1536    D = args[0]
1537    return CylinderVol(D/2.,[2.*L,])
1538   
1539def CylinderARVol(R,args):
1540    ''' Compute cylinder volume for radius & aspect ratio = L/D
1541    - numpy array friendly
1542    param float: R radius (A)
1543    param array args: [float AR]: =L/D=L/2R aspect ratio
1544    returns float:volume
1545    '''
1546    AR = args[0]
1547    return CylinderVol(R,[2.*R*AR,])
1548   
1549def UniSphereVol(R,args=()):
1550    ''' Compute volume of sphere
1551    - numpy array friendly
1552    param float R: sphere radius
1553    param array args: ignored
1554    returns float: volume
1555    '''
1556    return SphereVol(R)
1557   
1558def UniRodVol(R,args):
1559    ''' Compute cylinder volume for radius & length
1560    - numpy array friendly
1561    param float R: diameter (A)
1562    param array args: [float L]: length (A)
1563    returns float:volume (A^3)
1564    '''
1565    L = args[0]
1566    return CylinderVol(R,[L,])
1567   
1568def UniRodARVol(R,args):
1569    ''' Compute rod volume for radius & aspect ratio
1570    - numpy array friendly
1571    param float R: diameter (A)
1572    param array args: [float AR]: =L/D=L/2R aspect ratio
1573    returns float:volume (A^3)
1574    '''
1575    AR = args[0]
1576    return CylinderARVol(R,[AR,])
1577   
1578def UniDiskVol(R,args):
1579    ''' Compute disk volume for radius & thickness
1580    - numpy array friendly
1581    param float R: diameter (A)
1582    param array args: [float T]: thickness
1583    returns float:volume (A^3)
1584    '''
1585    T = args[0]
1586    return CylinderVol(R,[T,])
1587   
1588def UniTubeVol(R,args):
1589    ''' Compute tube volume for radius, length & wall thickness
1590    - numpy array friendly
1591    param float R: diameter (A)
1592    param array args: [float L,T]: tube length & wall thickness(A)
1593    returns float: volume (A^3) of tube wall
1594    '''
1595    L,T = args[:2]
1596    return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
1597   
1598################################################################################
1599#### Distribution functions & their cumulative fxns
1600################################################################################
1601
1602def LogNormalDist(x,pos,args):
1603    ''' Standard LogNormal distribution - numpy friendly on x axis
1604    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
1605    param float x: independent axis (can be numpy array)
1606    param float pos: location of distribution
1607    param float scale: width of distribution (m)
1608    param float shape: shape - (sigma of log(LogNormal))
1609    returns float: LogNormal distribution
1610    '''
1611    scale,shape = args
1612    return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
1613   
1614def GaussDist(x,pos,args):
1615    ''' Standard Normal distribution - numpy friendly on x axis
1616    param float x: independent axis (can be numpy array)
1617    param float pos: location of distribution
1618    param float scale: width of distribution (sigma)
1619    param float shape: not used
1620    returns float: Normal distribution
1621    '''
1622    scale = args[0]
1623    return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
1624   
1625def LSWDist(x,pos,args=[]):
1626    ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
1627    ref:
1628    param float x: independent axis (can be numpy array)
1629    param float pos: location of distribution
1630    param float scale: not used
1631    param float shape: not used
1632    returns float: LSW distribution   
1633    '''
1634    redX = x/pos
1635    result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
1636   
1637    return np.nan_to_num(result/pos)
1638   
1639def SchulzZimmDist(x,pos,args):
1640    ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
1641    ref: http://goldbook.iupac.org/S05502.html
1642    param float x: independent axis (can be numpy array)
1643    param float pos: location of distribution
1644    param float scale: width of distribution (sigma)
1645    param float shape: not used
1646    returns float: Schulz-Zimm distribution
1647    '''
1648    scale = args[0]
1649    b = (2.*pos/scale)**2
1650    a = b/pos
1651    if b<70:    #why bother?
1652        return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
1653    else:
1654        return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
1655           
1656def LogNormalCume(x,pos,args):
1657    ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
1658    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
1659    param float x: independent axis (can be numpy array)
1660    param float pos: location of distribution
1661    param float scale: width of distribution (sigma)
1662    param float shape: shape parameter
1663    returns float: LogNormal cumulative distribution
1664    '''
1665    scale,shape = args
1666    lredX = np.log((x-pos)/scale)
1667    return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
1668   
1669def GaussCume(x,pos,args):
1670    ''' Standard Normal cumulative distribution - numpy friendly on x axis
1671    param float x: independent axis (can be numpy array)
1672    param float pos: location of distribution
1673    param float scale: width of distribution (sigma)
1674    param float shape: not used
1675    returns float: Normal cumulative distribution
1676    '''
1677    scale = args[0]
1678    redX = (x-pos)/scale
1679    return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
1680   
1681def LSWCume(x,pos,args=[]):
1682    ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
1683    param float x: independent axis (can be numpy array)
1684    param float pos: location of distribution
1685    param float scale: not used
1686    param float shape: not used
1687    returns float: LSW cumulative distribution
1688    '''
1689    nP = 500
1690    xMin,xMax = [0.,2.*pos]
1691    X = np.linspace(xMin,xMax,nP,True)
1692    fxn = LSWDist(X,pos)
1693    mat = np.outer(np.ones(nP),fxn)
1694    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
1695    return np.interp(x,X,cume,0,1)
1696   
1697def SchulzZimmCume(x,pos,args):
1698    ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
1699    param float x: independent axis (can be numpy array)
1700    param float pos: location of distribution
1701    param float scale: width of distribution (sigma)
1702    param float shape: not used
1703    returns float: Normal distribution
1704    '''
1705    scale = args[0]
1706    nP = 500
1707    xMin = np.max([0.,pos-4.*scale])
1708    xMax = np.min([pos+4.*scale,1.e5])
1709    X = np.linspace(xMin,xMax,nP,True)
1710    fxn = LSWDist(X,pos)
1711    mat = np.outer(np.ones(nP),fxn)
1712    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
1713    return np.interp(x,X,cume,0,1)
1714   
1715    return []
1716   
1717         
1718################################################################################
1719##### SB-MaxEnt
1720################################################################################
1721
1722def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
1723    '''Calculates the response matrix :math:`G(Q,r)`
1724   
1725    :param float q: :math:`Q`
1726    :param float r: :math:`r`
1727    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
1728    :param function FFfxn: form factor function FF(q,r,args)
1729    :param function Volfxn: volume function Vol(r,args)
1730    :returns float: G(Q,r)
1731    '''
1732    FF = FFfxn(q,r,args)
1733    Vol = Volfxn(r,args)
1734    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
1735   
1736'''
1737sbmaxent
1738
1739Entropy maximization routine as described in the article
1740J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
1741("MNRAS": "Monthly Notices of the Royal Astronomical Society")
1742
1743:license: Copyright (c) 2013, UChicago Argonne, LLC
1744:license: This file is distributed subject to a Software License Agreement found
1745     in the file LICENSE that is included with this distribution.
1746
1747References:
1748
17491. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
17502. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
1751   Neutron Scattering Data Analysis, Rutherford
1752   Appleton Laboratory, UK, 1986; ed. MW Johnson,
1753   IOP Conference Series 81 (1986) 81 - 86, Institute
1754   of Physics, Bristol, UK.
17553. ID Culverwell and GP Clarke; Ibid. 87 - 96.
17564. JA Potton, GK Daniell, & BD Rainford,
1757   J APPL CRYST 21 (1988) 663 - 668.
17585. JA Potton, GJ Daniell, & BD Rainford,
1759   J APPL CRYST 21 (1988) 891 - 897.
1760
1761'''
1762
1763class MaxEntException(Exception): 
1764    '''Any exception from this module'''
1765    pass
1766
1767def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
1768    '''
1769    do the complete Maximum Entropy algorithm of Skilling and Bryan
1770   
1771    :param float datum[]:
1772    :param float sigma[]:
1773    :param float[][] G: transformation matrix
1774    :param float base[]:
1775    :param int IterMax:
1776    :param obj image_to_data: opus function (defaults to opus)
1777    :param obj data_to_image: tropus function (defaults to tropus)
1778   
1779    :returns float[]: :math:`f(r) dr`
1780    '''
1781       
1782    TEST_LIMIT        = 0.05                    # for convergence
1783    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
1784    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
1785    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
1786    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
1787   
1788    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
1789    MOVE_PASSES       = 0.001                   # convergence test in routine: move
1790
1791    def tropus (data, G):
1792        '''
1793        tropus: transform data-space -> solution-space:  [G] * data
1794       
1795        default definition, caller can use this definition or provide an alternative
1796       
1797        :param float[M] data: observations, ndarray of shape (M)
1798        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
1799        :returns float[N]: calculated image, ndarray of shape (N)
1800        '''
1801        return G.dot(data)
1802
1803    def opus (image, G):
1804        '''
1805        opus: transform solution-space -> data-space:  [G]^tr * image
1806       
1807        default definition, caller can use this definition or provide an alternative
1808       
1809        :param float[N] image: solution, ndarray of shape (N)
1810        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
1811        :returns float[M]: calculated data, ndarray of shape (M)
1812        '''
1813        return np.dot(G.T,image)    #G.transpose().dot(image)
1814
1815    def Dist(s2, beta):
1816        '''measure the distance of this possible solution'''
1817        w = 0
1818        n = beta.shape[0]
1819        for k in range(n):
1820            z = -sum(s2[k] * beta)
1821            w += beta[k] * z
1822        return w
1823   
1824    def ChiNow(ax, c1, c2, s1, s2):
1825        '''
1826        ChiNow
1827       
1828        :returns tuple: (ChiNow computation of ``w``, beta)
1829        '''
1830       
1831        bx = 1 - ax
1832        a =   bx * c2  -  ax * s2
1833        b = -(bx * c1  -  ax * s1)
1834   
1835        beta = ChoSol(a, b)
1836        w = 1.0
1837        for k in range(SEARCH_DIRECTIONS):
1838            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
1839        return w, beta
1840   
1841    def ChoSol(a, b):
1842        '''
1843        ChoSol: ? chop the solution vectors ?
1844       
1845        :returns: new vector beta
1846        '''
1847        n = b.shape[0]
1848        fl = np.zeros((n,n))
1849        bl = np.zeros_like(b)
1850       
1851        #print_arr("ChoSol: a", a)
1852        #print_vec("ChoSol: b", b)
1853   
1854        if (a[0][0] <= 0):
1855            msg = "ChoSol: a[0][0] = " 
1856            msg += str(a[0][0])
1857            msg += '  Value must be positive'
1858            raise MaxEntException(msg)
1859   
1860        # first, compute fl from a
1861        # note fl is a lower triangular matrix
1862        fl[0][0] = math.sqrt (a[0][0])
1863        for i in (1, 2):
1864            fl[i][0] = a[i][0] / fl[0][0]
1865            for j in range(1, i+1):
1866                z = 0.0
1867                for k in range(j):
1868                    z += fl[i][k] * fl[j][k]
1869                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
1870                z = a[i][j] - z
1871                if j == i:
1872                    y = math.sqrt(z)
1873                else:
1874                    y = z / fl[j][j]
1875                fl[i][j] = y
1876        #print_arr("ChoSol: fl", fl)
1877   
1878        # next, compute bl from fl and b
1879        bl[0] = b[0] / fl[0][0]
1880        for i in (1, 2):
1881            z = 0.0
1882            for k in range(i):
1883                z += fl[i][k] * bl[k]
1884                #print "\t", i, k, z
1885            bl[i] = (b[i] - z) / fl[i][i]
1886        #print_vec("ChoSol: bl", bl)
1887   
1888        # last, compute beta from bl and fl
1889        beta = np.empty((n))
1890        beta[-1] = bl[-1] / fl[-1][-1]
1891        for i in (1, 0):
1892            z = 0.0
1893            for k in range(i+1, n):
1894                z += fl[k][i] * beta[k]
1895                #print "\t\t", i, k, 'z=', z
1896            beta[i] = (bl[i] - z) / fl[i][i]
1897        #print_vec("ChoSol: beta", beta)
1898   
1899        return beta
1900
1901    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
1902        '''
1903        move beta one step closer towards the solution
1904        '''
1905        a_lower, a_upper = 0., 1.          # bracket  "a"
1906        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
1907        #print "MaxEntMove: cmin = %g" % cmin
1908        if cmin*chisq > chizer:
1909            ctarg = (1.0 + cmin)/2
1910        else:
1911            ctarg = chizer/chisq
1912        f_lower = cmin - ctarg
1913        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
1914        f_upper = c_upper - ctarg
1915   
1916        fx = 2*MOVE_PASSES      # just to start off
1917        loop = 1
1918        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
1919            a_new = (a_lower + a_upper) * 0.5           # search by bisection
1920            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
1921            fx = c_new - ctarg
1922            # tighten the search range for the next pass
1923            if f_lower*fx > 0:
1924                a_lower, f_lower = a_new, fx
1925            if f_upper*fx > 0:
1926                a_upper, f_upper = a_new, fx
1927            loop += 1
1928   
1929        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
1930            msg = "MaxEntMove: Loop counter = " 
1931            msg += str(MAX_MOVE_LOOPS)
1932            msg += '  No convergence in alpha chop'
1933            raise MaxEntException(msg)
1934   
1935        w = Dist (s2, beta);
1936        m = SEARCH_DIRECTIONS
1937        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
1938            for k in range(m):
1939                beta[k] *= math.sqrt (fSum/(blank*w))
1940        chtarg = ctarg * chisq
1941        return w, chtarg, loop, a_new, fx, beta
1942       
1943#MaxEnt_SB starts here   
1944
1945    if image_to_data == None:
1946        image_to_data = opus
1947    if data_to_image == None:
1948        data_to_image = tropus
1949    n   = len(base)
1950    npt = len(datum)
1951
1952    # Note that the order of subscripts for
1953    # "xi" and "eta" has been reversed from
1954    # the convention used in the FORTRAN version
1955    # to enable parts of them to be passed as
1956    # as vectors to "image_to_data" and "data_to_image".
1957    xi      = np.zeros((SEARCH_DIRECTIONS, n))
1958    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
1959    beta    = np.zeros((SEARCH_DIRECTIONS))
1960    # s1      = np.zeros((SEARCH_DIRECTIONS))
1961    # c1      = np.zeros((SEARCH_DIRECTIONS))
1962    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
1963    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
1964
1965    # TODO: replace blank (scalar) with base (vector)
1966    blank = sum(base) / len(base)   # use the average value of base
1967
1968    chizer, chtarg = npt*1.0, npt*1.0
1969    f = base * 1.0                              # starting distribution is base
1970
1971    fSum  = sum(f)                              # find the sum of the f-vector
1972    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
1973    chisq = sum(z*z)                            # Chi^2, SB eq. 4
1974
1975    for iter in range(IterMax):
1976        ox = -2 * z / sigma                        # gradient of Chi^2
1977
1978        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
1979        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
1980        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
1981        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
1982        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
1983
1984        a = 1.0
1985        b = 1.0 / cnorm
1986        if iter == 0:
1987            test = 0.0     # mismatch between entropy and ChiSquared gradients
1988        else:
1989            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
1990            a = 0.5 / (snorm * test)
1991            b *= 0.5 / test
1992        xi[0] = f * cgrad / cnorm
1993        xi[1] = f * (a * sgrad - b * cgrad)
1994
1995        eta[0] = image_to_data (xi[0], G);          # image --> data
1996        eta[1] = image_to_data (xi[1], G);          # image --> data
1997        ox = eta[1] / (sigma * sigma)
1998        xi[2] = data_to_image (ox, G);              # data --> image
1999        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
2000        xi[2] = f * xi[2] * a
2001        eta[2] = image_to_data (xi[2], G)           # image --> data
2002       
2003#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
2004#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
2005
2006        # prepare the search directions for the conjugate gradient technique
2007        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
2008        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
2009#         print_vec("MaxEnt: c1", c1)
2010#         print_vec("MaxEnt: s1", s1)
2011
2012        for k in range(SEARCH_DIRECTIONS):
2013            for l in range(k+1):
2014                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
2015                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
2016#         print_arr("MaxEnt: c2", c2)
2017#         print_arr("MaxEnt: s2", s2)
2018
2019        # reflect across the body diagonal
2020        for k, l in ((0,1), (0,2), (1,2)):
2021            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
2022            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
2023 
2024        beta[0] = -0.5 * c1[0] / c2[0][0]
2025        beta[1] = 0.0
2026        beta[2] = 0.0
2027        if (iter > 0):
2028            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
2029
2030        f_old = f.copy()    # preserve the last image
2031        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
2032       
2033        # As mentioned at the top of p.119,
2034        # need to protect against stray negative values.
2035        # In this case, set them to RESET_STRAYS * base[i]
2036        #f = f.clip(RESET_STRAYS * blank, f.max())
2037        for i in range(n):
2038            if f[i] <= 0.0:
2039                f[i] = RESET_STRAYS * base[i]
2040        df = f - f_old
2041        fSum = sum(f)
2042        fChange = sum(df)
2043
2044        # calculate the normalized entropy
2045        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
2046        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
2047        chisq = sum(z*z)                            # report this ChiSq
2048
2049        if report:
2050            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
2051            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
2052            if iter > 0:
2053                value = 100*( math.sqrt(chisq/chtarg)-1)
2054            else:
2055                value = 0
2056      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
2057            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
2058
2059        # See if we have finished our task.
2060        # do the hardest test first
2061        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
2062            print ' Convergence achieved.'
2063            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
2064    print ' No convergence! Try increasing Error multiplier.'
2065    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
2066
2067   
2068###############################################################################
2069#### IPG/TNNLS Routines
2070###############################################################################
2071
2072def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
2073    ''' An implementation of the Interior-Point Gradient method of
2074    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
2075    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
2076    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
2077    Problem addressed: Total Non-Negative Least Squares (TNNLS)
2078    :param float datum[]:
2079    :param float sigma[]:
2080    :param float[][] G: transformation matrix
2081    :param int IterMax:
2082    :param float Qvec: data positions for Power = 0-4
2083    :param float approach: 0.8 default fitting parameter
2084    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
2085   
2086    '''
2087    if Power < 0: 
2088        GmatE = G/sigma[:np.newaxis]
2089        IntE = datum/sigma
2090        pwr = 0
2091        QvecP = np.ones_like(datum)
2092    else:
2093        GmatE = G[:]
2094        IntE = datum[:]
2095        pwr = Power
2096        QvecP = Qvec**pwr
2097    Amat = GmatE*QvecP[:np.newaxis]
2098    AAmat = np.inner(Amat,Amat)
2099    Bvec = IntE*QvecP
2100    Xw = np.ones_like(Bins)*1.e-32
2101    calc = np.dot(G.T,Xw)
2102    nIter = 0
2103    err = 10.
2104    while (nIter<IterMax) and (err > 1.):
2105        #Step 1 in M&Z paper:
2106        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
2107        Dk = Xw/np.inner(AAmat,Xw)
2108        Pk = -Dk*Qk
2109        #Step 2 in M&Z paper:
2110        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
2111        alpWv = -Xw/Pk
2112        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
2113        #Step 3 in M&Z paper:
2114        shift = AkSt*Pk
2115        Xw += shift
2116        #done IPG; now results
2117        nIter += 1
2118        calc = np.dot(G.T,Xw)
2119        chisq = np.sum(((datum-calc)/sigma)**2)
2120        err = chisq/len(datum)
2121        if report:
2122            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
2123    return chisq,Xw,calc
2124
2125###############################################################################
2126#### SASD Utilities
2127###############################################################################
2128
2129def SetScale(Data,refData):
2130    Profile,Limits,Sample = Data
2131    refProfile,refLimits,refSample = refData
2132    x,y = Profile[:2]
2133    rx,ry = refProfile[:2]
2134    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
2135    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
2136    iBeg = np.searchsorted(x,Beg)
2137    iFin = np.searchsorted(x,Fin)
2138    sum = np.sum(y[iBeg:iFin])
2139    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
2140    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
2141   
2142def Bestimate(G,Rg,P):
2143    return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
2144   
2145###############################################################################
2146#### Size distribution
2147###############################################################################
2148
2149def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
2150    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
2151        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
2152        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
2153        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
2154        'Cylinder diam':[CylinderDFF,CylinderDVol]}
2155    Shape = data['Size']['Shape'][0]
2156    Parms = data['Size']['Shape'][1:]
2157    if data['Size']['logBins']:
2158        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
2159            data['Size']['Nbins']+1,True)/2.        #make radii
2160    else:
2161        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
2162            data['Size']['Nbins']+1,True)/2.        #make radii
2163    Dbins = np.diff(Bins)
2164    Bins = Bins[:-1]+Dbins/2.
2165    Contrast = Sample['Contrast'][1]
2166    Scale = Sample['Scale'][0]
2167    Sky = 10**data['Size']['MaxEnt']['Sky']
2168    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
2169    Back = data['Back']
2170    Q,Io,wt,Ic,Ib = Profile[:5]
2171    Qmin = Limits[1][0]
2172    Qmax = Limits[1][1]
2173    wtFactor = ProfDict['wtFactor']
2174    Ibeg = np.searchsorted(Q,Qmin)
2175    Ifin = np.searchsorted(Q,Qmax)
2176    BinMag = np.zeros_like(Bins)
2177    Ic[:] = 0.
2178    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
2179    if 'MaxEnt' == data['Size']['Method']:
2180        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
2181            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
2182            data['Size']['MaxEnt']['Niter'],report=True)
2183    elif 'IPG' == data['Size']['Method']:
2184        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
2185            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
2186            Power=data['Size']['IPG']['Power'],report=True)
2187    Ib[:] = Back[0]
2188    Ic[Ibeg:Ifin] += Back[0]
2189    print ' Final chi^2: %.3f'%(chisq)
2190    Vols = shapes[Shape][1](Bins,Parms)
2191    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
2192       
2193################################################################################
2194#### Modelling
2195################################################################################
2196
2197def ModelFit(Profile,ProfDict,Limits,Substances,Sample,Model):
2198    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
2199        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
2200        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
2201        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
2202        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
2203
2204    def GetModelParms():
2205        parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
2206            'PkInt','PkPos','PkSig','PkGam',]
2207        parmDict = {}
2208        varyList = []
2209        values = []
2210        levelTypes = []
2211        Back = Model['Back']
2212        if Back[1]:
2213            varyList += ['Back',]
2214            values.append(Back[0])
2215        parmDict['Back'] = Back[0]
2216        partData = Model['Particle']
2217        parmDict['Matrix density'] = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
2218        for i,level in enumerate(partData['Levels']):
2219            cid = str(i)+':'
2220            controls = level['Controls']
2221            Type = controls['DistType']
2222            levelTypes.append(Type)
2223            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
2224                if Type not in ['Monodosperse',]:
2225                    parmDict[cid+'NumPoints'] = controls['NumPoints']
2226                    parmDict[cid+'Cutoff'] = controls['Cutoff']
2227                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
2228                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
2229                parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
2230                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
2231                    if item in controls['FFargs']:
2232                        parmDict[cid+item] = controls['FFargs'][item][0]
2233                        if controls['FFargs'][item][1]:
2234                            varyList.append(cid+item)
2235                            values.append(controls['FFargs'][item][0])
2236            distDict = controls['DistType']
2237            for item in parmOrder:
2238                if item in level[distDict]:
2239                    parmDict[cid+item] = level[distDict][item][0]
2240                    if level[distDict][item][1]:
2241                        values.append(level[distDict][item][0])
2242                        varyList.append(cid+item)
2243        return levelTypes,parmDict,varyList,values
2244       
2245    def SetModelParms():
2246        print ' Refined parameters:'
2247        if 'Back' in varyList:
2248            Model['Back'][0] = parmDict['Back']
2249            print %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
2250        partData = Model['Particle']
2251        for i,level in enumerate(partData['Levels']):
2252            Type = level['Controls']['DistType']
2253            print ' Component %d: Type: %s:'%(i,Type)
2254            cid = str(i)+':'
2255            controls = level['Controls']
2256            Type = controls['DistType']
2257            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
2258                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
2259                    if cid+item in varyList:
2260                        controls['FFargs'][item][0] = parmDict[cid+item]
2261                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
2262            distDict = controls['DistType']
2263            for item in level[distDict]:
2264                if cid+item in varyList:
2265                    level[distDict][item][0] = parmDict[cid+item]
2266                    print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
2267                   
2268    def calcSASD(values,Q,Io,wt,levelTypes,parmDict,varyList):
2269        parmDict.update(zip(varyList,values))
2270        M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)-Io)
2271        return M
2272       
2273    def getSASD(Q,levelTypes,parmDict):
2274        Ic = np.ones_like(Q)
2275        Ic *= parmDict['Back']
2276        rhoMat = parmDict['Matrix density']
2277        for i,Type in enumerate(levelTypes):
2278            cid = str(i)+':'
2279            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
2280                FFfxn = parmDict[cid+'FormFact']
2281                Volfxn = parmDict[cid+'FFVolume']
2282                FFargs = []
2283                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
2284                    if item in parmDict: 
2285                        FFargs.append(parmDict[item])
2286                distDict = {}
2287                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
2288                    if item in parmDict:
2289                        distDict[item.split(':')[1]] = parmDict[item]
2290                rho = parmDict[cid+'XAnom density']
2291                contrast = rho**2-rhoMat**2
2292                rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
2293                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
2294                dist *= parmDict[cid+'Volume']
2295                Ic += np.dot(Gmat,dist)
2296            elif 'Unified' in Type:
2297                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
2298                    parmDict[cid+'P'],parmDict[cid+'Cutoff']
2299                Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
2300                Guin = G*np.exp(-(Q*Rg)**2/3)
2301                Porod = (B/Qstar**P)*np.exp(-(Q*Rgco)**2/3)
2302                Ic += Guin+Porod
2303            elif 'Porod' in Type:
2304                B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
2305                Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
2306                Ic += Porod
2307            elif 'Mono' in Type:
2308                FFfxn = parmDict[cid+'FormFact']
2309                Volfxn = parmDict[cid+'FFVolume']
2310                FFargs = []
2311                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
2312                    if item in parmDict: 
2313                        FFargs.append(parmDict[item])
2314                rho = parmDict[cid+'XAnom density']
2315                contrast = rho**2-rhoMat**2
2316                R = parmDict[cid+'Radius']
2317                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
2318                Ic += Gmat[0]*parmDict[cid+'Volume']
2319            elif 'Bragg' in distFxn:
2320                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
2321                    parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
2322        return Ic
2323       
2324    Q,Io,wt,Ic,Ib = Profile[:5]
2325    Qmin = Limits[1][0]
2326    Qmax = Limits[1][1]
2327    wtFactor = ProfDict['wtFactor']
2328    Ibeg = np.searchsorted(Q,Qmin)
2329    Ifin = np.searchsorted(Q,Qmax)
2330    Ic[:] = 0
2331    levelTypes,parmDict,varyList,values = GetModelParms()
2332    result = so.leastsq(calcSASD,values,full_output=True,   #ftol=Ftol,
2333        args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],levelTypes,parmDict,varyList))
2334    ncyc = int(result[2]['nfev']/2)
2335    chisq = np.sum(result[2]['fvec']**2)
2336    Rvals = {}
2337    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
2338    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
2339    parmDict.update(zip(varyList,result[0]))
2340    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
2341    try:
2342        sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
2343        sigDict = dict(zip(varyList,sig))
2344        print ' Results of small angle data modelling fit:'
2345        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
2346        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
2347        SetModelParms()
2348        covMatrix = result[1]*Rvals['GOF']
2349        return True,result,varyList,sig,Rvals,covMatrix
2350    except ValueError:
2351        return False,0,0,0,0,0
2352   
2353def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
2354    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
2355        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
2356        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
2357        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
2358        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
2359#    pdb.set_trace()
2360    partData = sasdData['Particle']
2361    rhoMat = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
2362    matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
2363    Scale = Sample['Scale']     #[value,flag]
2364    Back = sasdData['Back']
2365    Q,Io,wt,Ic,Ib = Profile[:5]
2366    Qmin = Limits[1][0]
2367    Qmax = Limits[1][1]
2368    wtFactor = ProfDict['wtFactor']
2369    Ibeg = np.searchsorted(Q,Qmin)
2370    Ifin = np.searchsorted(Q,Qmax)
2371    Ib[:] = Back[0]
2372    Ic[:] = 0
2373    Ic[Ibeg:Ifin] += Back[0]
2374    Rbins = []
2375    Dist = []
2376    for level in partData['Levels']:
2377        controls = level['Controls']
2378        distFxn = controls['DistType']
2379        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
2380            FFfxn = shapes[controls['FormFact']][0]
2381            Volfxn = shapes[controls['FormFact']][1]
2382            FFargs = []
2383            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
2384                if item in controls['FFargs']: 
2385                    FFargs.append(controls['FFargs'][item][0])
2386            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
2387            contrast = rho**2-rhoMat**2
2388            parmDict = level[controls['DistType']]
2389            distDict = {}
2390            for item in parmDict:
2391                distDict[item] = parmDict[item][0]
2392            rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
2393            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
2394            dist *= level[distFxn]['Volume'][0]
2395            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)
2396            Rbins.append(rBins)
2397            Dist.append(dist/(4.*dBins))
2398        elif 'Unified' in distFxn:
2399            parmDict = level[controls['DistType']]
2400            Rg,G,B,P,Rgco = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],    \
2401                parmDict['P'][0],parmDict['Cutoff'][0]
2402            Qstar = Q[Ibeg:Ifin]/(scsp.erf(Q[Ibeg:Ifin]*Rg/np.sqrt(6)))**3
2403            Guin = G*np.exp(-(Q[Ibeg:Ifin]*Rg)**2/3)
2404            Porod = (B/Qstar**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
2405            Ic[Ibeg:Ifin] += Guin+Porod
2406            Rbins.append([])
2407            Dist.append([])
2408        elif 'Porod' in distFxn:
2409            parmDict = level[controls['DistType']]
2410            B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
2411            Porod = (B/Q[Ibeg:Ifin]**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
2412            Ic[Ibeg:Ifin] += Porod
2413            Rbins.append([])
2414            Dist.append([])
2415        elif 'Mono' in distFxn:
2416            FFfxn = shapes[controls['FormFact']][0]
2417            Volfxn = shapes[controls['FormFact']][1]
2418            FFargs = []
2419            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
2420                if item in controls['FFargs']: 
2421                    FFargs.append(controls['FFargs'][item][0])
2422            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
2423            contrast = rho**2-rhoMat**2
2424            R = level[controls['DistType']]['Radius'][0]
2425            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
2426            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]
2427            Rbins.append([])
2428            Dist.append([])
2429        elif 'Bragg' in distFxn:
2430            parmDict = level[controls['DistType']]
2431            Ic[Ibeg:Ifin] += parmDict['PkInt'][0]*G2pwd.getPsVoigt(parmDict['PkPos'][0],
2432                parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
2433            Rbins.append([])
2434            Dist.append([])           
2435    sasdData['Size Calc'] = [Rbins,Dist]
2436   
2437def MakeDiamDist(DistName,nPoints,cutoff,distDict):
2438   
2439    if 'LogNormal' in DistName:
2440        distFxn = 'LogNormalDist'
2441        cumeFxn = 'LogNormalCume'
2442        pos = distDict['MinSize']
2443        args = [distDict['Mean'],distDict['StdDev']]
2444        step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))
2445        mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)
2446        minX = 1. #pos
2447        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
2448    elif 'Gauss' in DistName:
2449        distFxn = 'GaussDist'
2450        cumeFxn = 'GaussCume'
2451        pos = distDict['Mean']
2452        args = [distDict['StdDev']]
2453        step = 0.02*distDict['StdDev']
2454        mode = distDict['Mean']
2455        minX = np.max([mode-4.*distDict['StdDev'],1.])
2456        maxX = np.min([mode+4.*distDict['StdDev'],1.e5])
2457    elif 'LSW' in DistName:
2458        distFxn = 'LSWDist'
2459        cumeFxn = 'LSWCume'
2460        pos = distDict['Mean']
2461        args = []
2462        minX,maxX = [0.,2.*pos]
2463    elif 'Schulz' in DistName:
2464        distFxn = 'SchulzZimmDist'
2465        cumeFxn = 'SchulzZimmCume'
2466        pos = distDict['Mean']
2467        args = [distDict['StdDev']]
2468        minX = np.max([1.,pos-4.*distDict['StdDev']])
2469        maxX = np.min([pos+4.*distDict['StdDev'],1.e5])
2470    nP = 500
2471    Diam = np.logspace(0.,5.,nP,True)
2472    TCW = eval(cumeFxn+'(Diam,pos,args)')
2473    CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints+1,True)
2474    rBins = np.interp(CumeTgts,TCW,Diam,0,0)
2475    dBins = np.diff(rBins)
2476    rBins = rBins[:-1]+dBins/2.
2477    return rBins,dBins,eval(distFxn+'(rBins,pos,args)')
2478
2479################################################################################
2480#### MaxEnt testing stuff
2481################################################################################
2482
2483def print_vec(text, a):
2484    '''print the contents of a vector to the console'''
2485    n = a.shape[0]
2486    print "%s[ = (" % text,
2487    for i in range(n):
2488        s = " %g, " % a[i]
2489        print s,
2490    print ")"
2491
2492def print_arr(text, a):
2493    '''print the contents of an array to the console'''
2494    n, m = a.shape
2495    print "%s[][] = (" % text
2496    for i in range(n):
2497        print " (",
2498        for j in range(m):
2499            print " %g, " % a[i][j],
2500        print "),"
2501    print ")"
2502
2503def test_MaxEnt_SB(report=True):
2504    def readTextData(filename):
2505        '''return q, I, dI from a 3-column text file'''
2506        if not os.path.exists(filename):
2507            raise Exception("file not found: " + filename)
2508        buf = [line.split() for line in open(filename, 'r').readlines()]
2509        M = len(buf)
2510        buf = zip(*buf)         # transpose rows and columns
2511        q  = np.array(buf[0], dtype=np.float64)
2512        I  = np.array(buf[1], dtype=np.float64)
2513        dI = np.array(buf[2], dtype=np.float64)
2514        return q, I, dI
2515    print "MaxEnt_SB: "
2516    test_data_file = os.path.join( 'testinp', 'test.sas')
2517    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
2518    bkg   = 0.1     #   I = I - bkg
2519    dMin, dMax, nRadii = 25, 9000, 40
2520    defaultDistLevel = 1.0e-6
2521    IterMax = 40
2522    errFac = 1.05
2523   
2524    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
2525    dr   = r * (r[1]/r[0] - 1)          # step size
2526    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
2527    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
2528   
2529    qVec, I, dI = readTextData(test_data_file)
2530    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
2531   
2532    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
2533    if f_dr is None:
2534        print "no solution"
2535        return
2536   
2537    print "solution reached"
2538    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
2539        print '%10.4f %10.4f %12.4g'%(a,b,c)
2540
2541def tests():
2542    test_MaxEnt_SB(report=True)
2543
2544if __name__ == '__main__':
2545    tests()
2546
2547>>>>>>> .r1299
Note: See TracBrowser for help on using the repository browser.