diff --git a/include/evalhyd/detail/determinist/evaluator.hpp b/include/evalhyd/detail/determinist/evaluator.hpp
index 0fdab32c129ec69faa22bc9a0f4f93e6646fce1f..6295a20130e5cc34ad993185203cc8a9a34597ef 100644
--- a/include/evalhyd/detail/determinist/evaluator.hpp
+++ b/include/evalhyd/detail/determinist/evaluator.hpp
@@ -12,12 +12,13 @@ namespace evalhyd
 {
     namespace determinist
     {
+        template <class D2, class B2>
         class Evaluator
         {
         private:
             // members for input data
-            const xt::xtensor<double, 2>& q_obs;
-            const xt::xtensor<double, 2>& q_prd;
+            const D2& q_obs;
+            const D2& q_prd;
             xt::xtensor<bool, 3> t_msk;
             const std::vector<xt::xkeep_slice<int>>& b_exp;
 
@@ -43,9 +44,9 @@ namespace evalhyd
 
         public:
             // constructor method
-            Evaluator(const xt::xtensor<double, 2>& obs,
-                      const xt::xtensor<double, 2>& prd,
-                      const xt::xtensor<bool, 2>& msk,
+            Evaluator(const D2& obs,
+                      const D2& prd,
+                      const B2& msk,
                       const std::vector<xt::xkeep_slice<int>>& exp) :
                     q_obs{obs}, q_prd{prd}, b_exp{exp}
             {
@@ -111,7 +112,8 @@ namespace evalhyd
         // \assign mean_obs:
         //     Mean observed streamflow.
         //     shape: (subsets, samples, series, 1)
-        void Evaluator::calc_mean_obs()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_mean_obs()
         {
             mean_obs = xt::zeros<double>(mean_dims);
             for (int m = 0; m < n_msk; m++) {
@@ -137,7 +139,8 @@ namespace evalhyd
         // \assign mean_prd:
         //     Mean predicted streamflow.
         //     shape: (subsets, samples, series, 1)
-        void Evaluator::calc_mean_prd()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_mean_prd()
         {
             mean_prd = xt::zeros<double>(mean_dims);
             for (int m = 0; m < n_msk; m++) {
@@ -166,8 +169,8 @@ namespace evalhyd
         // \assign quad_err:
         //     Quadratic errors between observations and predictions.
         //     shape: (series, time)
-
-        void Evaluator::calc_quad_err()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_quad_err()
         {
             quad_err = xt::square(q_obs - q_prd);
         }
@@ -183,7 +186,8 @@ namespace evalhyd
         // \assign quad_obs:
         //     Quadratic errors between observations and mean observation.
         //     shape: (subsets, samples, series, time)
-        void Evaluator::calc_quad_obs()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_quad_obs()
         {
             quad_obs = xt::zeros<double>(inter_dims);
             for (int m = 0; m < n_msk; m++) {
@@ -210,7 +214,8 @@ namespace evalhyd
         // \assign quad_prd:
         //     Quadratic errors between predictions and mean prediction.
         //     shape: (subsets, samples, series, time)
-        void Evaluator::calc_quad_prd()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_quad_prd()
         {
             quad_prd = xt::zeros<double>(inter_dims);
             for (int m = 0; m < n_msk; m++) {
@@ -249,7 +254,8 @@ namespace evalhyd
         // \assign r_pearson:
         //     Pearson correlation coefficients.
         //     shape: (subsets, samples, series)
-        void Evaluator::calc_r_pearson()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_r_pearson()
         {
             // calculate error in timing and dynamics $r_{pearson}$
             // (Pearson's correlation coefficient)
@@ -298,7 +304,8 @@ namespace evalhyd
         // \assign alpha:
         //     Alphas, ratios of standard deviations.
         //     shape: (subsets, samples, series)
-        void Evaluator::calc_alpha()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_alpha()
         {
             // calculate error in spread of flow $alpha$
             alpha = xt::zeros<double>(inner_dims);
@@ -330,7 +337,8 @@ namespace evalhyd
         // \assign bias:
         //     Biases.
         //     shape: (subsets, samples, series)
-        void Evaluator::calc_bias()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_bias()
         {
             // calculate $bias$
             bias = xt::zeros<double>(inner_dims);
@@ -358,7 +366,8 @@ namespace evalhyd
         // \assign RMSE:
         //     Root-mean-square errors.
         //     shape: (series, subsets, samples)
-        void Evaluator::calc_RMSE()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_RMSE()
         {
             // compute RMSE
             RMSE = xt::zeros<double>(final_dims);
@@ -386,7 +395,8 @@ namespace evalhyd
         // \assign NSE:
         //     Nash-Sutcliffe efficiencies.
         //     shape: (series, subsets, samples)
-        void Evaluator::calc_NSE()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_NSE()
         {
             NSE = xt::zeros<double>(final_dims);
             for (int m = 0; m < n_msk; m++) {
@@ -424,7 +434,8 @@ namespace evalhyd
         // \assign KGE:
         //     Kling-Gupta efficiencies.
         //     shape: (series, subsets, samples)
-        void Evaluator::calc_KGE()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_KGE()
         {
             KGE = xt::zeros<double>(final_dims);
             for (int m = 0; m < n_msk; m++) {
@@ -459,7 +470,8 @@ namespace evalhyd
         // \assign KGEPRIME:
         //     Modified Kling-Gupta efficiencies.
         //     shape: (series, subsets, samples)
-        void Evaluator::calc_KGEPRIME()
+        template <class D2, class B2>
+        void Evaluator<D2, B2>::calc_KGEPRIME()
         {
             KGEPRIME = xt::zeros<double>(final_dims);
             for (int m = 0; m < n_msk; m++) {
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index 5dc85412c0d0bdb264e44e8aa787b7a8774428d9..98409964faf7c07c893f84c35690fbf5c323367b 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -310,7 +310,7 @@ namespace evalhyd
         }
 
         // instantiate determinist evaluator
-        eh::determinist::Evaluator evaluator(obs, prd, msk, exp);
+        eh::determinist::Evaluator<D2, B2> evaluator(obs, prd, msk, exp);
 
         // declare maps for memoisation purposes
         std::unordered_map<std::string, std::vector<std::string>> elt;