From 14118a1d5337385b7d397cdfb4f861691217703f Mon Sep 17 00:00:00 2001
From: Michael Uleysky <uleysky@poi.dvo.ru>
Date: Wed, 15 Mar 2023 13:20:19 +1000
Subject: [PATCH] Making averaging more generic

---
 include/NEMO.h         | 10 ++++-----
 include/basedata.h     | 43 ++++++++++++++++++++++++++++++++++++-
 include/simple2ddata.h | 48 ++++--------------------------------------
 3 files changed, 50 insertions(+), 51 deletions(-)

diff --git a/include/NEMO.h b/include/NEMO.h
index 2050d83..504096b 100644
--- a/include/NEMO.h
+++ b/include/NEMO.h
@@ -92,7 +92,7 @@ class NEMOData
   auto unit = f.A<MString>(name, "units");
   if(unit && unit.Get() == "m s-1") unitmul = 100.0;
 
-  Data data((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1), ye - yb + 1, Lonb(), Latb());
+  Data data((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1), ye - yb + 1, Lonb(), Latb(), 1.0 / 12.0, 1.0 / 12.0);
 
   if(xb < xe)
   {
@@ -154,17 +154,16 @@ class NEMOData
 
  Data ReadVar(const MString& name, const std::vector<size_t>& tindex) const
  {
-  Data out;
+  Averager<Data> out;
   if(tindex.size() == 0 || !isOk()) return out;
 
-  std::vector<size_t> count;
   for(size_t i = 0; i < tindex.size(); i++)
   {
    Data dat = ReadVar(name, tindex[i]);
    if(!dat) return Data();
-   out.Add(dat, count);
+   out.Add(dat);
   }
-  out.Div(count);
+  out.Div();
   return out;
  }
 
@@ -391,4 +390,3 @@ template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SAL>(const std
 template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::MLD>(const std::vector<size_t>& tindex) const { return ReadVar("mlotst", tindex); }
 template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SSH>(const std::vector<size_t>& tindex) const { return ReadVar("zos", tindex); }
 template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::W>(const std::vector<size_t>& tindex) const { return ReadVar("wo", tindex); }
-
diff --git a/include/basedata.h b/include/basedata.h
index 998dcb8..88fb9cf 100644
--- a/include/basedata.h
+++ b/include/basedata.h
@@ -24,6 +24,47 @@ class BaseData
 
  static real Fillval() { return fillval; }
 
- explicit operator bool() const { return data.size() != 0; }
+ explicit operator bool() const { return N() != 0; }
 };
 
+template<class Data> class Averager: public Data
+{
+ std::vector<size_t> count;
+
+ public:
+ void Add(const Data& d)
+ {
+  if(!d) return;
+  if(!*this)
+  {
+   *static_cast<Data*>(this) = d;
+   count.resize(Data::N());
+   for(size_t i = 0; i < Data::N(); i++) count[i] = (Data::V(i) == Data::Fillval()) ? 0 : 1;
+   return;
+  }
+
+  if(Data::N() != d.N()) return;
+  for(size_t i = 0; i < Data::N(); i++)
+   if(d(i) != Data::Fillval())
+   {
+    if(Data::V(i) == Data::Fillval())
+    {
+     Data::V(i) = d(i);
+     count[i]   = 1;
+    }
+    else
+    {
+     Data::V(i) += d(i);
+     count[i]++;
+    }
+   }
+ }
+
+ Averager& Div()
+ {
+  if(!*this) return *this;
+  for(size_t i = 0; i < Data::N(); i++)
+   if(count[i] != 0) Data::V(i) /= count[i];
+  return *this;
+ }
+};
diff --git a/include/simple2ddata.h b/include/simple2ddata.h
index a7a74e4..2c17c82 100644
--- a/include/simple2ddata.h
+++ b/include/simple2ddata.h
@@ -5,11 +5,12 @@ class Simple2DData: public BaseData
 {
  real   x0 = 0.0, y0 = 0.0;
  size_t nx = 0, ny = 0;
+ real   xstep = 0.0, ystep = 0.0;
 
  public:
  Simple2DData() = default;
 
- Simple2DData(size_t nx_, size_t ny_, real x0_, real y0_): BaseData(nx_ * ny_), x0(x0_), y0(y0_), nx(nx_), ny(ny_) {}
+ Simple2DData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_): BaseData(nx_ * ny_), x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_) {}
 
  const real& V(size_t i) const { return BaseData::V(i); }
  real&       V(size_t i) { return BaseData::V(i); }
@@ -26,50 +27,9 @@ class Simple2DData: public BaseData
  size_t Nx() const { return nx; }
  size_t Ny() const { return ny; }
 
- real Lon(size_t ix, size_t iy) const { return x0 + ix / 12.0; }
- real Lat(size_t ix, size_t iy) const { return y0 + iy / 12.0; }
+ real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return x0 + ix * xstep; }
+ real Lat([[maybe_unused]] size_t ix, size_t iy) const { return y0 + iy * ystep; }
 
  real Lon(size_t i) const { return Lon(i % nx, i / nx); }
  real Lat(size_t i) const { return Lat(i % nx, i / nx); }
-
- void InitCount(std::vector<size_t>& count) const
- {
-  count.resize(data.size());
-  for(size_t i = 0; i < N(); i++) count[i] = (V(i) == fillval) ? 0 : 1;
- }
-
- Simple2DData& Add(const Simple2DData& d, std::vector<size_t>& count)
- {
-  if(!d) return *this;
-  if(!*this)
-  {
-   *this = d;
-   InitCount(count);
-  }
-
-  if(N() != d.N() || N() != count.size()) return *this;
-  for(size_t i = 0; i < N(); i++)
-   if(d.V(i) != fillval)
-   {
-    if(V(i) == fillval)
-    {
-     V(i)     = d(i);
-     count[i] = 1;
-    }
-    else
-    {
-     V(i) += d(i);
-     count[i]++;
-    }
-   }
-  return *this;
- }
-
- Simple2DData& Div(const std::vector<size_t>& count)
- {
-  if(N() != count.size()) return *this;
-  for(size_t i = 0; i < N(); i++)
-   if(count[i] != 0) V(i) /= count[i];
-  return *this;
- }
 };