Detrending flags

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features
 
 project(
   EXP
-  VERSION "7.8.2"
+  VERSION "7.8.3"
   HOMEPAGE_URL https://github.com/EXP-code/EXP
   LANGUAGES C CXX Fortran)
 

expui/BiorthBasis.H

@@ -729,6 +729,7 @@ namespace BasisClasses
     int used;
 
     std::string cachename;
+    bool oldcache = false;
 
     using matT = std::vector<Eigen::MatrixXd>;
     using vecT = std::vector<Eigen::VectorXd>;

expui/BiorthBasis.cc

@@ -890,6 +890,7 @@ namespace BasisClasses
     "npca0",
     "nvtk",
     "cachename",
+    "oldcache",
     "eof_file",
     "override",
     "samplesz",
@@ -1061,6 +1062,7 @@ namespace BasisClasses
     ncylodd     = 9;
     ncylr       = 2000;
     cachename   = ".eof_cache_file";
+    oldcache    = false;
     Ignore      = false;
     deproject   = false;
 
@@ -1139,6 +1141,7 @@ namespace BasisClasses
       if (conf["ncylodd"   ])    ncylodd  = conf["ncylodd"   ].as<int>();
       if (conf["cachename" ])  cachename  = conf["cachename" ].as<std::string>();
       if (conf["eof_file"  ])  cachename  = conf["eof_file"  ].as<std::string>();
+      if (conf["oldcache"  ])   oldcache  = conf["oldcache"  ].as<bool>();
       if (conf["rnum"      ])       rnum  = conf["rnum"      ].as<int>();
       if (conf["pnum"      ])       pnum  = conf["pnum"      ].as<int>();
       if (conf["tnum"      ])       tnum  = conf["tnum"      ].as<int>();
@@ -1238,6 +1241,11 @@ namespace BasisClasses
     if (mlim>=0)  sl->set_mlim(mlim);
     if (EVEN_M)   sl->setEven(EVEN_M);
 
