Changeset 1244 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Mar 11, 2014 4:49:15 PM (8 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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIsasd.py

    r1237 r1244  
    5757   
    5858###############################################################################
    59 #### Particle form factors & volumes as class definitions
    60 ###############################################################################
    61 
    62 class SASDParticles(object):
    63     def __init__(self,name=None,parNames=None):
    64         self.name = name
    65         self.ParmNames = parmNames
    66        
    67     def FormFactor():
    68         FF = 0.
    69         return FF
    70        
    71     def Volume():
    72         Vol = 0.
    73         return Vol
    74        
    75 class Sphere(SASDParticles):
    76     def __init__(self,name=None,parmNames=None):
    77         self.Name = name
    78         if self.Name == None:
    79             self.Name = 'Sphere'
    80         self.ParmNames = parmNames
    81         if self.ParmNames == None:
    82             self.ParmNames = ['Radius',]
    83        
    84     def FormFactor(self,Q,Parms):
    85         ''' Compute hard sphere form factor - can use numpy arrays
    86         param float:Q Q value array (usually in A-1)
    87         param float:R sphere radius (Usually in A - must match Q-1 units)
    88         returns float: form factors as array as needed
    89         '''
    90         QR = Q*Parms[0]
    91         return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
    92        
    93     def Volume(self,Parms):
    94         ''' Compute volume of sphere
    95         - numpy array friendly
    96         param float:R sphere radius
    97         returns float: volume
    98         '''
    99         return (4./3.)*np.pi*Parms[0]**3
    100        
    101 class Spheroid(SASDParticles):
    102     def __init__(self,name=None,parmNames=None):
    103         self.Name = name
    104         if self.Name == None:
    105             self.Name = 'Spheroid'
    106         self.ParmNames = parmNames
    107         if self.ParmNames == None:
    108             self.ParmNames = ['Radius','Aspect ratio']
    109    
    110     def FormFactor(self,Q,Parms):
    111         ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
    112         - can use numpy arrays for R & AR; will return corresponding numpy array
    113         param float:Q Q value array (usually in A-1)
    114         param float R: radius along 2 axes of spheroid
    115         param float AR: aspect ratio so 3rd axis = R*AR
    116         returns float: form factors as array as needed
    117         '''
    118         R,AR = Parms
    119         NP = 50
    120         if 0.99 < AR < 1.01:
    121             return SphereFF(Q,R,0)
    122         else:
    123             cth = np.linspace(0,1.,NP)
    124             Rct = R*np.sqrt(1.+(AR**2-1.)*cth**2)
    125             return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
    126        
    127     def Volume(self,Parms):
    128         ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
    129         - numpy array friendly
    130         param float R: radius along 2 axes of spheroid
    131         param float AR: aspect ratio so radius of 3rd axis = R*AR
    132         returns float: volume
    133         '''
    134         R,AR = Parms
    135         return AR*(4./3.)*np.pi*R**3
    136        
    137 ###############################################################################
    13859#### Particle form factors
    13960###############################################################################
    14061
    141 def SphereFF(Q,R,dummy=0):
     62def SphereFF(Q,R,args=()):
    14263    ''' Compute hard sphere form factor - can use numpy arrays
    14364    param float:Q Q value array (usually in A-1)
     
    14566    returns float: form factors as array as needed
    14667    '''
    147     QR = Q*R
     68    QR = Q[:,np.newaxis]*R
    14869    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
    14970   
    150 def SpheroidFF(Q,R,AR):
     71def SpheroidFF(Q,R,args):
    15172    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
    15273    - can use numpy arrays for R & AR; will return corresponding numpy array
     
    15778    '''
    15879    NP = 50
     80    AR = args[0]
    15981    if 0.99 < AR < 1.01:
    16082        return SphereFF(Q,R,0)
     
    16486        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=1)/NP)
    16587           
    166 def CylinderFF(Q,R,L):
     88def CylinderFF(Q,R,args):
    16789    ''' Compute form factor for cylinders - can use numpy arrays
    16890    param float: Q Q value array (A-1)
     
    17193    returns float: form factor
    17294    '''
     95    L = args[0]
    17396    NP = 200
    17497    alp = np.linspace(0,np.pi/2.,NP)
     
    179102    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*(LBess*SBess)**2,axis=1)/NP)
    180103   
    181 def CylinderDFF(Q,L,D):
     104def CylinderDFF(Q,L,args):
    182105    ''' Compute form factor for cylinders - can use numpy arrays
    183106    param float: Q Q value array (A-1)
     
    186109    returns float: form factor
    187110    '''
     111    D = args[0]
    188112    return CylinderFF(Q,D/2.,L)   
    189113   
    190 def CylinderARFF(Q,R,AR):
     114def CylinderARFF(Q,R,args):
    191115    ''' Compute form factor for cylinders - can use numpy arrays
    192116    param float: Q Q value array (A-1)
     
    195119    returns float: form factor
    196120    '''
     121    AR = args[0]
    197122    return CylinderFF(Q,R,2.*R*AR)   
    198123   
    199 def UniSphereFF(Q,R,dummy=0):
     124def UniSphereFF(Q,R,args=0):
    200125    Rg = np.sqrt(3./5.)*R
    201126    B = 1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
     
    246171    FF += B1/QstV**4
    247172    return np.sqrt(FF)
    248    
    249    
    250    
    251    
    252173
    253174###############################################################################
     
    255176###############################################################################
    256177
    257 def SphereVol(R):
     178def SphereVol(R,arg=()):
    258179    ''' Compute volume of sphere
    259180    - numpy array friendly
     
    298219    '''
    299220    return CylinderVol(R,2.*R*AR)
    300        
    301    
     221   
     222def UniSphereVol(R,arg=()):
     223    ''' Compute volume of sphere
     224    - numpy array friendly
     225    param float:R sphere radius
     226    returns float: volume
     227    '''
     228    return SphereVol(R)
     229   
     230def UniRodVol(R,L):
     231    ''' Compute cylinder volume for radius & length
     232    - numpy array friendly
     233    param float: R diameter (A)
     234    param float: L length (A)
     235    returns float:volume (A^3)
     236    '''
     237    return CylinderVol(R,L)
     238   
     239def UniRodARVol(R,AR):
     240    return CylinderARVol(R,AR)
     241   
     242def UniDiskVol(R,T):
     243    return CylinderVol(R,T)
     244   
     245def UniTubeVol(R,L,T):
     246    ''' Compute tube volume for radius, length & wall thickness
     247    - numpy array friendly
     248    param float: R diameter (A)
     249    param float: L length (A)
     250    param float: T tube wall thickness (A)
     251    returns float: volume (A^3) of tube wall
     252    '''
     253    return CylinderVol(R,L)-CylinderVol(R-T,L)
     254         
     255################################################################################
     256##### SB-MaxEnt
     257################################################################################
     258
     259def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
     260    '''Calculates the response matrix :math:`G(Q,r)`
     261   
     262    :param float q: :math:`Q`
     263    :param float r: :math:`r`
     264    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
     265    :param function FFfxn: form factor function FF(q,r,args)
     266    :param function Volfxn: volume function Vol(r,args)
     267    :returns float: G(Q,r)
     268    '''
     269    FF = FFfxn(q,r,args)
     270    Vol = Volfxn(r,args)
     271    return 1.e-4*(contrast*Vol*FF**2)     #10^-20 vs 10^-24
     272   
     273'''
     274sbmaxent
     275
     276Entropy maximization routine as described in the article
     277J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
     278("MNRAS": "Monthly Notices of the Royal Astronomical Society")
     279
     280:license: Copyright (c) 2013, UChicago Argonne, LLC
     281:license: This file is distributed subject to a Software License Agreement found
     282     in the file LICENSE that is included with this distribution.
     283
     284References:
     285
     2861. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
     2872. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
     288   Neutron Scattering Data Analysis, Rutherford
     289   Appleton Laboratory, UK, 1986; ed. MW Johnson,
     290   IOP Conference Series 81 (1986) 81 - 86, Institute
     291   of Physics, Bristol, UK.
     2923. ID Culverwell and GP Clarke; Ibid. 87 - 96.
     2934. JA Potton, GK Daniell, & BD Rainford,
     294   J APPL CRYST 21 (1988) 663 - 668.
     2955. JA Potton, GJ Daniell, & BD Rainford,
     296   J APPL CRYST 21 (1988) 891 - 897.
     297
     298'''
     299
     300import os
     301import sys
     302import math
     303import numpy
     304
     305class MaxEntException(Exception):
     306    '''Any exception from this module'''
     307    pass
     308
     309def MaxEnt_SB(datum, sigma, base, IterMax, G, image_to_data=None, data_to_image=None, report=False):
     310    '''
     311    do the complete Maximum Entropy algorithm of Skilling and Bryan
     312   
     313    :param float datum[]:
     314    :param float sigma[]:
     315    :param float base[]:
     316    :param int IterMax:
     317    :param float[][] G: transformation matrix
     318    :param obj image_to_data: opus function (defaults to opus)
     319    :param obj data_to_image: tropus function (defaults to tropus)
     320   
     321    :returns float[]: :math:`f(r) dr`
     322    '''
     323       
     324    TEST_LIMIT        = 0.10                    # for convergence
     325    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
     326    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
     327    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
     328    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
     329   
     330    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
     331    MOVE_PASSES       = 0.001                   # convergence test in routine: move
     332
     333    def opus (data, G):
     334        '''
     335        opus: transform data-space -> solution-space:  [G] * data
     336       
     337        default definition, caller can use this definition or provide an alternative
     338       
     339        :param float[M] data: observations, ndarray of shape (M)
     340        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
     341        :returns float[N]: calculated image, ndarray of shape (N)
     342        '''
     343        return G.dot(data)
     344
     345    def tropus (image, G):
     346        '''
     347        tropus: transform solution-space -> data-space:  [G]^tr * image
     348       
     349        default definition, caller can use this definition or provide an alternative
     350       
     351        :param float[N] image: solution, ndarray of shape (N)
     352        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
     353        :returns float[M]: calculated data, ndarray of shape (M)
     354        '''
     355        return G.transpose().dot(image)
     356
     357    def Dist(s2, beta):
     358        '''measure the distance of this possible solution'''
     359        w = 0
     360        n = beta.shape[0]
     361        for k in range(n):
     362            z = -sum(s2[k] * beta)
     363            w += beta[k] * z
     364        return w
     365   
     366    def ChiNow(ax, c1, c2, s1, s2):
     367        '''
     368        ChiNow
     369       
     370        :returns tuple: (ChiNow computation of ``w``, beta)
     371        '''
     372       
     373        bx = 1 - ax
     374        a =   bx * c2  -  ax * s2
     375        b = -(bx * c1  -  ax * s1)
     376   
     377        beta = ChoSol(a, b)
     378        w = 1.0
     379        for k in range(SEARCH_DIRECTIONS):
     380            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
     381        return w, beta
     382   
     383    def ChoSol(a, b):
     384        '''
     385        ChoSol: ? chop the solution vectors ?
     386       
     387        :returns: new vector beta
     388        '''
     389        n = b.shape[0]
     390        fl = numpy.ndarray((n, n))*0
     391        bl = numpy.ndarray((n))*0
     392       
     393        #print_arr("ChoSol: a", a)
     394        #print_vec("ChoSol: b", b)
     395   
     396        if (a[0][0] <= 0):
     397            msg = "ChoSol: a[0][0] = "
     398            msg += str(a[0][0])
     399            msg += '  Value must be positive'
     400            raise MaxEntException(msg)
     401   
     402        # first, compute fl from a
     403        # note fl is a lower triangular matrix
     404        fl[0][0] = math.sqrt (a[0][0])
     405        for i in (1, 2):
     406            fl[i][0] = a[i][0] / fl[0][0]
     407            for j in range(1, i+1):
     408                z = 0.0
     409                for k in range(j):
     410                    z += fl[i][k] * fl[j][k]
     411                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
     412                z = a[i][j] - z
     413                if j == i:
     414                    y = math.sqrt(z)
     415                else:
     416                    y = z / fl[j][j]
     417                fl[i][j] = y
     418        #print_arr("ChoSol: fl", fl)
     419   
     420        # next, compute bl from fl and b
     421        bl[0] = b[0] / fl[0][0]
     422        for i in (1, 2):
     423            z = 0.0
     424            for k in range(i):
     425                z += fl[i][k] * bl[k]
     426                #print "\t", i, k, z
     427            bl[i] = (b[i] - z) / fl[i][i]
     428        #print_vec("ChoSol: bl", bl)
     429   
     430        # last, compute beta from bl and fl
     431        beta = numpy.ndarray((n))
     432        beta[-1] = bl[-1] / fl[-1][-1]
     433        for i in (1, 0):
     434            z = 0.0
     435            for k in range(i+1, n):
     436                z += fl[k][i] * beta[k]
     437                #print "\t\t", i, k, 'z=', z
     438            beta[i] = (bl[i] - z) / fl[i][i]
     439        #print_vec("ChoSol: beta", beta)
     440   
     441        return beta
     442
     443    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
     444        '''
     445        move beta one step closer towards the solution
     446        '''
     447        a_lower, a_upper = 0., 1.          # bracket  "a"
     448        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
     449        #print "MaxEntMove: cmin = %g" % cmin
     450        if cmin*chisq > chizer:
     451            ctarg = (1.0 + cmin)/2
     452        else:
     453            ctarg = chizer/chisq
     454        f_lower = cmin - ctarg
     455        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
     456        f_upper = c_upper - ctarg
     457   
     458        fx = 2*MOVE_PASSES      # just to start off
     459        loop = 1
     460        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
     461            a_new = (a_lower + a_upper) * 0.5           # search by bisection
     462            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
     463            fx = c_new - ctarg
     464            # tighten the search range for the next pass
     465            if f_lower*fx > 0:
     466                a_lower, f_lower = a_new, fx
     467            if f_upper*fx > 0:
     468                a_upper, f_upper = a_new, fx
     469            loop += 1
     470   
     471        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
     472            msg = "MaxEntMove: Loop counter = "
     473            msg += str(MAX_MOVE_LOOPS)
     474            msg += '  No convergence in alpha chop'
     475            raise MaxEntException(msg)
     476   
     477        w = Dist (s2, beta);
     478        m = SEARCH_DIRECTIONS
     479        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
     480            for k in range(m):
     481                beta[k] *= math.sqrt (fSum/(blank*w))
     482        chtarg = ctarg * chisq
     483        return w, chtarg, loop, a_new, fx, beta
     484#MaxEnt_SB starts here   
     485    if image_to_data == None:
     486        image_to_data = opus
     487    if data_to_image == None:
     488        data_to_image = tropus
     489    n   = len(base)
     490    npt = len(datum)
     491
     492    # Note that the order of subscripts for
     493    # "xi" and "eta" has been reversed from
     494    # the convention used in the FORTRAN version
     495    # to enable parts of them to be passed as
     496    # as vectors to "image_to_data" and "data_to_image".
     497    xi      = 0*numpy.ndarray((SEARCH_DIRECTIONS, n))
     498    eta     = 0*numpy.ndarray((SEARCH_DIRECTIONS, npt))
     499    beta    = 0*numpy.ndarray((SEARCH_DIRECTIONS))
     500    # s1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
     501    # c1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
     502    s2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
     503    c2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
     504
     505    # TODO: replace blank (scalar) with base (vector)
     506    blank = sum(base) / len(base)   # use the average value of base
     507
     508    chizer, chtarg = npt*1.0, npt*1.0
     509    f = base * 1.0                              # starting distribution is base
     510
     511    fSum  = sum(f)                              # find the sum of the f-vector
     512    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
     513    chisq = sum(z*z)                            # Chi^2, SB eq. 4
     514
     515    for iter in range(IterMax):
     516        ox = -2 * z / sigma                        # gradient of Chi^2
     517
     518        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
     519        sgrad = -numpy.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
     520        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
     521        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
     522        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
     523
     524        a = 1.0
     525        b = 1.0 / cnorm
     526        if iter == 0:
     527            test = 0.0     # mismatch between entropy and ChiSquared gradients
     528        else:
     529            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
     530            a = 0.5 / (snorm * test)
     531            b *= 0.5 / test
     532        xi[0] = f * cgrad / cnorm
     533        xi[1] = f * (a * sgrad - b * cgrad)
     534
     535        eta[0] = image_to_data (xi[0], G);          # image --> data
     536        eta[1] = image_to_data (xi[1], G);          # image --> data
     537        ox = eta[1] / (sigma * sigma)
     538        xi[2] = data_to_image (ox, G);              # data --> image
     539        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
     540        xi[2] = f * xi[2] * a
     541        eta[2] = image_to_data (xi[2], G)           # image --> data
     542       
     543#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
     544#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
     545
     546        # prepare the search directions for the conjugate gradient technique
     547        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
     548        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
     549#         print_vec("MaxEnt: c1", c1)
     550#         print_vec("MaxEnt: s1", s1)
     551
     552        for k in range(SEARCH_DIRECTIONS):
     553            for l in range(k+1):
     554                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
     555                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
     556#         print_arr("MaxEnt: c2", c2)
     557#         print_arr("MaxEnt: s2", s2)
     558
     559        # reflect across the body diagonal
     560        for k, l in ((0,1), (0,2), (1,2)):
     561            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
     562            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
     563 
     564        beta[0] = -0.5 * c1[0] / c2[0][0]
     565        beta[1] = 0.0
     566        beta[2] = 0.0
     567        if (iter > 0):
     568            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
     569
     570        f_old = f.copy()    # preserve the last image
     571        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
     572       
     573        # As mentioned at the top of p.119,
     574        # need to protect against stray negative values.
     575        # In this case, set them to RESET_STRAYS * base[i]
     576        #f = f.clip(RESET_STRAYS * blank, f.max())
     577        for i in range(n):
     578            if f[i] <= 0.0:
     579                f[i] = RESET_STRAYS * base[i]
     580        df = f - f_old
     581        fSum = sum(f)
     582        fChange = sum(df)
     583
     584        # calculate the normalized entropy
     585        S = -sum((f/fSum) * numpy.log(f/fSum))      # normalized entropy, S&B eq. 1
     586        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
     587        chisq = sum(z*z)                            # report this ChiSq
     588
     589        if report:
     590            print "%3d/%3d" % ((iter+1), IterMax)
     591            print " %5.2lf%% %8lg" % (100*test, S)
     592            if iter > 0:
     593                value = 100*( math.sqrt(chisq/chtarg)-1)
     594            else:
     595                value = 0
     596            print " %12.5lg %10.4lf" % ( math.sqrt (chtarg/npt), value )
     597            print "%12.6lg %8.2lf\n" % (fSum, 100*fChange/fSum)
     598
     599        # See if we have finished our task.
     600        # do the hardest test first
     601        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
     602            return f,image_to_data (f, G)     # solution FOUND returns here
     603   
     604    return f,image_to_data (f, G)       # no solution after IterMax iterations
     605
     606   
     607################################################################################
     608#### MaxEnt testing stuff
     609################################################################################
     610
     611def print_vec(text, a):
     612    '''print the contents of a vector to the console'''
     613    n = a.shape[0]
     614    print "%s[ = (" % text,
     615    for i in range(n):
     616        s = " %g, " % a[i]
     617        print s,
     618    print ")"
     619
     620def print_arr(text, a):
     621    '''print the contents of an array to the console'''
     622    n, m = a.shape
     623    print "%s[][] = (" % text
     624    for i in range(n):
     625        print " (",
     626        for j in range(m):
     627            print " %g, " % a[i][j],
     628        print "),"
     629    print ")"
     630
     631def test_MaxEnt_SB(report=True):
     632    def readTextData(filename):
     633        '''return q, I, dI from a 3-column text file'''
     634        if not os.path.exists(filename):
     635            raise Exception("file not found: " + filename)
     636        buf = [line.split() for line in open(filename, 'r').readlines()]
     637        M = len(buf)
     638        buf = zip(*buf)         # transpose rows and columns
     639        q  = numpy.array(buf[0], dtype=numpy.float64)
     640        I  = numpy.array(buf[1], dtype=numpy.float64)
     641        dI = numpy.array(buf[2], dtype=numpy.float64)
     642        return q, I, dI
     643    print "MaxEnt_SB: "
     644    test_data_file = os.path.join( 'testinp', 'test.sas')
     645    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
     646    bkg   = 0.1     #   I = I - bkg
     647    dMin, dMax, nRadii = 25, 9000, 40
     648    defaultDistLevel = 1.0e-6
     649    IterMax = 40
     650    errFac = 1.05
     651   
     652    r    = numpy.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
     653    dr   = r * (r[1]/r[0] - 1)          # step size
     654    f_dr = numpy.ndarray((nRadii)) * 0  # volume fraction histogram
     655    b    = numpy.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
     656   
     657    qVec, I, dI = readTextData(test_data_file)
     658    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
     659   
     660    f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
     661    if f_dr is None:
     662        print "no solution"
     663        return
     664   
     665    print "solution reached"
     666    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
     667        print '%10.4f %10.4f %12.4g'%(a,b,c)
     668
     669def tests():
     670    test_MaxEnt_SB(report=True)
     671
     672if __name__ == '__main__':
     673    tests()
    302674   
    303675###############################################################################
     
    306678
    307679def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
     680    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
     681        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
     682        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
     683        'Unified disk':[UniDiskFF,UniDiskVol]}
     684    Shape = data['Size']['Shape'][0]
     685    Parms = data['Size']['Shape'][1:]
    308686    if data['Size']['logBins']:
    309687        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
    310             data['Size']['Nbins']+1,True)
     688            data['Size']['Nbins']+1,True)/2.        #make radii
    311689    else:
    312690        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
    313             data['Size']['Nbins']+1,True)
     691            data['Size']['Nbins']+1,True)/2.        #make radii
    314692    Dbins = np.diff(Bins)
    315     BinMag = Dbins*np.ones_like(Dbins)
     693    Bins = Bins[:-1]+Dbins/2.
     694    BinsBack = np.ones_like(Bins)*1.e-6
     695    Back = data['Back']
     696    Q,Io,wt,Ic,Ib = Profile[:5]
     697    Qmin = Limits[1][0]
     698    Qmax = Limits[1][1]
     699    Contrast = Sample['Contrast'][1]
     700    Ibeg = np.searchsorted(Q,Qmin)
     701    Ifin = np.searchsorted(Q,Qmax)
     702    if Back[1]:
     703        Ib = Back[0]
     704        Ic[Ibeg:Ifin] = Back[0]
     705    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
     706    BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Io[Ibeg:Ifin]-Back[0],1./np.sqrt(wt[Ibeg:Ifin]),BinsBack,
     707        data['Size']['MaxEnt']['Niter'],Gmat)
     708    print BinMag.shape
     709    data['Size']['Distribution'] = [Bins,Dbins,BinMag]
    316710    print np.sum(BinMag)
    317711       
    318     print data['Size']
    319 #    print Limits
    320 #    print Substances
    321 #    print Sample
    322 #    print Profile
    323    
     712   
Note: See TracChangeset for help on using the changeset viewer.