Adaptive integration for truncated octahedron

[pkienzle]

WARNING Applies after #658. Change base to master before merging.

Modified the truncated octahedron function to use simple adaptive integration.

I reordered the axes for speed and accuracy, and adjusted the equations to minimize computation.

Note that the formula is suspicious since the limit as q -> 0 diverges. However setting t=1 (no truncation) and simplifying, the form matches that in 1, Appendix A, which also diverges. This cites 2, but I don't have access to verify the original derivation.

1. Wei-Ren Chen et al. "Scattering functions of Platonic solids".
   In: Journal of Applied Crystallography - J APPL CRYST 44 (June 2011).
   DOI: 10.1107/S0021889811011691

2. Stokes, A. R. & Wilson, A. J. C. (1942). Math. Proc. Camb. Philos.
   Soc. 38, 313–322. https://doi.org/10.1017/S0305004100021988

[marimperorclerc]

Today, I started updating the truncated octahedron model in order to have the convention t=0 for no truncature.

• figure octa-truncated.png is updated according to this new truncature convention in the \img folder on the adaptive-octahedron branch
• I also modified the octahedron-truncated.c c code. Is is only the exchange of t and tinv everywhere in the c code. Please double check it. I hope I didn't messed up somewhere ...

After, we will need to:

• change the name of the model to truncated_octahedron (not done yet)
• update the truncated_octahedron.py accordingly (description of the model and new name) (not done yet)

[sara-mokhtari]

Note: in the last commit, Marianne introduced tinv (inverse) that is defined as: tinv = 1 - t, t being the new truncation parameter defined as: 0 < t < 0.5 (0 being the octahedron, 0.5 being the cuboctahedron).

I have made some minors updates:

• I updated the truncation parameter name and interval in the .py file: from "t", "[0.5, 1]" to "truncation", "[0, 0.5]" (@marimperorclerc the tests failed simply because the parameter wasn't updated in the .py file)
• I updated the documentation as well
• I updated the name of the model: from "octahedron_truncated" to "truncated_octahedron"

Now I will look at the remaining issues (see TO DOs in the code), especially the instabilities when q tends to 0. From what I've observed, they're are quite significant for now.

[sara-mokhtari]

I must apologize because I forgot to ask something regarding the final version of the truncated_octahedron model before the code freeze and I'm not sure on how to process.

The initial truncated octahedron branch was merged a while ago and the code is not updated as we wish (especially the truncation parameter that is not consistent with the new tetrahedron model). The updated code is in the new branch called adaptative_octahedron, but in this branch, the code also contains a new adapative method that is still being studied. If I understood correctly, we can work on the adaptative method after the code freeze. So how should I update the model ?
By creating a new branch that will only contains the small changes for the consistency ?
The best would have been to change the code in the initial branch but it was merged already ...

Thank you in advance for the help.

@marimperorclerc @pkienzle @butlerpd

[butlerpd]

I have a couple of comments/questions I guess @sara-mokhtari - In general, if this was an update to a model that has been merged, you would normally make a new branch which should contain the 'incorrect" code and then you would alter that code. Is that what you did? or did you create a completely new model? If you did the 2nd than you are not fixing the old code at all. There will be two versions if we merge which would not be ideal in any case.?

Assuming your branch is changing the code for the existing truncated octahedron model then yes, you may want to make a new branch and change that model with only the changes that are not part of redoing the integration over orientations (leaving it as before) then making that a PR which can hopefully be merged before the code freeze. Depending on how the changes were committed and how comfortable you are with cherry picking, you could perhaps just cherry pick those commits to your new branch?

I would ask however whether there is not a world where the new integration cannot have a working, if not ideal, version ready for the code freeze?

To your question about how it should have been done, IMO I would always start with a clean branch if adding some fixes to the existing model. The other PR should just include the new integration code, but of course you need at least one model that will implement the new integration to test it and make sure it is useful. But I would have restricted all changes to the model in that branch to just what is needed to change the method of averaging. Others may suggest other approaches, but in my experience it is easiest to keep the pull requests "clean" with just the code needed for the stated purpose of that PR.

[sara-mokhtari]

Thank you a lot @butlerpd for your help !
So I guess it's pretty much what I had in mind. I will create a new branch and cherry pick the changes to update the already existing model. I hope renaming the model will note create issues though (octahedron_truncated to truncated octahedron).

Therefore, I will leave the adapatative-octahedron branch for its main purpose.
I can't really tell you if the new integration can be ready before the code freeze. I will try to work on it, but since Marianne is on vacations, I can't really promise anything ... So in conclusion: I will try to work on the new integration, if I can't work it out fast enough, I will make a PR for the other branch to at least have the model with the updated parameters and name. Hope this is okay, if not let me know !

[butlerpd]

From the purely technical point of view, changing the name of the model or of the parameters should not cause a problem right now because this model has never been in a released version. Once a model has been in a release version changing those names becomes difficult. The only other technical issue, specifically with the model name is that it is best of the file name be the same as the name attribute in the file. I believe there was a fix a while back to deal with times when that is not the case, but I'm not sure it has been well-tested? Moreover it should be best practice anyway IMO.

From the perspective of being acceptable for merging, the main issue with names is that they should follow the standard convention as much as possible to make it easier for users (and also future developers). I have not thought about the polyhedron models so cannot comment on truncated_octahedron vs octahedron_truncated. I assume that has been discussed already?

[sara-mokhtari]

Yes thank you @butlerpd, indeed the name of the file was discussed before, and it was chosen to put "truncated_octahedron" and "truncated_tetrahedron" for consistency with other models. The name attribute was also changed in the file, and the other changes were discussed as well !

[pkienzle]

The documentation says length for the truncated octahedron is from the center of the object. This is more like a radius than a length. For the cylinder model length goes from end to end.

I'm noting this as an inconsistency. Users may find radius less confusing.

[pkienzle]

Check against speed and accuracy targets:

$ python explore/check_adaptive.py truncated_octahedron                                    
== Speed and accuracy tests for all adaptive integration models ==
* target evaluation time is 2 s (running on a mac M2 chip)
*     q in [1e-5, 1] with 40 points per decade for 201 points total
*     warns if the adaptive model is 2x slower than a 76-point gaussian
* large models tested against 5000 point gaussian integration
*     q=[5e-4, 1e-3, 2e-3] with tol=1e-5 relative (measured q)
*     q=[0.01, 0.1] with tol=0.2 relative (slit resolution limits)
* small models tested against 5000 point gaussian integration
*     q in [1e-3, 1] with 1 points per decade
!!!! These tests run very slowly --- don't use as part of CI !!!!



=== small rods: a=20 b=40 c=200 ===


=== small disks: a=180 b=200 c=40 ===


=== small cubes: a=200 b=200 c=200 ===


=== big rods: a=1000 b=2000 c=200000 ===
truncated_octahedron background=0 sld=1 sld_solvent=0 length_a=500.0 b2a_ratio=2.0 c2a_ratio=200.0 truncation=0.25
! ** truncated_octahedron is slow: 2.3 s for 201 points in [1e-05, 1.0]


=== big disks: a=180000 b=200000 c=1000 ===
truncated_octahedron background=0 sld=1 sld_solvent=0 length_a=90000.0 b2a_ratio=1.1111111111111112 c2a_ratio=0.005555555555555556 truncation=0.25
  qvalue [1.00000000000000e-01]
  target [3.24338546013817e-02]
  actual [5.01627633867576e-02]
  relerr [5.46617384929039e-01]
! ** truncated_octahedron is slow: 3.3 s for 201 points in [1e-05, 1.0]


=== big cubes: a=200000 b=200000 c=200000 ===
truncated_octahedron background=0 sld=1 sld_solvent=0 length_a=100000.0 b2a_ratio=1.0 c2a_ratio=1.0 truncation=0.25
  qvalue [1.00000000000000e-02 1.00000000000000e-01]
  target [2.87810808892225e+00 3.11597219157298e-04]
  actual [7.96721938846658e-01 5.75473502932781e-05]
  relerr [7.23178590160245e-01 8.15314942640013e-01]
! ** truncated_octahedron is slow: 3.8 s for 201 points in [1e-05, 1.0]

[sara-mokhtari]

In the last commit (cee3e56), I changed "length_a", "length_b", "length_c" to "radius_a", "radius_b", "radius_c". In the documentation, I explained that they are are the distances from the center of the general octahedron to its 6 vertices and that they are equivalent to the circumradiuses of the general octahedron along the three directions.

I suppose that now the last thing to do is verify the adaptative method before the merge.

[sara-mokhtari]

I have compared the adaptative method to the non adaptative method (with different numbers of points: 20, 66, 150) and their computational time. I did it for two relatively big sizes (radius_a = 100 nm and 1000 nm).


@pkienzle @marimperorclerc

[pkienzle]

Regarding runtime speed for large shapes, I put a limiter on the outer loop so it doesn't exceed 76 points. Please rerun the test for speed and accuracy to make sure the error is still tolerable.

Regarding parameter renaming or changing parameter interpretation, it is worthwhile settling on something now. If we change the parameter name or interpretation later then we need to adapt the model parameters from old to new whenever we load an old model. Code for this will have to be added to sasmodels/convert.py here:

sasmodels/sasmodels/convert.py

Line 166 in 55c5b9c

def _hand_convert(name, oldpars, version=(3, 1, 2)):

Regarding low q error, it seems that the problems happen well below the qr range of interest to SAS, so no need to fix it. [Edit] Fixing the low q error might allow you to use single precision calculations running on the GPU.

[sara-mokhtari]

@pkienzle I ran the test again with the new version of the code (with the outer loop limited to 76 points). It's called "C kernel adaptative (new)" on the plot, the accuracy is still good and the computational time decreased a lot.