Adding truncated octahedron model (python and c) #687

sara-mokhtari · 2025-11-15T10:20:31Z

Description

Adding a new model for truncated octahedrons (python + C).

Fixes #684

How Has This Been Tested?

Unit tests worked.
Other tests were perfomed to validate the model: for truncated octahedrons with different sizes and truncatures, agreement was found with numerical calculation using Debye equation (Debye calculator: https://github.com/FrederikLizakJohansen/DebyeCalculator ). It shows good agreement between the two numerical techniques (Debye based on atomic positions and SasView model based on analytical expressions).

More tests should be made on small q. Indeed, we encouter issues when it comes to small q. More precise mathematical should be used (hospital rule, etc).

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)

…he Nan issue again)

marimperorclerc · 2025-11-15T14:20:54Z

This is a great start for us !!! :) About issues at small q, it is probably because no special treatment for singularities is included in the code. Singularities (i.e. division by zero) are happening for specific directions depending on the shape of the polyhedron (along vertices, edges and normal to the facets). So far, in what we have implemented, the code is running well because all quadrature points are not along any singularity direction. But a better solution would be a nice improvement, by taking into account singularity directions explicitely in the analytical expression of the form factor.

sasmodels/models/octahedron_truncated.c

wpotrzebowski · 2025-11-15T20:27:18Z

sasmodels/models/octahedron_truncated.c

+    	    const double qy = Qy * length_b;
+    	    const double qz = Qz * length_c;
+
+            const double AA = 1./((qy*qy-qz*qz)*(qy*qy-qx*qx))*((qy-qx)*sin(qy*(1.-t)-qx*t)+(qy+qx)*sin(qy*(1.-t)+qx*t))+


(qyqy - qzqz) probably needs to be taken care of. Otherwise we can end up in NaN/Inf scenario

Yes, we have this issue on how to handle singularities numerically.
To be discussed furthermore, in particular for the pure python model version of the truncated octahedron.
We will soon create another pull request for this one. :)
For this truncated octahedron model version written with both .py + .c code, my understanding is that because no points in the gauss-legendre quadrature correspond to a singularity direction, the code is running well without error messages.

sasmodels/models/octahedron_truncated.c

pkienzle · 2025-11-15T16:52:27Z

sasmodels/models/octahedron_truncated.c

+    //HERE: Octahedron formula
+    const double qx = qa * length_a;
+    const double qy = qb * length_b;
+    const double qz = qc * length_c;


My convention has been to use (qx, qy, qz) for the shape rotated relative to the beam and (qa, qb, qc) for the shape in canonical orientation with beam along the c axis. Furthermore, your qx variable is unitless, not 1/Å.

You may be able help the compiler using even more intermediate variables:

const double A = qa * length_a; const double Asq = square(A); const double At = A * t; const double Atc = A * (1 - t); ...

This makes the structure of the equations easier to see:

const double AA = ((B-A)*sin(Btc-At) + (B+A)*sin(Btc+At)) / (2*(Bsq-Csq)*(Bsq-Asq)) + ((C-A)*sin(Ctc-At) + (C+A)*sin(Ctc+At)) / (2*(Csq-Asq)*(Csq-Bsq)); ...

You could further move the factor of 1/2 from each of these equations to the AP equation:

const double AP = 3 * (AA + BB + CC) / (1 - 3*cube(1-t));

Dear Paul,
Thank you so much for your input and clarifications about notations for orientation !
We will improve with Sara this .c file to make it more readable and update it on the branch in sasmodels (adding_trOh_model_pythonC).
Note that we will create soon a new pull request with a pure python model for the same shape (truncated octahedron). :)

pkienzle · 2025-11-15T17:41:21Z

sasmodels/models/octahedron_truncated.c

