Changeset 1244 for trunk/GSASIIsasd.py

Timestamp:

Mar 11, 2014 4:49:15 PM (9 years ago)

Author:

vondreele

Message:

fix geometric correction in integrate - too many 1/cos(2-theta)
plot of size distribution from SASD
MaxEnt? size distribution in operation (some tuning/errors)

File:

• trunk/GSASIIsasd.py (modified) (12 diffs)

trunk/GSASIIsasd.py

@@ -57 +57 @@
     
 ###############################################################################
-#### Particle form factors & volumes as class definitions
-###############################################################################
-
-class SASDParticles(object):
-    def __init__(self,name=None,parNames=None):
-        self.name = name
-        self.ParmNames = parmNames
-        
-    def FormFactor():
-        FF = 0.
-        return FF
-        
-    def Volume():
-        Vol = 0.
-        return Vol
-        
-class Sphere(SASDParticles):
-    def __init__(self,name=None,parmNames=None):
-        self.Name = name
-        if self.Name == None:
-            self.Name = 'Sphere'
-        self.ParmNames = parmNames
-        if self.ParmNames == None:
-            self.ParmNames = ['Radius',]
-        
-    def FormFactor(self,Q,Parms):
-        ''' Compute hard sphere form factor - can use numpy arrays
-        param float:Q Q value array (usually in A-1)
-        param float:R sphere radius (Usually in A - must match Q-1 units)
-        returns float: form factors as array as needed
-        '''
-        QR = Q*Parms[0]
-        return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
-        
-    def Volume(self,Parms):
-        ''' Compute volume of sphere
-        - numpy array friendly
-        param float:R sphere radius
-        returns float: volume
-        '''
-        return (4./3.)*np.pi*Parms[0]**3
-        
-class Spheroid(SASDParticles):
-    def __init__(self,name=None,parmNames=None):
-        self.Name = name
-        if self.Name == None:
-            self.Name = 'Spheroid'
-        self.ParmNames = parmNames
-        if self.ParmNames == None:
-            self.ParmNames = ['Radius','Aspect ratio']
-    
-    def FormFactor(self,Q,Parms):
-        ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid) 
-        - can use numpy arrays for R & AR; will return corresponding numpy array
-        param float:Q Q value array (usually in A-1)
-        param float R: radius along 2 axes of spheroid
-        param float AR: aspect ratio so 3rd axis = R*AR
-        returns float: form factors as array as needed
-        '''
-        R,AR = Parms
-        NP = 50
-        if 0.99 < AR < 1.01:
-            return SphereFF(Q,R,0)
-        else:
-            cth = np.linspace(0,1.,NP)
-            Rct = R*np.sqrt(1.+(AR**2-1.)*cth**2)
-            return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
-        
-    def Volume(self,Parms):
-        ''' Compute volume of cylindrically symmetric ellipsoid (spheroid) 
-        - numpy array friendly
-        param float R: radius along 2 axes of spheroid
-        param float AR: aspect ratio so radius of 3rd axis = R*AR
-        returns float: volume
-        '''
-        R,AR = Parms
-        return AR*(4./3.)*np.pi*R**3
-        
-###############################################################################
 #### Particle form factors
 ###############################################################################
 
-def SphereFF(Q,R,dummy=0):
+def SphereFF(Q,R,args=()):
     ''' Compute hard sphere form factor - can use numpy arrays
     param float:Q Q value array (usually in A-1)
@@ -145 +66 @@
     returns float: form factors as array as needed
     '''
-    QR = Q*R
+    QR = Q[:,np.newaxis]*R
     return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
     
-def SpheroidFF(Q,R,AR):
+def SpheroidFF(Q,R,args):
     ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid) 
     - can use numpy arrays for R & AR; will return corresponding numpy array
@@ -157 +78 @@
     '''
     NP = 50
+    AR = args[0]
     if 0.99 < AR < 1.01:
         return SphereFF(Q,R,0)
@@ -164 +86 @@
         return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
             
-def CylinderFF(Q,R,L):
+def CylinderFF(Q,R,args):
     ''' Compute form factor for cylinders - can use numpy arrays
     param float: Q Q value array (A-1)
@@ -171 +93 @@
     returns float: form factor
     '''
+    L = args[0]
     NP = 200
     alp = np.linspace(0,np.pi/2.,NP)
@@ -179 +102 @@
     return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*(LBess*SBess)**2,axis=1)/NP)
     
