SasView / sasview

Code for the SasView application.
BSD 3-Clause "New" or "Revised" License
49 stars 41 forks source link

Should the low Q power law extrapolation for the invariant be removed #1495

Open butlerpd opened 4 years ago

butlerpd commented 4 years ago

There invariant perspective currently allows for calculating the invariant from the data with options to also estimate the low q and high contributions to the invariant from the curve at lower q than qmin and higher q than q max. This is done by fitting the last n points to a power law for high Q (some form of Porod) and a Guinier at low Q. Currently there is also the option to use a power law extrapolation to low Q. It is not at all clear why this low Q power law is an option. By definition a power law must diverge at q=0 meaning Q_star is infinite. Practically, it means that the experimental curve does not go to low enough Q to capture the full structure in which case the invariant is calculation is invalid. Only if one has a true Guinier region have you reached the full structural data. I realize for Geosciences this is often used even in the absence of a Guinier but in that case in principle they are measuring to very low Q not extrapolating and presumably assuming the the remaining low Q contribution is negligible?

If we remove the option, note that the unit tests based on data in sasinvariant/test/utest_use_cases.py have been removed while those in sasinvariant/test/utest_data_handling.py mostly remain but have notes in the doc strings about removing should we need to remove the feature.

smk78 commented 4 years ago

I would say that providing an extrapolation which does not comply with the mathematics of how the invariant is defined is just plain wrong and it should be removed!

If we want to offer an alternative to the Guinier, then I would suggest Vonk's function, as implemented in the original Corfunc:

C           Vonk model
            sigy=0.
            sigxsqu=0.
            sigxsquy=0.
            sigx4=0.

            DO j=start,chanend
C             Get sums
              temp=xdata(j)**2
              sigy=sigy+ydata(j)
              sigxsqu=sigxsqu+temp
              sigxsquy=sigxsquy+temp*ydata(j)
              sigx4=sigx4+temp*temp
            END DO

C           Calculate parameters
            npts=(chanend-start+1)
            h2=(sigxsqu*sigy-npts*sigxsquy)/(npts*sigx4-(sigxsqu**2))
            h1=(sigy+h2*sigxsqu)/npts

C           Get chi squ.
            chisqu=0.
            DO j=start,chanend
              chisqu=chisqu+(ydata(j)-h1+h2*xdata(j)*xdata(j))**2
            END DO
            chisqu=chisqu/npts

The references for this are actually at the bottom of the SasView Corfunc help file (although it is not actually implemented in SasView Corfunc; the references are there because they also discuss correlation functions).

In fact, the Vonk extrapolation could be implemented for both analyses.

butlerpd commented 4 years ago

I've not worked with the Vonk function model before. Looking at the code I can see invariant variations that looks like some sort of residual/uncertainty approach? .... but where is the extrapolation part? it looks like just the main integral?

Will have to look at that at some point in more detail I guess - thanks for the refs

smk78 commented 4 years ago

The functional form is in the last DO loop. So basically the difference is:

Vonk: I(Q) = h1-h2*Q^2

Guinier: I(Q) = A exp(BQ^2) ; where B is negative

lucas-wilkins commented 2 years ago

I'm not sure I follow this @butlerpd, you say

"By definition a power law must diverge at q=0 meaning Q_star is infinite."

If the power were necessarily less than -1 I would agree, but if its being fitted then its whatever fits. If it is negative 4, which is what the default fixed value is, then there's a problem.

So is the problem perhaps just the default value (which appears to be set to -4 because that's what the Porod law panel has)?

butlerpd commented 2 years ago

Doesn't q^-1 also diverge at q=0? I believe only q^0 (and higher exponents)? That said we need to talk to the porous materials people (geologists) to really understand the usage before suggesting anything specific so I would not do anything before understanding the use case properly.

lucas-wilkins commented 2 years ago

Oops, it seems I can't do basic integration, I should have said less than or equal to -3 (because the invariant is the integral of I q^2 dq). A fitted power of e.g. -2.9 should converge.

I do think there's something fishy about it - at first glance it doesn't seem to have any physical motivation, and it's quite a specific/constrained function to fit otherwise - but then again I am new to this, and there's quite likely something I'm missing.

On Mon, May 23, 2022 at 9:22 PM Paul Butler @.***> wrote:

Doesn't q^-1 also diverge at q=0? I believe only q^0 (and higher exponents)? That said we need to talk to the porous materials people (geologists) to really understand the usage before suggesting anything specific so I would not do anything before understanding the use case properly.

— Reply to this email directly, view it on GitHub https://github.com/SasView/sasview/issues/1495#issuecomment-1135102463, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACKU4SUHEVSCIGGNMT6PAXDVLPSGTANCNFSM4LS6CSQA . You are receiving this because you commented.Message ID: @.***>

--

Dr Lucas Wilkins +44 (0) 7505 915 726

Personal Website: http://www.lucaswilkins.com/ Alternate e-mail: @.***