Discrepancy in code for polymer_excl_volume (Trac #1073)

[smk78]

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.

Migrated from http://trac.sasview.org/ticket/1073

{
    "status": "closed",
    "changetime": "2018-03-04T17:29:02",
    "_ts": "2018-03-04 17:29:02.398708+00:00",
    "description": "Reported by pkienzle:\n\nJust 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?\n\nHere the code that Jae He added:\n{{{\n    $ git show -p ec1c5540\n    2010-10-04 ec1c5540 added polymerexclVolume model and its test Jae Cho\n    diff --git a/sansmodels/src/sans/models/PolymerExclVolume.py b/sansmodels/src/sans/models/PolymerExclVolume.py\n    ...\n    +        nu = 1.0 / mm\n    +    \n    +        Xx = x * x * rg * rg *(2.0 * nu + 1.0) * (2.0 * nu + 2.0) / 6.0\n    +        onu = 1.0 / nu\n    +        o2nu = 1.0 /(2.0 * nu)\n    +        Ps =(1.0 / (nu * pow(Xx,o2nu))) * (gamma(o2nu)*gammainc(o2nu,Xx) - \\\n    +                        1.0 / pow(Xx,o2nu) * gamma(onu)*gammainc(onu,Xx))\n    +\n    +        if x == 0:\n    +            Ps = 1.0\n    +\n    +        return (sc * Ps + bg);\n}}}\n\nHere's what is in Igor 7.12, !Analysis/Models/NewModels_2010/PolymerExcludVol_v40.ipf:\n{{{\n        nu=1/m\n        Xx=qval^2*Rg^2*(2*nu+1)*(2*nu+2)/6\n        onu=1.0/(nu)\n        o2nu=1.0/(2*nu)\n        Ps=(1/(nu*Xx^o2nu))*(gammaInc(o2nu,Xx,0)-(1/Xx^o2nu)*gammaInc(onu,Xx,0))\n        Debye=(2/Xx^2)*(exp(-Xx)-1+Xx)\n        \n        if(qval == 0)\n                Ps = 1\n        endif\n        \n        inten = I0*Ps + bgd \n\nIgor gammaInc(a, x, 0) =  gamma(a) - gammaInc(a, x) =  int_0^x e^-t t^(a-1) dt = scipy.special.gammainc(a, x).\n}}}\nNote: 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.\n",
    "reporter": "smk78",
    "cc": "",
    "resolution": "invalid",
    "workpackage": "SasView Bug Fixing",
    "time": "2018-02-14T11:32:40",
    "component": "sasmodels",
    "summary": "Discrepancy in code for polymer_excl_volume",
    "priority": "major",
    "keywords": "",
    "milestone": "sasmodels 1.0",
    "owner": "",
    "type": "defect"
}

[smk78]

Trac update at 2018/02/14 11:36:04:

• smk78 changed attachment from "" to "Hammouda1993PolymerModels.pdf"
• smk78 commented:

Hammouda reference for polymer_excl_volume model

[smk78]

Trac update at 2018/02/14 11:37:25: smk78 commented:

!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.

[butlerpd]

Trac update at 2018/02/14 14:03:45: butler changed _comment1 from:

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.5porod_exp) result= (porod_expupow (gamma(0.5porod_exp)gammainc(0.5porod_exp, usub) - upowgamma(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.5porod_exp) upow2 = power(usub, -porod_exp) result= (porod_expupow1 (gamma(0.5porod_exp)gammainc(0.5porod_exp, usub)upow2gamma(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 (1/2*nu for the first, and 1/nu for the second).

A second email after he sent the above:

• I see where the disagreement is from. Boualem’s paper from 1993 has the form that’s currently in SasView, but more recent publications like https://www.ncnr.nist.gov/staff/hammouda/publications/2017_hammouda_kim_j_mol_liq.pdf and mine https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf have the form that I’m talking about. So I guess it’s a matter of which is correct!

to:

1520184429808313

[sasview-bot]

Trac update at 2018/02/14 14:03:45:

• krzywon changed _comment0 from:

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.5porod_exp) result= (porod_expupow (gamma(0.5porod_exp)gammainc(0.5porod_exp, usub) - upowgamma(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.5porod_exp) upow2 = power(usub, -porod_exp) result= (porod_expupow1 (gamma(0.5porod_exp)gammainc(0.5porod_exp, usub)upow2gamma(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 (1/2*nu for the first, and 1/nu for the second).

A second email after hesent the above:

• I see where the disagreement is from. Boualem’s paper from 1993 has the form that’s currently in SasView, but more recent publications like https://www.ncnr.nist.gov/staff/hammouda/publications/2017_hammouda_kim_j_mol_liq.pdf and mine https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf have the form that I’m talking about. So I guess it’s a matter of which is correct!

to:

1518617078353696

• krzywon commented:

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 (1/2*nu for the first, and 1/nu for the second).

A second email after he sent the above:

• I see where the disagreement is from. Boualem’s paper from 1993 has the form that’s currently in !SasView, but more recent publications like https://www.ncnr.nist.gov/staff/hammouda/publications/2017_hammouda_kim_j_mol_liq.pdf and mine https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf have the form that I’m talking about. So I guess it’s a matter of which is correct!

[smk78]

Trac update at 2018/02/14 14:12:41:

• smk78 changed _comment0 from "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!!! Ergo, one of the Hammouda papers must be wrong." to "1518643484897285"
• smk78 commented:

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!!!

[butlerpd]

Trac update at 2018/02/14 19:30:46: butler changed _comment0 from "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." to "1520184500379176"

[sasview-bot]

Trac update at 2018/02/14 19:30:46: krzywon commented:

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.

[smk78]

Trac update at 2018/02/14 21:25:29: smk78 commented:

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

[sasview-bot]

Trac update at 2018/02/14 21:36:50:

• krzywon commented:

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

• krzywon changed resolution from "" to "invalid"
• krzywon changed status from "new" to "closed"

[smk78]

Trac update at 2018/02/15 12:29:17: smk78 commented:

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

[sasview-bot]

Trac update at 2018/02/15 12:59:50: Paul Kienzle pkienzle@nist.gov commented:

In changeset 547c6f0c4f0182017463ec1d667299383fe53598:

#!CommitTicketReference repository="sasmodels" revision="547c6f0c4f0182017463ec1d667299383fe53598"
polymer_excl_volume: update eq. to match paper. Refs #169

[sasview-bot]

Trac update at 2018/02/15 16:52:49: Paul Kienzle pkienzle@nist.gov commented:

In changeset aa900158f288ffb14cd35ce414d65dcc836750dd:

#!CommitTicketReference repository="sasmodels" revision="aa900158f288ffb14cd35ce414d65dcc836750dd"
polymer_excl_volume: note in the code that scipy gammainc is regularized, hence the gamma scale factor. Refs #169

[butlerpd]

Trac update at 2018/03/04 17:29:02:

• butler changed component from "SasView" to "sasmodels"
• butler changed milestone from "SasView 4.3.0" to "sasmodels 1.0"