+    // Cache override for old Eigen cache
+    //
+    if (oldcache) sl->AllowOldCache();
+
+
     // Attempt to read EOF cache
     //
     if (sl->read_cache() == 0) {

expui/expMSSA.H

@@ -84,8 +84,8 @@ namespace MSSA
     //! Number of components in the reconstruction
     int ncomp;
 
-    //! Normalization values
-    double totVar, totPow;
+    //! Normalization options
+    bool totVar, totPow;
 
     //! Toggle for detrending
     bool useMean;

expui/expMSSA.cc

@@ -1522,8 +1522,8 @@ namespace MSSA {
 	HighFive::Group recon = file.createGroup("reconstruction");
 
 	recon.createAttribute<int>   ("ncomp",  HighFive::DataSpace::From(ncomp) ).write(ncomp);
-	recon.createAttribute<double>("totVar", HighFive::DataSpace::From(totVar)).write(totVar);
-	recon.createAttribute<double>("totPow", HighFive::DataSpace::From(totVar)).write(totPow);
+	recon.createAttribute<bool>("totVar", HighFive::DataSpace::From(totVar)).write(totVar);
+	recon.createAttribute<bool>("totPow", HighFive::DataSpace::From(totVar)).write(totPow);
 
 	for (int n=0; n<keylist.size(); n++) {
 	  std::ostringstream scnt;
@@ -1961,13 +1961,22 @@ namespace MSSA {
     // Detrending style (totVar and totPow are only useful, so far, for
     // noise computation)
     //
-    if (params["totVar"])
-      type = TrendType::totVar;
-    else if (params["totPow"])
+    if (params["totVar"]) totVar = params["totVar"].as<bool>();
+    if (params["totPow"]) totPow = params["totPow"].as<bool>();
+
+    if (totPow==true) {
+      if (totVar==true) {
+        std::cerr << "expMSSA: both totVar and totPow are set to true."
+                  << "Using totPow." << std::endl;
+        totVar = false;
+      }
       type = TrendType::totPow;
-    else
+    } else if (totVar==true) {
+      type = TrendType::totVar;
+    } else {
+      // if nothing set go default
       type = TrendType::perChannel;
-
+    }
 
     // Set the SVD strategy for mSSA
     //

exputil/EmpCylSL.cc

@@ -65,7 +65,6 @@ double   EmpCylSL::RMIN            = 0.001;
 double   EmpCylSL::RMAX            = 20.0;
 double   EmpCylSL::HFAC            = 0.2;
 double   EmpCylSL::PPOW            = 4.0;
-bool     EmpCylSL::NewCache        = true;
 bool     EmpCylSL::NewCoefs        = true;
 
 
@@ -264,7 +263,15 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename)
       if (ver.compare(Version))
 	throw std::runtime_error("EmpCylSL: version mismatch");
     } else {
-      throw std::runtime_error("EmpCylSL: outdated cache");
+      if (allow_old_cache) {
+	if (myid==0)
+	  std::cout << "---- EmpCylSL::ReadH5Cache: "
+		    << "using a workaround for a HighFive HDF5 wrapper bug. "
+                    << "this workaround will be removed in EXP 7.9.0. "
+		    << "Please consider rebuilding your cache if possible!"
+		    << std::endl;
+      }
+      else throw std::runtime_error("EmpCylSL: outdated cache");
     }
 
     NMAX    = getH5<int>("nmaxfid", file);
@@ -7258,6 +7265,22 @@ void EmpCylSL::WriteH5Cache()
 	    << "wrote <" << cachefile << ">" << std::endl;
 }
 
+// Repack order for Eigen3 workaround from old version
+//
+Eigen::MatrixXd repack(const Eigen::MatrixXd& in)
+{
+  Eigen::MatrixXd tmp = in.transpose();
+  double * trn = tmp.data();
+  Eigen::MatrixXd ret(in.rows(), in.cols());
+  for (int j=0, c=0; j<in.cols(); j++) {
+    for (int i=0; i<in.rows(); i++) {
+      ret(i, j) = trn[c++];
+    }
+  }
+  return ret;
+};
+
+
 bool EmpCylSL::ReadH5Cache()
 {
   try {
@@ -7307,11 +7330,20 @@ bool EmpCylSL::ReadH5Cache()
     if (file.hasAttribute("Version")) {
       if (not checkStr(Version, "Version"))  return false;
     } else {
-      if (myid==0)
-	std::cout << "---- EmpCylSL::ReadH5Cache: "
-		  << "recomputing cache for HighFive API change"
-		  << std::endl;
-      return false;
+      if (allow_old_cache) {
+	if (myid==0)
+	  std::cout << "---- EmpCylSL::ReadH5Cache: "
+		    << "using a workaround for a HighFive HDF5 wrapper bug. "
+        << "this workaround will be removed in EXP 7.9.0. "
+		    << "Please consider rebuilding your cache if possible!" << std::endl;
+      }
+      else {
+	if (myid==0)
+	  std::cout << "---- EmpCylSL::ReadH5Cache: "
+		    << "recomputing cache for HighFive API change"
+		    << std::endl;
+	return false;
+      }
     }
 
     if (not checkStr(geometry, "geometry"))  return false;
@@ -7388,6 +7420,13 @@ bool EmpCylSL::ReadH5Cache()
 	rforceC[m][n] = order.getDataSet("rforceC").read<Eigen::MatrixXd>();
 	zforceC[m][n] = order.getDataSet("zforceC").read<Eigen::MatrixXd>();
 	densC  [m][n] = order.getDataSet("densC")  .read<Eigen::MatrixXd>();
+
+	if (allow_old_cache) {
+	  potC   [m][n] = repack(potC   [m][n]);
+	  rforceC[m][n] = repack(rforceC[m][n]);
+	  zforceC[m][n] = repack(zforceC[m][n]);
+	  densC  [m][n] = repack(densC  [m][n]);
+	}
       }
     }
 
@@ -7411,6 +7450,13 @@ bool EmpCylSL::ReadH5Cache()
 	rforceS[m][n] = order.getDataSet("rforceS").read<Eigen::MatrixXd>();
 	zforceS[m][n] = order.getDataSet("zforceS").read<Eigen::MatrixXd>();
 	densS  [m][n] = order.getDataSet("densS")  .read<Eigen::MatrixXd>();
+
+	if (allow_old_cache) {
+	  potS   [m][n] = repack(potS   [m][n]);
+	  rforceS[m][n] = repack(rforceS[m][n]);
+	  zforceS[m][n] = repack(zforceS[m][n]);
+	  densS  [m][n] = repack(densS  [m][n]);;
+	}
       }
     }
 

include/EmpCylSL.H

@@ -340,6 +340,10 @@ protected:
   //! MPI is active
   bool use_mpi;
 
+  //! Allow older cache files with the Eigen3 row-major bug to be read
+  //! with the new Eigen3 API
+  bool allow_old_cache = false;
+
 public:
 
   /*! Enum listing the possible selection algorithms for coefficient
@@ -447,9 +451,6 @@ public:
   //! Power exponent for EOF conditioning density function
   static double PPOW;
 
-  //! Use YAML header in cache file
-  static bool NewCache;
-
   //! Use YAML header in coefficient file
   static bool NewCoefs;
 
@@ -808,6 +809,11 @@ public:
   static std::map<std::string, std::string>
   cacheInfo(const std::string& cachefile, bool verbose=true);
 
+  /** Cache versioning override.  Allow read of old cache files
+      without triggering a cache recomputation. */
+  void AllowOldCache(bool allow=true) { allow_old_cache = allow; }
+
+
   //@{
   //! Get potential grid interpolated entries
   void getPotSC(int m, int j, double R, double z, double& pC, double& pS);

utils/ICs/cylcache.cc

@@ -471,13 +471,6 @@ main(int ac, char **av)
     }
   }
 
-  // Enable new YAML cache header
-  //
-  if (vm.count("newcache")) {
-    EmpCylSL::NewCache = true;
-  }
-
-
   // Convert mtype string to lower case
   //
   std::transform(mtype.begin(), mtype.end(), mtype.begin(),

utils/ICs/initial.cc

@@ -663,12 +663,6 @@ main(int ac, char **av)
     }
   }
 
-  // Enable new YAML cache header
-  //
-  if (vm.count("newcache")) {
-    EmpCylSL::NewCache = true;
-  }
-
   if (vm.count("spline")) {
     SphericalModelTable::linear = 0;
   } else {