feat: add NLL::project_with to do projecting and isolation in one step

This resets the `NLL` to its previous activated state after processing the weights

python/laddu/likelihoods/__init__.pyi

@@ -49,6 +49,9 @@ class NLL:
     def isolate(self, name: str | list[str]) -> None: ...
     def evaluate(self, parameters: list[float] | npt.NDArray[np.float64]) -> float: ...
     def project(self, parameters: list[float] | npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: ...
+    def project_with(
+        self, parameters: list[float] | npt.NDArray[np.float64], name: str | list[str]
+    ) -> npt.NDArray[np.float64]: ...
     def minimize(
         self,
         p0: list[float] | npt.NDArray[np.float64],

src/likelihoods.rs

@@ -193,6 +193,68 @@ impl NLL {
             .map(|(l, e)| e.weight * l.re / n_mc)
             .collect()
     }
+
+    /// Project the stored [`Expression`] over the events in the [`Dataset`] stored by the
+    /// [`Evaluator`] with the given values for free parameters to obtain weights for each
+    /// Monte-Carlo event. This method differs from the standard [`NLL::project`] in that it first
+    /// isolates the selected [`Amplitude`](`crate::amplitudes::Amplitude`)s by name, but returns
+    /// the [`NLL`] to its prior state after calculation.
+    ///
+    /// This method takes the real part of the given expression (discarding
+    /// the imaginary part entirely, which does not matter if expressions are coherent sums
+    /// wrapped in [`Expression::norm_sqr`]). Event weights are determined by the following
+    /// formula:
+    ///
+    /// ```math
+    /// \text{weight}(\vec{p}; e) = \text{weight}(e) \mathcal{L}(e) / N_{\text{MC}}
+    /// ```
+    #[cfg(feature = "rayon")]
+    pub fn project_with<T: AsRef<str>>(&self, parameters: &[Float], names: &[T]) -> Vec<Float> {
+        let current_active_data = self.data_evaluator.resources.read().active.clone();
+        let current_active_mc = self.mc_evaluator.resources.read().active.clone();
+        self.isolate_many(names);
+        let mc_result = self.mc_evaluator.evaluate(parameters);
+        let n_mc = self.mc_evaluator.dataset.len() as Float;
+        let res = mc_result
+            .par_iter()
+            .zip(self.mc_evaluator.dataset.par_iter())
+            .map(|(l, e)| e.weight * l.re / n_mc)
+            .collect();
+        self.data_evaluator.resources.write().active = current_active_data;
+        self.mc_evaluator.resources.write().active = current_active_mc;
+        res
+    }
+
+    /// Project the stored [`Expression`] over the events in the [`Dataset`] stored by the
+    /// [`Evaluator`] with the given values for free parameters to obtain weights for each
+    /// Monte-Carlo event. This method differs from the standard [`NLL::project`] in that it first
+    /// isolates the selected [`Amplitude`](`crate::amplitudes::Amplitude`)s by name, but returns
+    /// the [`NLL`] to its prior state after calculation.
+    ///
+    /// This method takes the real part of the given expression (discarding
+    /// the imaginary part entirely, which does not matter if expressions are coherent sums
+    /// wrapped in [`Expression::norm_sqr`]). Event weights are determined by the following
+    /// formula:
+    ///
+    /// ```math
+    /// \text{weight}(\vec{p}; e) = \text{weight}(e) \mathcal{L}(e) / N_{\text{MC}}
+    /// ```
+    #[cfg(not(feature = "rayon"))]
+    pub fn project_with<T: AsRef<str>>(&self, parameters: &[Float], names: &[T]) -> Vec<Float> {
+        let current_active_data = self.data_evaluator.resources.read().active.clone();
+        let current_active_mc = self.mc_evaluator.resources.read().active.clone();
+        self.isolate_many(names);
+        let mc_result = self.mc_evaluator.evaluate(parameters);
+        let n_mc = self.mc_evaluator.dataset.len() as Float;
+        let res = mc_result
+            .iter()
+            .zip(self.mc_evaluator.dataset.iter())
+            .map(|(l, e)| e.weight * l.re / n_mc)
+            .collect();
+        self.data_evaluator.resources.write().active = current_active_data;
+        self.mc_evaluator.resources.write().active = current_active_mc;
+        res
+    }
 }
 
 impl LikelihoodTerm for NLL {

src/python.rs

@@ -2016,6 +2016,53 @@ pub(crate) mod laddu {
             PyArray1::from_slice_bound(py, &self.0.project(&parameters))
         }
 
+        /// Project the model over the Monte Carlo dataset with the given parameter values, first
+        /// isolating the given terms by name. The NLL is then reset to its previous state of
+        /// activation.
+        ///
+        /// This is defined as
+        ///
+        /// .. math:: e_w(\vec{p}) = \frac{e_w}{N_{MC}} \mathcal{L}(e)
+        ///
+        /// Parameters
+        /// ----------
+        /// parameters : list of float
+        ///     The values to use for the free parameters
+        /// arg : str or list of str
+        ///     Names of Amplitudes to be isolated
+        ///
+        /// Returns
+        /// -------
+        /// result : array_like
+        ///     Weights for every Monte Carlo event which represent the fit to data
+        ///
+        /// Raises
+        /// ------
+        /// TypeError
+        ///     If `arg` is not a str or list of str
+        ///
+        fn project_with<'py>(
+            &self,
+            py: Python<'py>,
+            parameters: Vec<Float>,
+            arg: &Bound<'_, PyAny>,
+        ) -> PyResult<Bound<'py, PyArray1<Float>>> {
+            let names = if let Ok(string_arg) = arg.extract::<String>() {
+                vec![string_arg]
+            } else if let Ok(list_arg) = arg.downcast::<PyList>() {
+                let vec: Vec<String> = list_arg.extract()?;
+                vec
+            } else {
+                return Err(PyTypeError::new_err(
+                    "Argument must be either a string or a list of strings",
+                ));
+            };
+            Ok(PyArray1::from_slice_bound(
+                py,
+                &self.0.project_with(&parameters, &names),
+            ))
+        }
+
         /// Minimize the NLL with respect to the free parameters in the model
         ///
         /// This method "runs the fit". Given an initial position `p0` and optional `bounds`, this