Convolution #59
Convolution #59
Conversation
henrikjacobsenfys
added
[scope] enhancement
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## develop #59 +/- ##
===========================================
- Coverage 96.60% 95.89% -0.71%
===========================================
Files 12 14 +2
Lines 412 707 +295
===========================================
+ Hits 398 678 +280
- Misses 14 29 +15
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
|
||
| # Handle offset | ||
| if offset is None: | ||
| off = 0.0 |
There was a problem hiding this comment.
really bad name for offset - please consider changing to a meaningful name like _offset or rename the input offset to x_offset
There was a problem hiding this comment.
Good point! I've renamed it to offset_float, since it's the offset as a float.
| ) | ||
|
|
||
| # Handle temperature | ||
| T = None |
There was a problem hiding this comment.
Avoid single letter variable names unless used in tight loops/list comprehension.
There was a problem hiding this comment.
I don't know why it was even there - I put a decent amount of time into making sure I could pass Parameters or floats directly to _detailed_balance_factor, so I've just deleted it. Thanks!
| and isinstance(resolution_component, Gaussian) | ||
| ): | ||
| if isinstance(sample_component, Gaussian): | ||
| G, L = sample_component, resolution_component |
There was a problem hiding this comment.
Avoid single letter variables here as well, but in this case it's a bit more allowable, since all the usage is within several lines of code.
| ) -> List[DeltaFunction]: | ||
| """Return a list of DeltaFunction instances contained in `model`. | ||
|
|
||
| Args: |
There was a problem hiding this comment.
All other args: in docstring are lower case. Please be consistent :)
| if upsample_factor == 0: | ||
| # Check if the array is uniformly spaced. | ||
| x_diff = np.diff(x) | ||
| is_uniform = np.allclose(x_diff, x_diff[0]) | ||
| if not is_uniform: | ||
| raise ValueError( | ||
| "Input array `x` must be uniformly spaced if upsample_factor = 0." | ||
| ) | ||
| x_dense = x | ||
| else: | ||
| # Create an extended and upsampled x grid | ||
| x_min, x_max = x.min(), x.max() | ||
| span = x_max - x_min | ||
| extra = extension_factor * span | ||
| extended_min = x_min - extra | ||
| extended_max = x_max + extra | ||
| num_points = len(x) * upsample_factor | ||
| x_dense = np.linspace(extended_min, extended_max, num_points) |
There was a problem hiding this comment.
This is a pretty complex grid routine. Please consider extracting it to a helper function, where it can be properly unit tested.
|
|
||
| Args: | ||
| x : np.ndarray | ||
| 1D array of x values where the convolution is evaluated. |
There was a problem hiding this comment.
The docstring doesn't clearly explain:
- What physical quantities
xrepresents (energy transfer?) - Unit consistency requirements
- What happens if units don't match
Maybe add some examples? Just something to consider. Maybe imaging users know this by heart :)
>>> x = np.linspace(-10, 10, 1000) # Energy transfer in meV
>>> sample = Gaussian(center=0, width=1.0, area=1.0, unit="meV")
>>> resolution = Gaussian(center=0, width=0.5, area=1.0, unit="meV")
>>> result = convolution(x, sample, resolution)
| ) -> Tuple[bool, np.ndarray]: | ||
| """ | ||
| Attempt an analytic convolution for component pair (sample_component, resolution_component). | ||
| Returns (True, contribution) if handled, else (False, zeros). |
There was a problem hiding this comment.
You might want to consider improving the return description
Returns:
Tuple[bool, np.ndarray]:
- bool: True if analytical convolution was computed, False otherwise
- np.ndarray: The convolution result if computed, or zeros if not handled
| The evaluated Voigt profile values at x. | ||
| """ | ||
|
|
||
| return area * voigt_profile(x - center, g_width, l_width) |
There was a problem hiding this comment.
scipy's voigt_profile is normalized, so multiplying by area (which is the product of the two component areas) should work. However, this needs verification with unit tests. Please add some.
| if isinstance(resolution_component, DeltaFunction): | ||
| return True, resolution_component.area.value * sample_component.evaluate( | ||
| x - resolution_component.center.value - off | ||
| ) |
There was a problem hiding this comment.
Shouldn't this be x - off - sample_component.center.value for consistency? The same, but looks more proper maybe.
does offset shift the sample model, resolution model, or the convolution result? Would be good to mention that.
There was a problem hiding this comment.
You're right, offset needs to be explained better. It shifts the convolution result. The center should be the sum of the centers of the sample, resolution and offset. In this case, the sample center is included in the sample_component.evaluate() call, so it needs to be shifted by the other two terms (resolution and offset)
| if len(x_dense) % 2 == 0: | ||
| x_even_length_offset = -0.5 * dx | ||
| else: | ||
| x_even_length_offset = 0.0 |
There was a problem hiding this comment.
This looks like a subtle correction - can you add some comment on why even-length arrays need this correction? Is it correct for all symmetric/asymmetric cases?
There was a problem hiding this comment.
It is indeed subtle - I've added an explanation:
# The convolution of two arrays of length N is of length 2N-1. When using 'same' mode, only the central N points are kept,
# so the output has the same length as the input.
# However, if N is even, the center falls between two points, leading to a half-bin offset.
# For example, if N=4, the convolution has length 7, and when we select the 4 central points we either get
# indices [2,3,4,5] or [1,2,3,4], both of which are offset by 0.5*dx from the true center at index 3.5.
| @@ -0,0 +1,846 @@ | |||
| import numpy as np | |||
There was a problem hiding this comment.
Potentially needs additional tests for:
- Offset with detailed balance
- Non-uniform spacing detection accuracy
- Edge cases (empty models, zero widths)
- Voigt area conservation
- Very wide/narrow peaks triggering warnings
There was a problem hiding this comment.
Yeah, looking at the code, I should probably redo the whole test file. The vast majority of the tests are not actually unit tests. For example, I'm testing each case in _try_analytic_pair indirectly by testing _analytical_convolution.
|
Okay, the main file could use another round of reviews before I clean up the mess of unit tests :) |
Implements convolution of SampleModels and Components with each other, including detailed balancing. The code is based on the SampleModel branch, not develop.
I have very extensive tests of the convolutions. Many of them were written a long time ago, and so they can probably be improved. However, before doing that, a review of the main code would be very useful.