Slab force update

expui/BiorthBasis.H

@@ -1194,6 +1194,9 @@ namespace BasisClasses
     //! SLGridSlab mesh size
     int ngrid = 1000;
 
+    //! Coordinate map type for SLGridSlab
+    std::string cmap = "linear";
+
     //! Target model type for SLGridSlab
     std::string type = "isothermal";
 
@@ -1255,6 +1258,20 @@ namespace BasisClasses
     void computeAccel(double x, double y, double z,
 		      Eigen::Ref<Eigen::Vector3d> vstore);
 
+    //@{
+    //! Covariance structures.  First index is T, second is the
+    //! flattened 3-d k vector
+    std::vector<Eigen::VectorXcd> meanV;
+    std::vector<Eigen::MatrixXcd> covrV;
+    std::vector<Eigen::VectorXcd> dvarV;
+
+    void init_covariance();
+    void zero_covariance();
+    void writeCovarH5Params(HighFive::File& file);
+    int  sampT = 100;
+    bool diagcov = true;
+    //@}
+
   public:
 
     //! Constructor from YAML node
@@ -1311,6 +1328,58 @@ namespace BasisClasses
       return true;
     }
 
+    //! Write coefficient covariance data to an HDF5 file
+    virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0)
+    {
+      if (covarStore) {
+	SubsampleCovariance::CovarData elem;
+	std::get<0>(elem) = sampleCounts;
+	std::get<1>(elem) = sampleMasses;
+
+	int sampT = meanV.size();
+	int Nsize = meanV[0].size();
+
+	std::get<2>(elem) = SubsampleCovariance::CoefType(sampT, 1, Nsize);
+	std::get<3>(elem) = SubsampleCovariance::CovrType(sampT, 1, Nsize, Nsize);
+
+	for (int T=0; T<sampT; ++T) {
+	  for (int n1=0; n1<Nsize; n1++) {
+	    std::get<2>(elem)(T, 0, n1) = meanV[T](n1);
+	    for (int n2=0; n2<Nsize; n2++) {
+	      std::complex<double> cov = 0.0;
+	      if (diagcov) {
+		if (n1 == n2) cov = dvarV[T](n1);
+	      } else if (covar && T < covrV.size()) {
+		cov = covrV[T](n1, n2);
+	      }
+	      std::get<3>(elem)(T, 0, n1, n2) = cov;
+	    }
+	  }
+	}
+	covarStore->writeCoefCovariance(compname, runtag, elem, time);
+      } else {
+	throw std::runtime_error("Slab::writeCoefCovariance: covariance storage not initialized");
+      }
+    }
+
+
+    //! Enable coefficient covariance computation
+    virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100,
+				      bool ftype=false,
+				      bool covr_tot=false,
+				      bool covar_in=false)
+    {
+      pcavar = pcavar_in;
+      sampT  = sampT_in;
+      scovr  = covr_tot;
+      covar  = covar_in;
+      if (pcavar) {
+	init_covariance();
+	covarStore = std::make_shared<SubsampleCovariance>
+	  ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, scovr, covar);
+      }
+    }
+
   };
 
   /**

expui/BiorthBasis.cc

@@ -3416,11 +3416,14 @@ namespace BasisClasses
     "hslab",
     "zmax",
     "ngrid",
+    "cmap",
     "type",
     "knots",
     "verbose",
     "check",
-    "method"
+    "method",
+    "self_consistent",
+    "cachename"
   };
 
   Slab::Slab(const YAML::Node& CONF) : BiorthBasis(CONF, "slab")
@@ -3476,11 +3479,14 @@ namespace BasisClasses
       if (conf["hslab"])      hslab = conf["hslab"].as<double>();
       if (conf["zmax" ])      zmax  = conf["zmax" ].as<double>();
       if (conf["ngrid"])      ngrid = conf["ngrid"].as<int>();
+      if (conf["cmap" ])      cmap  = conf["cmap" ].as<std::string>();
       if (conf["type" ])      type  = conf["type" ].as<std::string>();
 
       if (conf["knots"])      knots = conf["knots"].as<int>();
 
       if (conf["check"])      check = conf["check"].as<bool>();
+
+      if (conf["cachename"])  cachename = conf["cachename"].as<std::string>();
     } 
     catch (YAML::Exception & error) {
       if (myid==0) std::cout << "Error parsing parameter stanza for <"
@@ -3493,6 +3499,19 @@ namespace BasisClasses
       throw std::runtime_error("Slab: error parsing YAML");
     }
 
+    // Check for non-null cache file name.  This must be specified
+    // to prevent recomputation and unexpected behavior.
+    //
+    if (cachename.size() == 0) {
+      throw std::runtime_error
+	("SlabSL requires a specified cachename in your YAML config\n"
+	 "for consistency with previous invocations and existing coefficient\n"
+	 "sets.  Please add explicitly add 'cachename: name' to your config\n"
+	 "with new 'name' for creating a basis or an existing 'name' for\n"
+	 "reading a previously generated basis cache\n");
+    }
+
+
     // Finally, make the basis
     //
     SLGridSlab::mpi  = 0;
@@ -3502,7 +3521,8 @@ namespace BasisClasses
 
     int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy;
 
-    ortho = std::make_shared<SLGridSlab>(nnmax, nmaxz, ngrid, zmax, type);
+    ortho = std::make_shared<SLGridSlab>(nnmax, nmaxz, ngrid, zmax, cachename,
+					 cmap, type);
 
     // Orthogonality sanity check
     //
@@ -3649,10 +3669,8 @@ namespace BasisClasses
 	  ortho->get_pot(zpot, z, iiy, iix);
 
 	for (int iz=0; iz<imz; iz++) {
-	                       // +--- density in orthogonal series
-                               // |    is 4.0*M_PI rho
-                               // v
-	  expcoef(ix, iy, iz) += -4.0*M_PI*mass*facx*facy*zpot[iz];
+	  // Biorthogonality gives a minus sign here
+	  expcoef(ix, iy, iz) += -mass*facx*facy*zpot[iz];
 	}
       }
     }
@@ -3974,6 +3992,57 @@ namespace BasisClasses
     "fullCovar"
   };
 
+  void Slab::init_covariance()
+  {
+    if (pcavar) {
+
+      int jmax = imx * imy * imz;
+
+      meanV.resize(sampT);
+      for (auto& v : meanV) {
+	v.resize(jmax);
+      }
+
+      if (covar) {
+	if (diagcov) {
+	  dvarV.resize(sampT);
+	  for (auto& v : dvarV) {
+	    v.resize(jmax);
+	  }
+	} else {
+	  covrV.resize(sampT);
+	  for (auto& v : covrV) {
+	    v.resize(jmax, jmax);
+	  }
+	}
+      } else {
+	covrV.clear();
+	dvarV.clear();
+      }
+
+      sampleCounts.resize(sampT);
+      sampleMasses.resize(sampT);
+
+      zero_covariance();
+    }
+  }
+
+
+  void Slab::zero_covariance()
+  {
+    for (int T=0; T<sampT; T++) {
+      meanV[T].setZero();
+      if (covar) {
+	if (diagcov) dvarV[T].setZero();
+	else         covrV[T].setZero();
+      }
+    }
+
+    sampleCounts.setZero();
+    sampleMasses.setZero();
+  }
+
+
   Cube::Cube(const YAML::Node& CONF) : BiorthBasis(CONF, "cube")
   {
     initialize();
@@ -4525,6 +4594,8 @@ namespace BasisClasses
       coef = std::make_shared<CoefClasses::CylStruct>();
     else if (name.compare("flatdisk") == 0)
       coef = std::make_shared<CoefClasses::CylStruct>();
+    else if (name.compare("slab") == 0)
+      coef = std::make_shared<CoefClasses::SlabStruct>();
     else if (name.compare("cube") == 0)
       coef = std::make_shared<CoefClasses::CubeStruct>();
     else {
@@ -4590,6 +4661,10 @@ namespace BasisClasses
       coefret = std::make_shared<CoefClasses::CylStruct>();
     else if (name.compare("flatdisk") == 0)
       coefret = std::make_shared<CoefClasses::CylStruct>();
+    else if (name.compare("slab") == 0)
+      coefret = std::make_shared<CoefClasses::SlabStruct>();
+    else if (name.compare("cube") == 0)
+      coefret = std::make_shared<CoefClasses::CubeStruct>();
     else {
       std::ostringstream sout;
       sout << "Basis::createCoefficients: basis <" << name << "> not recognized"
@@ -5217,6 +5292,16 @@ namespace BasisClasses
     file.createAttribute<double>("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax);
   }
 
+  void Slab::writeCovarH5Params(HighFive::File& file)
+  {
+    file.createAttribute<int>("nminx", HighFive::DataSpace::From(nminx)).write(nminx);
+    file.createAttribute<int>("nminy", HighFive::DataSpace::From(nminy)).write(nminy);
+    file.createAttribute<int>("nminz", HighFive::DataSpace::From(nminz)).write(nminz);
+    file.createAttribute<int>("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx);
+    file.createAttribute<int>("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy);
+    file.createAttribute<int>("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz);
+  }
+
   void Cube::writeCovarH5Params(HighFive::File& file)
   {
     file.createAttribute<int>("nminx", HighFive::DataSpace::From(nminx)).write(nminx);