NLL Terms

Cargo.toml

@@ -29,9 +29,10 @@ rand = "0.8.5"
 rayon = { version = "1.10.0", optional = true }
 pyo3 = { version = "0.22.5", optional = true, features = ["num-complex"] }
 numpy = { version = "0.22.0", optional = true, features = ["nalgebra"] }
-ganesh = "0.12.1"
+ganesh = "0.12.2"
 thiserror = "1.0.64"
 shellexpand = "3.1.0"
+accurate = "0.4.1"
 
 [dev-dependencies]
 approx = "0.5.1"
@@ -46,7 +47,7 @@ default = ["rayon", "python"]
 extension-module = ["pyo3/extension-module"]
 rayon = ["dep:rayon"]
 f32 = []
-python = ["dep:pyo3", "dep:numpy"]
+python = ["dep:pyo3", "dep:numpy", "extension-module"]
 
 [profile.release]
 lto = true

README.md

@@ -42,6 +42,7 @@
     - [Calculating a Likelihood](#calculating-a-likelihood)
   - [Python](#python)
     - [Fitting Data](#fitting-data)
+    - [Other Examples](#other-examples)
 - [Data Format](#data-format)
 - [Future Plans](#future-plans)
 - [Alternatives](#alternatives)
@@ -218,6 +219,24 @@ if __name__ == "__main__":
     main()
 ```
 This example would probably make the most sense for a binned fit, since there isn't actually any mass dependence in any of these amplitudes (so it will just plot the relative amount of each wave over the entire dataset).
+
+### Other Examples
+You can find other Python examples in the `python_examples` folder. They should each have a corresponding `requirements_[#].txt` file.
+
+#### Example 1
+
+The first example script uses data generated with [gen_amp](https://github.com/JeffersonLab/halld_sim/tree/962c1fffc29eb4801b146d0a7f1e9aecb417374a/src/programs/Simulation/gen_amp). These data consist of a data file with two resonances, an $`f_0(1500)`$ modeled as a Breit-Wigner with a mass of $`1506\text{ MeV}/c^2`$ and a width of $`112\text{ MeV}/c^2`$ and an $`f_2'(1525)`$, also modeled as a Breit-Wigner, with a mass of $`1517\text{ MeV}/c^2`$ and a width of $`86\text{ MeV}/c^2`$, as per the [PDG](https://pdg.lbl.gov/2020/tables/rpp2020-tab-mesons-light.pdf). These were generated to decay to pairs of $`K_S^0`$s and are produced via photoproduction off a proton target (as in the GlueX experiment). The beam photon is polarized with an angle of $`0`$ degrees relative to the production plane and a polarization magnitude of $`0.3519`$ (out of unity). The configuration file used to generate the corresponding data and Monte Carlo files can also be found in the `python_examples`, and the datasets contain $`100,000`$ data events and $`1,000,000`$ Monte Carlo events (generated with the `-f` argument to create a Monte Carlo file without resonances). The result of this fit can be seen in the following image (using the default 50 bins):
+
+<p align="center">
+  <img
+    width="800"
+    src="python_examples/example_1.svg"
+  />
+</p>
+
+> [!NOTE]
+> There appears to be an overall scale factor difference between `gen_amp` and `laddu`, so the raw values for the real and imaginary parts of each wave do not match even though the overall fit, the ratio of S- to D-wave, and relative phase between the waves are consistent.
+
 # Data Format
 The data format for `laddu` is a bit different from some of the alternatives like [`AmpTools`](https://github.com/mashephe/AmpTools). Since ROOT doesn't yet have bindings to Rust and projects to read ROOT files are still largely works in progress (although I hope to use [`oxyroot`](https://github.com/m-dupont/oxyroot) in the future when I can figure out a few bugs), the primary interface for data in `laddu` is Parquet files. These are easily accessible from almost any other language and they don't take up much more space than ROOT files. In the interest of future compatibility with any number of experimental setups, the data format consists of an arbitrary number of columns containing the four-momenta of each particle, the polarization vector of each particle (optional) and a single column for the weight. These columns all have standardized names. For example, the following columns would describe a dataset with four particles, the first of which is a polarized photon beam, as in the GlueX experiment:
 | Column name | Data Type | Interpretation |

benches/kmatrix_benchmark.rs

@@ -1,13 +1,16 @@
+use std::time::Duration;
+
 use criterion::{black_box, criterion_group, criterion_main, BatchSize, Criterion};
 use laddu::{
     amplitudes::{
         constant,
         kmatrix::{KopfKMatrixA0, KopfKMatrixA2, KopfKMatrixF0, KopfKMatrixF2},
         parameter,
         zlm::Zlm,
-        Manager, NLL,
+        Manager,
     },
     data::open,
+    likelihoods::{LikelihoodTerm, NLL},
     utils::{
         enums::{Frame, Sign},
         variables::{Angles, Mass, Polarization},
@@ -136,7 +139,7 @@ fn kmatrix_nll_benchmark(c: &mut Criterion) {
     let neg_re = (&s0n * z00n.real()).norm_sqr();
     let neg_im = (&s0n * z00n.imag()).norm_sqr();
     let model = pos_re + pos_im + neg_re + neg_im;
-    let nll = NLL::new(&manager, &ds_data, &ds_mc);
+    let nll = NLL::new(&manager, &ds_data, &ds_mc, &model);
     let mut rng = rand::thread_rng();
     let range = Uniform::new(-100.0, 100.0);
     c.bench_function("kmatrix benchmark (nll)", |b| {
@@ -147,11 +150,15 @@ fn kmatrix_nll_benchmark(c: &mut Criterion) {
                     .collect();
                 p
             },
-            |p| black_box(nll.evaluate(&model, &p)),
+            |p| black_box(nll.evaluate(&p)),
             BatchSize::SmallInput,
         )
     });
 }
 
-criterion_group!(benches, kmatrix_nll_benchmark);
+criterion_group! {
+    name = benches;
+    config = Criterion::default().measurement_time(Duration::from_secs(30)).sample_size(5000);
+    targets = kmatrix_nll_benchmark
+}
 criterion_main!(benches);

pyproject.toml

@@ -11,9 +11,21 @@ classifiers = [
   "Programming Language :: Python :: Implementation :: PyPy",
 ]
 dynamic = ["version"]
-dependencies = ["numpy"]
+dependencies = [
+  "numpy",
+  "docopt-ng",
+  "loguru",
+  "pandas",
+  "uproot",
+  "fastparquet",
+]
+
+[project.scripts]
+amptools-to-laddu = "laddu:convert.run"
+
 [project.optional-dependencies]
 tests = ["pytest"]
+
 [tool.maturin]
 python-source = "python"
 features = ["pyo3/extension-module"]

python/laddu/__init__.py

@@ -1,35 +1,40 @@
 from abc import ABCMeta, abstractmethod
 
-from laddu.amplitudes import NLL, Expression, Manager, Status, constant, parameter
+from laddu.amplitudes import Manager, constant, parameter
 from laddu.amplitudes.breit_wigner import BreitWigner
 from laddu.amplitudes.common import ComplexScalar, PolarComplexScalar, Scalar
 from laddu.amplitudes.ylm import Ylm
 from laddu.amplitudes.zlm import Zlm
+from laddu.convert import convert_from_amptools
 from laddu.data import BinnedDataset, Dataset, open, open_binned
+from laddu.likelihoods import NLL, LikelihoodManager, Status
 from laddu.utils.variables import Angles, CosTheta, Mass, Phi, PolAngle, Polarization, PolMagnitude
 from laddu.utils.vectors import Vector3, Vector4
 
-from . import amplitudes, data, utils
+from . import amplitudes, convert, data, likelihoods, utils
 from .laddu import version
 
 __version__ = version()
 
 
 class Observer(metaclass=ABCMeta):
     @abstractmethod
-    def callback(self, step: int, status: Status, expression: Expression) -> tuple[Status, Expression, bool]:
+    def callback(self, step: int, status: Status) -> tuple[Status, bool]:
         pass
 
 
 __all__ = [
     "__version__",
+    "convert",
+    "convert_from_amptools",
     "Dataset",
     "open",
     "BinnedDataset",
     "open_binned",
     "utils",
     "data",
     "amplitudes",
+    "likelihoods",
     "Vector3",
     "Vector4",
     "CosTheta",
@@ -40,8 +45,8 @@ def callback(self, step: int, status: Status, expression: Expression) -> tuple[S
     "Polarization",
     "Mass",
     "Manager",
+    "LikelihoodManager",
     "NLL",
-    "Expression",
     "Status",
     "Observer",
     "parameter",

python/laddu/__init__.pyi

@@ -1,24 +1,28 @@
 from abc import ABCMeta, abstractmethod
 
-from laddu.amplitudes import NLL, Expression, Manager, Status, constant, parameter
+from laddu.amplitudes import Expression, Manager, constant, parameter
 from laddu.amplitudes.breit_wigner import BreitWigner
 from laddu.amplitudes.common import ComplexScalar, PolarComplexScalar, Scalar
 from laddu.amplitudes.ylm import Ylm
 from laddu.amplitudes.zlm import Zlm
+from laddu.convert import convert_from_amptools
 from laddu.data import BinnedDataset, Dataset, open, open_binned
+from laddu.likelihoods import NLL, LikelihoodManager, Status
 from laddu.utils.variables import Angles, CosTheta, Mass, Phi, PolAngle, Polarization, PolMagnitude
 from laddu.utils.vectors import Vector3, Vector4
 
-from . import amplitudes, data, utils
+from . import amplitudes, convert, data, utils
 
 class Observer(metaclass=ABCMeta):
     @abstractmethod
-    def callback(self, step: int, status: Status, expression: Expression) -> tuple[Status, Expression, bool]: ...
+    def callback(self, step: int, status: Status) -> tuple[Status, bool]: ...
 
 __version__: str
 
 __all__ = [
     "__version__",
+    "convert",
+    "convert_from_amptools",
     "Dataset",
     "open",
     "BinnedDataset",
@@ -36,6 +40,7 @@ __all__ = [
     "Polarization",
     "Mass",
     "Manager",
+    "LikelihoodManager",
     "NLL",
     "Expression",
     "Status",

python/laddu/amplitudes/__init__.py

@@ -1,14 +1,11 @@
 from laddu.amplitudes import breit_wigner, common, kmatrix, ylm, zlm
 from laddu.laddu import (
-    NLL,
     Amplitude,
     AmplitudeID,
-    Bound,
     Evaluator,
     Expression,
     Manager,
     ParameterLike,
-    Status,
     constant,
     parameter,
 )
@@ -19,7 +16,6 @@
     "Amplitude",
     "Manager",
     "Evaluator",
-    "NLL",
     "ParameterLike",
     "parameter",
     "constant",
@@ -28,6 +24,4 @@
     "zlm",
     "breit_wigner",
     "kmatrix",
-    "Status",
-    "Bound",
 ]

python/laddu/amplitudes/__init__.pyi

@@ -1,5 +1,3 @@
-from typing import Literal
-
 import numpy as np
 import numpy.typing as npt
 
@@ -30,7 +28,7 @@ class Amplitude: ...
 class Manager:
     def __init__(self) -> None: ...
     def register(self, amplitude: Amplitude) -> AmplitudeID: ...
-    def load(self, dataset: Dataset) -> Evaluator: ...
+    def load(self, dataset: Dataset, expression: Expression) -> Evaluator: ...
 
 class Evaluator:
     parameters: list[str]
@@ -39,58 +37,14 @@ class Evaluator:
     def deactivate(self, name: str | list[str]) -> None: ...
     def deactivate_all(self) -> None: ...
     def isolate(self, name: str | list[str]) -> None: ...
-    def evaluate(
-        self, expression: Expression, parameters: list[float] | npt.NDArray[np.float64]
-    ) -> npt.NDArray[np.complex128]: ...
-
-class NLL:
-    parameters: list[str]
-    def __init__(self, manager: Manager, ds_data: Dataset, ds_mc: Dataset) -> None: ...
-    def activate(self, name: str | list[str]) -> None: ...
-    def activate_all(self) -> None: ...
-    def deactivate(self, name: str | list[str]) -> None: ...
-    def deactivate_all(self) -> None: ...
-    def isolate(self, name: str | list[str]) -> None: ...
-    def evaluate(self, expression: Expression, parameters: list[float] | npt.NDArray[np.float64]) -> float: ...
-    def project(
-        self, expression: Expression, parameters: list[float] | npt.NDArray[np.float64]
-    ) -> npt.NDArray[np.float64]: ...
-    def minimize(
-        self,
-        expression: Expression,
-        p0: list[float],
-        bounds: list[tuple[float | None, float | None]] | None = None,
-        method: Literal["lbfgsb", "nelder_mead"] = "lbfgsb",
-        max_steps: int = 4000,
-        debug: bool = False,  # noqa: FBT001, FBT002
-        verbose: bool = False,  # noqa: FBT001, FBT002
-        **kwargs,
-    ) -> Status: ...
-
-class Status:
-    x: npt.NDArray[np.float64]
-    err: npt.NDArray[np.float64] | None
-    x0: npt.NDArray[np.float64]
-    fx: float
-    cov: npt.NDArray[np.float64] | None
-    hess: npt.NDArray[np.float64] | None
-    message: str
-    converged: bool
-    bounds: list[Bound] | None
-    n_f_evals: int
-    n_g_evals: int
-
-class Bound:
-    lower: float
-    upper: float
+    def evaluate(self, parameters: list[float] | npt.NDArray[np.float64]) -> npt.NDArray[np.complex128]: ...
 
 __all__ = [
     "AmplitudeID",
     "Expression",
     "Amplitude",
     "Manager",
     "Evaluator",
-    "NLL",
     "ParameterLike",
     "parameter",
     "constant",
@@ -99,6 +53,4 @@ __all__ = [
     "zlm",
     "breit_wigner",
     "kmatrix",
-    "Status",
-    "Bound",
 ]