libROM  v1.0
Data-driven physical simulation library
BasisGenerator.cpp
1 
11 // Description: The abstract wrapper class for an abstract SVD algorithm and
12 // sampler. This class provides interfaces to each so that an
13 // application only needs to instantiate one concrete
14 // implementation of this class to control all aspects of basis
15 // vector generation.
16 
17 #include "BasisGenerator.h"
18 #include "svd/StaticSVD.h"
19 #include "svd/RandomizedSVD.h"
20 #include "svd/IncrementalSVDStandard.h"
21 #include "svd/IncrementalSVDFastUpdate.h"
22 #include "svd/IncrementalSVDBrand.h"
23 
24 #include <iomanip>
25 #include <fstream>
26 
27 namespace CAROM {
28 
29 BasisGenerator::BasisGenerator(
30  Options options,
31  bool incremental,
32  const std::string& basis_file_name,
33  Database::formats file_format) :
34  d_dim(options.dim),
35  d_incremental(incremental),
36  d_basis_writer(0),
37  d_basis_reader(0),
38  d_write_snapshots(options.write_snapshots)
39 {
40  CAROM_VERIFY(options.dim > 0);
41  CAROM_VERIFY(options.max_num_samples > 0);
42  CAROM_VERIFY(options.singular_value_tol >= 0);
43  CAROM_VERIFY(options.max_basis_dimension > 0);
44  if (incremental)
45  {
46  CAROM_VERIFY(options.linearity_tol > 0.0);
47  CAROM_VERIFY(options.initial_dt > 0.0);
48  CAROM_VERIFY(options.sampling_tol > 0.0);
49  CAROM_VERIFY(options.max_time_between_samples > 0.0);
50  CAROM_VERIFY(options.min_sampling_time_step_scale >= 0.0);
51  CAROM_VERIFY(options.sampling_time_step_scale >= 0.0);
52  CAROM_VERIFY(options.max_sampling_time_step_scale >= 0.0);
53  CAROM_VERIFY(options.min_sampling_time_step_scale <=
55  }
56 
57  if (!basis_file_name.empty()) {
58  d_basis_writer = new BasisWriter(this, basis_file_name, file_format);
59  }
60  d_update_right_SV = options.update_right_SV;
61  if (incremental)
62  {
63  d_tol = options.sampling_tol;
64  d_max_time_between_samples = options.max_time_between_samples;
65  d_min_sampling_time_step_scale = options.min_sampling_time_step_scale;
66  d_sampling_time_step_scale = options.sampling_time_step_scale;
67  d_max_sampling_time_step_scale = options.max_sampling_time_step_scale;
68  d_dt = options.initial_dt;
69  d_next_sample_time = 0.0;
70 
71  if (options.fast_update_brand) {
72  d_svd.reset(
74  options,
75  basis_file_name));
76  }
77  else if (options.fast_update) {
78  d_svd.reset(
80  options,
81  basis_file_name));
82  }
83  else {
84  d_svd.reset(
86  options,
87  basis_file_name));
88  }
89 
90  // Get the number of processors.
91  int mpi_init;
92  MPI_Initialized(&mpi_init);
93  if (mpi_init) {
94  MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs);
95  }
96  else {
97  d_num_procs = 1;
98  }
99  }
100  else
101  {
102  if (options.randomized) {
103  d_svd.reset(
104  new RandomizedSVD(
105  options));
106  }
107  else {
108  d_svd.reset(
109  new StaticSVD(
110  options));
111  }
112  }
113 }
114 
115 bool
117  double time)
118 {
119  CAROM_VERIFY(time >= 0.0);
120  if (d_incremental)
121  {
122  if(d_update_right_SV)
123  return true;
124  else
125  return time >= d_next_sample_time;
126  }
127 
128  return true;
129 }
130 
131 bool
133  double* u_in,
134  bool add_without_increase)
135 {
136  CAROM_VERIFY(u_in != 0);
137  CAROM_VERIFY(d_svd->getNumSamples() < d_svd->getMaxNumSamples());
138 
139  // Check that u_in is not non-zero.
140  Vector u_vec(u_in, getDim(), true);
141  if (u_vec.norm() == 0.0) {
142  printf("WARNING: BasisGenerator::takeSample skipped trivial sample.\n");
143  return false;
144  }
145 
146  return d_svd->takeSample(u_in, add_without_increase);
147 }
148 
149 void
150 BasisGenerator::loadSampleRange(const std::string& base_file_name,
151  const std::string& kind,
152  int col_min,
153  int col_max,
154  Database::formats db_format)
155 {
156  CAROM_ASSERT(!base_file_name.empty());
157  CAROM_VERIFY(kind == "basis" || kind == "snapshot");
158 
159  if (d_basis_reader) delete d_basis_reader;
160 
161  d_basis_reader = new BasisReader(base_file_name, db_format, d_dim);
162  const Matrix* mat;
163  const Vector* singular_vals;
164 
165  if (kind == "basis") {
167  singular_vals = d_basis_reader->getSingularValues();
168  }
169  else if (kind == "snapshot") {
171  }
172 
173  int num_rows = mat->numRows();
174  int num_cols = mat->numColumns();
175  if (col_min < 0) col_min = 0;
176  if (col_max > num_cols-1) col_max = num_cols-1;
177 
178  CAROM_VERIFY(col_max >= col_min);
179 
180  for (int j = col_min; j <= col_max; j++) {
181  double* u_in = new double[num_rows];
182  for (int i = 0; i < num_rows; i++) {
183  if (kind == "basis") {
184  u_in[i] = mat->item(i,j) * singular_vals->item(j);
185  }
186  else {
187  u_in[i] = mat->item(i,j);
188  }
189  }
190  d_svd->takeSample(u_in, false);
191  delete[] u_in;
192  }
193 }
194 
195 void
196 BasisGenerator::loadSamples(const std::string& base_file_name,
197  const std::string& kind,
198  int cutoff,
199  Database::formats db_format)
200 {
201  loadSampleRange(base_file_name, kind, 0, cutoff-1, db_format);
202 }
203 
204 double
206  double* u_in,
207  double* rhs_in,
208  double time)
209 {
210  CAROM_VERIFY(u_in != 0);
211  CAROM_VERIFY(rhs_in != 0);
212  CAROM_VERIFY(time >= 0);
213  if (d_incremental)
214  {
215 
216  // Check that u_in is not non-zero.
217  int dim = getDim();
218  Vector u_vec(u_in, dim, true);
219  if (u_vec.norm() == 0.0) {
220  return d_next_sample_time;
221  }
222 
223  // Get the current basis vectors.
224  const Matrix* basis = getSpatialBasis();
225 
226  // Compute l = basis' * u
227  Vector* l = basis->transposeMult(u_vec);
228 
229  // basisl = basis * l
230  Vector* basisl = basis->mult(l);
231 
232  // Compute u - basisl.
233  Vector* eta = u_vec.minus(basisl);
234 
235  delete l;
236  delete basisl;
237 
238  // Compute l = basis' * rhs
239  Vector rhs_vec(rhs_in, dim, true);
240  l = basis->transposeMult(rhs_vec);
241 
242  // basisl = basis * l
243  basisl = basis->mult(l);
244 
245  // Compute rhs - basisl.
246  Vector* eta_dot = rhs_vec.minus(basisl);
247 
248  delete l;
249  delete basisl;
250 
251  // Compute the l-inf norm of eta + d_dt*eta_dot.
252  double global_norm;
253  double local_norm = 0.0;
254  for (int i = 0; i < dim; ++i) {
255  double val = fabs(eta->item(i) + d_dt*eta_dot->item(i));
256  if (val > local_norm) {
257  local_norm = val;
258  }
259  }
260  delete eta;
261  delete eta_dot;
262  if (d_num_procs == 1) {
263  global_norm = local_norm;
264  }
265  else {
266  MPI_Allreduce(&local_norm,
267  &global_norm,
268  1,
269  MPI_DOUBLE,
270  MPI_MAX,
271  MPI_COMM_WORLD);
272  }
273 
274  // Compute dt from this norm.
275  double tmp = d_sampling_time_step_scale*sqrt(d_tol/global_norm);
276  if (tmp < d_min_sampling_time_step_scale) {
277  d_dt *= d_min_sampling_time_step_scale;
278  }
279  else if (tmp > d_max_sampling_time_step_scale) {
280  d_dt *= d_max_sampling_time_step_scale;
281  }
282  else {
283  d_dt *= tmp;
284  }
285  if (d_dt < 0) {
286  d_dt = 0.0;
287  }
288  else if (d_dt > d_max_time_between_samples) {
289  d_dt = d_max_time_between_samples;
290  }
291 
292  // Return next sample time.
293  d_next_sample_time = time + d_dt;
294  return d_next_sample_time;
295  }
296 
297  return time;
298 }
299 
300 void
301 BasisGenerator::resetDt(
302  double new_dt)
303 {
304  if (d_incremental)
305  {
306  d_dt = new_dt;
307  }
308 }
309 
310 void
312  const double energyFractionThreshold,
313  int & cutoff,
314  const std::string & cutoffOutputPath,
315  const int first_sv)
316 {
317  const int rom_dim = getSpatialBasis()->numColumns();
318  const Vector* sing_vals = getSingularValues();
319 
320  CAROM_VERIFY(rom_dim <= sing_vals->dim());
321 
322  double sum = 0.0;
323  for (int sv = first_sv; sv < sing_vals->dim(); ++sv) {
324  sum += (*sing_vals)(sv);
325  }
326 
327  int p = std::floor(-std::log10(energyFractionThreshold));
328  std::vector<double> energy_fractions(p);
329 
330  for (int i = 0; i < p; ++i) {
331  energy_fractions[i] = 1 - std::pow(10, -1 - i);
332  }
333 
334  cutoff = first_sv;
335  bool reached_cutoff = false;
336  double partialSum = 0.0;
337  int count = 0;
338 
339  std::ostream* output_stream;
340 
341  if (!cutoffOutputPath.empty()) {
342  output_stream = new std::ofstream(cutoffOutputPath);
343  } else {
344  output_stream = &std::cout;
345  }
346 
347  for (int sv = first_sv; sv < sing_vals->dim(); ++sv) {
348  partialSum += (*sing_vals)(sv);
349  for (int i = count; i < p; ++i)
350  {
351  if (partialSum / sum > 1.0 - std::pow(10, -1 - i))
352  {
353  *output_stream << "For energy fraction: 0.";
354  for (int j = 0; j < i+1; ++j) *output_stream << "9";
355  *output_stream << ", take first " << sv+1 << " of "
356  << sing_vals->dim() << " basis vectors" << std::endl;
357  count += 1;
358  }
359  else
360  {
361  break;
362  }
363  }
364  if (!reached_cutoff && partialSum / sum > 1.0 - energyFractionThreshold)
365  {
366  cutoff = sv+1;
367  reached_cutoff = true;
368  }
369  }
370 
371  if (!reached_cutoff) cutoff = sing_vals->dim();
372  *output_stream << std::fixed << std::setprecision(p+1);
373  *output_stream << "For energy fraction: " << 1.0 - energyFractionThreshold <<
374  ", take first "
375  << cutoff << " of " << sing_vals->dim() << " basis vectors" << std::endl;
376 
377  if (!cutoffOutputPath.empty()) {
378  static_cast<std::ofstream*>(output_stream)->close();
379  delete output_stream;
380  }
381 }
382 
384 {
385  if (d_basis_writer) {
386  delete d_basis_writer;
387  }
388  if (d_basis_reader) {
389  delete d_basis_reader;
390  }
391 }
392 
393 }
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::formats::HDF5)
Load previously saved sample (basis or state) within a column range.
BasisReader * d_basis_reader
Reader of basis vectors.
bool isNextSample(double time)
Returns true if it is time for the next svd sample.
BasisWriter * d_basis_writer
Writer of basis vectors.
const Matrix * getSpatialBasis()
Returns the basis vectors for the current time interval as a Matrix.
void loadSamples(const std::string &base_file_name, const std::string &kind="basis", int cut_off=1e9, Database::formats db_format=Database::formats::HDF5)
Load previously saved sample (basis or state).
boost::shared_ptr< SVD > d_svd
Pointer to the abstract SVD algorithm object.
~BasisGenerator()
Destructor.
double computeNextSampleTime(double *u_in, double *rhs_in, double time)
Computes next time an svd sample is needed.
void finalSummary(const double energyFractionThreshold, int &cutoff, const std::string &cutoffOutputPath="", const int first_sv=0)
Prints the summary of recommended numbers of basis vectors.
const Vector * getSingularValues()
Returns the singular values for the current time interval as a Vector.
bool takeSample(double *u_in, bool add_without_increase=false)
Sample the new state, u_in, at the given time.
Vector * getSingularValues()
Returns the singular values for the requested time.
Matrix * getSnapshotMatrix()
Returns the snapshot matrix for the requested time.
Matrix * getSpatialBasis()
Returns the spatial basis vectors as a Matrix.
Definition: BasisReader.cpp:76
formats
Implemented database file formats. Add to this enum each time a new database format is implemented.
Definition: Database.h:311
int numColumns() const
Returns the number of columns in the Matrix. This method will return the same value from each process...
Definition: Matrix.h:221
Matrix * transposeMult(const Matrix &other) const
Multiplies the transpose of this Matrix with other and returns the product, reference version.
Definition: Matrix.h:643
Matrix * mult(const Matrix &other) const
Multiplies this Matrix with other and returns the product, reference version.
Definition: Matrix.h:284
const double & item(int row, int col) const
Const Matrix member access. Matrix data is stored in row-major format.
Definition: Matrix.h:1015
int numRows() const
Returns the number of rows of the Matrix on this processor.
Definition: Matrix.h:195
bool fast_update_brand
If true use the exact Brand's algorithm for the incremental SVD.
Definition: Options.h:302
double max_time_between_samples
The upper bound on time between samples of the incremental SVD algorithm.
Definition: Options.h:290
int dim
The dimension of the system on this processor.
Definition: Options.h:215
double max_sampling_time_step_scale
The maximum overall scale factor to apply to time step of the incremental SVD algorithm.
Definition: Options.h:341
int max_num_samples
The maximum number of samples.
Definition: Options.h:220
double singular_value_tol
The singular value tolerance used in the SVD algorithm.
Definition: Options.h:240
double linearity_tol
The tolerance of the incremental SVD algorithm to determine whether or not a sample is linearly depen...
Definition: Options.h:271
double sampling_tol
The sampling control tolerance of the incremental SVD algorithm. Limits error in projection of soluti...
Definition: Options.h:284
double min_sampling_time_step_scale
The minimum overall scale factor to apply to time step of the incremental SVD algorithm.
Definition: Options.h:329
bool randomized
Whether to use the randomized SVD algorithm.
Definition: Options.h:252
bool update_right_SV
Whether to update the right singular values.
Definition: Options.h:225
bool fast_update
If true use the fast update version of the incremental SVD algorithm.
Definition: Options.h:296
double initial_dt
The initial simulation time step of the incremental SVD algorithm.
Definition: Options.h:277
int max_basis_dimension
The maximum dimension of the basis.
Definition: Options.h:235
double sampling_time_step_scale
The scaling factor to apply to sampling algorithm of the incremental SVD algorithm.
Definition: Options.h:335
double norm() const
Form the norm of this.
Definition: Vector.cpp:242
Vector * minus(const Vector &other) const
Subtracts other and this and returns the result, reference version.
Definition: Vector.h:538
const double & item(int i) const
Const Vector member access.
Definition: Vector.h:656
int dim() const
Returns the dimension of the Vector on this processor.
Definition: Vector.h:266