refactor bindings to avoid copies of 1D arrays in evalp

Since, by definition, `xt::pytensor` cannot be reshaped, and since
`py::array` cannot be manipulated in an intermediate step on the C++
side without compromising on it being accepted as an `xt::xexpression`
type afterwards, the only viable option seemed to be to add a Python
layer to the bindings so that the 1D numpy arrays can be reshaped into
2D numpy array views before calling the C++ extension.

As part of this refactoring, the default values for optional parameters
are also removed, and the function parameter types are specified
directly in the Python layer.

evalhyd/__init__.py

+from .evald import evald
+from .evalp import evalp

evalhyd/evald.py

+from typing import List, Dict
+from numpy import dtype
+from numpy.typing import NDArray
+import evalhyd.core
+
+
+def evald(q_obs: NDArray[dtype('float64')],
+          q_prd: NDArray[dtype('float64')],
+          metrics: List[str],
+          transform: str = None,
+          exponent: float = None,
+          epsilon: float = None,
+          t_msk: NDArray[dtype('bool')] = None,
+          m_cdt: NDArray[dtype('S32')] = None,
+          bootstrap: Dict[str, int] = None,
+          dts: List[str] = None) -> List[NDArray[dtype('float64')]]:
+    """Function to evaluate determinist streamflow predictions"""
+    # required arguments
+    kwargs = {
+        # convect 1D array into 2D array view
+        'q_obs': q_obs.reshape(1, q_obs.size) if q_obs.ndim == 1 else q_obs,
+        'q_prd': q_prd.reshape(1, q_prd.size) if q_prd.ndim == 1 else q_prd,
+        'metrics': metrics
+    }
+
+    # optional arguments
+    if transform is not None:
+        kwargs['transform'] = transform
+    if exponent is not None:
+        kwargs['exponent'] = exponent
+    if epsilon is not None:
+        kwargs['epsilon'] = epsilon
+    if t_msk is not None:
+        kwargs['t_msk'] = t_msk
+    if m_cdt is not None:
+        kwargs['m_cdt'] = m_cdt
+    if bootstrap is not None:
+        kwargs['bootstrap'] = bootstrap
+    if dts is not None:
+        kwargs['dts'] = dts
+
+    return evalhyd.core._evald(**kwargs)

evalhyd/evalp.py

+from typing import List, Dict
+from numpy import dtype
+from numpy.typing import NDArray
+import evalhyd.core
+
+
+def evalp(q_obs: NDArray[dtype('float64')],
+          q_prd: NDArray[dtype('float64')],
+          metrics: List[str],
+          q_thr: NDArray[dtype('float64')] = None,
+          t_msk: NDArray[dtype('bool')] = None,
+          m_cdt: NDArray[dtype('S32')] = None,
+          bootstrap: Dict[str, int] = None,
+          dts: List[str] = None) -> List[NDArray[dtype('float64')]]:
+    """Function to evaluate probabilist streamflow predictions"""
+    # required arguments
+    kwargs = {
+        'q_obs': q_obs,
+        'q_prd': q_prd,
+        'metrics': metrics
+    }
+
+    # optional arguments
+    if q_thr is not None:
+        kwargs['q_thr'] = q_thr
+    if t_msk is not None:
+        kwargs['t_msk'] = t_msk
+    if m_cdt is not None:
+        kwargs['m_cdt'] = m_cdt
+    if bootstrap is not None:
+        kwargs['bootstrap'] = bootstrap
+    if dts is not None:
+        kwargs['dts'] = dts
+
+    return evalhyd.core._evalp(**kwargs)

src/evalhyd-python.cpp → evalhyd/src/core.cpp

 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 #include <array>
+#include <optional>
 
 #define STRINGIFY(x) #x
 #define MACRO_STRINGIFY(x) STRINGIFY(x)
 
 #define FORCE_IMPORT_ARRAY
+#include <xtl/xoptional.hpp>
 #include <xtensor/xview.hpp>
 #include <xtensor-python/pytensor.hpp>
 
@@ -15,83 +17,96 @@
 namespace py = pybind11;
 using namespace py::literals;
 
-// reshape 1D tensors to 2D tensors
-auto evald_1d(
-    const xt::pytensor<double, 1>& q_obs,
-    const xt::pytensor<double, 1>& q_prd,
+auto evald(
+    const xt::pytensor<double, 2>& q_obs,
+    const xt::pytensor<double, 2>& q_prd,
     const std::vector<std::string>& metrics,
-    const std::string& transform,
-    const double exponent,
-    double epsilon,
+    std::optional<std::string> transform,
+    std::optional<double> exponent,
+    std::optional<double> epsilon,
     const xt::pytensor<bool, 2>& t_msk,
     const xt::pytensor<std::array<char, 32>, 1>& m_cdt,
-    const std::unordered_map<std::string, int>& bootstrap,
+    std::optional<std::unordered_map<std::string, int>> bootstrap,
     const std::vector<std::string>& dts
 )
 {
-    return evalhyd::evald<xt::pytensor<double, 2>, xt::pytensor<bool, 2>>(
-        xt::pytensor<double, 2>(xt::view(q_obs, xt::newaxis(), xt::all())),
-        xt::pytensor<double, 2>(xt::view(q_prd, xt::newaxis(), xt::all())),
+    return evalhyd::evald(
+        q_obs,
+        q_prd,
         metrics,
-        transform,
-        exponent,
-        epsilon,
+        (transform.has_value()) ? transform.value() : xtl::missing<std::string>(),
+        (exponent.has_value()) ? exponent.value() : xtl::missing<double>(),
+        (epsilon.has_value()) ? epsilon.value() : xtl::missing<double>(),
         t_msk,
         m_cdt,
-        bootstrap,
+        (bootstrap.has_value())
+        ? bootstrap.value()
+        : xtl::missing<std::unordered_map<std::string, int>>(),
+        dts
+    );
+}
+
+auto evalp(
+    const xt::pytensor<double, 2>& q_obs,
+    const xt::pytensor<double, 4>& q_prd,
+    const std::vector<std::string>& metrics,
+    const xt::pytensor<double, 2>& q_thr,
+    const xt::pytensor<bool, 4>& t_msk,
+    const xt::pytensor<std::array<char, 32>, 2>& m_cdt,
+    std::optional<std::unordered_map<std::string, int>> bootstrap,
+    const std::vector<std::string>& dts
+)
+{
+    return evalhyd::evalp(
+        q_obs,
+        q_prd,
+        metrics,
+        q_thr,
+        t_msk,
+        m_cdt,
+        (bootstrap.has_value())
+        ? bootstrap.value()
+        : xtl::missing<std::unordered_map<std::string, int>>(),
         dts
     );
 }
 
 // Python Module and Docstrings
-PYBIND11_MODULE(evalhyd, m)
+PYBIND11_MODULE(core, m)
 {
     xt::import_numpy();
 
-    m.doc() = "Utility for evaluation of streamflow predictions";
+    m.doc() = "Python bindings for the C++ core of evalhyd";
 
     // deterministic evaluation
     m.def(
-        "evald",
-        &evald_1d,
-        "Function to evaluate deterministic streamflow predictions (1D)",
-        py::arg("q_obs"), py::arg("q_prd"), py::arg("metrics"),
-        py::arg("transform") = "none",
-        py::arg("exponent") = 1,
-        py::arg("epsilon") = -9,
-        py::arg("t_msk") = xt::pytensor<bool, 2>({0}),
-        py::arg("m_cdt") = xt::pytensor<std::array<char, 32>, 1>({}),
-        py::arg("bootstrap") =
-            py::dict("n_samples"_a=-9, "len_sample"_a=-9, "summary"_a=0),
-        py::arg("dts") = py::list()
-    );
-
-    m.def(
-        "evald",
-        &evalhyd::evald<xt::pytensor<double, 2>, xt::pytensor<bool, 2>>,
-        "Function to evaluate deterministic streamflow predictions (2D)",
-        py::arg("q_obs"), py::arg("q_prd"), py::arg("metrics"),
-        py::arg("transform") = "none",
-        py::arg("exponent") = 1,
-        py::arg("epsilon") = -9,
+        "_evald",
+        &evald,
+        "Function to evaluate determinist streamflow predictions (2D)",
+        py::arg("q_obs"),
+        py::arg("q_prd"),
+        py::arg("metrics"),
+        py::arg("transform") = py::none(),
+        py::arg("exponent") = py::none(),
+        py::arg("epsilon") = py::none(),
         py::arg("t_msk") = xt::pytensor<bool, 2>({0}),
         py::arg("m_cdt") = xt::pytensor<std::array<char, 32>, 1>({}),
-        py::arg("bootstrap") =
-            py::dict("n_samples"_a=-9, "len_sample"_a=-9, "summary"_a=0),
+        py::arg("bootstrap") = py::none(),
         py::arg("dts") = py::list()
     );
 
     // probabilistic evaluation
     m.def(
-        "evalp",
-        &evalhyd::evalp<xt::pytensor<double, 2>, xt::pytensor<double, 4>, xt::pytensor<bool, 4>>,
-        "Function to evaluate probabilistic streamflow predictions",
-        py::arg("q_obs"), py::arg("q_prd"), py::arg("metrics"),
+        "_evalp",
+        &evalp,
+        "Function to evaluate probabilist streamflow predictions",
+        py::arg("q_obs"),
+        py::arg("q_prd"),
+        py::arg("metrics"),
         py::arg("q_thr") = xt::pytensor<double, 2>({0}),
         py::arg("t_msk") = xt::pytensor<bool, 4>({0}),
         py::arg("m_cdt") = xt::pytensor<std::array<char, 32>, 2>({0}),
-        py::arg("bootstrap") =
-            py::dict("n_samples"_a=-9, "len_sample"_a=-9, "summary"_a=0),
+        py::arg("bootstrap") = py::none(),
         py::arg("dts") = py::list()
     );
 

setup.py

@@ -61,8 +61,8 @@ for dep, version, url in deps:
 # configure Python extension
 ext_modules = [
     Pybind11Extension(
-        "evalhyd",
-        ['src/evalhyd-python.cpp'],
+        "evalhyd.core",
+        ['evalhyd/src/core.cpp'],
         include_dirs=[
             numpy.get_include(),
             os.path.join(sys.prefix, 'include'),
@@ -83,8 +83,9 @@ setup(
     url='https://gitlab.irstea.fr/hycar-hydro/evalhyd/evalhyd-python',
     description='Python bindings for EvalHyd',
     long_description='An evaluator for streamflow predictions.',
+    packages=["evalhyd"],
     ext_modules=ext_modules,
     cmdclass={'build_ext': build_ext},
-    extras_require={"tests": "numpy>=1.16"},
+    extras_require={'tests': 'numpy>=1.16'},
     zip_safe=False,
 )