-def CylinderDFF(Q,L,D):
+def CylinderDFF(Q,L,args):
     ''' Compute form factor for cylinders - can use numpy arrays
     param float: Q Q value array (A-1)
@@ -186 +109 @@
     returns float: form factor
     '''
+    D = args[0]
     return CylinderFF(Q,D/2.,L)    
     
-def CylinderARFF(Q,R,AR): 
+def CylinderARFF(Q,R,args): 
     ''' Compute form factor for cylinders - can use numpy arrays
     param float: Q Q value array (A-1)
@@ -195 +119 @@
     returns float: form factor
     '''
+    AR = args[0]
     return CylinderFF(Q,R,2.*R*AR)    
     
-def UniSphereFF(Q,R,dummy=0):
+def UniSphereFF(Q,R,args=0):
     Rg = np.sqrt(3./5.)*R
     B = 1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
@@ -246 +171 @@
     FF += B1/QstV**4
     return np.sqrt(FF)
-    
-    
-    
-    
 
 ###############################################################################
@@ -255 +176 @@
 ###############################################################################
 
-def SphereVol(R):
+def SphereVol(R,arg=()):
     ''' Compute volume of sphere
     - numpy array friendly
@@ -298 +219 @@
     '''
     return CylinderVol(R,2.*R*AR)