+    for(int i=0; i<GAUSS_N; i++) {
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
+        double sin_theta, cos_theta;
+        SINCOS(theta, sin_theta, cos_theta);


Use the following to save a little work in the inner loop:

const double qab = q * sin_theta; const double qc = q * cos_theta;

then in the inner loop:

const double qa = qab * cos_phi; const double qb = qab * sin_phi;

Intermediate variables related to qc can also be set in the outer loop.

We kept the initial code here, but added comment lines to explain the notations between rescaled components and initial ones.

pkienzle · 2025-11-15T18:00:51Z

sasmodels/models/octahedron_truncated.py

+    q_z = q\,\cos\theta
+
+θ is the angle between the scattering vector and the z axis of the octahedron (length_c).
+ϕ is the angle between the scattering vector (lying in the xy plane) and the x axis (length_a).


The octahedron uses the axes $(a, b, c)$, hence length_a, length_b and length_c, so I'm slightly bothered by the phrase "the z axis of the octahedron".

I prefer to think of this as an integral over $(q_a, q_b, q_c)$. However, we are integrating over the shape in canonical orientation, with $(q_x, q_y, q_z)$ identical to $(q_a, q_b, q_c)$, so thinking about it as an integral over $(q_x, q_y, q_z)$ is equally valid.

My concern is that the language suggests we are rotating the octahedron rather than the q vector. This is clearly not the case since the octahedron requires three rotation vectors (θ, φ, ψ) to specify its orientation.

So maybe:

θ is the angle between scattering vector and the $z$-axis, and φ is the rotation in the $xy$ plane. The octahedron is in its canonical orientation, with the $c$-axis aligned along $z$ and the $a$-axis aligned along $x$.

OK, we will modify and clarify this in the .py file.

We kept qx, qy, qz in the code but we specified that a is along x, b along y and c along z in the reference orientation.

pkienzle · 2025-11-15T19:19:53Z

sasmodels/models/octahedron_truncated.py

+.. math::
+    q_x = q\,\sin\theta\,\cos\phi, \qquad
+    q_y = q\,\sin\theta\,\sin\phi, \qquad
+    q_z = q\,\cos\theta


The math doesn't match the code. In the code you have q_x = q_a * length_a.

The code uses AA, BB and CC instead of AA, AB, and AC.

We changed into AA, BB, CC.

pkienzle · 2025-11-15T21:02:37Z

sasmodels/models/octahedron_truncated.py

+ϕ is the angle between the scattering vector (lying in the xy plane) and the x axis (length_a).
+
+The normalized form factor in 1D is obtained after averaging over all possible
+orientations and this part is implemented in the same way as in the rectangular prism model.


This is misleading since we are not rotating the shape through all possible orientations in the integral, but instead integrating over a sphere of radius q for a single orientation. This is equivalent to averaging over all possible orientations.

Yes, we will clarify this as well.

We changed the documentation in our last version taking into account what you said.

pkienzle · 2025-11-15T21:05:09Z

sasmodels/models/octahedron_truncated.c

+}
+
+static double
+Iq(double q,


Don't need Iq when you have Fq

Agree with this, but we keep here the initial code as Iq() is called most of the time.

Iq is not use if Fq is present.

The code in generate.py uses the have_Fq attribute to generate either a call to Fq or to Iq:

sasmodels/sasmodels/generate.py

Lines 1063 to 1070 in ed30b0f

if model_info.have_Fq:

pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)

call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars

clear_iq = "#undef CALL_FQ"

else:

pars = ",".join(["_q"] + model_refs)

call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars

clear_iq = "#undef CALL_IQ"

The code in kernel_iq.c checks if call Fq or call Iq is defined:

sasmodels/sasmodels/kernel_iq.c

Lines 499 to 571 in ed30b0f

#if defined(CALL_FQ) // COMPUTE_F1_F2 is true

// unoriented 1D returning <F> and <F^2>

// Note that F1 and F2 are returned from CALL_FQ by reference, and the

// user of the CALL_KERNEL macro below is assuming that F1 and F2 are defined.

double qk;

double F1, F2;

#define FETCH_Q() do { qk = q[q_index]; } while (0)

#define BUILD_ROTATION() do {} while(0)

#define APPLY_ROTATION() do {} while(0)

#define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table)

#elif defined(CALL_FQ_A)

// unoriented 2D return <F> and <F^2>

// Note that the CALL_FQ_A macro is computing _F1_slot and _F2_slot by

// reference then returning _F2_slot. We are calling them _F1_slot and

