Commit bca92a97 authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

avoid computing errors twice

by storing errors before computing quadratic errors
Showing with 100 additions and 65 deletions
+100 -65
......@@ -17,22 +17,16 @@ namespace evalhyd
{
/// Compute the Pearson correlation coefficient.
///
/// \param q_obs
/// Streamflow observations.
/// shape: (series, time)
/// \param q_prd
/// Streamflow predictions.
/// shape: (series, time)
/// \param mean_obs
/// Mean observed streamflow.
/// shape: (subsets, samples, series, 1)
/// \param mean_prd
/// Mean predicted streamflow.
/// shape: (subsets, samples, series, 1)
/// \param quad_obs
/// \param err_obs
/// Errors between observations and mean observation.
/// shape: (subsets, samples, series, time)
/// \param err_prd
/// Errors between predictions and mean prediction.
/// shape: (subsets, samples, series, time)
/// \param quad_err_obs
/// Quadratic errors between observations and mean observation.
/// shape: (subsets, samples, series, time)
/// \param quad_prd
/// \param quad_err_prd
/// Quadratic errors between predictions and mean prediction.
/// shape: (subsets, samples, series, time)
/// \param t_msk
......@@ -50,14 +44,11 @@ namespace evalhyd
/// \return
/// Pearson correlation coefficients.
/// shape: (subsets, samples, series)
template <class XD2>
inline xt::xtensor<double, 3> calc_r_pearson(
const XD2& q_obs,
const XD2& q_prd,
const xt::xtensor<double, 4>& mean_obs,
const xt::xtensor<double, 4>& mean_prd,
const xt::xtensor<double, 4>& quad_obs,
const xt::xtensor<double, 4>& quad_prd,
const xt::xtensor<double, 4>& err_obs,
const xt::xtensor<double, 4>& err_prd,
const xt::xtensor<double, 4>& quad_err_obs,
const xt::xtensor<double, 4>& quad_err_prd,
const xt::xtensor<bool, 3>& t_msk,
const std::vector<xt::xkeep_slice<int>>& b_exp,
std::size_t n_srs,
......@@ -72,26 +63,15 @@ namespace evalhyd
for (std::size_t m = 0; m < n_msk; m++)
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
auto prd_masked = xt::where(xt::view(t_msk, xt::all(), m),
q_prd, NAN);
auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
q_obs, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
{
auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
auto r_num = xt::nansum(
(prd - xt::view(mean_prd, m, e))
* (obs - xt::view(mean_obs, m, e)),
-1
);
auto prd = xt::view(err_prd, m, e, xt::all(), b_exp[e]);
auto obs = xt::view(err_obs, m, e, xt::all(), b_exp[e]);
auto r_num = xt::nansum(prd * obs, -1);
auto prd2 = xt::view(quad_prd, m, e, xt::all(), b_exp[e]);
auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
auto prd2 = xt::view(quad_err_prd, m, e, xt::all(), b_exp[e]);
auto obs2 = xt::view(quad_err_obs, m, e, xt::all(), b_exp[e]);
auto r_den = xt::sqrt(
xt::nansum(prd2, -1) * xt::nansum(obs2, -1)
);
......@@ -269,7 +249,7 @@ namespace evalhyd
/// \param quad_err
/// Quadratic errors between observations and predictions.
/// shape: (series, time)
/// \param quad_obs
/// \param quad_err_obs
/// Quadratic errors between observations and mean observation.
/// shape: (subsets, samples, series, time)
/// \param t_msk
......@@ -289,7 +269,7 @@ namespace evalhyd
/// shape: (series, subsets, samples)
inline xt::xtensor<double, 3> calc_NSE(
const xt::xtensor<double, 2>& quad_err,
const xt::xtensor<double, 4>& quad_obs,
const xt::xtensor<double, 4>& quad_err_obs,
const xt::xtensor<bool, 3>& t_msk,
const std::vector<xt::xkeep_slice<int>>& b_exp,
std::size_t n_srs,
......@@ -314,7 +294,7 @@ namespace evalhyd
auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]);
xt::xtensor<double, 1> f_num =
xt::nansum(err2, -1);
auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
auto obs2 = xt::view(quad_err_obs, m, e, xt::all(), b_exp[e]);
xt::xtensor<double, 1> f_den =
xt::nansum(obs2, -1);
......
......@@ -171,7 +171,7 @@ namespace evalhyd
return xt::square(err);
}
/// Compute the quadratic error between observations and mean observation.
/// Compute the error between observations and mean observation.
///
/// \param q_obs
/// Streamflow observations.
......@@ -191,10 +191,10 @@ namespace evalhyd
/// \param n_exp
/// Number of bootstrap samples.
/// \return
/// Quadratic errors between observations and mean observation.
/// Errors between observations and mean observation.
/// shape: (subsets, samples, series, time)
template <class XD2>
inline xt::xtensor<double, 4> calc_quad_obs(
inline xt::xtensor<double, 4> calc_err_obs(
const XD2& q_obs,
const xt::xtensor<double, 4>& mean_obs,
const xt::xtensor<bool, 3>& t_msk,
......@@ -204,7 +204,7 @@ namespace evalhyd
std::size_t n_exp
)
{
xt::xtensor<double, 4> quad_obs =
xt::xtensor<double, 4> err_obs =
xt::zeros<double>({n_msk, n_exp, n_srs, n_tim});
for (std::size_t m = 0; m < n_msk; m++)
......@@ -216,16 +216,31 @@ namespace evalhyd
for (std::size_t e = 0; e < n_exp; e++)
{
xt::view(quad_obs, m, e) = xt::square(
xt::view(err_obs, m, e) = (
obs_masked - xt::view(mean_obs, m, e)
);
}
}
return quad_obs;
return err_obs;
}
/// Compute the quadratic error between predictions and mean prediction.
/// Compute the quadratic error between observations and mean observation.
///
/// \param err_obs
/// Errors between observations and mean observation.
/// shape: (subsets, samples, series, time)
/// \return
/// Quadratic errors between observations and mean observation.
/// shape: (subsets, samples, series, time)
inline xt::xtensor<double, 4> calc_quad_err_obs(
const xt::xtensor<double, 4>& err_obs
)
{
return xt::square(err_obs);
}
/// Compute the error between predictions and mean prediction.
///
/// \param q_prd
/// Streamflow predictions.
......@@ -245,10 +260,10 @@ namespace evalhyd
/// \param n_exp
/// Number of bootstrap samples.
/// \return
/// Quadratic errors between predictions and mean prediction.
/// Errors between predictions and mean prediction.
/// shape: (subsets, samples, series, time)
template <class XD2>
inline xt::xtensor<double, 4> calc_quad_prd(
inline xt::xtensor<double, 4> calc_err_prd(
const XD2& q_prd,
const xt::xtensor<double, 4>& mean_prd,
const xt::xtensor<bool, 3>& t_msk,
......@@ -258,7 +273,7 @@ namespace evalhyd
std::size_t n_exp
)
{
xt::xtensor<double, 4> quad_prd =
xt::xtensor<double, 4> quad_err_prd =
xt::zeros<double>({n_msk, n_exp, n_srs, n_tim});
for (std::size_t m = 0; m < n_msk; m++)
......@@ -270,13 +285,28 @@ namespace evalhyd
for (std::size_t e = 0; e < n_exp; e++)
{
xt::view(quad_prd, m, e) = xt::square(
xt::view(quad_err_prd, m, e) = (
prd_masked - xt::view(mean_prd, m, e)
);
}
}
return quad_prd;
return quad_err_prd;
}
/// Compute the quadratic error between predictions and mean prediction.
///
/// \param err_obs
/// Errors between predictions and mean prediction.
/// shape: (subsets, samples, series, time)
/// \return
/// Quadratic errors between predictions and mean prediction.
/// shape: (subsets, samples, series, time)
inline xt::xtensor<double, 4> calc_quad_err_prd(
const xt::xtensor<double, 4>& err_prd
)
{
return xt::square(err_prd);
}
}
......
......@@ -42,8 +42,10 @@ namespace evalhyd
xtl::xoptional<xt::xtensor<double, 2>, bool> err;
xtl::xoptional<xt::xtensor<double, 2>, bool> abs_err;
xtl::xoptional<xt::xtensor<double, 2>, bool> quad_err;
xtl::xoptional<xt::xtensor<double, 4>, bool> quad_obs;
xtl::xoptional<xt::xtensor<double, 4>, bool> quad_prd;
xtl::xoptional<xt::xtensor<double, 4>, bool> err_obs;
xtl::xoptional<xt::xtensor<double, 4>, bool> quad_err_obs;
xtl::xoptional<xt::xtensor<double, 4>, bool> err_prd;
xtl::xoptional<xt::xtensor<double, 4>, bool> quad_err_prd;
xtl::xoptional<xt::xtensor<double, 3>, bool> r_pearson;
xtl::xoptional<xt::xtensor<double, 3>, bool> alpha;
xtl::xoptional<xt::xtensor<double, 3>, bool> gamma;
......@@ -116,28 +118,50 @@ namespace evalhyd
return quad_err.value();
};
xt::xtensor<double, 4> get_quad_obs()
xt::xtensor<double, 4> get_err_obs()
{
if (!quad_obs.has_value())
if (!err_obs.has_value())
{
quad_obs = elements::calc_quad_obs(
err_obs = elements::calc_err_obs(
q_obs, get_mean_obs(), t_msk,
n_srs, n_tim, n_msk, n_exp
);
}
return quad_obs.value();
return err_obs.value();
};
xt::xtensor<double, 4> get_quad_prd()
xt::xtensor<double, 4> get_quad_err_obs()
{
if (!quad_prd.has_value())
if (!quad_err_obs.has_value())
{
quad_prd = elements::calc_quad_prd(
quad_err_obs = elements::calc_quad_err_obs(
get_err_obs()
);
}
return quad_err_obs.value();
};
xt::xtensor<double, 4> get_err_prd()
{
if (!err_prd.has_value())
{
err_prd = elements::calc_err_prd(
q_prd, get_mean_prd(), t_msk,
n_srs, n_tim, n_msk, n_exp
);
}
return quad_prd.value();
return err_prd.value();
};
xt::xtensor<double, 4> get_quad_err_prd()
{
if (!quad_err_prd.has_value())
{
quad_err_prd = elements::calc_quad_err_prd(
get_err_prd()
);
}
return quad_err_prd.value();
};
xt::xtensor<double, 3> get_r_pearson()
......@@ -145,8 +169,9 @@ namespace evalhyd
if (!r_pearson.has_value())
{
r_pearson = elements::calc_r_pearson(
q_obs, q_prd, get_mean_obs(), get_mean_prd(),
get_quad_obs(), get_quad_prd(), t_msk, b_exp,
get_err_obs(), get_err_prd(),
get_quad_err_obs(), get_quad_err_prd(),
t_msk, b_exp,
n_srs, n_msk, n_exp
);
}
......@@ -272,7 +297,7 @@ namespace evalhyd
if (!NSE.has_value())
{
NSE = metrics::calc_NSE(
get_quad_err(), get_quad_obs(), t_msk, b_exp,
get_quad_err(), get_quad_err_obs(), t_msk, b_exp,
n_srs, n_msk, n_exp
);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment