Skip to content

Implement Perrichon2022 McStas TOSCA function as lookup table#43

Open
ajjackson wants to merge 10 commits intomainfrom
perrichon
Open

Implement Perrichon2022 McStas TOSCA function as lookup table#43
ajjackson wants to merge 10 commits intomainfrom
perrichon

Conversation

@ajjackson
Copy link
Member

@ajjackson ajjackson commented Mar 14, 2025

This scheme is a "scaled" lookup table:

  • At each energy-transfer a Gaussian fit was performed to determine a nominal width parameter "sigma". The relationship between sigma and energy-transfer is in turn described by a fitted polynomial.
  • Each row in the lookup table corresponds to a different energy-transfer.
  • The columns correspond to a zero-centered energy mesh (i.e. these are kernels) as a multiple of sigma.
  • This scaling keeps the overall widths in the table fairly consistent, allowing clean linear interpolation of the shape.
  • x-values are shifted and scaled appropriately in order to produce interpolated peak functions at their true locations and scales.

The table is quite compact but still a bit chunky to include in the .yaml files. Instead we store all the model parameters in a .npz file and point to this from the .yaml.

Some pre-processing was performed on the original data:

  • smoothing with a Scipy "smoothing spline"
  • truncation to 10σ with a smooth Tukey window function
  • normalisation to constant area

@ajjackson ajjackson changed the base branch from main to implement_kernels_convolutions March 14, 2025 22:02
@ajjackson
Copy link
Member Author

Results in usage are visually pretty similar to the AbINS model but we do see some divergence especially at high energy transfer.

plot

@ajjackson ajjackson changed the base branch from implement_kernels_convolutions to main April 16, 2025 09:49
@ajjackson ajjackson marked this pull request as ready for review April 16, 2025 10:10
Copy link
Collaborator

@RastislavTuranyi RastislavTuranyi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some comments, mostly to do with the most recent changes from #39, but there are two major things that stick out:

  1. No docstrings for get_kernel and get_peak
  2. No unit tests

Could you add these? That said, tho, not sure how viable unit tests are in this case, ig they'll just have to be reference values to test for regression.

Comment on lines 117 to 118
assert len(omega_q.shape) == 2 and omega_q.shape[1] == 1
energy = omega_q
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
assert len(omega_q.shape) == 2 and omega_q.shape[1] == 1
energy = omega_q
assert points.ndim == 2 and points.shape[1] == 1
energy = points

First, related to the discussion #39, is this check necessary? Secondly, are you sure you want to keep it as a 2D array?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The check seems a good idea for now, this should maybe become something more systematic once we really start using Q as well.

The 2-D array leads to scale factors with the right dimensionality later... Arguably "for the wrong reason"? The 2-D gets expanded to 3-D in the meshgrid step.

It might indeed be a bit clearer to reduce to 1 and than expand dimensions where needed.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I vaguely remember we discussed both of these things elsewhere and decided to do checks and to use Nx1 2D arrays as inputs? If that's the case, I think this can be ignored

omega_q: Float[np.ndarray, 'sample dimension=1'],
mesh: Float[np.ndarray, 'mesh'],
) -> Float[np.ndarray, 'sample mesh']:
shifted_meshes = [mesh - energy for energy in omega_q[:, 0]]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
shifted_meshes = [mesh - energy for energy in omega_q[:, 0]]
try:
points = points[:, 0]
except IndexError as e:
raise InvalidPointsError(
f'The provided array of points (shape={points.shape}) is not valid. The points '
f'array must be a Nx1 2D array where N is the number of energy transfers.'
) from e
shifted_meshes = [mesh - energy for energy in points]

@RastislavTuranyi RastislavTuranyi added the enhancement New feature or request label May 6, 2025
@RastislavTuranyi RastislavTuranyi linked an issue May 6, 2025 that may be closed by this pull request
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Add TOSCA model by A. Perrichon

2 participants