class rework Timstat2.cc

src/Timstat2.cc

@@ -263,155 +263,190 @@ rmsd_compute(size_t gridsize, double missval, const Varray<size_t> &nofvals, Var
   return nmiss;
 }
 
-void *
-Timstat2(void *process)
+class ModuleTimstat2
 {
   CdiDateTime vDateTime{};
 
-  cdo_initialize(process);
-
-  // clang-format off
-  cdo_operator_add("timcor",   FieldFunc_Cor,   5, nullptr);
-  cdo_operator_add("timcovar", FieldFunc_Covar, 3, nullptr);
-  cdo_operator_add("timrmsd",  FieldFunc_Rmsd,  1, nullptr);
-  // clang-format on
-
-  auto operatorID = cdo_operator_id();
-  auto operfunc = cdo_operator_f1(operatorID);
-  auto nwork = cdo_operator_f2(operatorID);
-  auto timeIsConst = (operfunc == FieldFunc_Rmsd);
+  CdoStreamID streamID1;
+  int taxisID1;
+  CdoStreamID streamID2;
 
-  operator_check_argc(0);
+  CdoStreamID streamID3;
+  int taxisID3;
 
-  auto streamID1 = cdo_open_read(0);
-  auto streamID2 = cdo_open_read(1);
+  int nrecs1;
 
-  auto vlistID1 = cdo_stream_inq_vlist(streamID1);
-  auto vlistID2 = cdo_stream_inq_vlist(streamID2);
-  auto vlistID3 = vlistDuplicate(vlistID1);
-
-  vlist_compare(vlistID1, vlistID2, CMP_ALL);
+  int operfunc;
 
   VarList varList1, varList2;
-  varListInit(varList1, vlistID1);
-  varListInit(varList2, vlistID2);
-
-  auto nvars = vlistNvars(vlistID1);
-  auto nrecs1 = vlistNrecs(vlistID1);
-  std::vector<int> recVarID(nrecs1), recLevelID(nrecs1);
-
-  auto taxisID1 = vlistInqTaxis(vlistID1);
-  // auto taxisID2 = vlistInqTaxis(vlistID2);
-  auto taxisID3 = taxisDuplicate(taxisID1);
-
-  if (timeIsConst)
-    for (int varID = 0; varID < nvars; ++varID) vlistDefVarTimetype(vlistID3, varID, TIME_CONSTANT);
-
+  std::vector<int> recVarID, recLevelID;
   Field field1, field2;
+  Varray4D<double> work;
+  Varray3D<size_t> nofvals;
 
-  vlistDefTaxis(vlistID3, taxisID3);
-  auto streamID3 = cdo_open_write(2);
-  cdo_def_vlist(streamID3, vlistID3);
-
-  Varray4D<double> work(nvars);
-  Varray3D<size_t> nofvals(nvars);
-
-  for (int varID = 0; varID < nvars; ++varID)
-    {
-      auto gridsize = varList1[varID].gridsize;
-      auto nlevels = varList1[varID].nlevels;
-
-      work[varID].resize(nlevels);
-      nofvals[varID].resize(nlevels);
-
-      for (int levelID = 0; levelID < nlevels; ++levelID)
-        {
-          nofvals[varID][levelID].resize(gridsize, 0);
-          work[varID][levelID].resize(nwork);
-          for (int i = 0; i < nwork; ++i) work[varID][levelID][i].resize(gridsize, 0.0);
-        }
-    }
-
-  int tsID = 0;
-  while (true)
-    {
-      auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
-      if (nrecs == 0) break;
-
-      vDateTime = taxisInqVdatetime(taxisID1);
-
-      auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
-      if (nrecs != nrecs2) cdo_warning("Input streams have different number of records!");
-
-      for (int recID = 0; recID < nrecs; ++recID)
-        {
-          int varID, levelID;
-          cdo_inq_record(streamID1, &varID, &levelID);
-          cdo_inq_record(streamID2, &varID, &levelID);
-
-          field1.init(varList1[varID]);
-          field2.init(varList2[varID]);
-
-          if (tsID == 0)
-            {
-              recVarID[recID] = varID;
-              recLevelID[recID] = levelID;
-            }
-
-          auto gridsize = varList1[varID].gridsize;
-
-          cdo_read_record(streamID1, field1);
-          cdo_read_record(streamID2, field2);
-
-          auto &rwork = work[varID][levelID];
-          auto &rnofvals = nofvals[varID][levelID];
-
-          if (operfunc == FieldFunc_Cor)
-            {
-              correlation_init(gridsize, field1, field2, rnofvals, rwork[0], rwork[1], rwork[2], rwork[3], rwork[4]);
-            }
-          else if (operfunc == FieldFunc_Covar)
-            {
-              covariance_init(gridsize, field1, field2, rnofvals, rwork[0], rwork[1], rwork[2]);
-            }
-          else if (operfunc == FieldFunc_Rmsd) { rmsd_init(gridsize, field1, field2, rnofvals, rwork[0]); }
-        }
-
-      tsID++;
-    }
-
-  tsID = 0;
-  taxisDefVdatetime(taxisID3, vDateTime);
-  cdo_def_timestep(streamID3, tsID);
+public:
+  void
+  init(void *process)
+  {
+    cdo_initialize(process);
 
-  for (int recID = 0; recID < nrecs1; ++recID)
-    {
-      auto varID = recVarID[recID];
-      auto levelID = recLevelID[recID];
-
-      auto gridsize = varList1[varID].gridsize;
-      auto missval = varList1[varID].missval;
-
-      auto &rwork = work[varID][levelID];
-      const auto &rnofvals = nofvals[varID][levelID];
-
-      size_t nmiss = 0;
-      if (operfunc == FieldFunc_Cor)
-        {
-          nmiss = correlation(gridsize, missval, rnofvals, rwork[0], rwork[1], rwork[2], rwork[3], rwork[4]);
-        }
-      else if (operfunc == FieldFunc_Covar) { nmiss = covariance(gridsize, missval, rnofvals, rwork[0], rwork[1], rwork[2]); }
-      else if (operfunc == FieldFunc_Rmsd) { nmiss = rmsd_compute(gridsize, missval, rnofvals, rwork[0]); }
-
-      cdo_def_record(streamID3, varID, levelID);
-      cdo_write_record(streamID3, rwork[0].data(), nmiss);
-    }
+    // clang-format off
+  cdo_operator_add("timcor",   FieldFunc_Cor,   5, nullptr);
+  cdo_operator_add("timcovar", FieldFunc_Covar, 3, nullptr);
+  cdo_operator_add("timrmsd",  FieldFunc_Rmsd,  1, nullptr);
+    // clang-format on
+
+    auto operatorID = cdo_operator_id();
+    operfunc = cdo_operator_f1(operatorID);
+    auto nwork = cdo_operator_f2(operatorID);
+    auto timeIsConst = (operfunc == FieldFunc_Rmsd);
+
+    operator_check_argc(0);
+
+    streamID1 = cdo_open_read(0);
+    streamID2 = cdo_open_read(1);
+
+    auto vlistID1 = cdo_stream_inq_vlist(streamID1);
+    auto vlistID2 = cdo_stream_inq_vlist(streamID2);
+    auto vlistID3 = vlistDuplicate(vlistID1);
+
+    vlist_compare(vlistID1, vlistID2, CMP_ALL);
+
+    varListInit(varList1, vlistID1);
+    varListInit(varList2, vlistID2);
+
+    auto nvars = vlistNvars(vlistID1);
+    nrecs1 = vlistNrecs(vlistID1);
+    recVarID = std::vector<int>(nrecs1);
+    recLevelID = std::vector<int>(nrecs1);
+
+    taxisID1 = vlistInqTaxis(vlistID1);
+    // auto taxisID2 = vlistInqTaxis(vlistID2);
+    taxisID3 = taxisDuplicate(taxisID1);
+
+    if (timeIsConst)
+      for (int varID = 0; varID < nvars; ++varID) vlistDefVarTimetype(vlistID3, varID, TIME_CONSTANT);
+
+    vlistDefTaxis(vlistID3, taxisID3);
+    streamID3 = cdo_open_write(2);
+    cdo_def_vlist(streamID3, vlistID3);
+
+    work = Varray4D<double>(nvars);
+    nofvals = Varray3D<size_t>(nvars);
+
+    for (int varID = 0; varID < nvars; ++varID)
+      {
+        auto gridsize = varList1[varID].gridsize;
+        auto nlevels = varList1[varID].nlevels;
+
+        work[varID].resize(nlevels);
+        nofvals[varID].resize(nlevels);
+
+        for (int levelID = 0; levelID < nlevels; ++levelID)
+          {
+            nofvals[varID][levelID].resize(gridsize, 0);
+            work[varID][levelID].resize(nwork);
+            for (int i = 0; i < nwork; ++i) work[varID][levelID][i].resize(gridsize, 0.0);
+          }
+      }
+  }
+  void
+  run()
+  {
+    int tsID = 0;
+    while (true)
+      {
+        auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
+        if (nrecs == 0) break;
+
+        vDateTime = taxisInqVdatetime(taxisID1);
+
+        auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
+        if (nrecs != nrecs2) cdo_warning("Input streams have different number of records!");
+
+        for (int recID = 0; recID < nrecs; ++recID)
+          {
+            int varID, levelID;
+            cdo_inq_record(streamID1, &varID, &levelID);
+            cdo_inq_record(streamID2, &varID, &levelID);
+
+            field1.init(varList1[varID]);
+            field2.init(varList2[varID]);
+
+            if (tsID == 0)
+              {
+                recVarID[recID] = varID;
+                recLevelID[recID] = levelID;
+              }
+
+            auto gridsize = varList1[varID].gridsize;
+
+            cdo_read_record(streamID1, field1);
+            cdo_read_record(streamID2, field2);
+
+            auto &rwork = work[varID][levelID];
+            auto &rnofvals = nofvals[varID][levelID];
+
+            if (operfunc == FieldFunc_Cor)
+              {
+                correlation_init(gridsize, field1, field2, rnofvals, rwork[0], rwork[1], rwork[2], rwork[3], rwork[4]);
+              }
+            else if (operfunc == FieldFunc_Covar)
+              {
+                covariance_init(gridsize, field1, field2, rnofvals, rwork[0], rwork[1], rwork[2]);
+              }
+            else if (operfunc == FieldFunc_Rmsd) { rmsd_init(gridsize, field1, field2, rnofvals, rwork[0]); }
+          }
+
+        tsID++;
+      }
+
+    tsID = 0;
+    taxisDefVdatetime(taxisID3, vDateTime);
+    cdo_def_timestep(streamID3, tsID);
+
+    for (int recID = 0; recID < nrecs1; ++recID)
+      {
+        auto varID = recVarID[recID];
+        auto levelID = recLevelID[recID];
+
+        auto gridsize = varList1[varID].gridsize;
+        auto missval = varList1[varID].missval;
+
+        auto &rwork = work[varID][levelID];
+        const auto &rnofvals = nofvals[varID][levelID];
+
+        size_t nmiss = 0;
+        if (operfunc == FieldFunc_Cor)
+          {
+            nmiss = correlation(gridsize, missval, rnofvals, rwork[0], rwork[1], rwork[2], rwork[3], rwork[4]);
+          }
+        else if (operfunc == FieldFunc_Covar) { nmiss = covariance(gridsize, missval, rnofvals, rwork[0], rwork[1], rwork[2]); }
+        else if (operfunc == FieldFunc_Rmsd) { nmiss = rmsd_compute(gridsize, missval, rnofvals, rwork[0]); }
+
+        cdo_def_record(streamID3, varID, levelID);
+        cdo_write_record(streamID3, rwork[0].data(), nmiss);
+      }
+  }
+  void
+  close()
+  {
+    cdo_stream_close(streamID3);
+    cdo_stream_close(streamID2);
+    cdo_stream_close(streamID1);
+
+    cdo_finish();
+  }
+};
 
-  cdo_stream_close(streamID3);
-  cdo_stream_close(streamID2);
-  cdo_stream_close(streamID1);
+void *
+Timstat2(void *process)
+{
+  ModuleTimstat2 timstat2;
 
-  cdo_finish();
+  timstat2.init(process);
+  timstat2.run();
+  timstat2.close();
 
   return nullptr;
 }