// _F2_slot here so they don't conflict with _F1 and _F2 in the macro

// expansion, or with the use of F2 = CALL_KERNEL() when it is used below.

double qx, qy;

double _F1_slot, _F2_slot;

#define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)

#define BUILD_ROTATION() do {} while(0)

#define APPLY_ROTATION() do {} while(0)

#define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),_F1_slot,_F2_slot,local_values.table)

#elif defined(CALL_IQ)

// unoriented 1D return <F^2>

double qk;

#define FETCH_Q() do { qk = q[q_index]; } while (0)

#define BUILD_ROTATION() do {} while(0)

#define APPLY_ROTATION() do {} while(0)

#define CALL_KERNEL() CALL_IQ(qk,local_values.table)

#elif defined(CALL_IQ_A)

// unoriented 2D

double qx, qy;

#define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)

#define BUILD_ROTATION() do {} while(0)

#define APPLY_ROTATION() do {} while(0)

#define CALL_KERNEL() CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table)

#elif defined(CALL_IQ_AC)

// oriented symmetric 2D

double qx, qy;

#define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)

double qa, qc;

QACRotation rotation;

// theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.

#define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);

#define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc)

#define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)

#elif defined(CALL_IQ_ABC)

// oriented asymmetric 2D

double qx, qy;

#define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)

double qa, qb, qc;

QABCRotation rotation;

// theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.

// psi and dpsi are only for IQ_ABC, so they are processed here.

const double psi = values[details->theta_par+4];

local_values.table.psi = 0.;

#define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)

#define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)

#define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)

#elif defined(CALL_IQ_XY)

// direct call to qx,qy calculator

double qx, qy;

#define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)

#define BUILD_ROTATION() do {} while(0)

#define APPLY_ROTATION() do {} while(0)

#define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table)

#endif

OK ! Thank you for this.

More precise documentation and some changes of variables (AA, BB, CC) (cf conversations).

Same changes : volume changed in form volume, adding more comments

Issue when writing s*= form_volume(length_a, b2a_ratio,c2a_ratio, t)

Added checks for NaN and infinity values in calculations.

sara-mokhtari · 2025-12-04T10:33:21Z

We updated the codes (both octahedron_truncated.c and octahedron_trucated.py) with more specific documentation and fixed some mistakes. We also added a new image that describes better the truncation parameter of the model.
Feel free to review it again. We explained some of our choices in the previous conversations. If you see remaining issues, let us know. Thank you !

pkienzle · 2025-12-04T15:05:08Z

sasmodels/models/octahedron_truncated.py

+    [{"background": 0, "scale": 1, "length_a": 100, "t": 1, "sld": 1., "sld_solvent": 0.},
+     0.01, 120.57218749910827],
+    [{"background": 0, "scale": 1, "length_a": 100, "t": 1, "sld": 1., "sld_solvent": 0.},
+     [0.01, 0.1], [120.57218749910827, 0.3741758252985113]],


Please note the source of the target values.

Make sure you are using an appropriate number of sigfigs on the target value. For example, if your target values come from a numerical integral only report the digits supported by the integration stopping conditions.

The sasmodels test code is only looking at the first five sigfigs of the target during the comparison. You may include more but for now we won't check them.

You could add a truncated octahedron to shape2sas if you wanted a further source of validation values.

To @pkienzle:

About validation, for information, here we are sharing the numerical validation we have performed by comparing with Debye-calculator for different truncated octahedrons:

This is not yet published but it is a strong validation of the analytical model as the calculation method based on atomic coordinates is completely different.

And thank you for mentioning shape2sas !

About the source of the target value, it is simply obtained by running the model in SASView. We will keep only 5 digits in the test block at the end of the python file.

Keep as many digits as you trust. Even though the current version only tests against 5 digits, some future version may want to test with higher precision.

[Optional] Perhaps you can run the Debye calculator multiple times to estimate the uncertainty of your test value.

[Optional] You could also use mpmath to calculate values using 500 bit arithmetic. The mpmath.quad function can be used to evaluate the numerical integral at a single q point. It may be slow but it only needs to be done once. If you do this, put the validation code into the sasmodels/explore directory.