-        
-    
+    
+def UniSphereVol(R,arg=()):
+    ''' Compute volume of sphere
+    - numpy array friendly
+    param float:R sphere radius
+    returns float: volume
+    '''
+    return SphereVol(R)
+    
+def UniRodVol(R,L):
+    ''' Compute cylinder volume for radius & length
+    - numpy array friendly
+    param float: R diameter (A)
+    param float: L length (A)
+    returns float:volume (A^3)
+    '''
+    return CylinderVol(R,L)
+    
+def UniRodARVol(R,AR):
+    return CylinderARVol(R,AR)
+    
+def UniDiskVol(R,T):
+    return CylinderVol(R,T)
+    
+def UniTubeVol(R,L,T):
+    ''' Compute tube volume for radius, length & wall thickness
+    - numpy array friendly
+    param float: R diameter (A)
+    param float: L length (A)
+    param float: T tube wall thickness (A)
+    returns float: volume (A^3) of tube wall
+    '''
+    return CylinderVol(R,L)-CylinderVol(R-T,L)
+          
+################################################################################
+##### SB-MaxEnt
+################################################################################
+
+def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
+    '''Calculates the response matrix :math:`G(Q,r)` 
+    
+    :param float q: :math:`Q`
+    :param float r: :math:`r`
+    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
+    :param function FFfxn: form factor function FF(q,r,args)
+    :param function Volfxn: volume function Vol(r,args)
+    :returns float: G(Q,r)
+    '''
+    FF = FFfxn(q,r,args)
+    Vol = Volfxn(r,args)
+    return 1.e-4*(contrast*Vol*FF**2)     #10^-20 vs 10^-24
+    
+'''
+sbmaxent
+
+Entropy maximization routine as described in the article
+J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
+("MNRAS": "Monthly Notices of the Royal Astronomical Society")
+
+:license: Copyright (c) 2013, UChicago Argonne, LLC
+:license: This file is distributed subject to a Software License Agreement found
+     in the file LICENSE that is included with this distribution. 
+
+References:
+
+1. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
+2. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
+   Neutron Scattering Data Analysis, Rutherford
+   Appleton Laboratory, UK, 1986; ed. MW Johnson,
+   IOP Conference Series 81 (1986) 81 - 86, Institute
+   of Physics, Bristol, UK.
+3. ID Culverwell and GP Clarke; Ibid. 87 - 96.
+4. JA Potton, GK Daniell, & BD Rainford,
+   J APPL CRYST 21 (1988) 663 - 668.
+5. JA Potton, GJ Daniell, & BD Rainford,
+   J APPL CRYST 21 (1988) 891 - 897.
+
+'''
+
+import os
+import sys
+import math
+import numpy
+
+class MaxEntException(Exception): 
+    '''Any exception from this module'''
+    pass
+
+def MaxEnt_SB(datum, sigma, base, IterMax, G, image_to_data=None, data_to_image=None, report=False):
+    '''
+    do the complete Maximum Entropy algorithm of Skilling and Bryan
+    
+    :param float datum[]:
+    :param float sigma[]:
+    :param float base[]:
+    :param int IterMax:
+    :param float[][] G: transformation matrix
+    :param obj image_to_data: opus function (defaults to opus)
+    :param obj data_to_image: tropus function (defaults to tropus)
+    
+    :returns float[]: :math:`f(r) dr`
+    '''
+        
+    TEST_LIMIT        = 0.10                    # for convergence
+    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
+    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
+    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
+    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
+    
+    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move, 
+    MOVE_PASSES       = 0.001                   # convergence test in routine: move
+
+    def opus (data, G):
+        '''
+        opus: transform data-space -> solution-space:  [G] * data
+        
+        default definition, caller can use this definition or provide an alternative
+        
+        :param float[M] data: observations, ndarray of shape (M)
+        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
+        :returns float[N]: calculated image, ndarray of shape (N)
+        '''
+        return G.dot(data)
+
+    def tropus (image, G):
+        '''
+        tropus: transform solution-space -> data-space:  [G]^tr * image
+        
+        default definition, caller can use this definition or provide an alternative
+        
+        :param float[N] image: solution, ndarray of shape (N)
+        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
+        :returns float[M]: calculated data, ndarray of shape (M)
+        '''
+        return G.transpose().dot(image)
+
+    def Dist(s2, beta):
+        '''measure the distance of this possible solution'''
+        w = 0
+        n = beta.shape[0]
+        for k in range(n):
+            z = -sum(s2[k] * beta)
+            w += beta[k] * z
+        return w
+    
+    def ChiNow(ax, c1, c2, s1, s2):
+        '''
+        ChiNow
+        
+        :returns tuple: (ChiNow computation of ``w``, beta)
+        '''
+        
+        bx = 1 - ax
+        a =   bx * c2  -  ax * s2
+        b = -(bx * c1  -  ax * s1)
+    
+        beta = ChoSol(a, b)
+        w = 1.0
+        for k in range(SEARCH_DIRECTIONS):
+            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
+        return w, beta
+    
+    def ChoSol(a, b):
+        '''
+        ChoSol: ? chop the solution vectors ?
+        
+        :returns: new vector beta
+        '''
+        n = b.shape[0]
+        fl = numpy.ndarray((n, n))*0
+        bl = numpy.ndarray((n))*0
+        
+        #print_arr("ChoSol: a", a)
+        #print_vec("ChoSol: b", b)
+    
+        if (a[0][0] <= 0):
+            msg = "ChoSol: a[0][0] = " 
+            msg += str(a[0][0])
+            msg += '  Value must be positive'
+            raise MaxEntException(msg)
+    
+        # first, compute fl from a
+        # note fl is a lower triangular matrix
+        fl[0][0] = math.sqrt (a[0][0])
+        for i in (1, 2):
+            fl[i][0] = a[i][0] / fl[0][0]
+            for j in range(1, i+1):
+                z = 0.0
+                for k in range(j):
+                    z += fl[i][k] * fl[j][k]
+                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
+                z = a[i][j] - z
+                if j == i:
+                    y = math.sqrt(z)
+                else:
+                    y = z / fl[j][j]
+                fl[i][j] = y
+        #print_arr("ChoSol: fl", fl)
+    
+        # next, compute bl from fl and b
+        bl[0] = b[0] / fl[0][0]
+        for i in (1, 2):
+            z = 0.0
+            for k in range(i):
+                z += fl[i][k] * bl[k]
+                #print "\t", i, k, z
+            bl[i] = (b[i] - z) / fl[i][i]
+        #print_vec("ChoSol: bl", bl)
+    
+        # last, compute beta from bl and fl
+        beta = numpy.ndarray((n))
+        beta[-1] = bl[-1] / fl[-1][-1]
+        for i in (1, 0):
+            z = 0.0
+            for k in range(i+1, n):
+                z += fl[k][i] * beta[k]
+                #print "\t\t", i, k, 'z=', z
+            beta[i] = (bl[i] - z) / fl[i][i]
+        #print_vec("ChoSol: beta", beta)
+    
+        return beta
+
+    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
+        '''
+        move beta one step closer towards the solution
+        '''
+        a_lower, a_upper = 0., 1.          # bracket  "a"
+        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
+        #print "MaxEntMove: cmin = %g" % cmin
+        if cmin*chisq > chizer:
+            ctarg = (1.0 + cmin)/2
+        else:
+            ctarg = chizer/chisq
+        f_lower = cmin - ctarg
+        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
+        f_upper = c_upper - ctarg
+    
+        fx = 2*MOVE_PASSES      # just to start off
+        loop = 1
+        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
+            a_new = (a_lower + a_upper) * 0.5           # search by bisection
+            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
+            fx = c_new - ctarg
+            # tighten the search range for the next pass
+            if f_lower*fx > 0:
+                a_lower, f_lower = a_new, fx
+            if f_upper*fx > 0:
+                a_upper, f_upper = a_new, fx
+            loop += 1
+    
+        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
+            msg = "MaxEntMove: Loop counter = " 
+            msg += str(MAX_MOVE_LOOPS)
+            msg += '  No convergence in alpha chop'
+            raise MaxEntException(msg)
+    
+        w = Dist (s2, beta);
+        m = SEARCH_DIRECTIONS
+        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
+            for k in range(m):
+                beta[k] *= math.sqrt (fSum/(blank*w))
+        chtarg = ctarg * chisq
+        return w, chtarg, loop, a_new, fx, beta
+#MaxEnt_SB starts here    
+    if image_to_data == None:
+        image_to_data = opus
+    if data_to_image == None:
+        data_to_image = tropus
+    n   = len(base)
+    npt = len(datum)
+
+    # Note that the order of subscripts for
+    # "xi" and "eta" has been reversed from
+    # the convention used in the FORTRAN version
+    # to enable parts of them to be passed as
+    # as vectors to "image_to_data" and "data_to_image".
+    xi      = 0*numpy.ndarray((SEARCH_DIRECTIONS, n))
+    eta     = 0*numpy.ndarray((SEARCH_DIRECTIONS, npt))
+    beta    = 0*numpy.ndarray((SEARCH_DIRECTIONS))
+    # s1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
+    # c1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
+    s2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
+    c2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
+
+    # TODO: replace blank (scalar) with base (vector)
+    blank = sum(base) / len(base)   # use the average value of base
+
+    chizer, chtarg = npt*1.0, npt*1.0
+    f = base * 1.0                              # starting distribution is base
+
+    fSum  = sum(f)                              # find the sum of the f-vector
+    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
+    chisq = sum(z*z)                            # Chi^2, SB eq. 4
+
+    for iter in range(IterMax):
+        ox = -2 * z / sigma                        # gradient of Chi^2
+
+        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
+        sgrad = -numpy.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
+        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
+        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
+        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST 
+
+        a = 1.0
+        b = 1.0 / cnorm
+        if iter == 0:
+            test = 0.0     # mismatch between entropy and ChiSquared gradients
+        else:
+            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
+            a = 0.5 / (snorm * test)
+            b *= 0.5 / test
+        xi[0] = f * cgrad / cnorm
+        xi[1] = f * (a * sgrad - b * cgrad)
+
+        eta[0] = image_to_data (xi[0], G);          # image --> data
+        eta[1] = image_to_data (xi[1], G);          # image --> data
+        ox = eta[1] / (sigma * sigma)
+        xi[2] = data_to_image (ox, G);              # data --> image
+        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
+        xi[2] = f * xi[2] * a
+        eta[2] = image_to_data (xi[2], G)           # image --> data
+        
+#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
+#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
+
+        # prepare the search directions for the conjugate gradient technique
+        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
+        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
+#         print_vec("MaxEnt: c1", c1)
+#         print_vec("MaxEnt: s1", s1)
+
+        for k in range(SEARCH_DIRECTIONS):
+            for l in range(k+1):
+                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
+                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
+#         print_arr("MaxEnt: c2", c2)
+#         print_arr("MaxEnt: s2", s2)
+
+        # reflect across the body diagonal
+        for k, l in ((0,1), (0,2), (1,2)):
+            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
+            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
+ 
+        beta[0] = -0.5 * c1[0] / c2[0][0]
+        beta[1] = 0.0
+        beta[2] = 0.0
+        if (iter > 0):
+            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
+
+        f_old = f.copy()    # preserve the last image
+        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
+        
+        # As mentioned at the top of p.119,
+        # need to protect against stray negative values.
+        # In this case, set them to RESET_STRAYS * base[i]
+        #f = f.clip(RESET_STRAYS * blank, f.max())
+        for i in range(n):
+            if f[i] <= 0.0:
+                f[i] = RESET_STRAYS * base[i]
+        df = f - f_old
+        fSum = sum(f)
+        fChange = sum(df)
+
+        # calculate the normalized entropy
+        S = -sum((f/fSum) * numpy.log(f/fSum))      # normalized entropy, S&B eq. 1
+        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
+        chisq = sum(z*z)                            # report this ChiSq
+
+        if report:
+            print "%3d/%3d" % ((iter+1), IterMax)
+            print " %5.2lf%% %8lg" % (100*test, S)
+            if iter > 0:
+                value = 100*( math.sqrt(chisq/chtarg)-1)
+            else:
+                value = 0
+            print " %12.5lg %10.4lf" % ( math.sqrt (chtarg/npt), value )
+            print "%12.6lg %8.2lf\n" % (fSum, 100*fChange/fSum)
+
+        # See if we have finished our task.
+        # do the hardest test first
+        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
+            return f,image_to_data (f, G)     # solution FOUND returns here
+    
+    return f,image_to_data (f, G)       # no solution after IterMax iterations
+
+    
+################################################################################
+#### MaxEnt testing stuff
+################################################################################
+
+def print_vec(text, a):
+    '''print the contents of a vector to the console'''
+    n = a.shape[0]
+    print "%s[ = (" % text,
+    for i in range(n):
+        s = " %g, " % a[i]
+        print s,
+    print ")"
+
+def print_arr(text, a):
+    '''print the contents of an array to the console'''
+    n, m = a.shape
+    print "%s[][] = (" % text
+    for i in range(n):
+        print " (",
+        for j in range(m):
+            print " %g, " % a[i][j],
+        print "),"
+    print ")"
+
+def test_MaxEnt_SB(report=True):
+    def readTextData(filename):
+        '''return q, I, dI from a 3-column text file'''
+        if not os.path.exists(filename):
+            raise Exception("file not found: " + filename)
+        buf = [line.split() for line in open(filename, 'r').readlines()]
+        M = len(buf)
+        buf = zip(*buf)         # transpose rows and columns
+        q  = numpy.array(buf[0], dtype=numpy.float64)
+        I  = numpy.array(buf[1], dtype=numpy.float64)
+        dI = numpy.array(buf[2], dtype=numpy.float64)
+        return q, I, dI
+    print "MaxEnt_SB: "
+    test_data_file = os.path.join( 'testinp', 'test.sas')
+    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
+    bkg   = 0.1     #   I = I - bkg
+    dMin, dMax, nRadii = 25, 9000, 40
+    defaultDistLevel = 1.0e-6
+    IterMax = 40
+    errFac = 1.05
+    
+    r    = numpy.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
+    dr   = r * (r[1]/r[0] - 1)          # step size
+    f_dr = numpy.ndarray((nRadii)) * 0  # volume fraction histogram
+    b    = numpy.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
+    
+    qVec, I, dI = readTextData(test_data_file)
+    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
+    
+    f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
+    if f_dr is None:
+        print "no solution"
+        return
+    
+    print "solution reached"
+    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
+        print '%10.4f %10.4f %12.4g'%(a,b,c)
+
+def tests():
+    test_MaxEnt_SB(report=True)
+
+if __name__ == '__main__':
+    tests()
     
 ###############################################################################
