Add fully-eigen interface to IOscCalc, forwards to vector interface b…

…y default. Change Eigen return type to Array from Matrix. It seems to me these are more a list of values than something truly vector-y.
chenel · Oct 29, 2020 · 284e963 · 284e963
1 parent 22fd646
commit 284e963
Show file tree

Hide file tree

Showing 6 changed files with 18 additions and 8 deletions.
diff --git a/OscLib/IOscCalc.cxx b/OscLib/IOscCalc.cxx
@@ -23,9 +23,8 @@ namespace{
 
 namespace osc
 {
-  template<class T>
-  Eigen::Matrix<T,Eigen::Dynamic,1>
-  _IOscCalc<T>::P(int flavBefore, int flavAfter, const std::vector<double> &E)
+  template<class T> Eigen::Array<T, Eigen::Dynamic, 1> _IOscCalc<T>::
+  P(int flavBefore, int flavAfter, const std::vector<double>& E)
   {
     Eigen::Matrix<T,Eigen::Dynamic,1> ret(E.size());
     for(auto i = 0u; i < E.size(); i++) {
@@ -34,6 +33,12 @@ namespace osc
     return ret.array().isNaN().select(0, ret);
   }
 
+  template<class T> Eigen::Array<T, Eigen::Dynamic, 1> _IOscCalc<T>::
+  P(int flavBefore, int flavAfter, const Eigen::ArrayXd& E)
+  {
+    return P(flavBefore, flavAfter, std::vector<double>(&E[0], &E[0]+E.size()));
+  }
+
   //---------------------------------------------------------------------------
   template<class T> _IOscCalcAdjustable<T>::~_IOscCalcAdjustable()
   {

diff --git a/OscLib/IOscCalc.h b/OscLib/IOscCalc.h
@@ -26,7 +26,12 @@ namespace osc
 
     /// E in GeV; flavors as PDG codes (so, neg==>antinu)
     virtual T P(int flavBefore, int flavAfter, double E) = 0;
-    virtual Eigen::Matrix<T,Eigen::Dynamic,1> P(int flavBefore, int flavAfter, const std::vector<double> &E);
+    /// Default implementation forwards to non-vector version using a simple
+    /// loop. Override if your calculator has a more efficient implementation.
+    virtual Eigen::Array<T, Eigen::Dynamic, 1> P(int flavBefore, int flavAfter, const std::vector<double>& E);
+    /// Default implementation forawrds to vector<double> version. Override if
+    /// your calculator has a more efficient implementation.
+    virtual Eigen::Array<T, Eigen::Dynamic, 1> P(int flavBefore, int flavAfter, const Eigen::ArrayXd& E);
 
     /// \brief Use to check two calculators are in the same state
     ///

diff --git a/OscLib/OscCalcDMP.cxx b/OscLib/OscCalcDMP.cxx
@@ -81,7 +81,7 @@ namespace osc
 
 
   template<typename T>
-  Matrix<T, Dynamic, 1> _OscCalcDMP<T>::P(int flavBefore, int flavAfter, const std::vector<double> &E)
+  Array<T, Dynamic, 1> _OscCalcDMP<T>::P(int flavBefore, int flavAfter, const std::vector<double> &E)
   {
     if (fCache.energies.size() != (size_t) fCache.probabilities.cols() &&
         fCache.energies.size() != 0)

diff --git a/OscLib/OscCalcDMP.h b/OscLib/OscCalcDMP.h
@@ -47,7 +47,7 @@ namespace osc
     T P(int flavBefore, int flavAfter, double E, bool fast_and_loose);
 
     Array<T,Dynamic,Dynamic> P() {return this->fCache.probabilities;}
-    Matrix<T,Dynamic,1>      P(int flavBefore, int flavAfter, const std::vector<double> &E) override;
+    Array<T,Dynamic,1>      P(int flavBefore, int flavAfter, const std::vector<double> &E) override;
 
     void SetL     (double L       ) override {SaveLastParams(); this->fL      = L;}
     void SetRho   (double rho     ) override {SaveLastParams(); this->fRho    = rho;}

diff --git a/OscLib/OscCalcPMNSOptEigen.cxx b/OscLib/OscCalcPMNSOptEigen.cxx
@@ -93,7 +93,7 @@ namespace osc {
   }
 
 
-  Eigen::VectorXd
+  Eigen::ArrayXd
   OscCalcPMNSOptEigen::P(int flavBefore, int flavAfter, const std::vector<double>& E)
   {
     // Are there probabilities cached and can we use them?

diff --git a/OscLib/OscCalcPMNSOptEigen.h b/OscLib/OscCalcPMNSOptEigen.h
@@ -86,7 +86,7 @@ namespace osc
 
     double P(int flavBefore, int flavAfter, double E) override;
     double P(int flavBefore, int flavAfter, double E, bool fast_and_loose);
-    Eigen::VectorXd P(int flavBefore, int flavAfter, 
+    Eigen::ArrayXd P(int flavBefore, int flavAfter,
 		     const std::vector<double>& E) override;