5 digits are OK with us.

To @pkienzle and @wpotrzebowski:

Sara and I finished to review the two model files. :)
Let us know if you see something else to do ?

And we may merge this pull request to the main branch afterwards ?

Then we will start working on other polyhedron models using new branches and new issue numbers. These models will be pure python models with isotropic average computed using Fibonnaci sphere quadrature and the Lebedev quadrature.
Maybe we will start with the nanoprism models ... or just one simple shape to start with.
Like again the truncated octahedron, but with a pure python model version.
But it will wait after the Xmas vacation time. :)

fix _{ratio} in math blocks

Trying "text" in the math equations, probably didn't work before cause i forgot an empty line between math and the expression

sara-mokhtari and others added 25 commits November 13, 2025 14:28

Add truncated octahedron (python+C) model

11ba48c

Add truncated octahedron (python+C) model updated (image added)

71974fb

Add truncated octahedron (python+C) model updated (PEP8 convention)

e08dac9

Add truncated octahedron (python+C) model updated again (PEP8)

4c49098

Add truncated octahedron (python+C) model updated (test failed)

53a2fd5

Add truncated octahedron (python+C) model updated (fixing tests)

b0c14de

Add truncated octahedron (python+C) model updated (trying to remove t…

fb79959

…he Nan issue again)

Add truncated octahedron (python+C) model updated (NaN issue and test)

fa35edd

Add truncated octahedron (python+C) model updated (rst issue)

8653749

Add truncated octahedron (python+C) model updated (rst issue and ruff)

eec5414

Add truncated octahedron (python+C) model updated

9f29bda

Add truncated octahedron (python+C) model updated

81c6d17

Add truncated octahedron (python+C) model updated (rst issue and ruff)

f01a7c1

Add truncated octahedron (python+C) model updated (rst issue and ruff)

54ed97c

Add truncated octahedron (python+C) model updated

b132d27

Add truncated octahedron (python+C) model updated

6c31722

Add truncated octahedron (python+C) model updated

fcfe74e

Add truncated octahedron (python+C) model updated

6807855

Add truncated octahedron (python+C) model updated

557209d

Update octahedron_truncated.py

76adf42

Update octahedron_truncated.py

f11bc68

Update octahedron_truncated.py

969da09

Update octahedron_truncated.py

20321f6

changed documentation

b871ab1

changed documentation

1834d80

marimperorclerc self-assigned this Nov 15, 2025

add new category by changing category in octohedron model

679262c

wpotrzebowski reviewed Nov 15, 2025

View reviewed changes

sasmodels/models/octahedron_truncated.c Show resolved Hide resolved

wpotrzebowski reviewed Nov 15, 2025

View reviewed changes

sasmodels/models/octahedron_truncated.c Show resolved Hide resolved

wpotrzebowski reviewed Nov 15, 2025

View reviewed changes

sasmodels/models/octahedron_truncated.c Outdated Show resolved Hide resolved

pkienzle reviewed Nov 15, 2025

View reviewed changes

marimperorclerc and others added 7 commits November 21, 2025 15:04

update of octa-truncated.png image file

5cac0ab

Merge branch 'master' into adding_trOh_model_pythonC

ff9442c

Changing the documentation in the python code

6bf1b94

More precise documentation and some changes of variables (AA, BB, CC) (cf conversations).

Updating the c code

12a2efb

Same changes : volume changed in form volume, adding more comments

Changing octahedron_truncated c and py code: changing documentation

47b02ba

Changing octahedron_truncated.c

42eea23

Issue when writing s*= form_volume(length_a, b2a_ratio,c2a_ratio, t)

Handle NaN and infinity in octahedron calculations

de8df47

Added checks for NaN and infinity values in calculations.

pkienzle approved these changes Dec 4, 2025

View reviewed changes

sara-mokhtari and others added 6 commits December 10, 2025 10:07

Fix test values (sig)

7b6bc2c

Fix LaTeX formatting in octahedron_truncated.py

0666f47

Update octahedron_truncated.py

63e5f69

fix _{ratio} in math blocks

latex doc (math expressions) in octahedron_truncated.py

dd32c89