@@ -306 +678 @@
 
 def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
+    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
+        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
+        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
+        'Unified disk':[UniDiskFF,UniDiskVol]}
+    Shape = data['Size']['Shape'][0]
+    Parms = data['Size']['Shape'][1:]
     if data['Size']['logBins']:
         Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
-            data['Size']['Nbins']+1,True)
+            data['Size']['Nbins']+1,True)/2.        #make radii
     else:
         Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
-            data['Size']['Nbins']+1,True)
+            data['Size']['Nbins']+1,True)/2.        #make radii
     Dbins = np.diff(Bins)
-    BinMag = Dbins*np.ones_like(Dbins)
+    Bins = Bins[:-1]+Dbins/2.
+    BinsBack = np.ones_like(Bins)*1.e-6
+    Back = data['Back']
+    Q,Io,wt,Ic,Ib = Profile[:5]
+    Qmin = Limits[1][0]
+    Qmax = Limits[1][1]
+    Contrast = Sample['Contrast'][1]
+    Ibeg = np.searchsorted(Q,Qmin)
+    Ifin = np.searchsorted(Q,Qmax)
+    if Back[1]:
+        Ib = Back[0]
+        Ic[Ibeg:Ifin] = Back[0]
+    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
+    BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Io[Ibeg:Ifin]-Back[0],1./np.sqrt(wt[Ibeg:Ifin]),BinsBack,
+        data['Size']['MaxEnt']['Niter'],Gmat)
+    print BinMag.shape
+    data['Size']['Distribution'] = [Bins,Dbins,BinMag]
     print np.sum(BinMag)
         
-    print data['Size']
-#    print Limits
-#    print Substances
-#    print Sample
-#    print Profile
-    
+