Modify spectra handling

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,

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()