Trying "text" in the math equations, probably didn't work before cause i forgot an empty line between math and the expression

Fix LaTeX formatting in octahedron_truncated.py (again)

5afaff5

latex doc last fix

7a444c7

	if model_info.have_Fq:
	pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
	call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars
	clear_iq = "#undef CALL_FQ"
	else:
	pars = ",".join(["_q"] + model_refs)
	call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
	clear_iq = "#undef CALL_IQ"

	#if defined(CALL_FQ) // COMPUTE_F1_F2 is true
	// unoriented 1D returning <F> and <F^2>
	// Note that F1 and F2 are returned from CALL_FQ by reference, and the
	// user of the CALL_KERNEL macro below is assuming that F1 and F2 are defined.
	double qk;
	double F1, F2;
	#define FETCH_Q() do { qk = q[q_index]; } while (0)
	#define BUILD_ROTATION() do {} while(0)
	#define APPLY_ROTATION() do {} while(0)
	#define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table)

	#elif defined(CALL_FQ_A)
	// unoriented 2D return <F> and <F^2>
	// Note that the CALL_FQ_A macro is computing _F1_slot and _F2_slot by
	// reference then returning _F2_slot. We are calling them _F1_slot and
	// _F2_slot here so they don't conflict with _F1 and _F2 in the macro
	// expansion, or with the use of F2 = CALL_KERNEL() when it is used below.
	double qx, qy;
	double _F1_slot, _F2_slot;
	#define FETCH_Q() do { qx = q[2q_index]; qy = q[2q_index+1]; } while (0)
	#define BUILD_ROTATION() do {} while(0)
	#define APPLY_ROTATION() do {} while(0)
	#define CALL_KERNEL() CALL_FQ_A(sqrt(qxqx+qyqy),_F1_slot,_F2_slot,local_values.table)

	#elif defined(CALL_IQ)
	// unoriented 1D return <F^2>
	double qk;
	#define FETCH_Q() do { qk = q[q_index]; } while (0)
	#define BUILD_ROTATION() do {} while(0)
	#define APPLY_ROTATION() do {} while(0)
	#define CALL_KERNEL() CALL_IQ(qk,local_values.table)

	#elif defined(CALL_IQ_A)
	// unoriented 2D
	double qx, qy;
	#define FETCH_Q() do { qx = q[2q_index]; qy = q[2q_index+1]; } while (0)
	#define BUILD_ROTATION() do {} while(0)
	#define APPLY_ROTATION() do {} while(0)
	#define CALL_KERNEL() CALL_IQ_A(sqrt(qxqx+qyqy), local_values.table)

	#elif defined(CALL_IQ_AC)
	// oriented symmetric 2D
	double qx, qy;
	#define FETCH_Q() do { qx = q[2q_index]; qy = q[2q_index+1]; } while (0)
	double qa, qc;
	QACRotation rotation;
	// theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
	#define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
	#define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc)
	#define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)

	#elif defined(CALL_IQ_ABC)
	// oriented asymmetric 2D
	double qx, qy;
	#define FETCH_Q() do { qx = q[2q_index]; qy = q[2q_index+1]; } while (0)
	double qa, qb, qc;
	QABCRotation rotation;
	// theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
	// psi and dpsi are only for IQ_ABC, so they are processed here.
	const double psi = values[details->theta_par+4];
	local_values.table.psi = 0.;
	#define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
	#define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
	#define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)

	#elif defined(CALL_IQ_XY)
	// direct call to qx,qy calculator
	double qx, qy;
	#define FETCH_Q() do { qx = q[2q_index]; qy = q[2q_index+1]; } while (0)
	#define BUILD_ROTATION() do {} while(0)
	#define APPLY_ROTATION() do {} while(0)
	#define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table)
	#endif

Adding truncated octahedron model (python and c) #687

Are you sure you want to change the base?

Adding truncated octahedron model (python and c) #687

Uh oh!

Conversation

sara-mokhtari commented Nov 15, 2025

Description

How Has This Been Tested?

Review Checklist:

Uh oh!

marimperorclerc commented Nov 15, 2025

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

sara-mokhtari commented Dec 4, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants