bift.py result on BSA test data shows "saw" on top of data

[maaeli]

When running bift.py on the BSA test data, there is a clear saw pattern with wavelength 2Δr on top of the density. Increasing MC sampling reduces the amplitude, but not sufficiently. Taking the averages of two subsequent points removes the pattern:

Also see issue #18 .

[maaeli]

This effect should be quantifiable by looking at the last value of the abs(FFT) of the density

[kif]

Could you tell me which branch you are using ?

[maaeli]

I am on kif/freesas/master, the last commit is the following:

commit 8e4f4d7 (HEAD -> master, origin/master)
Author: Jerome Kieffer [email protected]
Date: Thu Apr 30 18:43:40 2020 +0200

increase version number as bift algorithm changed

[kif]

I have been mostly hacking on the Guinier region code recently , with reverse engineering of Atsas' code, but this lead also some improvement in bift. Sorry I did not merge the PR as I am not happy with the quality indicator yet.
All the recent development has been done in the fit_autorg branch:
https://github.com/kif/freesas/tree/fit_autorg

[maaeli]

Unfortunately, I get the same phenomenon with v0.8.3

The investigation continues...

[maaeli]

Apparently smooth_density introduces a discontinuity in a smooth function:

import numpy
from freesas._bift import distribution_parabola, distribution_sphere, smooth_density
ps = numpy.asarray(distribution_sphere(60.1, 10.6, 100))
smooth = numpy.zeros(101, numpy.float64)
smooth_density(ps,smooth)
plt.figure()
plt.plot(ps,label="sphere")
plt.plot(smooth,label="smooth")
plt.legend()

This seems quite suspicious.

@kif: Does the original paper specify how to smooth? I no longer have access to it...

[kif]

Hi Martha, nice catch.
The publication states:

[kif]

I confirm you that the bug comes from this part as I don't have the same results as you for the smoothed curve: I see no such discontinuity:

After 10 iteration, it looks like:

Could you tell me which compiler and version of cython you are using ? because the code is very simple (even if I noticed some oddities in the order of the operations. This what Hansen initially used in his Fortran code:

c******************************
c     Smoothness prior
c******************************
      DO 64 J=2,NTOT-1
      M(J)=SUMFACT*M(J)
      m(j)=(f(j-1)+f(j+1))/2
   64 CONTINUE
      m(1)=f(2)/2
      m(ntot)=m(ntot-1)/2
      if(answer5.eq.'C') m(1)=f(2)
c*******************************
c     smoothness end
c*******************************

Beside the difference of indexing, it is translated to Cython like:

cpdef inline void smooth_density(double[::1] raw,
                                 double[::1] smooth) nogil:
        """This function applies the smoothing of the density plot
        
        :param raw: raw density, called *f* in eq.19 
        :param smooth: smoothed density, called *m* in eq.19
        
        The smoothing is performed in place
        """
        cdef:
            int k, npt
        #assert raw.shape[0] == smooth.shape[0]
        npt = raw.shape[0] - 1
        # This enforces the boundary values to be null
        smooth[0] = raw[0] = 0.0        
        smooth[npt] = raw[npt] = 0.0        
        smooth[1] = raw[2] * 0.5
        for k in range(2, npt-1):
            smooth[k] = 0.5 * (raw[k-1] + raw[k+1])
        smooth[npt-1] = smooth[npt-2] * 0.5 # is it p or f on the RHS? does this enforce a smoother tail ?

[maaeli]

Hi,

my python version is 0.29.20 and python 3.8.
And GCC --version says this: Apple LLVM version 10.0.1 (clang-1001.0.46.4)

However, I after changing the code and reverting back to the original, I also cannot reproduce the above problem and more. Not that the "new" result is correct:

ps = numpy.ones(101, numpy.float64)
smooth = numpy.zeros(101, numpy.float64)
smooth_density(ps,smooth)
print(smooth[98],smooth[99],smooth[100])
plt.figure()
plt.plot(ps,label="start")
plt.plot(smooth,label="smooth")
print(smooth[1])
plt.legend()

gives the following output:
0.0 0.0 0.0

Now, as I understand it (from the code), smooth[99] should 0.5 and smooth[98] should be 1.0???

Cheers,
Martha

PS:
The "magic" was that after npt = raw.shape[0] - 1 I inserted

with gil:
    print(npt)

[maaeli]

By the way, I would suggest to split the current smooth function into 2 functions:

• ensureEdgesZero(f_r) which sets the first and last element of f_r to 0
• smooth(raw, smooth) which does not change raw but ensures that smooth[0] = raw[0] and smooth[-1] = raw[-1]

[kif]

Hi Martha,
I tried with my mac, python3.8, compiler based on LLVM8 (sorry old mac) and cython version 0.29.20 and 3.0-beta. All give the proper answer.
I am now upgrading MacOS, and hopefully not break my mac.

By any chance, can you run run python setup.py build --force-cython I suspect there are former files laying arround which could interfere. Sorry the build environment in freesas is not as hardened as the one of silx.

[kif]

I upgraded my macbook to the latest version to be able to upgrade the compiler:

(py38) mac13:freesas kieffer$ python --version
Python 3.8.0
(py38) mac13:freesas kieffer$ cython --version
Cython version 0.29.20
(py38) mac13:freesas kieffer$ gcc --version
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/4.2.1
Apple clang version 11.0.3 (clang-1103.0.32.62)
Target: x86_64-apple-darwin19.5.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

And after rebuild, the code behaves as expected:

In [1]: import numpy 
   ...: from freesas._bift import distribution_parabola, distribution_sphere, smooth_density 
   ...: ps = numpy.ones(101, numpy.float64) 
   ...: smooth = numpy.zeros(101, numpy.float64) 
   ...: smooth_density(ps,smooth) 
   ...: print(smooth)                                                                                                                                                                         
[0.  0.5 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
 1.  1.  1.  1.  1.  1.  1.  1.  1.  0.5 0. ]

[maaeli]

Hi,

the answer was: Use clang11 instead of clang10 ;)

gcc --version
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/4.2.1
Apple clang version 11.0.0 (clang-1100.0.33.17)
Target: x86_64-apple-darwin18.7.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

Things look reasonable now!

Thanks for your help,
Martha