nse.hpp

#include <vector>
#include <xtensor/xexpression.hpp>
#include <xtensor/xmath.hpp>

template <class A>
A nse(const xt::xexpression<A>& sim,
      const xt::xexpression<A>& obs,
      int axis = 0);

template <class A>
A nse(const xt::xexpression<A>& sim,
      const xt::xexpression<A>& obs,
      int axis)
{
    const A& q_sim = sim.derived_cast();
    const A& q_obs = obs.derived_cast();

    A q_avg = xt::mean(q_obs, xt::keep_dims);

    A f_num = xt::sum(xt::square(q_obs - q_sim), std::vector<int> { axis }, xt::keep_dims);
    A f_den = xt::sum(xt::square(q_obs - q_avg), xt::keep_dims);

    return 1 - (f_num / f_den);
}