From 3881614033ba878e63baa8288a66fe0748fb5d11 Mon Sep 17 00:00:00 2001 From: "Siu Wun \"Tony\" Cheung" Date: Mon, 11 Mar 2024 09:13:44 -0700 Subject: [PATCH] Refactor examples (#253) * Add example results * Count DMD prediction timer only at final time * Move final summary for basis generator into library * Stylize * Add finalSummary to combine_samples * Enable capabilities in combine_samples.cpp * Astyle * Remove spaces * Run astyle * Change to 1-base * Change comments * Change to energy fraction threshold * Pass the string by reference * Astyle * Modify brief * Run astyle --- examples/dmd/dg_advection.cpp | 8 +- examples/dmd/dg_euler.cpp | 8 +- examples/dmd/heat_conduction.cpp | 8 +- examples/dmd/heat_conduction_dmdc.cpp | 8 +- examples/dmd/nonlinear_elasticity.cpp | 28 +++-- examples/dmd/wave_equation.cpp | 6 +- examples/misc/combine_samples.cpp | 29 ++++- examples/prom/dg_advection_global_rom.cpp | 36 +++---- examples/prom/mixed_nonlinear_diffusion.cpp | 51 +-------- .../prom/nonlinear_elasticity_global_rom.cpp | 65 +---------- lib/linalg/BasisGenerator.cpp | 102 ++++++++++++++++-- lib/linalg/BasisGenerator.h | 40 ++++++- 12 files changed, 203 insertions(+), 186 deletions(-) diff --git a/examples/dmd/dg_advection.cpp b/examples/dmd/dg_advection.cpp index 9e1d8b040..6bb905ac0 100644 --- a/examples/dmd/dg_advection.cpp +++ b/examples/dmd/dg_advection.cpp @@ -710,9 +710,7 @@ int main(int argc, char *argv[]) Vector true_solution_u(U->Size()); true_solution_u = U->GetData(); - dmd_prediction_timer.Start(); - - // 14. Predict the state at t_final using DMD. + // 14. Predict using DMD. if (myid == 0) { std::cout << "Predicting solution using DMD" << std::endl; @@ -755,9 +753,9 @@ int main(int argc, char *argv[]) } } - dmd_prediction_timer.Stop(); - + dmd_prediction_timer.Start(); result_u = dmd_U.predict(t_final); + dmd_prediction_timer.Stop(); // 15. Calculate the relative error between the DMD final solution and the true solution. Vector dmd_solution_u(result_u->getData(), result_u->dim()); diff --git a/examples/dmd/dg_euler.cpp b/examples/dmd/dg_euler.cpp index ae4f06885..30b596193 100644 --- a/examples/dmd/dg_euler.cpp +++ b/examples/dmd/dg_euler.cpp @@ -574,9 +574,7 @@ int main(int argc, char *argv[]) Vector true_solution_e(u_block.GetBlock(3).Size()); true_solution_e = u_block.GetBlock(3).GetData(); - dmd_prediction_timer.Start(); - - // 14. Predict the state at t_final using DMD. + // 14. Predict using DMD. if (mpi.WorldRank() == 0) { std::cout << "Predicting density, momentum, and energy using DMD" << std::endl; @@ -639,12 +637,12 @@ int main(int argc, char *argv[]) } } - dmd_prediction_timer.Stop(); - + dmd_prediction_timer.Start(); result_dens = dmd_dens->predict(t_final); result_x_mom = dmd_x_mom->predict(t_final); result_y_mom = dmd_y_mom->predict(t_final); result_e = dmd_e->predict(t_final); + dmd_prediction_timer.Stop(); // 15. Calculate the relative error between the DMD final solution and the true solution. Vector dmd_solution_dens(result_dens->getData(), result_dens->dim()); diff --git a/examples/dmd/heat_conduction.cpp b/examples/dmd/heat_conduction.cpp index baea9b2c7..3b229a87d 100644 --- a/examples/dmd/heat_conduction.cpp +++ b/examples/dmd/heat_conduction.cpp @@ -498,9 +498,7 @@ int main(int argc, char *argv[]) Vector true_solution_u(u.Size()); true_solution_u = u.GetData(); - dmd_prediction_timer.Start(); - - // 14. Predict the state at t_final using DMD. + // 14. Predict using DMD. if (myid == 0) { std::cout << "Predicting temperature using DMD" << std::endl; @@ -540,9 +538,9 @@ int main(int argc, char *argv[]) } } - dmd_prediction_timer.Stop(); - + dmd_prediction_timer.Start(); result_u = dmd_u->predict(t_final); + dmd_prediction_timer.Stop(); // 15. Calculate the relative error between the DMD final solution and the true solution. Vector dmd_solution_u(result_u->getData(), result_u->dim()); diff --git a/examples/dmd/heat_conduction_dmdc.cpp b/examples/dmd/heat_conduction_dmdc.cpp index c1fa4d56e..91347c89a 100644 --- a/examples/dmd/heat_conduction_dmdc.cpp +++ b/examples/dmd/heat_conduction_dmdc.cpp @@ -534,9 +534,7 @@ int main(int argc, char *argv[]) Vector true_solution_u(u.Size()); true_solution_u = u.GetData(); - dmd_prediction_timer.Start(); - - // 14. Predict the state at t_final using DMDc. + // 14. Predict using DMDc. if (myid == 0) { std::cout << "Predicting temperature using DMDc" << std::endl; @@ -576,9 +574,9 @@ int main(int argc, char *argv[]) } } - dmd_prediction_timer.Stop(); - + dmd_prediction_timer.Start(); result_u = dmd_u->predict(t_final); + dmd_prediction_timer.Stop(); // 15. Calculate the relative error between the DMDc final solution and the true solution. Vector dmd_solution_u(result_u->getData(), result_u->dim()); diff --git a/examples/dmd/nonlinear_elasticity.cpp b/examples/dmd/nonlinear_elasticity.cpp index 2c093c818..fce0eecb1 100644 --- a/examples/dmd/nonlinear_elasticity.cpp +++ b/examples/dmd/nonlinear_elasticity.cpp @@ -515,6 +515,12 @@ int main(int argc, char *argv[]) if (last_step || (ti % windowNumSamples) == 0) { + // Calculate DMD modes + if (myid == 0 && rdim != -1 && ef != -1) + { + std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl; + } + if (rdim != -1) { if (myid == 0) @@ -608,10 +614,10 @@ int main(int argc, char *argv[]) w_gf.Save(ee_ofs); } - // 13. Calculate the DMD modes. - if (myid == 0 && rdim != -1 && ef != -1) + // 13. Predict using DMD. + if (myid == 0) { - std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl; + std::cout << "Predicting position and velocity using DMD" << std::endl; } Vector true_solution_x(x_gf.GetTrueVector().Size()); @@ -620,14 +626,6 @@ int main(int argc, char *argv[]) Vector true_solution_v(v_gf.GetTrueVector().Size()); true_solution_v = v_gf.GetTrueVector(); - dmd_prediction_timer.Start(); - - // 14. Predict the state at t_final using DMD. - if (myid == 0) - { - std::cout << "Predicting position and velocity using DMD" << std::endl; - } - curr_window = 0; CAROM::Vector* result_x = dmd_x[curr_window]->predict(ts[0]); CAROM::Vector* result_v = dmd_v[curr_window]->predict(ts[0]); @@ -687,12 +685,12 @@ int main(int argc, char *argv[]) } } - dmd_prediction_timer.Stop(); - + dmd_prediction_timer.Start(); result_x = dmd_x[curr_window]->predict(t_final); result_v = dmd_v[curr_window]->predict(t_final); + dmd_prediction_timer.Stop(); - // 15. Calculate the relative error between the DMD final solution and the true solution. + // 14. Calculate the relative error between the DMD final solution and the true solution. Vector dmd_solution_x(result_x->getData(), result_x->dim()); Vector diff_x(true_solution_x.Size()); subtract(dmd_solution_x, true_solution_x, diff_x); @@ -722,7 +720,7 @@ int main(int argc, char *argv[]) dmd_prediction_timer.RealTime()); } - // 16. Free the used memory. + // 15. Free the used memory. delete ode_solver; delete pmesh; delete result_x; diff --git a/examples/dmd/wave_equation.cpp b/examples/dmd/wave_equation.cpp index 52326e3e5..aa0f800cb 100644 --- a/examples/dmd/wave_equation.cpp +++ b/examples/dmd/wave_equation.cpp @@ -525,8 +525,7 @@ int main(int argc, char *argv[]) dudt_gf.Save(osol); } - // 10. Predict the state at t_final using DMD. - dmd_prediction_timer.Start(); + // 10. Predict using DMD. cout << "Predicting temperature using DMD" << endl; CAROM::Vector* result_u = nullptr; VisItDataCollection dmd_visit_dc("DMD_Wave_Equation", mesh); @@ -564,8 +563,9 @@ int main(int argc, char *argv[]) } } } - dmd_prediction_timer.Stop(); + dmd_prediction_timer.Start(); result_u = dmd_u[curr_window]->predict(t_final); + dmd_prediction_timer.Stop(); // 11. Calculate the relative error between the DMD final solution and the true solution. Vector dmd_solution_u(result_u->getData(), result_u->dim()); diff --git a/examples/misc/combine_samples.cpp b/examples/misc/combine_samples.cpp index f191368f6..22429a29c 100644 --- a/examples/misc/combine_samples.cpp +++ b/examples/misc/combine_samples.cpp @@ -53,15 +53,19 @@ int main(int argc, char* argv[]) std::vector sample_names; int snaps = 0; int dim = 0; + int col_min = 1; + int col_max = 1e9; bool subtract_mean = false; bool subtract_offset = false; + std::string generator_filename = "total"; std::string kind = "snapshot"; std::string offset_file; bool offset_arg = false; if (argc >= 2) { - for (int i = 1; i < argc; i++) { + int i = 1; + while (i < argc) { if (!strcmp(argv[i], "basis") || !strcmp(argv[i], "-b")) { if (rank==0) std::cout << "Argument " << i << " identified as basis or -b" << std::endl; @@ -80,9 +84,25 @@ int main(int argc, char* argv[]) if (rank==0) std::cout << "Will subtract mean" << std::endl; subtract_mean = true; } + else if (!strcmp(argv[i], "file") || !strcmp(argv[i], "-f")) { + i += 1; + generator_filename = argv[i]; + if (rank==0) std::cout << "Output prefix = " << generator_filename << std::endl; + } + else if (!strcmp(argv[i], "col_min") || !strcmp(argv[i], "-cmin")) { + i += 1; + col_min = std::stoi(argv[i]); + if (rank==0) std::cout << "First column to read = " << col_min << std::endl; + } + else if (!strcmp(argv[i], "col_max") || !strcmp(argv[i], "-cmax")) { + i += 1; + col_max = std::stoi(argv[i]); + if (rank==0) std::cout << "Last column to read = " << col_max << std::endl; + } else { sample_names.push_back(argv[i]); } + i += 1; } } else { @@ -112,7 +132,6 @@ int main(int argc, char* argv[]) CAROM_VERIFY((snaps > 0) && (dim > 0)); /*-- Load data from input files --*/ - std::string generator_filename = "total"; std::unique_ptr static_basis_generator; static_basis_generator.reset(new CAROM::BasisGenerator( CAROM::Options(dim, snaps).setMaxBasisDimension(snaps), false, @@ -120,19 +139,22 @@ int main(int argc, char* argv[]) if (rank==0) std::cout << "Loading data from " << kind << std::endl; for(const auto& sample_name: sample_names) { - static_basis_generator->loadSamples(sample_name, kind); + static_basis_generator->loadSampleRange(sample_name, kind, col_min-1, + col_max-1); } if (rank==0) std::cout << "Saving data uploaded as a snapshot matrix" << std::endl; static_basis_generator->writeSnapshot(); + int rdim; if (!subtract_mean && !subtract_offset) { /*-- Compute SVD and save file --*/ if (rank==0) std::cout << "Computing SVD" << std::endl; int rom_dim = static_basis_generator->getSpatialBasis()->numColumns(); if (rank==0) std::cout << "U ROM Dimension: " << rom_dim << std::endl; static_basis_generator->endSamples(); + if (rank==0) static_basis_generator->finalSummary(1e-8, rdim); } else { /*-- load data from hdf5 file to find the mean and subtract it --*/ @@ -188,6 +210,7 @@ int main(int argc, char* argv[]) int rom_dim = static_basis_generator2->getSpatialBasis()->numColumns(); if (rank==0) std::cout << "U ROM Dimension: " << rom_dim << std::endl; static_basis_generator2->endSamples(); + if (rank==0) static_basis_generator2->finalSummary(1e-8, rdim); static_basis_generator2 = nullptr; } diff --git a/examples/prom/dg_advection_global_rom.cpp b/examples/prom/dg_advection_global_rom.cpp index e438766e8..18dc48ea0 100644 --- a/examples/prom/dg_advection_global_rom.cpp +++ b/examples/prom/dg_advection_global_rom.cpp @@ -10,14 +10,21 @@ // libROM MFEM Example: DG Advection (adapted from ex9p.cpp) // +// ================================================================================= +// // Compile with: make dg_advection_global_rom // -// For ROM (reproductive case): +// Sample runs and results for reproductive case: // dg_advection_global_rom -offline // dg_advection_global_rom -merge -ns 1 // dg_advection_global_rom -online // -// For ROM (global rom): +// Output: +// Relative l2 error of ROM solution 0.0005641839849343575 +// +// ================================================================================= +// +// Sample runs and results for global ROM: // Offline phase: dg_advection_global_rom -offline -ff 1.0 -id 0 // dg_advection_global_rom -offline -ff 1.1 -id 1 // dg_advection_global_rom -offline -ff 1.2 -id 2 @@ -28,19 +35,10 @@ // // Online phase: dg_advection_global_rom -online -ff 1.15 // -// Sample runs: -// mpirun -np 4 dg_advection_global_rom -p 0 -dt 0.005 -// mpirun -np 4 dg_advection_global_rom -p 0 -dt 0.01 -// mpirun -np 4 dg_advection_global_rom -p 1 -dt 0.005 -tf 9 -// mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.002 -tf 9 -// mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.02 -s 13 -tf 9 -// mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.004 -tf 9 -// mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.005 -tf 9 -// mpirun -np 4 dg_advection_global_rom -p 3 -rp 2 -dt 0.0025 -tf 9 -vs 20 -// mpirun -np 4 dg_advection_global_rom -p 0 -o 2 -rp 1 -dt 0.01 -tf 8 -// mpirun -np 4 dg_advection_global_rom -p 0 -rs 2 -dt 0.005 -tf 2 -// mpirun -np 4 dg_advection_global_rom -p 0 -rs 1 -o 2 -tf 2 -// mpirun -np 3 dg_advection_global_rom -p 1 -rs 1 -rp 0 -dt 0.005 -tf 0.5 +// Output: +// Relative l2 error of ROM solution 0.0004333183604809453 +// +// ================================================================================= // // Description: This example code solves the time-dependent advection equation // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and @@ -128,10 +126,7 @@ class AIR_prec : public Solver delete AIR_solver; AIR_solver = new HypreBoomerAMG(A_s); AIR_solver->SetAdvectiveOptions(1, "", "FA"); - AIR_solver->SetPrintLevel( - 0); // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine - // this mesh further in parallel to increase the resolution. Once the - // parallel mesh is defined, the serial mesh can be deleted. + AIR_solver->SetPrintLevel(0); AIR_solver->SetMaxLevels(50); } @@ -464,9 +459,6 @@ int main(int argc, char *argv[]) { cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; } - delete mesh; // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine - // this mesh further in parallel to increase the resolution. Once the - // parallel mesh is defined, the serial mesh can be deleted. return 3; } diff --git a/examples/prom/mixed_nonlinear_diffusion.cpp b/examples/prom/mixed_nonlinear_diffusion.cpp index 47f01f284..236e0ca2f 100644 --- a/examples/prom/mixed_nonlinear_diffusion.cpp +++ b/examples/prom/mixed_nonlinear_diffusion.cpp @@ -427,54 +427,6 @@ CAROM::Matrix* GetFirstColumns(const int N, const CAROM::Matrix* A) } // TODO: move this to the library? -void BasisGeneratorFinalSummary(CAROM::BasisGenerator* bg, - const double energyFraction, int & cutoff, const std::string cutoffOutputPath) -{ - const int rom_dim = bg->getSpatialBasis()->numColumns(); - const CAROM::Vector* sing_vals = bg->getSingularValues(); - - MFEM_VERIFY(rom_dim <= sing_vals->dim(), ""); - - double sum = 0.0; - for (int sv = 0; sv < sing_vals->dim(); ++sv) { - sum += (*sing_vals)(sv); - } - - vector energy_fractions = {0.9999, 0.999, 0.99, 0.9}; - bool reached_cutoff = false; - - ofstream outfile(cutoffOutputPath); - - double partialSum = 0.0; - for (int sv = 0; sv < sing_vals->dim(); ++sv) { - partialSum += (*sing_vals)(sv); - for (int i = energy_fractions.size() - 1; i >= 0; i--) - { - if (partialSum / sum > energy_fractions[i]) - { - outfile << "For energy fraction: " << energy_fractions[i] << ", take first " - << sv+1 << " of " << sing_vals->dim() << " basis vectors" << endl; - energy_fractions.pop_back(); - } - else - { - break; - } - } - - if (!reached_cutoff && partialSum / sum > energyFraction) - { - cutoff = sv+1; - reached_cutoff = true; - } - } - - if (!reached_cutoff) cutoff = sing_vals->dim(); - outfile << "Take first " << cutoff << " of " << sing_vals->dim() << - " basis vectors" << endl; - outfile.close(); -} - void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots, std::string name) { @@ -496,9 +448,10 @@ void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots, generator.endSamples(); // save the merged basis file int cutoff = 0; - BasisGeneratorFinalSummary(&generator, 0.9999, cutoff, "mergedSV_" + name); + generator.finalSummary(1e-4, cutoff, "mergedSV_" + name); } +// TODO: move this to the library? const CAROM::Matrix* GetSnapshotMatrix(const int dimFOM, const int nparam, const int max_num_snapshots, std::string name) { diff --git a/examples/prom/nonlinear_elasticity_global_rom.cpp b/examples/prom/nonlinear_elasticity_global_rom.cpp index 677aff4ee..37003869d 100644 --- a/examples/prom/nonlinear_elasticity_global_rom.cpp +++ b/examples/prom/nonlinear_elasticity_global_rom.cpp @@ -307,68 +307,6 @@ CAROM::Matrix *GetFirstColumns(const int N, const CAROM::Matrix *A) } // TODO: move this to the library? -void BasisGeneratorFinalSummary(CAROM::BasisGenerator *bg, - const double energyFraction, int &cutoff, - const std::string cutoffOutputPath) -{ - const int rom_dim = bg->getSpatialBasis()->numColumns(); - const CAROM::Vector *sing_vals = bg->getSingularValues(); - - MFEM_VERIFY(rom_dim <= sing_vals->dim(), ""); - - double sum = 0.0; - for (int sv = 0; sv < sing_vals->dim(); ++sv) - { - sum += (*sing_vals)(sv); - } - - vector energy_fractions = {0.9999999, 0.999999, 0.99999, 0.9999, 0.999, 0.99, 0.9}; - bool reached_cutoff = false; - - ofstream outfile(cutoffOutputPath); - - double partialSum = 0.0; - stringstream prec; - int ctr = 1; - char buffer[100]; - for (int sv = 0; sv < sing_vals->dim(); ++sv) - { - partialSum += (*sing_vals)(sv); - for (int i = energy_fractions.size() - 1; i >= 0; i--) - { - if (partialSum / sum > energy_fractions[i]) - { - // Format string - prec.str(std::string()); - prec << "%." << ctr << "f"; - sprintf(buffer, prec.str().c_str(), energy_fractions[i]); - ctr++; - - // Output string - outfile << "For energy fraction: " << string(buffer) << ", take first " - << sv + 1 << " of " << sing_vals->dim() << " basis vectors" << endl; - energy_fractions.pop_back(); - } - else - { - break; - } - } - - if (!reached_cutoff && partialSum / sum > energyFraction) - { - cutoff = sv + 1; - reached_cutoff = true; - } - } - - if (!reached_cutoff) - cutoff = sing_vals->dim(); - outfile << "Take first " << cutoff << " of " << sing_vals->dim() << - " basis vectors" << endl; - outfile.close(); -} - void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots, std::string name) { @@ -390,8 +328,7 @@ void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots, generator.endSamples(); // save the merged basis file int cutoff = 0; - BasisGeneratorFinalSummary(&generator, 0.9999999, cutoff, - "mergedSV_" + name + ".txt"); + generator.finalSummary(1e-7, cutoff, "mergedSV_" + name + ".txt"); } const CAROM::Matrix *GetSnapshotMatrix(const int dimFOM, const int nparam, diff --git a/lib/linalg/BasisGenerator.cpp b/lib/linalg/BasisGenerator.cpp index 945c34bda..46b9bda9f 100644 --- a/lib/linalg/BasisGenerator.cpp +++ b/lib/linalg/BasisGenerator.cpp @@ -21,6 +21,9 @@ #include "svd/IncrementalSVDFastUpdate.h" #include "svd/IncrementalSVDBrand.h" +#include +#include + namespace CAROM { BasisGenerator::BasisGenerator( @@ -143,10 +146,11 @@ BasisGenerator::takeSample( } void -BasisGenerator::loadSamples(const std::string& base_file_name, - const std::string& kind, - int cut_off, - Database::formats db_format) +BasisGenerator::loadSampleRange(const std::string& base_file_name, + const std::string& kind, + int col_min, + int col_max, + Database::formats db_format) { CAROM_ASSERT(!base_file_name.empty()); CAROM_VERIFY(kind == "basis" || kind == "snapshot"); @@ -167,10 +171,12 @@ BasisGenerator::loadSamples(const std::string& base_file_name, int num_rows = mat->numRows(); int num_cols = mat->numColumns(); - int max_cols = num_cols; - if (cut_off < num_cols) max_cols = cut_off; + if (col_min < 0) col_min = 0; + if (col_max > num_cols-1) col_max = num_cols-1; + + CAROM_VERIFY(col_max >= col_min); - for (int j = 0; j < max_cols; j++) { + for (int j = col_min; j <= col_max; j++) { double* u_in = new double[num_rows]; for (int i = 0; i < num_rows; i++) { if (kind == "basis") { @@ -185,6 +191,15 @@ BasisGenerator::loadSamples(const std::string& base_file_name, } } +void +BasisGenerator::loadSamples(const std::string& base_file_name, + const std::string& kind, + int cutoff, + Database::formats db_format) +{ + loadSampleRange(base_file_name, kind, 0, cutoff-1, db_format); +} + double BasisGenerator::computeNextSampleTime( double* u_in, @@ -291,6 +306,79 @@ BasisGenerator::resetDt( } } +void +BasisGenerator::finalSummary( + const double energyFractionThreshold, + int & cutoff, + const std::string & cutoffOutputPath, + const int first_sv) +{ + const int rom_dim = getSpatialBasis()->numColumns(); + const Vector* sing_vals = getSingularValues(); + + CAROM_VERIFY(rom_dim <= sing_vals->dim()); + + double sum = 0.0; + for (int sv = first_sv; sv < sing_vals->dim(); ++sv) { + sum += (*sing_vals)(sv); + } + + int p = std::floor(-std::log10(energyFractionThreshold)); + std::vector energy_fractions(p); + + for (int i = 0; i < p; ++i) { + energy_fractions[i] = 1 - std::pow(10, -1 - i); + } + + cutoff = first_sv; + bool reached_cutoff = false; + double partialSum = 0.0; + int count = 0; + + std::ostream* output_stream; + + if (!cutoffOutputPath.empty()) { + output_stream = new std::ofstream(cutoffOutputPath); + } else { + output_stream = &std::cout; + } + + for (int sv = first_sv; sv < sing_vals->dim(); ++sv) { + partialSum += (*sing_vals)(sv); + for (int i = count; i < p; ++i) + { + if (partialSum / sum > 1.0 - std::pow(10, -1 - i)) + { + *output_stream << "For energy fraction: 0."; + for (int j = 0; j < i+1; ++j) *output_stream << "9"; + *output_stream << ", take first " << sv+1 << " of " + << sing_vals->dim() << " basis vectors" << std::endl; + count += 1; + } + else + { + break; + } + } + if (!reached_cutoff && partialSum / sum > 1.0 - energyFractionThreshold) + { + cutoff = sv+1; + reached_cutoff = true; + } + } + + if (!reached_cutoff) cutoff = sing_vals->dim(); + *output_stream << std::fixed << std::setprecision(p+1); + *output_stream << "For energy fraction: " << 1.0 - energyFractionThreshold << + ", take first " + << cutoff << " of " << sing_vals->dim() << " basis vectors" << std::endl; + + if (!cutoffOutputPath.empty()) { + static_cast(output_stream)->close(); + delete output_stream; + } +} + BasisGenerator::~BasisGenerator() { if (d_basis_writer) { diff --git a/lib/linalg/BasisGenerator.h b/lib/linalg/BasisGenerator.h index 7ba7b02c5..ec6f7898d 100644 --- a/lib/linalg/BasisGenerator.h +++ b/lib/linalg/BasisGenerator.h @@ -138,20 +138,39 @@ class BasisGenerator } } + /** + * @brief Load previously saved sample (basis or state) + * within a column range. + * + * @param[in] base_file_name The base part of the name of the files + * holding the basis/snapshot vectors. + * @param[in] kind A string equal to "basis" or "snapshot", representing + * which kind of data to load. + * @param[in] col_min The first basis/snapshot vector to read. + * @param[in] col_max The last basis/snapshot vector to read. + * @param[in] db_format Format of the file to read. + */ + void + loadSampleRange(const std::string& base_file_name, + const std::string& kind = "basis", + int col_min = 0, + int col_max = 1e9, + Database::formats db_format = Database::HDF5); + /** * @brief Load previously saved sample (basis or state). * * @param[in] base_file_name The base part of the name of the files - * holding the basis / snapshot vectors. + * holding the basis/snapshot vectors. * @param[in] kind A string equal to "basis" or "snapshot", representing * which kind of data to load. - * @param[in] cut_off The maximum number of bases or snapshots to read. + * @param[in] cutoff The maximum number of basis/snapshot vectors to read. * @param[in] db_format Format of the file to read. */ void loadSamples(const std::string& base_file_name, const std::string& kind = "basis", - int cut_off = 1e9, + int cutoff = 1e9, Database::formats db_format = Database::HDF5); /** @@ -228,6 +247,21 @@ class BasisGenerator return d_svd->getNumSamples(); } + /** + * @brief Prints the summary of recommended numbers of basis vectors. + * + * @param[in] energyFractionThreshold Energy Fraction threshold + * (energy fraction = 1.0 - energyFractionThreshold). + * @param[in] cutoff Number of basis vectors selected. + * @param[in] cutoffOutputPath Path of the summary file. + * @param[in] first_sv First singular vector in the calculaton of energy. + */ + void finalSummary( + const double energyFractionThreshold, + int & cutoff, + const std::string & cutoffOutputPath = "", + const int first_sv = 0); + protected: /** * @brief Writer of basis vectors.