Adding a pure python truncated octahedron model with Fibonacci orientation averaging

[sara-mokhtari]

Description

Adding a pure python model for truncated octahedrons with orientation averaging using Fibonacci quadrature.

Fixes #686

How Has This Been Tested?

Unit tests worked.
This model was compared with the model using C code (octahedron_truncated.c) and shows good agreement.
In general, the octahedron model was compared to an another numerical calculation : for truncated octahedrons with different sizes and truncatures, agreement was found with the Debye equation (Debye calculator: https://github.com/FrederikLizakJohansen/DebyeCalculator). The Debye equation is based on atomic positions while SasView model is based on analytical expressions.

More tests should be made on small q. Indeed, like the previous model (octahedron_truncated.c), we encouter issues when it comes to small q (q<10^-6 Angs). More precise mathematical should be used (hospital rule, etc).
Note : the fibonacci quadrature code (sasmodels/quadratures/fibonacci.py) was added to a new repository called "quadratures" because it could be also useful for other models.

Review Checklist:

[if using the editor, use [x] in place of [ ] to check a box]

Documentation (check at least one)

• There is nothing that needs documenting
• Documentation changes are in this PR
• There is an issue open for the documentation (link?)

Installers

• There is a chance this will affect the installers, if so
  • Windows installer (GH artifact) has been tested (installed and worked)
  • MacOSX installer (GH artifact) has been tested (installed and worked)
  • Wheels installer (GH artifact) has been tested (installed and worked)

Licensing (untick if necessary)

• The introduced changes comply with SasView license (BSD 3-Clause)

[pkienzle]

Having different versions with different integration schemes is potentially confusing. If one is clearly better than the other we should be using that.

Controlling the accuracy of the model with a model parameter is unusual. This should never be a fitting parameter. If we want to support this then we need a new parameter type ("constant" or "fixed") to make sure that it doesn't show up as fittable in the GUI.

If fibonacci integration works better than 2D gaussian integration on the sphere then we should be applying it to all our shapes without rotational symmetry (see #658).

Support for Fq and beta is not yet available for pure python models (see #570).

Support for Iqac and Iqabc (orientation and jitter applied before calling kernel) is not yet implemented for pure python models, nor is support for magnetism.

My preference would be to leave this in the model marketplace until sasmodels has complete support for python models, or move the implementation of fibonacci integration to the C model.

I'm concerned about the singularities in the function. The usual trick of switching to a Taylor expansion around the singularity won't work here because the calculation is spread across the polyhedron vertices. Maybe a scheme which collects the bad edges and triangles for later processing would work, substituting zero instead of ≈1e18 in the summation, and correcting the integral for singularities when done. This applies both to the C and the python models regardless of integration scheme.

[pkienzle]

Setting denom to 1e-18 will lead to terms on the order of 1e18. Even when these eventually cancel they will effectively erase the rest of the integral during summation, leading to random values.

You can see if this is a problem using fibonacci lattice point k, setting c2a_ratio such that qny = qnz. Solving, c2a_ratio = b2a_ratio*qb/qc = b2a_ratio*q_unit[k, 1]/q_unit[k, 2].

[marimperorclerc]

Maybe we could remove this eps_den method ?
And make sure that for all Fibonacci points qny=qnz never occurs for any rational value of c2a_ratio ? Taking advantage that gold number is irrational ?
Not sure about this, but it could be tested.

[pkienzle]

Memory required is 3 * 8 bytes * len(q) * npoints_fibonacci, which is 120 MB for 5000 fibonacci points and 1000 q points. Probably okay.

[pkienzle]

You can make an Fq version using f1, f2 = sum(amplitude), sum(intensity)

[marimperorclerc]

Yes, but usually, the orientation average over the amplitudes won't be used.
But we can calculate it just in case.

[pkienzle]

This partially address the eps=1e-18 problem above by turning it into eps=1e-12. I'm expecting that you will have 5 or 6 bits of precision in your answer if you happen to hit a singularity.

[marimperorclerc]

Thank you for all the suggestions. :)
This is the first model with Fibonacci integration we would like to integrate in SASView.
We will also create another branch very soon for the nanoprism model we have already coded.
It is also using Fibonacci integration (pure python model) and nanoprisms will be new shapes, which are useful for inorganic nanoparticles.

Having different versions with different integration schemes is potentially confusing. If one is clearly better than the other we should be using that.

Yes, sure.
Here this second model for the same shape (truncated octahedron) is built in order to test the performances of the new Fibonacci integration scheme for a pure python model.
After the comparision with the c code model will be completely discussed together, we might decide what is the best model to keep in the SASView model list !

Controlling the accuracy of the model with a model parameter is unusual. This should never be a fitting parameter. If we want to support this then we need a new parameter type ("constant" or "fixed") to make sure that it doesn't show up as fittable in the GUI.

Yes, it is a good suggestion to have this possibility to add a new parameter type to avoid 'fiting' with it.
In this specific case, it it convenient for us to keep the number of Fibonacci points in the parameters list in order to adjust easily its value depending of q-range and type of scattering curve.
We are still in the process to determine a good default value for it, something between 300 or 500 points. And for that, we will probably need to test different form factors as well (see the nanoprisms model mentionned before).

If fibonacci integration works better than 2D gaussian integration on the sphere then we should be applying it to all our shapes without rotational symmetry (see #658).

Yes, this is the idea, and it is why we are comparing the two models for truncated octahedron.

Support for Fq and beta is not yet available for pure python models (see #570).
Support for Iqac and Iqabc (orientation and jitter applied before calling kernel) is not yet implemented for pure python models, nor is support for magnetism.

Would be great to have Fqabc in pure python model as well. What I have in mind here is a core/shell model for two polyhedron shapes. It is something we do have a lot using seed-mediated growth.
One shape is internal (seed) inside the whole particle. And for that, the combination of the two Fqabc is necessary (with a translation and relative rotation between the core particle and the whole particle) ... but this we may discuss later on ! ;)

My preference would be to leave this in the model marketplace until sasmodels has complete support for python models, or move the implementation of fibonacci integration to the C model.

These two options are fine with us.
Again, at this stage, we need to have the two models in the same branch to compare them smoothly.
Implementation of Fibonacci in C model would be difficult for us to implement, as we are much more confortable with python.

A last idea: maybe usage of numba package could help decreasing the calculation time for pure python even more ?

I'm concerned about the singularities in the function. The usual trick of switching to a Taylor expansion around the singularity won't work here because the calculation is spread across the polyhedron vertices. Maybe a scheme which collects the bad edges and triangles for later processing would work, substituting zero instead of ≈1e18 in the summation, and correcting the integral for singularities when done. This applies both to the C and the python models regardless of integration scheme.

Yes, we are aware of this issue. It is present for any polyhedron shape, and there are some solutions in the litterature to treat them (but not implemented here, it was too complicated for us to do that until now).
But I can provide the references and we can discuss about them.
In the framework of this model, if we can make sure that singularities are occuring only in the plateau region, at very small q values, I believe that it will be good enough. It is what we get from our own experience when using this model and the nanoprism model.

Best wishes,
Marianne

[pkienzle]

A last idea: maybe usage of numba package could help decreasing the calculation time for pure python even more ?

Yes, we could use numba and/or torch for models, but they use the pure python interface. Oriented 2D, magnetism and model reparameterization are only available for C models.

Python models do support polydispersity by looping a parameter like radius over (radius ± Δradius), but we can't do this for orientational dispersion. At the pole (θ=0°) the scattering from φ+Δφ is the same as φ, so ranging over φ±Δφ would have no effect. Instead we first apply dispersion to the shape, rotating the (a, b, c) axes by (Δθ, Δφ, Δψ) [pitch-yaw-roll] then we rotate by (θ, φ, ψ) to set the orientation. We then call Iqabc with the equivalent q vector for the particle in canonical orientation, with the c-axis aligned along the beam. The code for this is in kernel_iq.c for C models, but not yet in kernelpy.py for python models.

[marimperorclerc]

A last idea: maybe usage of numba package could help decreasing the calculation time for pure python even more ?

Yes, we could use numba and/or torch for models, but they use the pure python interface. Oriented 2D, magnetism and model reparameterization are only available for C models.

Python models do support polydispersity by looping a parameter like radius over (radius ± Δradius), but we can't do this for orientational dispersion. At the pole (θ=0°) the scattering from φ+Δφ is the same as φ, so ranging over φ±Δφ would have no effect. Instead we first apply dispersion to the shape, rotating the (a, b, c) axes by (Δθ, Δφ, Δψ) [pitch-yaw-roll] then we rotate by (θ, φ, ψ) to set the orientation. We then call Iqabc with the equivalent q vector for the particle in canonical orientation, with the c-axis aligned along the beam. The code for this is in kernel_iq.c for C models, but not yet in kernelpy.py for python models.

Thank you, that's great that you are mentioning here oriented 2D images.
For some other experimental data, I'm currently working with the oriented cylinder SASView model for fiting 2D images of nanorods which are oriented perpendicular to the beam by a magnetic or an electric field.
I'm adding two polydispersity functions for phi and theta, and also apply the appropriate rotation to get the c axis (long axis of the particles) perpendicular to the field.
But I have a doubt about how polydispersity is applied, because the orientation is not the canonical one, with the long axis along the beam.
We would like to add the possibility to model the orientation with a single polydispersity function as the system is uniaxial around the applied field.

Could we discuss this in a separate location in the github ?
Sorry, I'm not yet familar with the different ways for categorizing the different topics.
Is there already an appropiate issue to discuss this point ???
Or should I initiate a new one ?

[pkienzle]

I added ticket #695 for orientation support with pure python models.

Orientation dispersion is described in the docs here

Note that the gaussian distribution only supports low dispersion levels. For large Δθ, Δφ, Δψ you need to switch to a cyclic gaussian, or the equivalent Maier-Saupe distribution. These are not available in the standard distributions, but are provided as examples of plugin distributions in sasmodels repository.