diff --git a/src/dodal/common/general_maths/absorber_geometry.py b/src/dodal/common/general_maths/absorber_geometry.py new file mode 100644 index 00000000000..6a45b5b5ad0 --- /dev/null +++ b/src/dodal/common/general_maths/absorber_geometry.py @@ -0,0 +1,130 @@ +from collections.abc import Callable +from typing import Annotated, Literal, Self + +from pydantic import ( + BaseModel, + Field, + PrivateAttr, + StrictFloat, + StrictInt, + field_validator, + model_validator, +) + +from .arithmetic_conversions import ( + convert_cm_to_mm, + convert_microns_to_cm, + convert_mm_to_cm, + get_straight_line_y, +) + +SupportedThicknessUnits = Annotated[ + Literal["cm", "mm", "um", "micron"], "Supported thickness units" +] + + +class FoilGeometry(BaseModel): + """Represents flat shape of a specific foil absorber. + + Initialised with a foil thickness quantity in a supported thickness unit and, + where necessary, converting this to a cm thickness. + N.B. The missing "air" slot in a foil wheel is represented as a canonical / system "OUT" position, + of the absorber - rather than as a zero thickness foil. + + Args: + unit (str): one element from the set of supported choices - "cm", "mm", "um" or "micron" + numerical_value (float): non-zero positive - physical thickness of the foil in the given units + """ + + # Class-level private dictionary - defined once only when the class is loaded + _convertors: dict[str, Callable[[float], float]] = { + "cm": lambda t: t, + "mm": convert_mm_to_cm, + "um": convert_microns_to_cm, + "micron": convert_microns_to_cm, + } + + unit: str + numerical_value: Annotated[float, Field(gt=0.0)] + _foil_thickness: float = PrivateAttr(default=0.0) + + @model_validator(mode="after") + def _calculate_thickness(self) -> "FoilGeometry": + # Look up directly from the conversion dictionary + self._foil_thickness = self._convertors[self.unit](self.numerical_value) + return self + + def get_thickness_cm(self) -> float: + return self._foil_thickness + + +class WedgeGeometry(BaseModel): + """Represents a semi infinite triangular prism wedge shape, growing out indefinitely from a tip. + Initialised with cotangent of the wedge angle and offset between motor axis zero and wedge tip. + N.B. absorber thickness are best worked out in cm (science), motor positions prefer mm (controls). + Values used for inputs come from beamline scientist calibrations on real absorber. + + Args: + tip_mm (float): unrestricted value, point on motor axis scale where modelled taper reaches zero thickness. + taper_cotangent (float): taper angle cotangent, sign depends on motor direction convention, bounded to keep wedge angle thin. + """ + + tip_mm: StrictFloat | StrictInt + taper_cotangent: StrictFloat | StrictInt + _taper_gradient: float = PrivateAttr(default=0.0) + + @field_validator("taper_cotangent") + @classmethod + def validate_taper_is_thin(cls, value: float) -> float: + bound = 5.0 # round number cotangent limit for thin wedge + if -bound < value < bound: + msg: str = f"taper_contangent must correspond to thin wedge, of magnitude {bound} or larger" + raise ValueError(msg) + return value + + @model_validator(mode="after") + def _calculate_gradient(self) -> Self: + self._taper_gradient = 1.0 / self.taper_cotangent + return self + + def motor_position_mm_for_thickness_cm( + self, thickness_cm: Annotated[float, Field(ge=0.0)] + ) -> float: + """Hypothetical motor position at which the wedge thickness would take on the requested value. + N.B. absorber thickness are best worked out in cm (science), motor positions prefer mm (controls). + + Args: + thickness_cm: (float) non-negative float, zero permitted - the required target thickness in cm + + Returns: + motor position in mm where a mathematical model wedge puts the target thickness in the x-ray beam + """ + _line_offset = ( + self.tip_mm + ) # Given that thickness (line x) is zero at this motor position (line y) + _line_gradient = ( + self.taper_cotangent + ) # small changes in thickness lead to large changes in position + _x = convert_cm_to_mm(thickness_cm) # do the line calculations in mm + return get_straight_line_y(_line_offset, _line_gradient, _x) + + def thickness_cm_at_motor_position_mm(self, motor_position_mm: float) -> float: + """Thickness of modelled wedge a given motor position in most useful units. + N.B. The wedge geometry is not bounded by lateral position limits, but does reject, + negative thickness, + unlikely large taper angles (beyond about 11.31 degrees). + + Args: + motor_position_mm: (float) unrestricted float, zero permitted - the required target thickness + + Returns: + position in mm on the motor axis scale where a mathematical wedge would have the target thickness + """ + _line_offset = ( + -self.tip_mm * self._taper_gradient + ) # not immediately obvious - but after some algebra + _line_gradient = self._taper_gradient + _thickness_mm = get_straight_line_y( + _line_offset, _line_gradient, motor_position_mm + ) + return convert_mm_to_cm(_thickness_mm) diff --git a/src/dodal/common/general_maths/arithmetic_conversions.py b/src/dodal/common/general_maths/arithmetic_conversions.py index d0465337938..47e8ba96be0 100644 --- a/src/dodal/common/general_maths/arithmetic_conversions.py +++ b/src/dodal/common/general_maths/arithmetic_conversions.py @@ -1,5 +1,6 @@ """Provides functions to convert between common units for attenuation.""" +from numpy.polynomial import polynomial from pydantic import StrictFloat, validate_call @@ -105,3 +106,21 @@ def convert_ev_to_kev(energy_ev: StrictFloat) -> float: energy_ev (float): a value of kilo-electron volts """ return energy_ev * 1e-3 + + +@validate_call +def get_straight_line_y( + line_offset: StrictFloat, line_gradient: StrictFloat, x: StrictFloat +) -> float: + """Convenient internal conversion method which delegates to numpy for a 1st order polynomial (line). + + The traditional y = mx + c to compute y. + + Args: + line_offset: the offset, where the modelled line intercepts the abscissa. + line_gradient: the gradient of the modelled line + x: the input coordinate of the point on the line, from which y is computed. + """ + _coefficients = [line_offset, line_gradient] + _y = polynomial.polyval(x, _coefficients) + return float(_y) diff --git a/src/dodal/common/general_maths/material_absorption_maths.py b/src/dodal/common/general_maths/material_absorption_maths.py index 33741d6ded7..fa3ff1f6e37 100644 --- a/src/dodal/common/general_maths/material_absorption_maths.py +++ b/src/dodal/common/general_maths/material_absorption_maths.py @@ -1,8 +1,11 @@ +from collections.abc import Callable from typing import Annotated +from numpy.polynomial import polynomial from pydantic import ( Field, StrictFloat, + StrictInt, validate_call, ) @@ -12,18 +15,108 @@ ) +class AbsorptionCalculator: + """Interface for calculating absorption (due to mass attenuation) per cm as a function of x-ray energy in keV. + Reports back a natural log based absorption factor as function of x-ray energy. + """ + + def __init__(self, _calculation: Callable[[float], float]): + self._calculate = _calculation + + def absorption_coefficient_per_cm( + self, energy_kev: Annotated[StrictFloat | StrictInt, Field(gt=0)] + ) -> float: + """Logarithmic contribution to x-ray absorption per cm (depth into absorber) at a particular photon energy. + + Notes: + 1) The cm as a typical unit length scale is a de facto standard in the field. + 2a) The absorption coefficient is in factors of e. + 2b) Barnett absorption units are only used once this coefficient is combined with depth in cm. + 3a) Negative contributions to a sum have to be permitted so the output, + of a specific calculator instance, can be a float of either sign. + 3b) Sub-classes may implicitly restrict this expectation. + Real materials the final absorption will be positive. + + Args: + energy_kev (StrictFloat|StrictInt): The individual energy per photon in kilo-electronvolts + + Returns: + (float): Model adjustment of attenuation per cm of absorber material depth. + """ + return self._calculate(energy_kev) + + +class CompoundAbsorptionCalculator(AbsorptionCalculator): + """Advanced physics model for mass attenuation per cm as a function of x-ray energy in keV. + Applies effects from several absorption term calculations. + + Attributes: + contributers (list[AbsorptionCalculator]) + """ + + def __init__(self, contributions: list[AbsorptionCalculator]): + super().__init__( + # lambda k is energy in keV + lambda k: sum(c.absorption_coefficient_per_cm(k) for c in contributions) + ) + + +class PolynomialAbsorptionCorrection(AbsorptionCalculator): + """Corrective model for mass attenuation per cm as a function of x-ray energy in keV. + Provides terms for correcting a baseline mass attenuation. + + Attributes: + corrective_terms (float): Polynomial coefficients to correct the baseline modelled absorption per cm + """ + + def __init__(self, coefficients_per_cm: list[float]): + + def _calculate_correction(energy_kev: float) -> float: + correction = polynomial.polyval(energy_kev, coefficients_per_cm) + return float(correction) # numpy did not specify float as the return type + + super().__init__(_calculate_correction) + + +class SingleRollOffAbsorptionCalculator(AbsorptionCalculator): + """Simplest physics model for mass attenuation per cm as a function of x-ray energy in keV. + Typically appropriate where one element dominates the absorption, + and energies passed in as calculation arguments are above the elements absorption resonance edge. + + Attributes: + material_factor_per_cm (StrictFloat): Positive non-zero material constant + (hypothetical trend offset equivalent to an absorption per cm at 1 keV) + roll_off: negative exponent of energy dependence above the resonant edge. + """ + + def __init__( + self, + material_factor_per_cm: Annotated[StrictFloat | StrictInt, Field(gt=0)], + roll_off: Annotated[StrictFloat | StrictInt, Field(lt=0)], + ): + # lambda k is energy in keV + super().__init__( + lambda k: photon_mass_attenuation_per_unit_length( + k, material_factor_per_cm, roll_off + ) + ) + + @validate_call def photon_mass_attenuation_per_unit_length( - energy_kev: Annotated[StrictFloat, Field(gt=0)], - photon_absorption_factor_per_unit_length: Annotated[StrictFloat, Field(gt=0)], - energy_dependence_exponent: Annotated[StrictFloat, Field(lt=0)], + energy_kev: Annotated[StrictFloat | StrictInt, Field(gt=0)], + photon_absorption_factor_per_unit_length: Annotated[ + StrictFloat | StrictInt, Field(gt=0) + ], + energy_dependence_exponent: Annotated[StrictFloat | StrictInt, Field(lt=0)], ) -> float: """Calculates mass attenuation per unit length. + See for example: https://en.wikipedia.org/wiki/Mass_attenuation_coefficient. + Args: energy_kev (StrictFloat): X-ray energy in keV (positive values only). - photon_absorption_factor_per_unit_length (StrictFloat): Factors of e, - in X-ray flux reduction, per unit depth of absorber (positive values only). + photon_absorption_factor_per_unit_length (StrictFloat): Logarithmic constant, scaled in factors of e. energy_dependence_exponent (StrictFloat): Roll off of absorption factor (negative values only). Returns: @@ -36,8 +129,8 @@ def photon_mass_attenuation_per_unit_length( @validate_call def attenuation_at_depth_cm( - depth_cm: Annotated[StrictFloat, Field(ge=0)], - absorption_coefficient_per_cm: Annotated[StrictFloat, Field(ge=0)], + depth_cm: Annotated[StrictFloat | StrictInt, Field(ge=0)], + absorption_coefficient_per_cm: Annotated[StrictFloat | StrictInt, Field(gt=0)], ) -> float: """Calculates attenuation in Barnett units, where 1000 Bn equivalent to 1/e, 0Bn to 1 and 2000 Bn to 1/(e^2). @@ -61,8 +154,8 @@ def attenuation_at_depth_cm( @validate_call def thickness_cm_required_to_attenuate( - target_attenuation_bn: Annotated[StrictFloat, Field(ge=0)], - absorption_coefficient_per_cm: Annotated[StrictFloat, Field(gt=1.0e-8)], + target_attenuation_bn: Annotated[StrictFloat | StrictInt, Field(ge=0)], + absorption_coefficient_per_cm: Annotated[StrictFloat | StrictInt, Field(gt=1.0e-8)], ) -> float: """Calculates material depth in cm. @@ -72,7 +165,6 @@ def thickness_cm_required_to_attenuate( absorption_coefficient_per_cm (StrictFloat): Factors of e per cm, reduction in flux. N.B. This coefficient is positive (and greater than a lower bound for realism). - Raises: ValueError: if attenuation is below zero, or absorption is less than the lower bound set at the round magnitude 1.0e-8, diff --git a/src/dodal/common/general_maths/transmission_interconversion.py b/src/dodal/common/general_maths/transmission_interconversion.py index 8495176e5bd..a4338d19e7c 100644 --- a/src/dodal/common/general_maths/transmission_interconversion.py +++ b/src/dodal/common/general_maths/transmission_interconversion.py @@ -1,6 +1,7 @@ import math +from typing import Annotated -from pydantic import StrictFloat, validate_call +from pydantic import Field, StrictFloat, validate_call _CANONICAL_BARNETT_CONVERSION = -1.0e3 _REVERSE_BARNETT_CONVERSION = -1.0e-3 @@ -8,8 +9,7 @@ @validate_call def attenuation_from_natural_log_of_transmission(ln_t: StrictFloat) -> float: - """Converts from natural log of transmission fraction into Barnett attenuation units - . + """Converts from natural log of transmission fraction into Barnett attenuation units. Args: ln_t (StrictFloat): natural log of transmission fraction @@ -21,7 +21,9 @@ def attenuation_from_natural_log_of_transmission(ln_t: StrictFloat) -> float: @validate_call -def attenuation_from_transmission(transmission_as_fraction: StrictFloat) -> float: +def attenuation_from_transmission( + transmission_as_fraction: Annotated[StrictFloat, Field(ge=0, le=1)], +) -> float: """Converts from transmission fraction into Barnett attenuation units. Args: @@ -35,25 +37,28 @@ def attenuation_from_transmission(transmission_as_fraction: StrictFloat) -> floa @validate_call -def natural_log_of_transmission_from_attenuation(attenuation_bn: StrictFloat) -> float: - """Converts from Barnett attenuation units into natural log of transmission fraction - . +def natural_log_of_transmission_from_attenuation( + attenuation_bn: Annotated[StrictFloat, Field(ge=0)], +) -> float: + """Converts from Barnett attenuation units into natural log of transmission fraction. Args: - attenuation_bn (StrictFloat): Barnett attenuation units + attenuation_bn (StrictFloat): Barnett attenuation units - zero or a positive number Returns: - (float): natural log of transmission fraction + (float): natural log of transmission fraction - zero or, more usually, a negative number """ return _REVERSE_BARNETT_CONVERSION * attenuation_bn @validate_call -def transmission_from_attenutation(attenuation_bn: StrictFloat) -> float: +def transmission_from_attenutation( + attenuation_bn: Annotated[StrictFloat, Field(ge=0)], +) -> float: """Converts from Barnett attenuation units into transmission fraction. Args: - attenuation_bn (StrictFloat): Barnett attenuation units + attenuation_bn (StrictFloat): Barnett attenuation units - zero or a positive number Returns: (float): transmission fraction diff --git a/tests/common/general_maths/test_absorber_geometry.py b/tests/common/general_maths/test_absorber_geometry.py new file mode 100644 index 00000000000..95e31527a2f --- /dev/null +++ b/tests/common/general_maths/test_absorber_geometry.py @@ -0,0 +1,180 @@ +import math + +import pydantic +import pytest + +from dodal.common.general_maths.absorber_geometry import FoilGeometry, WedgeGeometry + +# happy path + + +@pytest.mark.parametrize( + "distance_unit, numerical_thickness_in_specified_unit, expected_thickness_cm", + [ + ("cm", 0.193, 0.193), + ("um", 150.0, 0.015), + ("mm", 4.2, 0.42), + ("micron", 812.5, 0.08125), + ], +) +def test_foil_absorber_reports_thickness_in_cm_when_validly_specified( + distance_unit, numerical_thickness_in_specified_unit, expected_thickness_cm +): + flat_absorber = FoilGeometry( + unit=distance_unit, numerical_value=numerical_thickness_in_specified_unit + ) + assert flat_absorber.get_thickness_cm() == pytest.approx(expected_thickness_cm) + + +@pytest.mark.parametrize( + "tip, taper_cotangent, probed_position_mm, expected_thickness_cm", + [ + (-5.0, 10.0, -5.0, 0.0), + (2.0, -10.0, -15.2, 0.172), + (-12.8, 9.7, 80.6, 0.9628866), + (12.8, 11.1, 22.4, 0.08648649), + ], +) +def test_wedge_geometry_reports_correct_thickness_in_cm( + tip, taper_cotangent, probed_position_mm, expected_thickness_cm +): + wedge_absorber = WedgeGeometry(tip_mm=tip, taper_cotangent=taper_cotangent) + assert wedge_absorber.thickness_cm_at_motor_position_mm( + probed_position_mm + ) == pytest.approx(expected_thickness_cm) + + +@pytest.mark.parametrize( + "tip, taper_cotangent, required_thickness_cm, expected_motor_position_mm", + [ + (-5.0, 10.0, 0.04, -1.0), + (2.0, -10.0, 0.172, -15.2), + (-12.8, 9.7, 0.9628866, 80.6), + (12.8, 11.1, 0.08648649, 22.4), + ], +) +def test_wedge_geometry_reports_correct_motor_position_mm_to_achieve_requested_thickness( + tip, taper_cotangent, required_thickness_cm, expected_motor_position_mm +): + wedge_absorber = WedgeGeometry(tip_mm=tip, taper_cotangent=taper_cotangent) + assert wedge_absorber.motor_position_mm_for_thickness_cm( + required_thickness_cm + ) == pytest.approx(expected_motor_position_mm) + + +# inauspicious path + + +@pytest.mark.parametrize( + "bad_input", + [ + "", + "k", + math.cos, + None, + [], + {}, + object(), + False, + ], +) +def test_foil_absorber_raises_error_when_given_invalid_input_for_numerical_thickness( + bad_input, +): + with pytest.raises(pydantic.ValidationError): + FoilGeometry(unit="cm", numerical_value=bad_input) + + +@pytest.mark.parametrize( + "unsupported_unit", + [ + "", + "mile", + "km", + "furlong", + "yard", + "ft", + "lb", + "s", + "Bq", + "A", + "eV", + ], +) +def test_foil_absorber_raises_error_when_asked_to_use_unsupported_distance_units( + unsupported_unit, +): + with pytest.raises(KeyError): + FoilGeometry(unit=unsupported_unit, numerical_value=1.1) + + +@pytest.mark.parametrize( + "bad_input", + [ + -9, + 4.7, + math.sin, + None, + [], + {}, + object(), + True, + ], +) +def test_foil_absorber_raises_error_when_given_invalid_distance_units(bad_input): + with pytest.raises(pydantic.ValidationError): + FoilGeometry(unit=bad_input, numerical_value=1.1) + + +@pytest.mark.parametrize( + "taper_cotangent", + [ + 0.0, + 1.0, + -2.1, + 0.03, + 4.58, + 1.16, + ], +) +def test_wedge_absorber_raises_error_when_given_chunky_wedge_angle(taper_cotangent): + with pytest.raises(pydantic.ValidationError): + WedgeGeometry(tip_mm=4.8, taper_cotangent=taper_cotangent) + + +@pytest.mark.parametrize( + "bad_input", + [ + "", + "k", + math.cos, + None, + [], + {}, + object(), + True, + ], +) +def test_wedge_absorber_raises_error_when_given_invalid_input_for_taper(bad_input): + with pytest.raises(pydantic.ValidationError): + WedgeGeometry(tip_mm=3.17, taper_cotangent=bad_input) + + +@pytest.mark.parametrize( + "bad_input", + [ + "", + "k", + math.cos, + None, + [], + {}, + object(), + True, + ], +) +def test_wedge_absorber_raises_error_when_given_invalid_input_for_tip_position( + bad_input, +): + with pytest.raises(pydantic.ValidationError): + WedgeGeometry(tip_mm=bad_input, taper_cotangent=-7.6) diff --git a/tests/common/general_maths/test_arithmetic_conversions.py b/tests/common/general_maths/test_arithmetic_conversions.py index d9ebaba1e95..60c03c4896e 100644 --- a/tests/common/general_maths/test_arithmetic_conversions.py +++ b/tests/common/general_maths/test_arithmetic_conversions.py @@ -13,6 +13,7 @@ convert_mm_to_cm, convert_mm_to_microns, convert_percentage_to_factor, + get_straight_line_y, ) from .operator_inversion_pairing import OperatorInversionPairing @@ -59,6 +60,29 @@ def test_conversion_from_microns_to_centimetres(input, result): assert convert_microns_to_cm(input) == pytest.approx(result) +@pytest.mark.parametrize( + "c,m,x,expected_y", + [ + (1.5, 1.0, 2.0, 3.5), + (-1.5, 1.0, 10.0, 8.5), + (3.5, -4.2, 7.2, -26.74), + (-0.8, -0.14, -9.1, 0.474), + (-2.4, 1.7, 2.8, 2.36), + (0.0, 15.8, 0.0, 0.0), + (0.0, -9.14, 0.0, 0.0), + (0.0, 3.2, 11.6, 37.12), + (0.0, -3.2, 11.6, -37.12), + (0.0, 3.2, -11.6, -37.12), + (0.1, -3.2, -11.6, 37.22), + (5.2, 0.0, -8.64, 5.2), + ], +) +def test_straight_line_conversion(c, m, x, expected_y): + assert get_straight_line_y(line_offset=c, line_gradient=m, x=x) == pytest.approx( + expected_y + ) + + # Circular "sanity check" tests, exercise pairs of reciprocating functions # proving the result of applying a function and its inverse results in the original value @@ -107,9 +131,11 @@ def test_reciprocal_function_pairs_nest_consistent_with_identity( # The inauspicuous path + + @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), False], ) def test_convert_microns_to_cm_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): @@ -118,7 +144,7 @@ def test_convert_microns_to_cm_raises_error_with_bad_input(bad_input): @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), False], ) def test_convert_ev_to_kev_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): @@ -127,7 +153,7 @@ def test_convert_ev_to_kev_raises_error_with_bad_input(bad_input): @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), False], ) def test_convert_microns_to_mm_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): @@ -136,7 +162,7 @@ def test_convert_microns_to_mm_raises_error_with_bad_input(bad_input): @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), False], ) def test_convert_mm_to_microns_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): @@ -145,7 +171,7 @@ def test_convert_mm_to_microns_raises_error_with_bad_input(bad_input): @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), True], ) def test_convert_factor_to_percentage_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): @@ -154,7 +180,7 @@ def test_convert_factor_to_percentage_raises_error_with_bad_input(bad_input): @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), False], ) def test_convert_percentage_to_factor_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): @@ -163,7 +189,7 @@ def test_convert_percentage_to_factor_raises_error_with_bad_input(bad_input): @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), False], ) def test_convert_mm_to_cm_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): @@ -172,8 +198,41 @@ def test_convert_mm_to_cm_raises_error_with_bad_input(bad_input): @pytest.mark.parametrize( "bad_input", - ["a", [], None, math.sin, object(), False], + ["", "a", [], None, math.sin, object(), True], ) def test_convert_cm_to_mm_raises_error_with_bad_input(bad_input): with pytest.raises(pydantic.ValidationError): convert_cm_to_mm(bad_input) + + +@pytest.mark.parametrize( + "bad_input", + ["", "a", [], None, math.log, object(), False], +) +def test_straight_line_calculator_raises_error_with_bad_offset(bad_input): + _probe_x = 11.1 + _line_gradient = -2.07 + with pytest.raises(pydantic.ValidationError): + get_straight_line_y(bad_input, _line_gradient, _probe_x) + + +@pytest.mark.parametrize( + "bad_input", + ["", "a", [], None, math.sin, object(), True], +) +def test_straight_line_calculator_raises_error_with_bad_gradient(bad_input): + _probe_x = -11.1 + _line_offset = 14.2 + with pytest.raises(pydantic.ValidationError): + get_straight_line_y(_line_offset, bad_input, _probe_x) + + +@pytest.mark.parametrize( + "bad_input", + ["", "a", [], None, math.tan, object(), False], +) +def test_straight_line_calculator_raises_error_with_bad_x_value(bad_input): + _line_offset = 0.1 + _line_gradient = 5.4 + with pytest.raises(pydantic.ValidationError): + get_straight_line_y(_line_offset, _line_gradient, bad_input) diff --git a/tests/common/general_maths/test_material_absorption_maths.py b/tests/common/general_maths/test_material_absorption_maths.py index 9d34d047a3b..400afcced8a 100644 --- a/tests/common/general_maths/test_material_absorption_maths.py +++ b/tests/common/general_maths/test_material_absorption_maths.py @@ -1,24 +1,101 @@ import math +from typing import cast +from unittest.mock import MagicMock, patch import pytest from pydantic import ValidationError from dodal.common.general_maths.material_absorption_maths import ( + AbsorptionCalculator, + CompoundAbsorptionCalculator, + PolynomialAbsorptionCorrection, + SingleRollOffAbsorptionCalculator, attenuation_at_depth_cm, photon_mass_attenuation_per_unit_length, thickness_cm_required_to_attenuate, ) - # happy path + + +# N.B. the mathematical validity of results from a simple absorption calculator +# are not tested here - since that would duplicate the validatity tests of numerical calculations +# already checked against the photon_mass_attenuation function +@patch( + "dodal.common.general_maths.material_absorption_maths.photon_mass_attenuation_per_unit_length" +) +def test_simple_absorption_calculator_invokes_photon_mass_attenuation_calculator( + mock_pma_calc: MagicMock, +): + _material_factor = 3.1415926 + _roll_off = -2.71828 + calculator = SingleRollOffAbsorptionCalculator(_material_factor, _roll_off) + _x_ray_kev = 15.6302 + calculator.absorption_coefficient_per_cm(_x_ray_kev) + mock_pma_calc.assert_called_once_with(_x_ray_kev, _material_factor, _roll_off) + + +# N.B. the mathematical validity of results from a compound absorption calculator +# are not tested here - since that would duplicate the validatity tests of numerical calculations +# already checked against the contributing calculators within +def test_compound_absorption_calculator_invokes_contributing_calculators(): + mock_calcs: list[MagicMock] = [ + MagicMock(spec=AbsorptionCalculator), + MagicMock(spec=AbsorptionCalculator), + ] + calcs = cast(list[AbsorptionCalculator], mock_calcs) + compound_calculator = CompoundAbsorptionCalculator(calcs) + _x_ray_kev = 12.6039 + compound_calculator.absorption_coefficient_per_cm(_x_ray_kev) + for mocked_calculator in mock_calcs: + mocked_calculator.absorption_coefficient_per_cm.assert_called_once_with( + _x_ray_kev + ) + + +@pytest.mark.parametrize( + "energy_kev, polynomial_coefficients, expected_absorption_coefficient", + [ + ( + 8.3328, + [50.9211, -23.6148, 4.2138, -0.3814, 1.867e-2, -4.709e-4, 4.796e-6], + -1.24317951, + ), + ( + 12.3984, + [43.4, -22.8148, 4.138, -0.3714, 1.967e-2, -5.709e-4, 3.796e-6], + 0.11254471963, + ), + ( + 11.9187, + [50.9211, -23.6148, 4.2138, -0.3814, 1.867e-2, -4.709e-4, 4.796e-6], + -0.45286288, + ), + ], +) +def test_corrective_polynomial_absorption_calculator( + energy_kev, polynomial_coefficients, expected_absorption_coefficient +): + corrective_calculator = PolynomialAbsorptionCorrection(polynomial_coefficients) + assert corrective_calculator.absorption_coefficient_per_cm( + energy_kev + ) == pytest.approx(expected_absorption_coefficient) + + @pytest.mark.parametrize( "energy_kev, photon_absorption_factor_per_unit_length, energy_dependence_exponent, " "result", [ - (5.042, 1.98e2, -2.717, 2.44170544), # Arbitrary Energy - (8.3328, 2.5706e3, -2.83, 6.3708311), # Nickel - (11.9187, 1.48e3, -2.93, 1.03970725), # Gold-Three - (25.514, 6.48e3, -2.41, 2.63778077), # Silver + (5.042, 1.98e2, -2.717, 2.44170544), # At an arbitrary energy 5.042 keV + (8.3328, 2.5706e3, -2.83, 6.3708311), # Nickel energy 8.3328 keV + (11.9187, 1.48e3, -2.93, 1.03970725), # Gold-Three energy 11.9187 keV + (25.514, 6.48e3, -2.41, 2.63778077), # Silver energy 25.514 keV + ( + 20, + 1201, + -2, + 3.0025, + ), # 20 keV: Test that integers, for input args, are not rejected ], ) def test_photon_mass_attenuation_per_unit_length( @@ -84,6 +161,8 @@ def test_attenuation_at_depth_cm(depth_cm, absorption_coefficient_per_cm, result # inauspicious path + + @pytest.mark.parametrize("negative_energy", [-0.1, -12.3, -9739.43]) def test_photon_mass_attenuation_per_unit_length_errors_with_negative_energy( negative_energy, @@ -127,7 +206,7 @@ def test_photon_mass_attenuation_per_unit_length_errors_with_invalid_exponent( invalid_roll_off, ): _energy_kev = 13.819 - _absorption_coefficient = 24.8 + _absorption_coefficient = 2.148 with pytest.raises(ValidationError): photon_mass_attenuation_per_unit_length( _energy_kev, _absorption_coefficient, invalid_roll_off @@ -194,7 +273,7 @@ def test_attenuation_at_depth_raises_error_with_invalid_depth(invalid_depth): @pytest.mark.parametrize( - "invalid_absorption", ["a", [], None, math.sin, object(), False] + "invalid_absorption", ["a", [], None, math.tan, object(), False] ) def test_attenuation_at_depth_raises_error_with_invalid_attenuation(invalid_absorption): _depth_cm = 0.425 diff --git a/tests/common/general_maths/test_transmission_interconversion.py b/tests/common/general_maths/test_transmission_interconversion.py index e891f60757b..0d5b937a19a 100644 --- a/tests/common/general_maths/test_transmission_interconversion.py +++ b/tests/common/general_maths/test_transmission_interconversion.py @@ -1,4 +1,5 @@ import math +from collections.abc import Callable import pytest from pydantic import ValidationError @@ -10,6 +11,8 @@ transmission_from_attenutation, ) +from .operator_inversion_pairing import OperatorInversionPairing + @pytest.mark.parametrize( "attenuation_bn,result", @@ -80,7 +83,7 @@ def test_attenuation_from_transmission(transmission_as_fraction, result): @pytest.mark.parametrize( - "ln_t,result", + "ln_t, result", [ (-1, 1000), # tests negative unity log of transmission is 1000 (canonical) (0, 0), # tests natural log of transparency is zero (canonical) @@ -101,24 +104,50 @@ def test_attenuation_from_natural_log_of_transmission(ln_t, result): # Circular tests (all numbers here arbitrary) -@pytest.mark.parametrize("input", [1.0, 10.0, 100.0]) -def test_circular_attenuation_from_log_and_back(input): - assert attenuation_from_natural_log_of_transmission( - natural_log_of_transmission_from_attenuation(input) - ) == pytest.approx(input) - assert natural_log_of_transmission_from_attenuation( - attenuation_from_natural_log_of_transmission(input) - ) == pytest.approx(input) - - -@pytest.mark.parametrize("input", [1.0, 10.0, 100.0]) -def test_circular_attenuation_from_transmission_and_back(input): - assert attenuation_from_transmission( - transmission_from_attenutation(input) - ) == pytest.approx(input) - assert transmission_from_attenutation( - attenuation_from_transmission(input) - ) == pytest.approx(input) +# with different ranges of valid input - the argument values must differ for each direction +@pytest.mark.parametrize( + "f, g, numerical_args", + [ + ( + attenuation_from_natural_log_of_transmission, + natural_log_of_transmission_from_attenuation, + [-0.075, -1.2, -6.3], + ), + ( + natural_log_of_transmission_from_attenuation, + attenuation_from_natural_log_of_transmission, + [0.04, 0.91, 2.02, 5.7], + ), + ], +) +def test_attenuation_log_transmssion_conversions_reciprocate_as_expected( + f: Callable[[float], float], + g: Callable[[float], float], + numerical_args: list[float], +): + op_pair = OperatorInversionPairing(f, g) + for x in numerical_args: + assert op_pair.composed_operator_is_consistent_with_identity_operator(x) + + +@pytest.mark.parametrize( + "f, g, numerical_args", + [ + ( + attenuation_from_transmission, + transmission_from_attenutation, + [0.92, 0.68, 0.147, 0.00013], + ) + ], +) +def test_transmission_attenuation_conversions_pair_off_to_form_identity_operation( + f: Callable[[float], float], + g: Callable[[float], float], + numerical_args: list[float], +): + for op_pair in [OperatorInversionPairing(f, g), OperatorInversionPairing(g, f)]: + for x in numerical_args: + assert op_pair.composed_operator_is_consistent_with_identity_operator(x) # inauspicious: