tweaks

src/mr_util/infra/mav.h

@@ -105,13 +105,6 @@ class fmav_info
     stride_t str;
     size_t sz;
 
-    static size_t prod(const shape_t &shape)
-      {
-      size_t res=1;
-      for (auto sz: shape)
-        res*=sz;
-      return res;
-      }
     static stride_t shape2stride(const shape_t &shp)
       {
       auto ndim = shp.size();
@@ -128,7 +121,7 @@ class fmav_info
 
   public:
     fmav_info(const shape_t &shape_, const stride_t &stride_)
-      : shp(shape_), str(stride_), sz(prod(shp))
+      : shp(shape_), str(stride_), sz(reduce(shp.begin(),shp.end(),size_t(1),multiplies<>()))
       {
       MR_assert(shp.size()>0, "at least 1D required");
       MR_assert(shp.size()==str.size(), "dimensions mismatch");
@@ -175,13 +168,6 @@ template<size_t ndim> class mav_info
     stride_t str;
     size_t sz;
 
-    static size_t prod(const shape_t &shape)
-      {
-      size_t res=1;
-      for (auto sz: shape)
-        res*=sz;
-      return res;
-      }
     static stride_t shape2stride(const shape_t &shp)
       {
       stride_t res;
@@ -197,7 +183,7 @@ template<size_t ndim> class mav_info
 
   public:
     mav_info(const shape_t &shape_, const stride_t &stride_)
-      : shp(shape_), str(stride_), sz(prod(shp)) {}
+      : shp(shape_), str(stride_), sz(reduce(shp.begin(),shp.end(),size_t(1),multiplies<>())) {}
     mav_info(const shape_t &shape_)
       : mav_info(shape_, shape2stride(shape_)) {}
     size_t size() const { return sz; }

src/mr_util/math/es_kernel.h

@@ -42,14 +42,13 @@ class ES_Kernel
   {
   private:
     double beta;
-    float fbeta;
     int p;
     vector<double> x, wgt, psi;
     size_t supp;
 
   public:
     ES_Kernel(size_t supp_, double ofactor, size_t nthreads)
-      : beta(get_beta(supp_,ofactor)*supp_), fbeta(float(beta)),
+      : beta(get_beta(supp_,ofactor)*supp_),
         p(int(1.5*supp_+2)), supp(supp_)
       {
       GL_Integrator integ(2*p,nthreads);
@@ -62,10 +61,12 @@ class ES_Kernel
     ES_Kernel(size_t supp_, size_t nthreads)
       : ES_Kernel(supp_, 2., nthreads){}
 
-    double operator()(double v) const
-      { return (v*v>1.) ? 0. : exp(beta*(std::sqrt(1.-v*v)-1.)); }
-    float operator()(float v) const
-      { return (v*v>1.f) ? 0.f : exp(fbeta*(std::sqrt(1.f-v*v)-1.f)); }
+    template<typename T> T operator()(T v) const
+      {
+      auto tmp = (1-v)*(1+v);
+      auto tmp2 = tmp>=0;
+      return tmp2*exp(T(beta)*(sqrt(tmp*tmp2)-1));
+      }
     /* Compute correction factors for the ES gridding kernel
        This implementation follows eqs. (3.8) to (3.10) of Barnett et al. 2018 */
     double corfac(double v) const