Opened 7 years ago

Closed 7 years ago

Last modified 7 years ago

#1073 closed defect (invalid)

Discrepancy in code for polymer_excl_volume

Reported by: smk78 Owned by:
Priority: major Milestone: sasmodels 1.0
Component: sasmodels Keywords:
Cc: Work Package: SasView Bug Fixing

Description

Reported by pkienzle:

Just went through the code for polymer_excl_volume. The result should match older versions of SasView, but this does not match the documentation, nor does it match Igor. For some reason, when Jae He added the model he scaled the gammainc term by gamma. Maybe Paul or Mat has some memory of why?

Here the code that Jae He added:

    $ git show -p ec1c5540
    2010-10-04 ec1c5540 added polymerexclVolume model and its test Jae Cho
    diff --git a/sansmodels/src/sans/models/PolymerExclVolume.py b/sansmodels/src/sans/models/PolymerExclVolume.py
    ...
    +        nu = 1.0 / mm
    +    
    +        Xx = x * x * rg * rg *(2.0 * nu + 1.0) * (2.0 * nu + 2.0) / 6.0
    +        onu = 1.0 / nu
    +        o2nu = 1.0 /(2.0 * nu)
    +        Ps =(1.0 / (nu * pow(Xx,o2nu))) * (gamma(o2nu)*gammainc(o2nu,Xx) - \
    +                        1.0 / pow(Xx,o2nu) * gamma(onu)*gammainc(onu,Xx))
    +
    +        if x == 0:
    +            Ps = 1.0
    +
    +        return (sc * Ps + bg);

Here's what is in Igor 7.12, Analysis/Models/NewModels_2010/PolymerExcludVol_v40.ipf:

        nu=1/m
        Xx=qval^2*Rg^2*(2*nu+1)*(2*nu+2)/6
        onu=1.0/(nu)
        o2nu=1.0/(2*nu)
        Ps=(1/(nu*Xx^o2nu))*(gammaInc(o2nu,Xx,0)-(1/Xx^o2nu)*gammaInc(onu,Xx,0))
        Debye=(2/Xx^2)*(exp(-Xx)-1+Xx)
        
        if(qval == 0)
                Ps = 1
        endif
        
        inten = I0*Ps + bgd 

Igor gammaInc(a, x, 0) =  gamma(a) - gammaInc(a, x) =  int_0^x e^-t t^(a-1) dt = scipy.special.gammainc(a, x).

Note: The current version looks different from Jae He's because I updated it during one of my cleanup passes. I'm pretty sure that I didn't change the result, but probably should revert it so that the equations in the code better match the equations in the documentation.

Change History (11)

Changed 7 years ago by smk78

Hammouda reference for polymer_excl_volume model

comment:1 Changed 7 years ago by smk78

SasView documentation for polymer_excl_volume model references:

H Benoit, Comptes Rendus, 245 (1957) 2244-2247

B Hammouda, SANS from Homogeneous Polymer Mixtures - A Unified Overview, Advances in Polym. Sci. 106(1993) 87-133

Second of these is attached to this ticket.

comment:2 Changed 7 years ago by krzywon

From Mike Hore:

  • The reason that gammainc is scaled by gamma in the SasView function is that gammainc in scipy is normalized by gamma, so that is 100% correct to do.
  • Looking at the function in the version of SasView (4.1.2) that I have installed currently, the issue is in the following block of code:
    usub = (q*rg)**2 * (2.0/porod_exp + 1.0) * (2.0/porod_exp + 2.0)/6.0
    with errstate(divide='ignore', invalid='ignore'):
        upow = power(usub, -0.5*porod_exp)
        result= (porod_exp*upow *
                (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp, usub) -
                upow*gamma(porod_exp)*gammainc(porod_exp, usub)))
    
  • I think it should be something like the following instead:
    usub = (q*rg)**2 * (2.0/porod_exp + 1.0) * (2.0/porod_exp + 2.0)/6.0
    with errstate(divide='ignore', invalid='ignore'):
        upow1 = power(usub, -0.5*porod_exp)
        upow2 = power(usub, -porod_exp)
        result= (porod_exp*upow1 *
                (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp,
                usub)upow2*gamma(porod_exp)*gammainc(porod_exp, usub)))
    
  • Because U is raised to a different power for the first and second terms in the correct version of the function (½*nu for the first, and 1/nu for the second).

A second email after he sent the above:

Last edited 7 years ago by butler (previous) (diff)

comment:3 Changed 7 years ago by smk78

Interestingly, https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf references the polymer with excluded volume form factor with exactly the same two references that are in the SasView documentation!!!

Last edited 7 years ago by smk78 (previous) (diff)

comment:4 Changed 7 years ago by krzywon

I just spoke with Boualem and we went through the math together. In the 1993 paper, [1/vX^(1/2v)] has been factored out. When multiplied through, the resulting equation is identical to the equation in the 2013 and 2017 papers. The calculation currently within SasView matches that in the 1993 paper so I think the model might be correct.

Last edited 7 years ago by butler (previous) (diff)

comment:5 Changed 7 years ago by smk78

So can we close this ticket now, or is there still an issue?

comment:6 Changed 7 years ago by krzywon

  • Resolution set to invalid
  • Status changed from new to closed

I contacted Mike Hore and he agrees, the mathematics is correct. Yes, we can close this ticket.

comment:7 Changed 7 years ago by smk78

I have pushed some updates to the documentation to Master to reflect what has been learnt from this exchange.

comment:8 Changed 7 years ago by Paul Kienzle <pkienzle@…>

In 547c6f0c4f0182017463ec1d667299383fe53598/sasmodels:

polymer_excl_volume: update eq. to match paper. Refs #1073

comment:9 Changed 7 years ago by Paul Kienzle <pkienzle@…>

In aa900158f288ffb14cd35ce414d65dcc836750dd/sasmodels:

polymer_excl_volume: note in the code that scipy gammainc is regularized, hence the gamma scale factor. Refs #1073

comment:10 Changed 7 years ago by butler

  • Component changed from SasView to sasmodels
  • Milestone changed from SasView 4.3.0 to sasmodels 1.0
Note: See TracTickets for help on using tickets.