diff --git a/docs/tutorials/notebooks/Tutorial_extending_fluopy.ipynb b/docs/tutorials/notebooks/Tutorial_extending_fluopy.ipynb index ec2598f..115a5c4 100644 --- a/docs/tutorials/notebooks/Tutorial_extending_fluopy.ipynb +++ b/docs/tutorials/notebooks/Tutorial_extending_fluopy.ipynb @@ -16,6 +16,20 @@ "Here we provide some hints on extending fluoropy with new fluorophore and transition data." ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "788a279a", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import numpy as np\n", + "\n", + "import fluopy" + ] + }, { "cell_type": "markdown", "id": "69c0ebd0-b041-4676-b14b-a51a47334bef", @@ -46,10 +60,75 @@ "metadata": {}, "source": [ "## Adding a fluorophore\n", - "1. Create a new instance of FluorophoreData in fluodata.py\n", - "2. Overwrite the known constants**\n", - "3. Create a Fluorophore instance with a name matching the variable in fluodata.py\n", - "4. Provide datafiles containing absorption and emission spectra* in a folder in fluorophore_spectra and keep the naming convention (e.g., absorption of S0 should be named absorption_S0.csv). The name of the folder should be provided as a str in fluo_data.py (data_files=str)." + "\n", + "A fluorophore that is not included with Fluopy can be defined in a notebook.\n", + "\n", + "1. Create an S0 absorption Spectrum from arrays or a CSV file. Also create an\n", + " emission Spectrum when bandpass filtering or energy-transfer calculations are\n", + " required.\n", + "2. Create a FluorophoreData object containing the spectra and photophysical constants.\n", + "3. Pass the FluorophoreData object to Fluorophore\n", + "4. Use FluorophoreSystem.load_transitions() to derive transitions automatically." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c44c57d3", + "metadata": {}, + "outputs": [], + "source": [ + "wavelengths = np.array([600, 620, 640, 660, 680])\n", + "\n", + "emission = fluopy.Spectrum.from_arrays(\n", + " wavelengths=wavelengths,\n", + " values=[0.0, 0.2, 1.0, 0.6, 0.1],\n", + ")\n", + "absorption_s0 = fluopy.Spectrum.from_arrays(\n", + " wavelengths=wavelengths,\n", + " values=[10000, 40000, 80000, 30000, 5000],\n", + ")\n", + "\n", + "custom_data = fluopy.FluorophoreData(\n", + " QUANTUM_YIELD=0.6,\n", + " FLUORESCENCE_LIFETIME=3e-9,\n", + " emission_spectrum=emission,\n", + " absorption_spectra={\"s0\": absorption_s0},\n", + ")\n", + "\n", + "custom_fluorophore = fluopy.Fluorophore(\n", + " name=\"custom\",\n", + " position=[0, 0],\n", + " constants=custom_data,\n", + ")\n", + "\n", + "fluorophore_system = fluopy.FluorophoreSystem(\n", + " fluorophores=[custom_fluorophore],\n", + ")\n", + "transitions = fluorophore_system.load_transitions(\n", + " wavelength=640,\n", + " energy_transfer=False,\n", + " dstorm=False,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "f42ff35f", + "metadata": {}, + "outputs": [], + "source": [ + "spectrum_dir = Path(fluopy.__file__).parent / \"fluorophore_spectra\" / \"atto643_data\"\n", + "\n", + "custom_data = fluopy.FluorophoreData(\n", + " QUANTUM_YIELD=0.6,\n", + " FLUORESCENCE_LIFETIME=3e-9,\n", + " emission_spectrum=fluopy.Spectrum.from_csv(spectrum_dir / \"emission.csv\"),\n", + " absorption_spectra={\n", + " \"s0\": fluopy.Spectrum.from_csv(spectrum_dir / \"absorption_s0.csv\"),\n", + " },\n", + ")" ] }, { @@ -92,7 +171,15 @@ }, "source": [ "### Spectra\n", - "The absorption spectra should contain absolute extinction coefficient values for each wavelength (nm) [200, 201, ..., 1000]. The emission spectrum should contain relative values where 1 corresponds to maximum emission, for each wavelength (nm) [200, 201, ..., 1000]." + "A spectrum consists of one-dimensional wavelength and value arrays. Wavelengths\n", + "are given in nm and must be strictly increasing. Spectrum values must be non-negative.\n", + "\n", + "Absorption spectra contain absolute molar extinction coefficients. Emission\n", + "spectra may contain relative intensities because they are normalized where\n", + "required.\n", + "\n", + "CSV files use 'Wavelengths' and 'y' as the default column names. Alternative\n", + "column names can be passed to Spectrum.from_csv()." ] }, { @@ -107,7 +194,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "fluopy (3.12.10)", "language": "python", "name": "python3" }, @@ -121,7 +208,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.4" + "version": "3.12.10" } }, "nbformat": 4, diff --git a/src/fluopy/emissions.py b/src/fluopy/emissions.py index 74fc10d..42f5ca4 100644 --- a/src/fluopy/emissions.py +++ b/src/fluopy/emissions.py @@ -14,7 +14,7 @@ from scipy.stats import binom, gamma, norm, poisson from . import figure as fi -from .fluorophores import Fluorophore +from .fluo_data import Spectrum from .simulation import ( Simulation, eval_floating_point_precision_error, @@ -67,8 +67,8 @@ def __init__( For possible input values, see https://pandas.pydata.org/docs/user_guide/timeseries.html -> Offset aliases. bandpass - The lowest and highest emission wavelength to be passed by the bandpass - filter. Requires emission spectrum data when specified. + The lowest and highest wavelength in nm passed by the bandpass filter. + Requires emission spectrum data when specified. seed A seed to initialize the BitGenerator. """ @@ -330,19 +330,20 @@ def get_emission_indices( rng = np.random.default_rng(seed) processed = [] collect_emission_indices = [] - data_dir = Path(__file__).parent / "fluorophore_spectra" for ( fluorophore ) in simulation.transition_set.fluorophore_system.fluorophores: - if fluorophore.constants is None: + constants = fluorophore.constants + if constants is None or constants.emission_spectrum is None: raise ValueError( "bandpass not None but emission data not available for " f"this kind of fluorophore: {fluorophore.name}" ) if fluorophore.name not in processed: p_passed = get_p_filter( - data_dir=data_dir, fluorophore=fluorophore, bandpass=bandpass + emission_spectrum=constants.emission_spectrum, + bandpass=bandpass, ) p_not_passed = 1 - p_passed sub_df = simulation.transition_set.transition_df.loc[ @@ -755,19 +756,16 @@ def load(cls, path: str | Path, name_extension: str = "") -> Emissions: def get_p_filter( - data_dir: str | Path, - fluorophore: Fluorophore, + emission_spectrum: Spectrum, bandpass: tuple[float, float], ) -> float: """ - Get the probability of a photon emitted by fluorophore passing the bandpass filter. + Get the fraction of an emission spectrum passing the bandpass filter. Parameters ---------- - data_dir - The directory of data files of fluorophores. - fluorophore - Contains attributes of a fluorophore. + emission_spectrum + Emission spectrum to which the bandpass filter is applied. bandpass The lowest and highest emission wavelength to be passed by the bandpass filter. @@ -777,26 +775,24 @@ def get_p_filter( The probability of a photon passing the bandpass filter. """ - if bandpass[0] < 200 or bandpass[0] > 1000: - raise ValueError("The lower bandpass limit has to be between 200 and 1000 nm.") - if bandpass[1] < 200 or bandpass[1] > 1000: - raise ValueError("The upper bandpass limit has to be between 200 and 1000 nm.") - if bandpass[0] >= bandpass[1]: + lower, upper = bandpass + + if not np.isfinite(lower) or not np.isfinite(upper): + raise ValueError("bandpass limits must be finite.") + if lower >= upper: raise ValueError( "The lower bandpass limit has to be smaller than the upper limit." ) - emission_data = pd.read_csv( - Path(data_dir) / fluorophore.constants.data_files / "emission.csv" - ) + total_emission = emission_spectrum.integral() + if total_emission == 0: + raise ValueError("emission spectrum has zero total intensity.") - minimum_wavelength = 200 - - emissions = emission_data["y"] - bandpass_low = bandpass[0] - minimum_wavelength - bandpass_high = bandpass_low + (bandpass[1] - bandpass[0]) - rel_emission = emissions[bandpass_low:bandpass_high] / emissions.sum() - p_passed = rel_emission.sum() + passed_emission = emission_spectrum.integral( + lower=lower, + upper=upper, + ) + p_passed = passed_emission / total_emission return p_passed @@ -828,17 +824,16 @@ def get_emitting_transition_ids( emitting_transition_ids = {} if bandpass is not None: processed = [] - data_dir = Path(__file__).parent / "fluorophore_spectra" for fluorophore in transition_set.fluorophore_system.fluorophores: - if fluorophore.constants is None: + constants = fluorophore.constants + if constants is None or constants.emission_spectrum is None: raise ValueError( "bandpass not None but emission data not available for " f"this kind of fluorophore: {fluorophore.name}" ) if fluorophore.name not in processed: p_passed = get_p_filter( - data_dir=data_dir, - fluorophore=fluorophore, + emission_spectrum=constants.emission_spectrum, bandpass=bandpass, ) sub_df = transition_set.transition_df.loc[fluorophore.name] @@ -849,6 +844,7 @@ def get_emitting_transition_ids( ].index.to_numpy() for emitting_transition_id in emitting_transition_ids_f: emitting_transition_ids[emitting_transition_id] = p_passed + processed.append(fluorophore.name) else: df = transition_set.combined_state_transitions_df emitting_transition_ids_ = df.loc[df["photon"]].index.to_numpy() diff --git a/src/fluopy/fluo_data.py b/src/fluopy/fluo_data.py index dd93762..80a7696 100644 --- a/src/fluopy/fluo_data.py +++ b/src/fluopy/fluo_data.py @@ -4,10 +4,222 @@ This module provides a dataclass container to hold photophysical constants. """ -from dataclasses import dataclass +from dataclasses import dataclass, field +from os import PathLike from pathlib import Path +from typing import Self -__all__: list[str] = ["FluorophoreData", "cy5_dna", "atto643"] +import numpy as np +import numpy.typing as npt +import pandas as pd + +__all__: list[str] = ["Spectrum", "FluorophoreData", "cy5_dna", "atto643"] + + +@dataclass +class Spectrum: + """ + Contains wavelength-dependent spectral data. + + Spectrum values can be provided directly as array-like objects or loaded from a csv + file. Wavelengths are given in nm and must be strictly increasing. + + Attributes + ---------- + wavelengths : 1-D array_like + The wavelength values in nm. + values : 1-D array_like + Spectrum values corresponding to wavelengths. + """ + + wavelengths: npt.ArrayLike + values: npt.ArrayLike + + @classmethod + def from_arrays( + cls, + wavelengths: npt.ArrayLike, + values: npt.ArrayLike, + ) -> Self: + """ + Create a spectrum from wavelength and value arrays. + + Parameters + ---------- + wavlengths + Wavelengths in nm. + values + Spectrum values corresponding to wavelenghts. + + Returns + ------- + Spectrum + Spectrum object containing copies of the input arrays. + """ + return cls(wavelengths=wavelengths, values=values) + + @classmethod + def from_csv( + cls, + path: str | PathLike[str], + wavelength_column: str = "Wavelengths", + value_column: str = "y", + ) -> Self: + """ + Create a spectrum from a CSV file. + + Parameters + ---------- + path + Path to the CSV file. + wavelength_column + Name of the column containing wavelengths in nm. + value_column + Name of the column containing spetrum values. + + Returns + ------- + Spectrum + Spectrum object loaded from the CSV file. + """ + data = pd.read_csv(path) + + missing_columns = { + wavelength_column, + value_column, + }.difference(data.columns) + + if missing_columns: + missing = ", ".join(sorted(missing_columns)) + raise ValueError(f"spectrum CSV is missing columns: {missing}.") + + return cls( + wavelengths=data[wavelength_column].to_numpy(), + values=data[value_column].to_numpy(), + ) + + def at(self, wavelength: float) -> float: + """ + Return the spectrum value at a wavelength. + + Values between given wavelengths are linearily interpolated. Extrapolation + outside the spectrum range is not supported. + + Parameters + ---------- + wavelength + Wavelength in nm. + + Returns + ------- + value + Spectrum value at the specified wavelength. + """ + if not np.isfinite(wavelength): + raise ValueError("wavelength must be finite.") + + minimum = self.wavelengths[0] + maximum = self.wavelengths[-1] + + if wavelength < minimum or wavelength > maximum: + raise ValueError( + f"wavelength {wavelength} nm is outside the spectrum range " + f"{minimum}-{maximum} nm." + ) + + value = float(np.interp(wavelength, self.wavelengths, self.values)) + + return value + + def integral( + self, + lower: float | None = None, + upper: float | None = None, + ) -> float: + """ + Integrate the spectrum over a wavelength interval. + + Integration limits outside the available spectrum are clipped to the spectrum + range. An interval without overlap has an integral of zero. + + Parameters + ---------- + lower + Lower integration limit in nm. If None, use the lowest available wavelength. + upper + Upper integration limit in nm. If None, use the highest available + wavelength. + + Returns + ------- + integral + Trapezoidal integral of the spectrum. + """ + if lower is None: + lower = float(self.wavelengths[0]) + if upper is None: + upper = float(self.wavelengths[-1]) + + if not np.isfinite(lower) or not np.isfinite(upper): + raise ValueError("integration limits must be finite.") + if lower >= upper: + raise ValueError( + "the lower integration limit must be smaller than the upper limit." + ) + + lower = max(lower, float(self.wavelengths[0])) + upper = min(upper, float(self.wavelengths[-1])) + + if lower >= upper: + return 0.0 + + inside = (self.wavelengths > lower) & (self.wavelengths < upper) + wavelengths = np.concatenate(([lower], self.wavelengths[inside], [upper])) + values = np.interp( + wavelengths, + self.wavelengths, + self.values, + ) + integral = float(np.trapezoid(values, wavelengths)) + + return integral + + def __post_init__(self) -> None: + self.wavelengths = np.asarray(self.wavelengths, dtype=float).copy() + self.values = np.asarray(self.values, dtype=float).copy() + + if self.wavelengths.ndim != 1: + raise ValueError("spectrum wavelengths must be one-dimensional.") + if self.values.ndim != 1: + raise ValueError("spectrum values must be one-dimensional.") + if self.wavelengths.size != self.values.size: + raise ValueError( + "spectrum wavelengths and values must have the same length." + ) + if self.wavelengths.size < 2: + raise ValueError("a spectrum must contain at least two data points.") + if not np.all(np.isfinite(self.wavelengths)): + raise ValueError("spectrum wavelengths must be finite.") + if not np.all(np.isfinite(self.values)): + raise ValueError("spectrum values must be finite.") + if not np.all(np.diff(self.wavelengths) > 0): + raise ValueError("spectrum wavelengths must be strictly increasing.") + if np.any(self.values < 0): + raise ValueError("spectrum values must be non-negative.") + + +def _load_bundled_spectra( + directory_name: str, +) -> tuple[Spectrum, dict[str, Spectrum]]: + directory = Path(__file__).parent / "fluorophore_spectra" / directory_name + + emission_spectrum = Spectrum.from_csv(directory / "emission.csv") + absorption_spectra = { + path.stem.removeprefix("absorption_"): Spectrum.from_csv(path) + for path in sorted(directory.glob("absorption_*.csv")) + } + + return emission_spectrum, absorption_spectra @dataclass @@ -18,11 +230,13 @@ class FluorophoreData: Attributes ---------- - data_files : str | Path | None - The name of the folder containing the spectra data files. The folder should be - located in src/fluopy/fluorophore_spectra. Needed to infer excitation rate - and energy transfer rates. If None, no automatic inference of rates will be - performed. + emission_spectrum : Spectrum | None + Emission spectrum used for bandpass filtering and as the donor spectrum in + energy-transfer calculations. + absorption_spectra : dict[str, Spectrum] + Absorption spectra indexed by lowercase acceptor-state names, for example + 's0', 't1', 'cis' or 'off'. The S0 spectrum is also used to infer the excitation + rate. QUANTUM_YIELD : float The fluorescence quantum yield of the fluorophore. Should be between 0 and 1. FLUORESCENCE_LIFETIME : float @@ -41,12 +255,12 @@ class FluorophoreData: PHOTOBLEACH_T1_RATE : float The photobleaching rate from T1 to B in 1/s. CROSS_SECTION_WAVELENGTH : int | None - The wavelength in nm at which absorption cross sections are defined. The - standard excitation of S0 is handled via data_files (entire absorption - spectrum), but for other transitions (e.g., cis absorption to define - photoinduced back-isomerization), a single cross section should be provided. - The cross_section_wavelength is used to check whether the provided cross - sections are given for the same wavelength as a specified wavelength. + The wavelength in nm at which individual absorption cross sections are defined. + Standard excitation from S0 is calculated using the S0 absorption spectrum in + absorption_spectra. For other transitions, such as photoinduced + back-isomerization from cis, an individual cross section can be provided. + CROSS_SECTION_WAVELENGTH is used to check whether these cross sections + correspond to the specified excitation wavelength. DSTORM_PET_T_RATE_MOL : float The concentration-dependent PET rate that targets T1 in 1/(M*s). DSTORM_PET_S_RATE_MOL : float @@ -84,7 +298,8 @@ class FluorophoreData: """ # spectra - data_files: str | Path | None = None + emission_spectrum: Spectrum | None = None + absorption_spectra: dict[str, Spectrum] = field(default_factory=dict) # general QUANTUM_YIELD: float = 0 @@ -118,9 +333,30 @@ class FluorophoreData: H2O_ATTACK_T: float = 0 BACK_REACTION: float = 0 + def __post_init__(self) -> None: + if self.emission_spectrum is not None and not isinstance( + self.emission_spectrum, Spectrum + ): + raise TypeError("emission_spectrum must be a Spectrum or None.") + + for state, spectrum in self.absorption_spectra.items(): + if not isinstance(state, str) or not state: + raise TypeError("absorption_spectra keys must be non-empty strings.") + if not isinstance(spectrum, Spectrum): + raise TypeError( + f"absorption spectrum for state {state!r} must be a Spectrum." + ) + + +_cy5_emission, _cy5_absorption = _load_bundled_spectra("cy5_data") +_atto643_emission, _atto643_absorption = _load_bundled_spectra("atto643_data") +_testfluo_1_emission, _testfluo_1_absorption = _load_bundled_spectra("testing_data_1") +_testfluo_2_emission, _testfluo_2_absorption = _load_bundled_spectra("testing_data_2") + cy5_dna = FluorophoreData( - data_files="cy5_data", + emission_spectrum=_cy5_emission, + absorption_spectra=_cy5_absorption, QUANTUM_YIELD=0.27, FLUORESCENCE_LIFETIME=1.7e-9, ISC_ST_RATE=8.3e5, @@ -149,7 +385,8 @@ class FluorophoreData: atto643 = FluorophoreData( - data_files="atto643_data", + emission_spectrum=_atto643_emission, + absorption_spectra=_atto643_absorption, QUANTUM_YIELD=0.6, FLUORESCENCE_LIFETIME=3e-9, S1_QUENCH_RATE=0, # to be updated @@ -165,7 +402,8 @@ class FluorophoreData: testfluo_1 = FluorophoreData( - data_files="testing_data_1", + emission_spectrum=_testfluo_1_emission, + absorption_spectra=_testfluo_1_absorption, QUANTUM_YIELD=0.27, FLUORESCENCE_LIFETIME=1e-9, ISC_ST_RATE=8.3e5, @@ -184,7 +422,8 @@ class FluorophoreData: testfluo_2 = FluorophoreData( - data_files="testing_data_2", + emission_spectrum=_testfluo_2_emission, + absorption_spectra=_testfluo_2_absorption, QUANTUM_YIELD=0.6, FLUORESCENCE_LIFETIME=3e-9, S1_QUENCH_RATE=0, diff --git a/src/fluopy/formulas.py b/src/fluopy/formulas.py index a8cb199..a3e4f42 100644 --- a/src/fluopy/formulas.py +++ b/src/fluopy/formulas.py @@ -311,9 +311,13 @@ def calculate_spectral_overlap_integral( if donor.size != acceptor.size or donor.size != wavelengths.size: raise ValueError("donor, acceptor and wavelengths have to be of the same size.") - donor = donor / np.trapezoid(donor) # normalize spectrum to area of 1 - not_integrated = donor * acceptor * wavelengths**4 - spectral_overlap_integral = np.trapezoid(not_integrated) + donor_area = np.trapezoid(donor, x=wavelengths) + if donor_area <= 0: + raise ValueError("donor emission spectrum must have positive area.") + + normalized_donor = donor / donor_area + integrand = normalized_donor * acceptor * wavelengths**4 + spectral_overlap_integral = np.trapezoid(integrand, x=wavelengths) return spectral_overlap_integral diff --git a/src/fluopy/miscellaneous.py b/src/fluopy/miscellaneous.py index dc35b47..c0dea4f 100644 --- a/src/fluopy/miscellaneous.py +++ b/src/fluopy/miscellaneous.py @@ -14,6 +14,7 @@ import numpy as np import numpy.typing as npt import pandas as pd +from PIL import Image, ImageOps if TYPE_CHECKING: from matplotlib.axes import Axes as mplAxes @@ -318,14 +319,11 @@ def compute_tight_bbox(fig, pad_inches: float = 0.0): def crop_to_content_with_padding( - in_file, out_file, - dpi: int = 300, - pad_inches: float = 2/72, - threshold: int = 255 + in_file, out_file, dpi: int = 300, pad_inches: float = 2 / 72, threshold: int = 255 ) -> None: """ Crops the image to the content and adds padding, then saves the image. - + Parameters ---------- in_file diff --git a/src/fluopy/transitions.py b/src/fluopy/transitions.py index d28a46d..6c0749a 100644 --- a/src/fluopy/transitions.py +++ b/src/fluopy/transitions.py @@ -11,7 +11,6 @@ from dataclasses import asdict, dataclass, field from enum import Enum from itertools import product -from pathlib import Path from typing import TYPE_CHECKING, Self import numpy as np @@ -21,11 +20,11 @@ from . import formulas as fo from . import network as net +from .fluo_data import FluorophoreData if TYPE_CHECKING: from matplotlib.axes import Axes as mplAxes - from fluopy.fluo_data import FluorophoreData from fluopy.fluorophores import Fluorophore, FluorophoreSystem @@ -1042,7 +1041,7 @@ def derive_energy_transfer_transitions( """ Derive energy transfer transitions based on the experimental conditions and the fluorophore-combinations to be mimicked. The type of energy transfer is determined - via the data file names. + by the state names in acceptor_data.absorption_spectra. Parameters ---------- @@ -1075,21 +1074,24 @@ def derive_energy_transfer_transitions( transitions : list[Transition] Contains energy transfer transitions of type Transition. """ - data_dir = Path(__file__).parent / "fluorophore_spectra" - donor_emission = pd.read_csv(data_dir / donor_data.data_files / "emission.csv") - acceptor_files = sorted(Path(data_dir / acceptor_data.data_files).iterdir()) - acceptor_abs_files = [ - data_file - for data_file in acceptor_files - if data_file.name.startswith("absorption") - ] + donor_emission = donor_data.emission_spectrum + if donor_emission is None: + raise ValueError( + "cannot derive energy-transfer transitions without a donor " + "emission spectrum." + ) + + acceptor_absorptions = acceptor_data.absorption_spectra + if not acceptor_absorptions: + raise ValueError( + "cannot derive energy-transfer transitions without acceptor " + "absorption spectra." + ) emission_rate = fo.calculate_emission_rate( quantum_yield=donor_data.QUANTUM_YIELD, fluorescence_lifetime=donor_data.FLUORESCENCE_LIFETIME, ) - minimum, maximum = 200, 1000 - wavelengths_of_interest = np.arange(minimum, maximum + 1, 1, dtype=float) which_et = { "s0": [(TransitionType.FRET, 1)], @@ -1164,24 +1166,69 @@ def derive_energy_transfer_transitions( ) transitions = [] - for acceptor_abs_file in acceptor_abs_files: - acceptor_abs = pd.read_csv( - Path(data_dir) / acceptor_data.data_files / acceptor_abs_file + for acceptor_state, acceptor_absorption in sorted(acceptor_absorptions.items()): + minimum = max( + donor_emission.wavelengths[0], + acceptor_absorption.wavelengths[0], ) - - J = fo.calculate_spectral_overlap_integral( - donor=donor_emission["y"], - acceptor=acceptor_abs["y"], - wavelengths=wavelengths_of_interest, + maximum = min( + donor_emission.wavelengths[-1], + acceptor_absorption.wavelengths[-1], ) + + if minimum >= maximum: + spectral_overlap_integral = 0.0 + else: + donor_inside = donor_emission.wavelengths[ + (donor_emission.wavelengths > minimum) + & (donor_emission.wavelengths < maximum) + ] + acceptor_inside = acceptor_absorption.wavelengths[ + (acceptor_absorption.wavelengths > minimum) + & (acceptor_absorption.wavelengths < maximum) + ] + wavelengths = np.unique( + np.concatenate( + ( + [minimum], + donor_inside, + acceptor_inside, + [maximum], + ) + ) + ) + + donor_values = np.interp( + wavelengths, + donor_emission.wavelengths, + donor_emission.values, + ) + acceptor_values = np.interp( + wavelengths, + acceptor_absorption.wavelengths, + acceptor_absorption.values, + ) + + spectral_overlap_integral = fo.calculate_spectral_overlap_integral( + donor=donor_values, + acceptor=acceptor_values, + wavelengths=wavelengths, + ) + rate = fo.calculate_fret_rate( distance=distance, emission_rate=emission_rate, - spectral_overlap_integral=J, + spectral_overlap_integral=spectral_overlap_integral, dipole_orientation_factor=dipole_orientation_factor, refractive_index=refractive_index, ) - acceptor_state = acceptor_abs_file.name.split("_")[1].split(".")[0] + + if acceptor_state not in which_et_new: + raise ValueError( + f"energy transfer to acceptor state {acceptor_state!r} " + "is not supported." + ) + if exclude is not None and acceptor_state in exclude: continue for transition_type, factor in which_et_new[acceptor_state]: @@ -1246,18 +1293,18 @@ def derive_transitions( ) _, _, frequency = fo.convert_wavenumber_wavelength_frequency(wavelength=wavelength) photon_flux = fo.calculate_photon_flux(irradiance=irradiance, frequency=frequency) - path_absorption = ( - Path(__file__).parent - / "fluorophore_spectra" - / fd.data_files - / "absorption_s0.csv" - ) + + if "s0" not in fd.absorption_spectra: + raise ValueError( + "cannot derive excitation transition without an S0 absorption spectrum." + ) + + absorption_spectrum = fd.absorption_spectra["s0"] + if fluorophore_ids is None: fluorophore_ids = [0] - dataframe_absorption = pd.read_csv(filepath_or_buffer=path_absorption, index_col=0) - - extinction_coefficient = dataframe_absorption.loc[int(wavelength), "y"] + extinction_coefficient = absorption_spectrum.at(wavelength) excitation_rate = fo.calculate_excitation_rate( photon_flux=photon_flux, extinction_coefficient=extinction_coefficient diff --git a/tests/test_emissions.py b/tests/test_emissions.py index bbcb5a6..6596d6c 100644 --- a/tests/test_emissions.py +++ b/tests/test_emissions.py @@ -7,51 +7,95 @@ import pytest from fluopy import emissions as em -from fluopy import fluorophores as fl +from fluopy import fluo_data as fd @pytest.mark.parametrize( "bandpass, expected", [ - [(650, 700), 0.685702066268812], - [(100, 750), "ValueError1"], - [(200, 1001), "ValueError2"], - [(450, 400), "ValueError3"], + [(650, 700), 0.6820037131347214], + [(450, 400), "ValueError"], [(200, 1000), 1.0], ], ) def test_get_p_filter(bandpass, expected): - data_dir = Path(__file__).parents[1] / "src" / "fluopy" / "fluorophore_spectra" - fluorophore = fl.Fluorophore(name="testfluo_1", position=[0, 0]) - if expected == "ValueError1": - with pytest.raises( - ValueError, - match="The lower bandpass limit has to be between 200 and 1000 nm.", - ): - em.get_p_filter( - data_dir=data_dir, fluorophore=fluorophore, bandpass=bandpass - ) - elif expected == "ValueError2": - with pytest.raises( - ValueError, - match="The upper bandpass limit has to be between 200 and 1000 nm.", - ): - em.get_p_filter( - data_dir=data_dir, fluorophore=fluorophore, bandpass=bandpass - ) - elif expected == "ValueError3": + emission_spectrum = fd.testfluo_1.emission_spectrum + assert emission_spectrum is not None + if expected == "ValueError": with pytest.raises( ValueError, - match="The lower bandpass limit has to be smaller than the upper limit.", + match=("The lower bandpass limit has to be smaller than the upper limit."), ): - em.get_p_filter( - data_dir=data_dir, fluorophore=fluorophore, bandpass=bandpass + p_passed = em.get_p_filter( + emission_spectrum=emission_spectrum, + bandpass=bandpass, ) else: p_passed = em.get_p_filter( - data_dir=data_dir, fluorophore=fluorophore, bandpass=bandpass + emission_spectrum=emission_spectrum, + bandpass=bandpass, + ) + assert p_passed == pytest.approx(expected) + + +@pytest.mark.parametrize("bandpass", [(np.nan, 700), (650, np.inf)]) +def test_get_p_filter_non_finite_bandpass(bandpass): + emission_spectrum = fd.Spectrum( + wavelengths=[500, 600], + values=[0, 1], + ) + with pytest.raises( + ValueError, + match="bandpass limits must be finite.", + ): + em.get_p_filter( + emission_spectrum=emission_spectrum, + bandpass=bandpass, + ) + + +def test_get_p_filter_with_in_memory_spectrum(): + emission_spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[0, 1, 0], + ) + + p_passed = em.get_p_filter( + emission_spectrum=emission_spectrum, + bandpass=(505, 515), + ) + + assert p_passed == pytest.approx(0.75) + + +def test_get_p_filter_zero_emission_spectrum(): + emission_spectrum = fd.Spectrum( + wavelengths=[500, 600], + values=[0, 0], + ) + + with pytest.raises( + ValueError, + match="emission spectrum has zero total intensity.", + ): + em.get_p_filter( + emission_spectrum=emission_spectrum, + bandpass=(500, 600), ) - assert p_passed == expected + + +def test_get_p_filter_without_spectral_overlap(): + emission_spectrum = fd.Spectrum( + wavelengths=[500, 600], + values=[0, 1], + ) + + p_passed = em.get_p_filter( + emission_spectrum=emission_spectrum, + bandpass=(700, 800), + ) + + assert p_passed == 0 @pytest.mark.parametrize( @@ -61,15 +105,15 @@ def test_get_p_filter(bandpass, expected): [ (650, 700), { - 4: 0.685702066268812, - 5: 0.685702066268812, - 6: 0.685702066268812, - 7: 0.685702066268812, - 38: 0.5859441764799607, - 39: 0.5859441764799607, - 40: 0.5859441764799607, - 41: 0.5859441764799607, - 42: 0.5859441764799607, + 4: 0.6820037131347214, + 5: 0.6820037131347214, + 6: 0.6820037131347214, + 7: 0.6820037131347214, + 38: 0.5847564420110373, + 39: 0.5847564420110373, + 40: 0.5847564420110373, + 41: 0.5847564420110373, + 42: 0.5847564420110373, }, ], ], @@ -144,9 +188,9 @@ def test_emissions_simulate(tr_set_1f_bl): frames=10, store_time_points=True, ) - assert emis.event_time_points.size == 205 + assert emis.event_time_points.size == 204 exp_event_time_series = pd.Series( - np.array([0, 80, 0, 0, 0, 0, 0, 16, 52, 7, 50], dtype=np.int64), + np.array([0, 80, 0, 0, 0, 0, 0, 16, 51, 7, 50], dtype=np.int64), index=np.linspace(0, 0.001, 11), ) pd.testing.assert_series_equal(emis.event_time_series, exp_event_time_series) @@ -207,15 +251,15 @@ def test_emissions_tcspc_parameters(tr_set_bl_et_2f_diff): store_time_points=True, ) emitting_transition_ids = { - 4: 0.685702066268812, - 5: 0.685702066268812, - 6: 0.685702066268812, - 7: 0.685702066268812, - 38: 0.5859441764799607, - 39: 0.5859441764799607, - 40: 0.5859441764799607, - 41: 0.5859441764799607, - 42: 0.5859441764799607, + 4: 0.6820037131347214, + 5: 0.6820037131347214, + 6: 0.6820037131347214, + 7: 0.6820037131347214, + 38: 0.5847564420110373, + 39: 0.5847564420110373, + 40: 0.5847564420110373, + 41: 0.5847564420110373, + 42: 0.5847564420110373, } args, kwargs = mock_tcspc.call_args np.testing.assert_array_equal( diff --git a/tests/test_fluo_data.py b/tests/test_fluo_data.py index 15c084e..d13f7c6 100644 --- a/tests/test_fluo_data.py +++ b/tests/test_fluo_data.py @@ -1,17 +1,312 @@ +from pathlib import Path + +import numpy as np +import pytest + from fluopy import fluo_data as fd +def test_init_spectrum(): + spectrum = fd.Spectrum( + wavelengths=[600, 610, 620], + values=[0.1, 1.0, 0.4], + ) + + np.testing.assert_array_equal( + spectrum.wavelengths, + np.array([600.0, 610.0, 620.0]), + ) + np.testing.assert_array_equal( + spectrum.values, + np.array([0.1, 1.0, 0.4]), + ) + assert spectrum.wavelengths.dtype == np.float64 + assert spectrum.values.dtype == np.float64 + + +def test_spectrum_copies_input_arrays(): + wavelengths = np.array([600.0, 610.0]) + values = np.array([0.2, 0.8]) + + spectrum = fd.Spectrum(wavelengths=wavelengths, values=values) + wavelengths[0] = 500 + values[0] = 1 + + assert spectrum.wavelengths[0] == 600 + assert spectrum.values[0] == 0.2 + + +@pytest.mark.parametrize( + "wavelengths, values, message", + [ + ( + [[600, 610], [620, 630]], + [0.1, 0.2, 0.3, 0.4], + "spectrum wavelengths must be one-dimensional.", + ), + ( + [600, 610, 620, 630], + [[0.1, 0.2], [0.3, 0.4]], + "spectrum values must be one-dimensional.", + ), + ( + [600, 610], + [0.1, 0.2, 0.3], + "spectrum wavelengths and values must have the same length.", + ), + ( + [600], + [0.1], + "a spectrum must contain at least two data points.", + ), + ( + [600, np.nan], + [0.1, 0.2], + "spectrum wavelengths must be finite.", + ), + ( + [600, 610], + [0.1, np.inf], + "spectrum values must be finite.", + ), + ( + [600, 620, 610], + [0.1, 0.2, 0.3], + "spectrum wavelengths must be strictly increasing.", + ), + ( + [600, 600, 610], + [0.1, 0.2, 0.3], + "spectrum wavelengths must be strictly increasing.", + ), + ( + [600, 610], + [0.1, -0.2], + "spectrum values must be non-negative.", + ), + ], +) +def test_spectrum_errors(wavelengths, values, message): + with pytest.raises(ValueError, match=message): + fd.Spectrum(wavelengths=wavelengths, values=values) + + def test_init_FluorophoreData(): - fluophore_data = fd.FluorophoreData() - assert fluophore_data.QUANTUM_YIELD == 0 + fluorophore_data = fd.FluorophoreData() + assert fluorophore_data.QUANTUM_YIELD == 0 + assert fluorophore_data.emission_spectrum is None + assert fluorophore_data.absorption_spectra == {} + + +def test_fluorophore_data_with_spectra(): + emission = fd.Spectrum( + wavelengths=[600, 610, 620], + values=[0.1, 1.0, 0.4], + ) + absorption = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[1000, 2000, 500], + ) + + fluorophore_data = fd.FluorophoreData( + QUANTUM_YIELD=0.7, + FLUORESCENCE_LIFETIME=3e-9, + emission_spectrum=emission, + absorption_spectra={"s0": absorption}, + ) + + assert fluorophore_data.emission_spectrum is emission + assert fluorophore_data.absorption_spectra["s0"] is absorption + + +def test_fluorophore_data_emission_spectrum_error(): + with pytest.raises( + TypeError, + match="emission_spectrum must be a Spectrum or None.", + ): + fd.FluorophoreData(emission_spectrum=[0.1, 0.2]) + + +def test_fluorophore_data_absorption_spectrum_error(): + with pytest.raises( + TypeError, + match="absorption spectrum for state 's0' must be a Spectrum.", + ): + fd.FluorophoreData(absorption_spectra={"s0": [0.1, 0.2]}) + + +def test_spectrum_from_arrays(): + spectrum = fd.Spectrum.from_arrays( + wavelengths=[500, 510, 520], + values=[0.1, 0.8, 0.2], + ) + + np.testing.assert_array_equal( + spectrum.wavelengths, + np.array([500.0, 510.0, 520.0]), + ) + np.testing.assert_array_equal( + spectrum.values, + np.array([0.1, 0.8, 0.2]), + ) + + +def test_spectrum_from_csv(tmp_path): + path = tmp_path / "spectrum.csv" + path.write_text("Wavelengths,y\n500,0.1\n510,0.8\n520,0.2\n", encoding="utf-8") + + spectrum = fd.Spectrum.from_csv(path) + + np.testing.assert_array_equal( + spectrum.wavelengths, + np.array([500.0, 510.0, 520.0]), + ) + np.testing.assert_array_equal( + spectrum.values, + np.array([0.1, 0.8, 0.2]), + ) + + +def test_spectrum_from_csv_custom_columns(tmp_path): + path = tmp_path / "spectrum.csv" + path.write_text("wavelength,intensity\n500,0.1\n510,0.8\n", encoding="utf-8") + + spectrum = fd.Spectrum.from_csv( + path, + wavelength_column="wavelength", + value_column="intensity", + ) + + np.testing.assert_array_equal( + spectrum.wavelengths, + np.array([500.0, 510.0]), + ) + np.testing.assert_array_equal( + spectrum.values, + np.array([0.1, 0.8]), + ) + + +def test_spectrum_from_csv_missing_column(tmp_path): + path = tmp_path / "spectrum.csv" + path.write_text("wavelength,intensity\n500,0.1\n510,0.8\n", encoding="utf-8") + + with pytest.raises( + ValueError, + match="spectrum CSV is missing columns: Wavelengths, y.", + ): + fd.Spectrum.from_csv(path) + + +def test_spectrum_from_existing_csv(): + data_dir = Path(fd.__file__).parent / "fluorophore_spectra" / "testing_data_1" + + spectrum = fd.Spectrum.from_csv(data_dir / "emission.csv") + + assert spectrum.wavelengths[0] == 200 + assert spectrum.wavelengths[-1] == 1000 + assert spectrum.wavelengths.size == spectrum.values.size + + +def test_spectrum_at_existing_wavelength(): + spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[1000, 2000, 500], + ) + + assert spectrum.at(510) == 2000 + + +def test_spectrum_at_interpolated_wavelength(): + spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[1000, 2000, 500], + ) + + assert spectrum.at(505) == 1500 + + +def test_spectrum_at_boundaries(): + spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[1000, 2000, 500], + ) + + assert spectrum.at(500) == 1000 + assert spectrum.at(520) == 500 + + +@pytest.mark.parametrize("wavelength", [499, 521]) +def test_spectrum_at_outside_range(wavelength): + spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[1000, 2000, 500], + ) + + with pytest.raises( + ValueError, + match="is outside the spectrum range", + ): + spectrum.at(wavelength) + + +def test_spectrum_at_non_finite_wavelength(): + spectrum = fd.Spectrum( + wavelengths=[500, 510], + values=[1000, 2000], + ) + + with pytest.raises(ValueError, match="wavelength must be finite."): + spectrum.at(np.nan) + + +def test_spectrum_integral(): + spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[0, 1, 0], + ) + + assert spectrum.integral() == pytest.approx(10) + + +def test_spectrum_partial_integral(): + spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[0, 1, 0], + ) + + assert spectrum.integral(505, 515) == pytest.approx(7.5) + + +def test_spectrum_integral_clips_to_range(): + spectrum = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[0, 1, 0], + ) + + assert spectrum.integral(400, 600) == pytest.approx(10) + assert spectrum.integral(400, 450) == 0 + + +def test_spectrum_integral_limits_error(): + spectrum = fd.Spectrum( + wavelengths=[500, 510], + values=[0, 1], + ) + + with pytest.raises( + ValueError, + match=("the lower integration limit must be smaller than the upper limit."), + ): + spectrum.integral(510, 500) def test_init_cy5_dna(): - fluophore_data = fd.cy5_dna + fluorophore_data = fd.cy5_dna # print(fluophore_data.__doc__) - assert fluophore_data.QUANTUM_YIELD == 0.27 + assert fluorophore_data.QUANTUM_YIELD == 0.27 def test_init_atto643(): - fluophore_data = fd.atto643 - assert fluophore_data.QUANTUM_YIELD == 0.6 + fluorophore_data = fd.atto643 + assert fluorophore_data.QUANTUM_YIELD == 0.6 diff --git a/tests/test_fluorophores.py b/tests/test_fluorophores.py index 075e7a1..9872754 100644 --- a/tests/test_fluorophores.py +++ b/tests/test_fluorophores.py @@ -3,8 +3,10 @@ import numpy as np import pytest +from fluopy import emissions as em from fluopy import fluorophores as fl -from fluopy.fluo_data import FluorophoreData, testfluo_1, testfluo_2 +from fluopy import transitions as tr +from fluopy.fluo_data import FluorophoreData, Spectrum, testfluo_1, testfluo_2 @pytest.mark.parametrize( @@ -272,3 +274,141 @@ def test_construct_fluorophores(name, distance, count, expected, caplog): for fluorophore, position in zip(fluorophores, expected): assert fluorophore.name == name np.testing.assert_allclose(fluorophore.position, position, rtol=1e-5) + + +def test_custom_fluorophore_automatic_transitions_and_bandpass(): + fluorophore_data = FluorophoreData( + QUANTUM_YIELD=0.5, + FLUORESCENCE_LIFETIME=2e-9, + emission_spectrum=Spectrum( + wavelengths=[500, 510, 520], + values=[0, 1, 0], + ), + absorption_spectra={ + "s0": Spectrum( + wavelengths=[500, 510, 520], + values=[1000, 2000, 1000], + ) + }, + ) + fluorophore = fl.Fluorophore( + name="custom", + position=[0, 0], + constants=fluorophore_data, + ) + fluorophore_system = fl.FluorophoreSystem(fluorophores=[fluorophore]) + + transitions = fluorophore_system.load_transitions( + wavelength=510, + irradiance=1, + energy_transfer=False, + dstorm=False, + ) + + assert list(transitions) == ["custom"] + + excitation = next( + transition + for transition in transitions["custom"] + if transition.transition_type is tr.TransitionType.EXCITATION + ) + emission = next( + transition + for transition in transitions["custom"] + if transition.transition_type is tr.TransitionType.FLUORESCENT_EMISSION + ) + + assert excitation.rate > 0 + assert emission.rate > 0 + + transition_set = tr.TransitionSet( + transitions=transitions, + fluorophore_system=fluorophore_system, + ) + emitting_transition_ids = em.get_emitting_transition_ids( + bandpass=(505, 515), + transition_set=transition_set, + ) + + assert emitting_transition_ids + assert all( + probability == pytest.approx(0.75) + for probability in emitting_transition_ids.values() + ) + + +def test_custom_fluorophores_automatic_energy_transfer(): + donor_data = FluorophoreData( + QUANTUM_YIELD=0.5, + FLUORESCENCE_LIFETIME=2e-9, + emission_spectrum=Spectrum( + wavelengths=[500, 510, 520], + values=[0, 1, 0], + ), + absorption_spectra={ + "s0": Spectrum( + wavelengths=[505, 515], + values=[1000, 2000], + ) + }, + ) + acceptor_data = FluorophoreData( + QUANTUM_YIELD=0.6, + FLUORESCENCE_LIFETIME=3e-9, + emission_spectrum=Spectrum( + wavelengths=[505, 515, 525], + values=[0, 1, 0], + ), + absorption_spectra={ + "s0": Spectrum( + wavelengths=[500, 510, 520], + values=[500, 2000, 500], + ) + }, + ) + + donor = fl.Fluorophore( + name="custom_donor", + position=[0, 0], + constants=donor_data, + ) + acceptor = fl.Fluorophore( + name="custom_acceptor", + position=[5, 0], + constants=acceptor_data, + ) + fluorophore_system = fl.FluorophoreSystem(fluorophores=[donor, acceptor]) + + transitions = fluorophore_system.load_transitions( + wavelength=510, + irradiance=1, + energy_transfer=True, + dstorm=False, + ) + + forward_key = "D: custom_donor, A: custom_acceptor, dist: 5.0" + reverse_key = "D: custom_acceptor, A: custom_donor, dist: 5.0" + + assert forward_key in transitions + assert reverse_key in transitions + + forward_fret = next( + transition + for transition in transitions[forward_key] + if transition.transition_type is tr.TransitionType.FRET + ) + reverse_fret = next( + transition + for transition in transitions[reverse_key] + if transition.transition_type is tr.TransitionType.FRET + ) + + assert forward_fret.rate > 0 + assert reverse_fret.rate > 0 + + transition_set = tr.TransitionSet( + transitions=transitions, + fluorophore_system=fluorophore_system, + ) + + assert not transition_set.transition_df.empty diff --git a/tests/test_formulas.py b/tests/test_formulas.py index 914856a..72edeea 100644 --- a/tests/test_formulas.py +++ b/tests/test_formulas.py @@ -180,6 +180,28 @@ def test_calculate_spectral_overlap_integral(donor, acceptor, wavelengths, expec np.testing.assert_allclose(result, expected) +def test_calculate_spectral_overlap_integral_nonuniform_grid(): + result = fo.calculate_spectral_overlap_integral( + donor=[1, 1, 1], + acceptor=[1, 1, 1], + wavelengths=[1, 2, 4], + ) + + assert result == pytest.approx(93.5) + + +def test_calculate_spectral_overlap_integral_zero_donor(): + with pytest.raises( + ValueError, + match="donor emission spectrum must have positive area.", + ): + fo.calculate_spectral_overlap_integral( + donor=[0, 0, 0], + acceptor=[1, 1, 1], + wavelengths=[500, 510, 520], + ) + + def test_calculate_fret_rate(): result = fo.calculate_fret_rate( distance=10, diff --git a/tests/test_transitions.py b/tests/test_transitions.py index 23e93a9..41674c3 100644 --- a/tests/test_transitions.py +++ b/tests/test_transitions.py @@ -5,6 +5,8 @@ import pandas as pd import pytest +from fluopy import fluo_data as fd +from fluopy import formulas as fo from fluopy import transitions as tr @@ -679,6 +681,98 @@ def test_derive_energy_transfer_transitions( assert assert_rate == 0.5 * assert_total_rate +def test_derive_energy_transfer_with_in_memory_spectra(): + donor_emission = fd.Spectrum( + wavelengths=[500, 510, 520], + values=[0, 1, 0], + ) + acceptor_absorption = fd.Spectrum( + wavelengths=[505, 515], + values=[1000, 1000], + ) + + donor_data = fd.FluorophoreData( + QUANTUM_YIELD=0.5, + FLUORESCENCE_LIFETIME=2e-9, + emission_spectrum=donor_emission, + ) + acceptor_data = fd.FluorophoreData( + absorption_spectra={"s0": acceptor_absorption}, + ) + + transitions = tr.derive_energy_transfer_transitions( + donor_data=donor_data, + acceptor_data=acceptor_data, + fluorophore_ids=[(0, 1)], + dipole_orientation_factor=2 / 3, + distance=5, + refractive_index=1.33, + ) + + assert len(transitions) == 1 + assert transitions[0].transition_type is tr.TransitionType.FRET + assert transitions[0].rate > 0 + + +def test_derive_energy_transfer_without_spectral_overlap(): + donor_data = fd.FluorophoreData( + QUANTUM_YIELD=0.5, + FLUORESCENCE_LIFETIME=2e-9, + emission_spectrum=fd.Spectrum( + wavelengths=[500, 510], + values=[0, 1], + ), + ) + acceptor_data = fd.FluorophoreData( + absorption_spectra={ + "s0": fd.Spectrum( + wavelengths=[600, 610], + values=[1000, 2000], + ) + }, + ) + + transitions = tr.derive_energy_transfer_transitions( + donor_data=donor_data, + acceptor_data=acceptor_data, + fluorophore_ids=[(0, 1)], + dipole_orientation_factor=2 / 3, + distance=5, + refractive_index=1.33, + ) + + assert len(transitions) == 1 + assert transitions[0].rate == 0 + + +def test_derive_energy_transfer_without_donor_emission(): + with pytest.raises( + ValueError, + match=( + "cannot derive energy-transfer transitions without " + "a donor emission spectrum." + ), + ): + tr.derive_energy_transfer_transitions( + donor_data=fd.FluorophoreData( + QUANTUM_YIELD=0.5, + FLUORESCENCE_LIFETIME=2e-9, + ), + acceptor_data=fd.FluorophoreData( + absorption_spectra={ + "s0": fd.Spectrum( + wavelengths=[500, 510], + values=[1000, 2000], + ) + } + ), + fluorophore_ids=[(0, 1)], + dipole_orientation_factor=2 / 3, + distance=5, + refractive_index=1.33, + ) + + @pytest.mark.parametrize( "irradiance, bleaching, dstorm, summarize", [[0, False, False, False], [1, False, False, True], [1, True, True, False]], @@ -736,6 +830,61 @@ def test_derive_transitions(irradiance, bleaching, dstorm, summarize, request): assert transition.abbreviation not in summarize_checker +def test_derive_transitions_with_in_memory_absorption(): + absorption = fd.Spectrum( + wavelengths=[600, 650], + values=[1000, 2000], + ) + fluorophore_data = fd.FluorophoreData( + QUANTUM_YIELD=0.5, + FLUORESCENCE_LIFETIME=2e-9, + absorption_spectra={"s0": absorption}, + ) + + transitions = tr.derive_transitions( + fluorophore_data=fluorophore_data, + wavelength=625, + irradiance=1, + dstorm=False, + ) + + excitation = next( + transition + for transition in transitions + if transition.transition_type is tr.TransitionType.EXCITATION + ) + + _, _, frequency = fo.convert_wavenumber_wavelength_frequency(wavelength=625) + photon_flux = fo.calculate_photon_flux( + irradiance=1, + frequency=frequency, + ) + expected = fo.calculate_excitation_rate( + photon_flux=photon_flux, + extinction_coefficient=1500, + ) + + assert excitation.rate == pytest.approx(expected) + + +def test_derive_transitions_without_absorption_error(): + fluorophore_data = fd.FluorophoreData( + QUANTUM_YIELD=0.5, + FLUORESCENCE_LIFETIME=2e-9, + ) + + with pytest.raises( + ValueError, + match=( + "cannot derive excitation transition without an S0 absorption spectrum." + ), + ): + tr.derive_transitions( + fluorophore_data=fluorophore_data, + dstorm=False, + ) + + def test_interpolate_data(): data = pd.DataFrame( {"Wavelengths": [617, 619, 620, 621], "y": [0.5, 0.6, 0.7, 0.6]}