libROM  v1.0
Data-driven physical simulation library
GreedySampler.cpp
1 
11 // Description: This class greedily selects parameter points
12 // for the construction of a ROM database.
13 
14 #include "GreedySampler.h"
15 #include "utils/HDFDatabase.h"
16 #include "mpi.h"
17 #include <cmath>
18 #include <algorithm>
19 #include <limits.h>
20 #include <fstream>
21 
22 namespace CAROM {
23 
24 GreedyErrorIndicatorPoint
25 createGreedyErrorIndicatorPoint(const std::shared_ptr<Vector> & point,
26  const std::shared_ptr<Vector> & localROM)
27 {
28  GreedyErrorIndicatorPoint result{point, localROM};
29  return result;
30 }
31 
32 Vector
33 getNearestPoint(const std::vector<Vector> & param_points, Vector & point)
34 {
35 
36  int closest_point_index = getNearestPointIndex(param_points, point);
37  return param_points[closest_point_index];
38 }
39 
40 int
41 getNearestPointIndex(const std::vector<Vector> & param_points, Vector & point)
42 {
43 
44  double closest_dist_to_points = INT_MAX;
45  int closest_point_index = -1;
46 
47  for (int i = 0; i < param_points.size(); i++)
48  {
49  Vector diff;
50  point.minus(param_points[i], diff);
51  double dist = diff.norm();
52  if (dist < closest_dist_to_points)
53  {
54  closest_dist_to_points = dist;
55  closest_point_index = i;
56  }
57  }
58 
59  return closest_point_index;
60 }
61 
62 double
63 getNearestPoint(std::vector<double> & param_points, double point)
64 {
65 
66  int closest_point_index = getNearestPointIndex(param_points, point);
67  return param_points[closest_point_index];
68 }
69 
70 int
71 getNearestPointIndex(const std::vector<double> & param_points, double point)
72 {
73 
74  double closest_dist_to_points = INT_MAX;
75  int closest_point_index = -1;
76 
77  for (int i = 0; i < param_points.size(); i++)
78  {
79  double dist = std::abs(point - param_points[i]);
80  if (dist < closest_dist_to_points)
81  {
82  closest_dist_to_points = dist;
83  closest_point_index = i;
84  }
85  }
86 
87  return closest_point_index;
88 }
89 
90 GreedySampler::GreedySampler(
91  const std::vector<Vector> & parameter_points,
92  bool check_local_rom,
93  double relative_error_tolerance,
94  double alpha,
95  double max_clamp,
96  int subset_size,
97  int convergence_subset_size,
98  std::string output_log_path,
99  std::string warm_start_file_name,
100  bool use_centroid,
101  int random_seed,
102  bool debug_algorithm)
103 {
104 
105  d_num_parameter_points = parameter_points.size();
106  CAROM_VERIFY(d_num_parameter_points >= 1);
107  d_parameter_points = parameter_points;
108 
109  constructObject(check_local_rom, relative_error_tolerance, alpha, max_clamp,
110  subset_size, convergence_subset_size, output_log_path, use_centroid,
111  random_seed, debug_algorithm);
112 }
113 
114 GreedySampler::GreedySampler(
115  const std::vector<double> & parameter_points,
116  bool check_local_rom,
117  double relative_error_tolerance,
118  double alpha,
119  double max_clamp,
120  int subset_size,
121  int convergence_subset_size,
122  std::string output_log_path,
123  std::string warm_start_file_name,
124  bool use_centroid,
125  int random_seed,
126  bool debug_algorithm)
127 {
128 
129  d_num_parameter_points = parameter_points.size();
130 
131  std::vector<Vector> parameter_points_vec;
132  for (int i = 0; i < parameter_points.size(); i++)
133  {
134  Vector vec(1, false);
135  vec.item(0) = parameter_points[i];
136  parameter_points_vec.push_back(vec);
137  }
138  CAROM_VERIFY(d_num_parameter_points >= 1);
139  d_parameter_points = parameter_points_vec;
140 
141  constructObject(check_local_rom, relative_error_tolerance, alpha, max_clamp,
142  subset_size, convergence_subset_size, output_log_path, use_centroid,
143  random_seed, debug_algorithm);
144 }
145 
146 GreedySampler::GreedySampler(
147  const Vector & param_space_min,
148  const Vector & param_space_max,
149  int num_parameter_points,
150  bool check_local_rom,
151  double relative_error_tolerance,
152  double alpha,
153  double max_clamp,
154  int subset_size,
155  int convergence_subset_size,
156  std::string output_log_path,
157  std::string warm_start_file_name,
158  bool use_centroid,
159  int random_seed,
160  bool debug_algorithm)
161 {
162 
163  d_min_param_point = param_space_min;
164  d_max_param_point = param_space_max;
165  d_num_parameter_points = num_parameter_points;
166 
167  constructObject(check_local_rom, relative_error_tolerance, alpha, max_clamp,
168  subset_size,
169  convergence_subset_size, output_log_path, use_centroid, random_seed,
170  debug_algorithm);
171 }
172 
173 GreedySampler::GreedySampler(
174  double param_space_min,
175  double param_space_max,
176  int num_parameter_points,
177  bool check_local_rom,
178  double relative_error_tolerance,
179  double alpha,
180  double max_clamp,
181  int subset_size,
182  int convergence_subset_size,
183  std::string output_log_path,
184  std::string warm_start_file_name,
185  bool use_centroid,
186  int random_seed,
187  bool debug_algorithm)
188 {
189 
190  Vector param_space_min_vec(1, false);
191  param_space_min_vec.item(0) = param_space_min;
192  Vector param_space_max_vec(1, false);
193  param_space_max_vec.item(0) = param_space_max;
194 
195  d_min_param_point = param_space_min_vec;
196  d_max_param_point = param_space_max_vec;
197  d_num_parameter_points = num_parameter_points;
198 
199  constructObject(check_local_rom, relative_error_tolerance, alpha, max_clamp,
200  subset_size,
201  convergence_subset_size, output_log_path, use_centroid, random_seed,
202  debug_algorithm);
203 }
204 
205 GreedySampler::GreedySampler(
206  std::string base_file_name,
207  std::string output_log_path)
208 {
209  CAROM_ASSERT(!base_file_name.empty());
210 
211  int mpi_init;
212  MPI_Initialized(&mpi_init);
213  if (mpi_init) {
214  MPI_Comm_rank(MPI_COMM_WORLD, &d_rank);
215  }
216  else {
217  d_rank = 0;
218  }
219 
220  d_output_log_path = output_log_path;
221 
222  load(base_file_name);
223 }
224 
225 void
226 GreedySampler::addDatabaseFromFile(
227  std::string const& warm_start_file_name)
228 {
229  char tmp[100];
230  std::string full_file_name = warm_start_file_name;
231  HDFDatabase database;
232  database.open(full_file_name, "r");
233 
234  sprintf(tmp, "num_parameter_sampled_indices");
235  int num_parameter_sampled_indices;
236  database.getInteger(tmp, num_parameter_sampled_indices);
237  if (num_parameter_sampled_indices > 0)
238  {
239  int temp_parameter_sampled_indices[num_parameter_sampled_indices];
240  sprintf(tmp, "parameter_sampled_indices");
241  database.getIntegerArray(tmp, &temp_parameter_sampled_indices[0],
242  num_parameter_sampled_indices);
243  for (int i = 0; i < num_parameter_sampled_indices; i++)
244  {
245  std::string vec_path = warm_start_file_name + "_" + std::to_string(i);
246  Vector point;
247  point.local_read(vec_path, 0);
248  for (int j = 0; j < d_parameter_points.size(); j++)
249  {
250  Vector diff;
251  point.minus(d_parameter_points[j], diff);
252  double dist = diff.norm();
253  if (dist > 1e-12)
254  {
255  d_parameter_points.push_back(point);
256  d_parameter_sampled_indices.insert(d_parameter_points.size() - 1);
257  }
258  }
259  }
260  }
261 }
262 
263 void
264 GreedySampler::load(
265  std::string base_file_name)
266 {
267  char tmp[100];
268  std::string full_file_name = base_file_name;
269  HDFDatabase database;
270  database.open(full_file_name, "r");
271 
272  sprintf(tmp, "num_parameter_points");
273  int num_parameter_points;
274  database.getInteger(tmp, num_parameter_points);
275  for (int i = 0; i < num_parameter_points; i++)
276  {
277  std::string vec_path = base_file_name + "_" + std::to_string(i);
278  Vector point;
279  point.local_read(vec_path, 0);
280  d_parameter_points.push_back(point);
281  }
282 
283  std::string vec_path = base_file_name + "_min_point";
284  d_min_param_point.local_read(vec_path, 0);
285 
286  vec_path = base_file_name + "_max_point";
287  d_max_param_point.local_read(vec_path, 0);
288 
289  sprintf(tmp, "num_parameter_sampled_indices");
290  int num_parameter_sampled_indices;
291  database.getInteger(tmp, num_parameter_sampled_indices);
292  if (num_parameter_sampled_indices > 0)
293  {
294  int temp_parameter_sampled_indices[num_parameter_sampled_indices];
295  sprintf(tmp, "parameter_sampled_indices");
296  database.getIntegerArray(tmp, &temp_parameter_sampled_indices[0],
297  num_parameter_sampled_indices);
298  for (int i = 0; i < num_parameter_sampled_indices; i++)
299  {
300  d_parameter_sampled_indices.insert(temp_parameter_sampled_indices[i]);
301  }
302  }
303 
304  int bool_int_temp;
305  sprintf(tmp, "procedure_completed");
306  database.getInteger(tmp, bool_int_temp);
307  d_procedure_completed = (bool) bool_int_temp;
308 
309  if (!d_procedure_completed)
310  {
311  sprintf(tmp, "max_error");
312  database.getDouble(tmp, d_max_error);
313  sprintf(tmp, "curr_relative_error");
314  database.getDouble(tmp, d_curr_relative_error);
315  sprintf(tmp, "error_indicator_tol");
316  database.getDouble(tmp, d_error_indicator_tol);
317  sprintf(tmp, "relative_error_tol");
318  database.getDouble(tmp, d_relative_error_tol);
319  sprintf(tmp, "alpha");
320  database.getDouble(tmp, d_alpha);
321  sprintf(tmp, "max_clamp");
322  database.getDouble(tmp, d_max_clamp);
323  sprintf(tmp, "max_num_parameter_points");
324  database.getInteger(tmp, d_num_parameter_points);
325  sprintf(tmp, "subset_size");
326  database.getInteger(tmp, d_subset_size);
327  sprintf(tmp, "convergence_subset_size");
328  database.getInteger(tmp, d_convergence_subset_size);
329  sprintf(tmp, "next_point_to_sample");
330  database.getInteger(tmp, d_next_point_to_sample);
331  sprintf(tmp, "next_point_requiring_error_indicator");
332  database.getInteger(tmp, d_next_point_requiring_error_indicator);
333  sprintf(tmp, "check_local_rom");
334  database.getInteger(tmp, bool_int_temp);
335  d_check_local_rom = (bool) bool_int_temp;
336  sprintf(tmp, "use_centroid");
337  database.getInteger(tmp, bool_int_temp);
338  d_use_centroid = (bool) bool_int_temp;
339  sprintf(tmp, "iteration_started");
340  database.getInteger(tmp, bool_int_temp);
341  d_iteration_started = (bool) bool_int_temp;
342  sprintf(tmp, "convergence_started");
343  database.getInteger(tmp, bool_int_temp);
344  d_convergence_started = (bool) bool_int_temp;
345  sprintf(tmp, "next_parameter_point_computed");
346  database.getInteger(tmp, bool_int_temp);
347  d_next_parameter_point_computed = (bool) bool_int_temp;
348  sprintf(tmp, "point_requiring_error_indicator_computed");
349  database.getInteger(tmp, bool_int_temp);
350  d_point_requiring_error_indicator_computed = (bool) bool_int_temp;
351  sprintf(tmp, "subset_created");
352  database.getInteger(tmp, bool_int_temp);
353  d_subset_created = (bool) bool_int_temp;
354  sprintf(tmp, "debug_algorithm");
355  database.getInteger(tmp, bool_int_temp);
356  d_debug_algorithm = (bool) bool_int_temp;
357  sprintf(tmp, "counter");
358  database.getInteger(tmp, d_counter);
359  sprintf(tmp, "subset_counter");
360  database.getInteger(tmp, d_subset_counter);
361  sprintf(tmp, "random_seed");
362  database.getInteger(tmp, d_random_seed);
363 
364  sprintf(tmp, "parameter_point_random_indices");
365  d_parameter_point_random_indices.resize(d_parameter_points.size());
366  database.getIntegerArray(tmp, &d_parameter_point_random_indices[0],
367  d_parameter_points.size());
368 
369  sprintf(tmp, "parameter_point_errors");
370  d_parameter_point_errors.resize(d_parameter_points.size());
371  database.getDoubleArray(tmp, &d_parameter_point_errors[0],
372  d_parameter_points.size());
373 
374  sprintf(tmp, "parameter_point_local_rom");
375  d_parameter_point_local_rom.resize(d_parameter_points.size());
376  database.getIntegerArray(tmp, &d_parameter_point_local_rom[0],
377  d_parameter_points.size());
378 
379  for (int i = 0; i < d_convergence_subset_size; i++)
380  {
381  std::string vec_path = base_file_name + "_conv_" + std::to_string(i);
382  Vector point;
383  point.local_read(vec_path, 0);
384  d_convergence_points.push_back(point);
385  }
386  }
387  database.close();
388 
389  rng.seed(d_random_seed + d_parameter_sampled_indices.size());
390 }
391 
392 void
393 GreedySampler::checkParameterPointInput()
394 {
395  CAROM_VERIFY(d_min_param_point.dim() == d_max_param_point.dim());
396  CAROM_VERIFY(d_num_parameter_points >= 1);
397 
398  bool isGreater = false;
399  for (int i = 0; i < d_min_param_point.dim(); i++)
400  {
401  if (d_num_parameter_points == 1)
402  {
403  if (d_max_param_point.item(i) >= d_min_param_point.item(i))
404  {
405  isGreater = true;
406  break;
407  }
408  }
409  else
410  {
411  if (d_max_param_point.item(i) > d_min_param_point.item(i))
412  {
413  isGreater = true;
414  break;
415  }
416  }
417  }
418  CAROM_VERIFY(isGreater);
419 
420  if (d_rank == 0)
421  {
422  std::string str;
423  str += "Total number of sample points: " + std::to_string(
424  d_num_parameter_points) + "\n";
425  str += "Parameter space minimum: [ ";
426  for (int i = 0 ; i < d_min_param_point.dim(); i++)
427  {
428  str += std::to_string(d_min_param_point.item(i)) + " ";
429  }
430  str += "]\n";
431  str += "Parameter space maximum: [ ";
432  for (int i = 0 ; i < d_max_param_point.dim(); i++)
433  {
434  str += std::to_string(d_max_param_point.item(i)) + " ";
435  }
436  str += "]\n";
437  agnosticPrint(str);
438  }
439 }
440 
441 void
442 GreedySampler::constructObject(
443  bool check_local_rom,
444  double relative_error_tolerance,
445  double alpha,
446  double max_clamp,
447  int subset_size,
448  int convergence_subset_size,
449  std::string output_log_path,
450  bool use_centroid,
451  int random_seed,
452  bool debug_algorithm)
453 {
454  CAROM_VERIFY(relative_error_tolerance > 0.0);
455  CAROM_VERIFY(alpha >= 1.0);
456  CAROM_VERIFY(max_clamp >= 1.0);
457  CAROM_VERIFY(subset_size > 0);
458  CAROM_VERIFY(convergence_subset_size > 0);
459  CAROM_VERIFY(subset_size < convergence_subset_size);
460  CAROM_VERIFY(random_seed > 0);
461 
462  int mpi_init;
463  MPI_Initialized(&mpi_init);
464  if (mpi_init) {
465  MPI_Comm_rank(MPI_COMM_WORLD, &d_rank);
466  }
467  else {
468  d_rank = 0;
469  }
470 
471  d_check_local_rom = check_local_rom;
472  d_error_indicator_tol = 0.0;
473  d_relative_error_tol = relative_error_tolerance;
474  d_alpha = alpha;
475  d_max_clamp = max_clamp;
476  d_subset_size = subset_size;
477  d_convergence_subset_size = convergence_subset_size;
478  d_output_log_path = output_log_path;
479  d_use_centroid = use_centroid;
480  d_max_error = 0;
481  d_next_point_to_sample = -1;
482  d_next_point_requiring_error_indicator = -1;
483  d_next_parameter_point_computed = false;
484  d_point_requiring_error_indicator_computed = false;
485  d_iteration_started = false;
486  d_convergence_started = false;
487  d_subset_created = false;
488  d_procedure_completed = false;
489  d_subset_counter = 0;
490  d_counter = -1;
491  d_debug_algorithm = debug_algorithm;
492  d_random_seed = random_seed;
493 
494  rng.seed(d_random_seed);
495 
496  if (d_rank == 0)
497  {
498  std::string str;
499  str += "Relative error tolerance: " + std::to_string(d_relative_error_tol) +
500  "\n";
501  str += "Alpha constant: " + std::to_string(d_alpha) + "\n";
502  str += "Max clamp constant: " + std::to_string(d_max_clamp) + "\n";
503  str += "Iteration subset size: " + std::to_string(d_subset_size) + "\n";
504  str += "Convergence subset size: " + std::to_string(d_convergence_subset_size) +
505  "\n";
506  agnosticPrint(str);
507  }
508 }
509 
510 void
511 GreedySampler::initializeParameterPoints()
512 {
513  CAROM_VERIFY(d_parameter_points.size() > 0);
514  CAROM_VERIFY(d_subset_size <= d_parameter_points.size());
515  CAROM_VERIFY(d_convergence_subset_size <= d_parameter_points.size());
516 
517  for (int i = 0; i < d_parameter_points.size() - 1; i++) {
518  CAROM_VERIFY(d_parameter_points[i].dim() == d_parameter_points[i + 1].dim());
519  }
520 
521  for (int i = 0 ; i < d_parameter_points.size(); i++) {
522  d_parameter_point_errors.push_back(INT_MAX);
523  d_parameter_point_local_rom.push_back(-1);
524  d_parameter_point_random_indices.push_back(i);
525  }
526 
527  generateConvergenceSubset();
528 }
529 
530 std::shared_ptr<Vector>
531 GreedySampler::getNextParameterPoint()
532 {
533  if (isComplete())
534  {
535  return std::shared_ptr<Vector>(nullptr);
536  }
537  if (d_iteration_started)
538  {
539  return std::shared_ptr<Vector>(nullptr);
540  }
541  if (d_next_parameter_point_computed)
542  {
543  Vector* result = new Vector(d_parameter_points[d_next_point_to_sample]);
544  return std::shared_ptr<Vector>(result);
545  }
546 
547  if (d_parameter_sampled_indices.size() == 0)
548  {
549  d_next_point_to_sample = getCenterPoint(d_parameter_points, d_use_centroid);
550  }
551 
552  d_max_error = 0;
553 
554  int curr_point_to_sample = d_next_point_to_sample;
555 
556  d_convergence_started = false;
557  d_subset_created = false;
558  d_subset_counter = 0;
559  d_counter = -1;
560 
561  auto search = d_parameter_sampled_indices.find(curr_point_to_sample);
562  CAROM_VERIFY(search == d_parameter_sampled_indices.end());
563 
564  d_parameter_sampled_indices.insert(curr_point_to_sample);
565 
566  if (d_rank == 0)
567  {
568  std::string str;
569  str += "\nPoint sampled at [ ";
570  for (int i = 0 ; i < d_parameter_points[curr_point_to_sample].dim(); i++)
571  {
572  str += std::to_string(d_parameter_points[curr_point_to_sample].item(i)) + " ";
573  }
574  str += "]\n";
575  agnosticPrint(str);
576  }
577 
578  d_next_parameter_point_computed = true;
579 
580  Vector* result = new Vector(d_parameter_points[curr_point_to_sample]);
581  return std::shared_ptr<Vector>(result);
582 }
583 
584 GreedyErrorIndicatorPoint
585 GreedySampler::getNextPointRequiringRelativeError()
586 {
587  if (isComplete())
588  {
589  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
590  }
591  if (!d_next_parameter_point_computed)
592  {
593  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
594  }
595 
596  Vector* result1 = new Vector(d_parameter_points[d_next_point_to_sample]);
597  Vector* result2 = nullptr;
598 
599  if (d_parameter_sampled_indices.size() == 1)
600  {
601  result2 = new Vector(d_parameter_points[getNearestROMIndexToParameterPoint(
602  d_next_point_to_sample, false)]);
603  }
604  else
605  {
606  result2 = new Vector(d_parameter_points[getNearestROMIndexToParameterPoint(
607  d_next_point_to_sample, true)]);
608  }
609 
610  return createGreedyErrorIndicatorPoint(std::shared_ptr<Vector>(result1),
611  std::shared_ptr<Vector>(result2));
612 }
613 
614 GreedyErrorIndicatorPoint
615 GreedySampler::getNextPointRequiringErrorIndicator()
616 {
617  if (isComplete())
618  {
619  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
620  }
621  if (!d_iteration_started)
622  {
623  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
624  }
625 
626  if (d_convergence_started)
627  {
628  return getNextConvergencePointRequiringErrorIndicator();
629  }
630  else
631  {
632  return getNextSubsetPointRequiringErrorIndicator();
633  }
634 }
635 
636 GreedyErrorIndicatorPoint
637 GreedySampler::getNextSubsetPointRequiringErrorIndicator()
638 {
639  if (d_point_requiring_error_indicator_computed)
640  {
641  Vector* result1 = new Vector(
642  d_parameter_points[d_next_point_requiring_error_indicator]);
643  Vector* result2 = new Vector(
644  d_parameter_points[getNearestROMIndexToParameterPoint(
645  d_next_point_requiring_error_indicator, false)]);
646  return createGreedyErrorIndicatorPoint(std::shared_ptr<Vector>(result1),
647  std::shared_ptr<Vector>(result2));
648  }
649  if (d_subset_counter == d_subset_size)
650  {
651  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
652  }
653 
654  if(!d_subset_created)
655  {
656  // generate random shuffle
657  if (!d_debug_algorithm)
658  {
659  std::shuffle(d_parameter_point_random_indices.begin(),
660  d_parameter_point_random_indices.end(), rng);
661  }
662  d_subset_created = true;
663  }
664 
665  d_next_point_requiring_error_indicator = -1;
666 
667  while (d_counter < (int) d_parameter_points.size() - 1)
668  {
669  d_counter++;
670  if (d_subset_counter == d_subset_size)
671  {
672  break;
673  }
674  auto search = d_parameter_sampled_indices.find(
675  d_parameter_point_random_indices[d_counter]);
676  if (search == d_parameter_sampled_indices.end())
677  {
678  d_subset_counter++;
679  double curr_error =
680  d_parameter_point_errors[d_parameter_point_random_indices[d_counter]];
681  if (curr_error > d_max_error)
682  {
683  // if we have already computed this error indicator at the same local rom, the error indicator will not improve
684  // no need to calculate the error indicator again
685  if (d_parameter_point_local_rom[d_parameter_point_random_indices[d_counter]] ==
686  getNearestROMIndexToParameterPoint(d_parameter_point_random_indices[d_counter],
687  false))
688  {
689  d_max_error = curr_error;
690  d_next_point_to_sample = d_parameter_point_random_indices[d_counter];
691  if (d_rank == 0)
692  {
693  std::string str;
694  str += "Error indicator at [ ";
695  for (int i = 0 ;
696  i < d_parameter_points[d_parameter_point_random_indices[d_counter]].dim(); i++)
697  {
698  str += std::to_string(
699  d_parameter_points[d_parameter_point_random_indices[d_counter]].item(i)) + " ";
700  }
701  str += "] skipped.\n";
702  str += "Error indicator " + std::to_string(curr_error) +
703  " already computed at the same local ROM.\n";
704  agnosticPrint(str);
705  }
706  }
707  else
708  {
709  d_next_point_requiring_error_indicator =
710  d_parameter_point_random_indices[d_counter];
711  d_point_requiring_error_indicator_computed = true;
712  Vector* result1 = new Vector(
713  d_parameter_points[d_next_point_requiring_error_indicator]);
714  Vector* result2 = new Vector(
715  d_parameter_points[getNearestROMIndexToParameterPoint(
716  d_next_point_requiring_error_indicator, false)]);
717  return createGreedyErrorIndicatorPoint(std::shared_ptr<Vector>(result1),
718  std::shared_ptr<Vector>(result2));
719  }
720  }
721  else
722  {
723  if (d_rank == 0)
724  {
725  std::string str;
726  str += "Error indicator at [ ";
727  for (int i = 0 ;
728  i < d_parameter_points[d_parameter_point_random_indices[d_counter]].dim(); i++)
729  {
730  str += std::to_string(
731  d_parameter_points[d_parameter_point_random_indices[d_counter]].item(i)) + " ";
732  }
733  str += "] skipped.\n";
734  str += "Error indicator " + std::to_string(curr_error) +
735  " is less than current max error " + std::to_string(d_max_error) + "\n";
736  agnosticPrint(str);
737  }
738  }
739  }
740  }
741  if (d_next_point_requiring_error_indicator == -1)
742  {
743  if (d_rank == 0)
744  {
745  std::string str;
746  str += "Ran out of points to calculate error indicator for in this iteration.\n";
747  agnosticPrint(str);
748  }
749  if (d_max_error < d_error_indicator_tol)
750  {
751  startConvergence();
752  }
753  else
754  {
755  printErrorIndicatorToleranceNotMet();
756  d_iteration_started = false;
757  }
758  }
759 
760  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
761 }
762 
763 GreedyErrorIndicatorPoint
764 GreedySampler::getNextConvergencePointRequiringErrorIndicator()
765 {
766  if (d_point_requiring_error_indicator_computed)
767  {
768  Vector* result1 = new Vector(
769  d_convergence_points[d_next_point_requiring_error_indicator]);
770  std::shared_ptr<Vector> result2 = getNearestROM(
771  d_convergence_points[d_next_point_requiring_error_indicator]);
772  return createGreedyErrorIndicatorPoint(std::shared_ptr<Vector>(result1),
773  result2);
774  }
775  if (d_counter == d_convergence_subset_size)
776  {
777  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
778  }
779 
780  d_next_point_requiring_error_indicator = -1;
781 
782  //get next point requiring error indicator
783  if (d_counter < (int) d_convergence_points.size())
784  {
785  d_next_point_requiring_error_indicator = d_counter;
786  d_point_requiring_error_indicator_computed = true;
787  Vector* result1 = new Vector(
788  d_convergence_points[d_next_point_requiring_error_indicator]);
789  std::shared_ptr<Vector> result2 = getNearestROM(
790  d_convergence_points[d_next_point_requiring_error_indicator]);
791  return createGreedyErrorIndicatorPoint(std::shared_ptr<Vector>(result1),
792  result2);
793  }
794 
795  if (d_next_point_requiring_error_indicator == -1)
796  {
797  printConvergenceAchieved();
798  d_procedure_completed = true;
799  }
800 
801  return createGreedyErrorIndicatorPoint(nullptr, nullptr);
802 }
803 
804 void
805 GreedySampler::printSamplingType(std::string sampling_type)
806 {
807  if (d_rank == 0)
808  {
809  std::string str;
810  str += "Sampling type: " + sampling_type + "\n";
811  agnosticPrint(str);
812  }
813 }
814 
815 void
816 GreedySampler::printConvergenceAchieved()
817 {
818  if (d_rank == 0)
819  {
820  std::string str;
821  str += "Convergence achieved.\n";
822  agnosticPrint(str);
823 
824  str = "\nSampled Parameter Points\n";
825  std::vector<std::pair<double, int>> first_dim_of_sampled_points;
826  for (auto itr = d_parameter_sampled_indices.begin();
827  itr != d_parameter_sampled_indices.end(); ++itr) {
828  first_dim_of_sampled_points.push_back(std::make_pair(
829  d_parameter_points[*itr].item(0), *itr));
830  }
831  sort(first_dim_of_sampled_points.begin(), first_dim_of_sampled_points.end());
832  for (int i = 0; i < first_dim_of_sampled_points.size(); i++)
833  {
834  str += "[ ";
835  for (int j = 0 ;
836  j < d_parameter_points[first_dim_of_sampled_points[i].second].dim(); j++)
837  {
838  str += std::to_string(
839  d_parameter_points[first_dim_of_sampled_points[i].second].item(j)) + " ";
840  }
841  str += "]\n";
842  }
843  agnosticPrint(str);
844  }
845 }
846 
847 void
848 GreedySampler::setPointRelativeError(double error)
849 {
850  CAROM_VERIFY(error >= 0);
851  CAROM_VERIFY(d_next_parameter_point_computed);
852 
853  if (!std::isfinite(error))
854  {
855  error = INT_MAX;
856  }
857 
858  if (d_rank == 0)
859  {
860  if (d_output_log_path == "")
861  {
862  std::string str;
863  str += "Relative error computed at [ ";
864  for (int i = 0 ; i < d_parameter_points[d_next_point_to_sample].dim(); i++)
865  {
866  str += std::to_string(d_parameter_points[d_next_point_to_sample].item(i)) + " ";
867  }
868  str += "]\n";
869  str += "Relative error: " + std::to_string(error) + "\n";
870  agnosticPrint(str);
871  }
872  }
873 
874  d_curr_relative_error = error;
875  d_next_parameter_point_computed = false;
876  d_iteration_started = true;
877 
878  double old_error_indicator_tol = d_error_indicator_tol;
879 
880  if (d_parameter_sampled_indices.size() > 1)
881  {
882  double max1 = d_alpha * d_error_indicator_tol;
883  double min1 = d_max_clamp * d_parameter_point_errors[d_next_point_to_sample];
884  double min2 = d_relative_error_tol *
885  d_parameter_point_errors[d_next_point_to_sample] / d_curr_relative_error;
886 
887  if (d_curr_relative_error <= d_relative_error_tol)
888  {
889  if (d_rank == 0)
890  {
891  std::string str;
892  str += "The relative error was smaller than the relative error tolerance. The error indicator tolerance must increase.\n";
893  str += "Alpha constant * relative error tolerance: " + std::to_string(
894  max1) + "\n";
895  str += "The minimum value the error indicator tolerance can take is: " +
896  std::to_string(max1) + "\n";
897  str += "Max clamp constant * current max error indicator: " + std::to_string(
898  min1) + "\n";
899  str += "Relative error tolerance * current max error indicator / current relative error: "
900  + std::to_string(min2) + "\n";
901  str += "The maximum value the error indicator tolerance can take is the minimum of the previous two values: "
902  + std::to_string(std::min(min1, min2)) + "\n";
903  agnosticPrint(str);
904  }
905  d_error_indicator_tol = std::max(max1, std::min(min1, min2));
906  }
907  else
908  {
909  double min1 = d_relative_error_tol *
910  d_parameter_point_errors[d_next_point_to_sample] / d_curr_relative_error;
911 
912  if (d_rank == 0)
913  {
914  std::string str;
915  str += "The relative error was larger than the relative error tolerance. The error indicator tolerance must decrease.\n";
916  str += "Current error indicator tolerance: " + std::to_string(
917  d_error_indicator_tol) + "\n";
918  str += "Relative error tolerance * current max error indicator / current relative error: "
919  + std::to_string(min1) + "\n";
920  str += "The minimum value the error indicator tolerance can take is: " +
921  std::to_string(std::min(d_error_indicator_tol, min1)) + "\n";
922  agnosticPrint(str);
923  }
924  d_error_indicator_tol = std::min(d_error_indicator_tol, min1);
925  }
926  if (d_rank == 0)
927  {
928  if (old_error_indicator_tol != d_error_indicator_tol)
929  {
930  std::string str;
931  str += "Error indicator tolerance was adaptively changed from " +
932  std::to_string(old_error_indicator_tol) + " to " + std::to_string(
933  d_error_indicator_tol) + "\n";
934  agnosticPrint(str);
935  }
936  }
937  }
938 
939  d_parameter_point_errors[d_next_point_to_sample] = 0;
940  d_parameter_point_local_rom[d_next_point_to_sample] = d_next_point_to_sample;
941 
942  if (d_parameter_sampled_indices.size() > 1
943  && d_curr_relative_error <= d_relative_error_tol)
944  {
945  startConvergence();
946  }
947  else
948  {
949  if (d_check_local_rom || d_parameter_sampled_indices.size() == 1)
950  {
951  d_next_point_requiring_error_indicator = d_next_point_to_sample;
952  d_point_requiring_error_indicator_computed = true;
953  }
954  else
955  {
956  // Precompute next error indicator point
957  // This will allow us to figure out if the greedy algorithm has terminated
958  // early without needing an extra call to the error indicator function.
959  getNextPointRequiringErrorIndicator();
960  }
961  }
962 }
963 
964 void
965 GreedySampler::setPointErrorIndicator(double error, int vec_size)
966 {
967  CAROM_VERIFY(error >= 0);
968  CAROM_VERIFY(d_point_requiring_error_indicator_computed);
969 
970  double proc_errors = pow(error, 2);
971  int total_vec_size = vec_size;
972  CAROM_VERIFY(MPI_Allreduce(MPI_IN_PLACE,
973  &proc_errors,
974  1,
975  MPI_DOUBLE,
976  MPI_SUM,
977  MPI_COMM_WORLD) == MPI_SUCCESS);
978  CAROM_VERIFY(MPI_Allreduce(MPI_IN_PLACE,
979  &total_vec_size,
980  1,
981  MPI_INT,
982  MPI_SUM,
983  MPI_COMM_WORLD) == MPI_SUCCESS);
984  proc_errors = sqrt(proc_errors);
985  proc_errors /= total_vec_size;
986 
987  if (!std::isfinite(proc_errors))
988  {
989  proc_errors = INT_MAX;
990  }
991 
992  if (d_convergence_started)
993  {
994  setConvergenceErrorIndicator(proc_errors);
995  }
996  else
997  {
998  setSubsetErrorIndicator(proc_errors);
999  }
1000 }
1001 
1002 void
1003 GreedySampler::printErrorIndicator(const Vector & errorIndicatorPoint,
1004  double proc_errors)
1005 {
1006  if (d_rank == 0)
1007  {
1008  std::string str;
1009  str += "Error indicator computed at [ ";
1010  for (int i = 0 ; i < errorIndicatorPoint.dim(); i++)
1011  {
1012  str += std::to_string(errorIndicatorPoint.item(i)) + " ";
1013  }
1014  str += "]\n";
1015  str += "Error indicator: " + std::to_string(proc_errors) + "\n";
1016  agnosticPrint(str);
1017  }
1018 }
1019 
1020 void
1021 GreedySampler::agnosticPrint(std::string str)
1022 {
1023  if (d_output_log_path == "")
1024  {
1025  std::cout << str;
1026  }
1027  else
1028  {
1029  std::ofstream database_history;
1030  database_history.open(d_output_log_path, std::ios::app);
1031  database_history << str;
1032  database_history.close();
1033  }
1034 }
1035 
1036 void
1037 GreedySampler::printErrorIndicatorToleranceNotMet()
1038 {
1039  if (d_rank == 0)
1040  {
1041  std::string str;
1042  str += "Error indicator tolerance " + std::to_string(d_error_indicator_tol) +
1043  " not met.\n";
1044  agnosticPrint(str);
1045  }
1046 }
1047 
1048 void
1049 GreedySampler::setSubsetErrorIndicator(double proc_errors)
1050 {
1051  if (d_check_local_rom || d_parameter_sampled_indices.size() == 1)
1052  {
1053  auto search = d_parameter_sampled_indices.find(
1054  d_next_point_requiring_error_indicator);
1055  if (search != d_parameter_sampled_indices.end())
1056  {
1057  d_parameter_point_errors[d_next_point_requiring_error_indicator] = proc_errors;
1058  d_parameter_point_local_rom[d_next_point_requiring_error_indicator] =
1059  d_next_point_requiring_error_indicator;
1060 
1061  double old_error_indicator_tol = d_error_indicator_tol;
1062  if (d_parameter_sampled_indices.size() == 1)
1063  {
1064  d_error_indicator_tol = std::max(d_error_indicator_tol, proc_errors);
1065  }
1066  if (d_rank == 0)
1067  {
1068  std::string str;
1069  str += "Local ROM error indicator computed at [ ";
1070  for (int i = 0 ;
1071  i < d_parameter_points[d_next_point_requiring_error_indicator].dim(); i++)
1072  {
1073  str += std::to_string(
1074  d_parameter_points[d_next_point_requiring_error_indicator].item(i)) + " ";
1075  }
1076  str += "]\n";
1077  str += "Local ROM error indicator (tolerance unchecked): " + std::to_string(
1078  proc_errors) + "\n";
1079  if (old_error_indicator_tol != d_error_indicator_tol)
1080  {
1081  str += "Error indicator at the local ROM was higher than the previous tolerance.\n";
1082  str += "The error indicator tolerance should always be at least the error indicator at the local ROM.\n";
1083  str += "Error indicator tolerance was adaptively changed from " +
1084  std::to_string(old_error_indicator_tol) + " to " + std::to_string(
1085  d_error_indicator_tol) + "\n";
1086  }
1087  agnosticPrint(str);
1088  }
1089 
1090  d_point_requiring_error_indicator_computed = false;
1091 
1092  // Precompute next error indicator point
1093  // This will allow us to figure out if the greedy algorithm has terminated
1094  // early without needing an extra call to the error indicator function.
1095  getNextPointRequiringErrorIndicator();
1096 
1097  return;
1098  }
1099  }
1100 
1101  if (proc_errors <
1102  d_parameter_point_errors[d_parameter_point_random_indices[d_counter]])
1103  {
1104  d_parameter_point_errors[d_parameter_point_random_indices[d_counter]] =
1105  proc_errors;
1106  d_parameter_point_local_rom[d_parameter_point_random_indices[d_counter]] =
1107  getNearestROMIndexToParameterPoint(d_parameter_point_random_indices[d_counter],
1108  false);
1109  }
1110 
1111  printErrorIndicator(
1112  d_parameter_points[d_parameter_point_random_indices[d_counter]], proc_errors);
1113 
1114  d_point_requiring_error_indicator_computed = false;
1115 
1116  if (proc_errors > d_max_error)
1117  {
1118  d_max_error = proc_errors;
1119  d_next_point_to_sample = d_parameter_point_random_indices[d_counter];
1120  }
1121 
1122  if (d_subset_counter == d_subset_size
1123  || d_counter == (int) d_parameter_points.size() - 1)
1124  {
1125  if (d_max_error < d_error_indicator_tol)
1126  {
1127  startConvergence();
1128  }
1129  else
1130  {
1131  d_iteration_started = false;
1132  printErrorIndicatorToleranceNotMet();
1133  }
1134  }
1135 
1136  // Precompute next error indicator point
1137  // This will allow us to figure out if the greedy algorithm has terminated
1138  // early without needing an extra call to the error indicator function.
1139  getNextPointRequiringErrorIndicator();
1140 
1141  return;
1142 }
1143 
1144 void
1145 GreedySampler::setConvergenceErrorIndicator(double proc_errors)
1146 {
1147  printErrorIndicator(d_convergence_points[d_counter], proc_errors);
1148 
1149  d_point_requiring_error_indicator_computed = false;
1150 
1151  if (proc_errors > d_max_error)
1152  {
1153  d_max_error = proc_errors;
1154  }
1155 
1156  if (proc_errors >= d_error_indicator_tol)
1157  {
1158  d_iteration_started = false;
1159  printErrorIndicatorToleranceNotMet();
1160  getNextParameterPointAfterConvergenceFailure();
1161  generateConvergenceSubset();
1162  }
1163  else
1164  {
1165  d_counter++;
1166  }
1167 
1168  if (d_counter == d_convergence_subset_size)
1169  {
1170  printConvergenceAchieved();
1171  d_procedure_completed = true;
1172  }
1173 
1174  // Precompute next error indicator point
1175  // This will allow us to figure out if the greedy algorithm has terminated
1176  // early without needing an extra call to the error indicator function.
1177  getNextPointRequiringErrorIndicator();
1178 }
1179 
1180 void
1181 GreedySampler::generateConvergenceSubset()
1182 {
1183  d_convergence_points.clear();
1184  d_convergence_points = generateRandomPoints(d_convergence_subset_size);
1185 }
1186 
1187 void
1188 GreedySampler::startConvergence()
1189 {
1190  d_convergence_started = true;
1191  d_max_error = 0;
1192  d_counter = 0;
1193  d_subset_counter = 0;
1194 
1195  if (d_rank == 0)
1196  {
1197  std::string str;
1198  str += "Error indicator tolerance " + std::to_string(d_error_indicator_tol) +
1199  " met. Computing convergence.\n";
1200  agnosticPrint(str);
1201  }
1202 
1203  // Precompute next error indicator point
1204  // This will allow us to figure out if the greedy algorithm has terminated
1205  // early without needing an extra call to the error indicator function.
1206  getNextPointRequiringErrorIndicator();
1207 }
1208 
1209 std::vector<Vector>
1210 GreedySampler::generateRandomPoints(int num_points)
1211 {
1212  std::vector<Vector> random_points;
1213 
1214  std::vector<std::uniform_real_distribution<double>> unif;
1215  for (int i = 0; i < d_min_param_point.dim(); i++)
1216  {
1217  unif.push_back(std::uniform_real_distribution<double>(d_min_param_point.item(i),
1218  d_max_param_point.item(i)));
1219  }
1220 
1221  for (int i = 0; i < num_points; i++)
1222  {
1223  Vector point(d_min_param_point.dim(), false);
1224  for (int j = 0; j < point.dim(); j++)
1225  {
1226  point.item(j) = unif[j](rng);
1227  }
1228  random_points.push_back(point);
1229  }
1230  return random_points;
1231 }
1232 
1233 std::shared_ptr<Vector>
1234 GreedySampler::getNearestROM(const Vector & point)
1235 {
1236  CAROM_VERIFY(point.dim() == d_parameter_points[0].dim());
1237 
1238  double closest_dist_to_points = INT_MAX;
1239  int closest_point_index = -1;
1240 
1241  for (auto itr = d_parameter_sampled_indices.begin();
1242  itr != d_parameter_sampled_indices.end(); ++itr) {
1243  Vector diff;
1244  point.minus(d_parameter_points[*itr], diff);
1245  double dist = diff.norm();
1246  if (dist < closest_dist_to_points)
1247  {
1248  closest_dist_to_points = dist;
1249  closest_point_index = *itr;
1250  }
1251  }
1252 
1253  if (closest_point_index == -1)
1254  {
1255  return std::shared_ptr<Vector>(nullptr);
1256  }
1257 
1258  Vector* result = new Vector(d_parameter_points[closest_point_index]);
1259  return std::shared_ptr<Vector>(result);
1260 }
1261 
1262 int
1263 GreedySampler::getNearestNonSampledPoint(const Vector & point)
1264 {
1265  double closest_dist_to_points = INT_MAX;
1266  int closest_point_index = -1;
1267 
1268  for (int i = 0; i < d_parameter_points.size(); i++)
1269  {
1270  auto search = d_parameter_sampled_indices.find(i);
1271  if (search == d_parameter_sampled_indices.end())
1272  {
1273  Vector diff;
1274  point.minus(d_parameter_points[i], diff);
1275  double dist = diff.norm();
1276  if (dist < closest_dist_to_points)
1277  {
1278  closest_dist_to_points = dist;
1279  closest_point_index = i;
1280  }
1281  }
1282  }
1283 
1284  return closest_point_index;
1285 }
1286 
1287 int
1288 GreedySampler::getNearestROMIndexToParameterPoint(int index, bool ignore_self)
1289 {
1290 
1291  CAROM_VERIFY(index >= 0 && index < d_parameter_points.size());
1292 
1293  double closest_dist_to_points = INT_MAX;
1294  int closest_point_index = -1;
1295 
1296  for (auto itr = d_parameter_sampled_indices.begin();
1297  itr != d_parameter_sampled_indices.end(); ++itr) {
1298  if (index != *itr)
1299  {
1300  Vector diff;
1301  d_parameter_points[index].minus(d_parameter_points[*itr], diff);
1302  double dist = diff.norm();
1303  if (dist < closest_dist_to_points ||
1304  (dist == closest_dist_to_points && *itr == d_parameter_point_local_rom[index]))
1305  {
1306  closest_dist_to_points = dist;
1307  closest_point_index = *itr;
1308  }
1309  }
1310  else if (!ignore_self)
1311  {
1312  closest_dist_to_points = 0;
1313  closest_point_index = *itr;
1314  break;
1315  }
1316  }
1317 
1318  return closest_point_index;
1319 }
1320 
1321 std::vector<Vector>
1322 GreedySampler::getParameterPointDomain()
1323 {
1324  return d_parameter_points;
1325 }
1326 
1327 std::vector<Vector>
1328 GreedySampler::getSampledParameterPoints()
1329 {
1330  std::vector<Vector> sampled_points;
1331  for (auto itr = d_parameter_sampled_indices.begin();
1332  itr != d_parameter_sampled_indices.end(); ++itr) {
1333  sampled_points.push_back(d_parameter_points[*itr]);
1334  }
1335  return sampled_points;
1336 }
1337 
1338 void
1339 GreedySampler::save(std::string base_file_name)
1340 {
1341  CAROM_ASSERT(!base_file_name.empty());
1342 
1343  if (d_rank == 0)
1344  {
1345  char tmp[100];
1346  std::string full_file_name = base_file_name;
1347  HDFDatabase database;
1348  database.create(full_file_name);
1349 
1350  sprintf(tmp, "num_parameter_points");
1351  database.putInteger(tmp, d_parameter_points.size());
1352  for (int i = 0; i < d_parameter_points.size(); i++)
1353  {
1354  std::string vec_path = base_file_name + "_" + std::to_string(i);
1355  d_parameter_points[i].write(vec_path);
1356  }
1357  for (int i = 0; i < d_convergence_points.size(); i++)
1358  {
1359  std::string vec_path = base_file_name + "_conv_" + std::to_string(i);
1360  d_convergence_points[i].write(vec_path);
1361  }
1362 
1363  std::string vec_path = base_file_name + "_min_point";
1364  d_min_param_point.write(vec_path);
1365 
1366  vec_path = base_file_name + "_max_point";
1367  d_max_param_point.write(vec_path);
1368 
1369  sprintf(tmp, "num_parameter_sampled_indices");
1370  database.putInteger(tmp, d_parameter_sampled_indices.size());
1371  if (d_parameter_sampled_indices.size() > 0)
1372  {
1373  sprintf(tmp, "parameter_sampled_indices");
1374  std::vector<int> d_parameter_sampled_indices_vec(
1375  d_parameter_sampled_indices.begin(), d_parameter_sampled_indices.end());
1376  database.putIntegerArray(tmp, &d_parameter_sampled_indices_vec[0],
1377  d_parameter_sampled_indices.size());
1378  }
1379 
1380  sprintf(tmp, "procedure_completed");
1381  database.putInteger(tmp, d_procedure_completed);
1382 
1383  if (!d_procedure_completed)
1384  {
1385  sprintf(tmp, "max_error");
1386  database.putDouble(tmp, d_max_error);
1387  sprintf(tmp, "curr_relative_error");
1388  database.putDouble(tmp, d_curr_relative_error);
1389  sprintf(tmp, "error_indicator_tol");
1390  database.putDouble(tmp, d_error_indicator_tol);
1391  sprintf(tmp, "relative_error_tol");
1392  database.putDouble(tmp, d_relative_error_tol);
1393  sprintf(tmp, "alpha");
1394  database.putDouble(tmp, d_alpha);
1395  sprintf(tmp, "max_clamp");
1396  database.putDouble(tmp, d_max_clamp);
1397  sprintf(tmp, "max_num_parameter_points");
1398  database.putInteger(tmp, d_num_parameter_points);
1399  sprintf(tmp, "subset_size");
1400  database.putInteger(tmp, d_subset_size);
1401  sprintf(tmp, "convergence_subset_size");
1402  database.putInteger(tmp, d_convergence_subset_size);
1403  sprintf(tmp, "next_point_to_sample");
1404  database.putInteger(tmp, d_next_point_to_sample);
1405  sprintf(tmp, "next_point_requiring_error_indicator");
1406  database.putInteger(tmp, d_next_point_requiring_error_indicator);
1407  sprintf(tmp, "check_local_rom");
1408  database.putInteger(tmp, d_check_local_rom);
1409  sprintf(tmp, "use_centroid");
1410  database.putInteger(tmp, d_use_centroid);
1411  sprintf(tmp, "iteration_started");
1412  database.putInteger(tmp, d_iteration_started);
1413  sprintf(tmp, "convergence_started");
1414  database.putInteger(tmp, d_convergence_started);
1415  sprintf(tmp, "next_parameter_point_computed");
1416  database.putInteger(tmp, d_next_parameter_point_computed);
1417  sprintf(tmp, "point_requiring_error_indicator_computed");
1418  database.putInteger(tmp, d_point_requiring_error_indicator_computed);
1419  sprintf(tmp, "subset_created");
1420  database.putInteger(tmp, d_subset_created);
1421  sprintf(tmp, "debug_algorithm");
1422  database.putInteger(tmp, d_debug_algorithm);
1423  sprintf(tmp, "counter");
1424  database.putInteger(tmp, d_counter);
1425  sprintf(tmp, "subset_counter");
1426  database.putInteger(tmp, d_subset_counter);
1427  sprintf(tmp, "random_seed");
1428  database.putInteger(tmp, d_random_seed);
1429 
1430  sprintf(tmp, "parameter_point_random_indices");
1431  database.putIntegerArray(tmp, &d_parameter_point_random_indices[0],
1432  d_parameter_point_random_indices.size());
1433  sprintf(tmp, "parameter_point_errors");
1434  database.putDoubleArray(tmp, &d_parameter_point_errors[0],
1435  d_parameter_point_errors.size());
1436  sprintf(tmp, "parameter_point_local_rom");
1437  database.putIntegerArray(tmp, &d_parameter_point_local_rom[0],
1438  d_parameter_point_local_rom.size());
1439  }
1440  database.close();
1441  }
1442  MPI_Barrier(MPI_COMM_WORLD);
1443 }
1444 
1445 bool
1446 GreedySampler::isComplete()
1447 {
1448  if (!d_procedure_completed
1449  && d_parameter_sampled_indices.size() == d_num_parameter_points)
1450  {
1451  if (d_rank == 0)
1452  {
1453  std::string str;
1454  str += "Maximum number of parameter points reached. Stopping now.\n";
1455  agnosticPrint(str);
1456  }
1457  d_procedure_completed = true;
1458  }
1459  return d_procedure_completed;
1460 }
1461 
1462 GreedySampler::~GreedySampler()
1463 {
1464 }
1465 
1466 }
CAROM::Database::getInteger
void getInteger(const std::string &key, int &data)
Reads an integer associated with the supplied key from the currently open database file.
Definition: Database.h:181
CAROM::Database::putDouble
void putDouble(const std::string &key, double data)
Writes a double associated with the supplied key to currently open database file.
Definition: Database.h:128
CAROM::Database::putInteger
void putInteger(const std::string &key, int data)
Writes an integer associated with the supplied key to currently open database file.
Definition: Database.h:94
CAROM::Database::getDouble
void getDouble(const std::string &key, double &data)
Reads a double associated with the supplied key from the currently open database file.
Definition: Database.h:215
CAROM::GreedySampler::printErrorIndicator
void printErrorIndicator(const Vector &errorIndicatorPoint, double proc_errors)
Print the error indicator.
Definition: GreedySampler.cpp:1003
CAROM::GreedySampler::GreedySampler
GreedySampler(const std::vector< Vector > &parameter_points, bool check_local_rom, double relative_error_tolerance, double alpha, double max_clamp, int subset_size, int convergence_subset_size, std::string output_log_path="", std::string warm_start_file_name="", bool use_centroid=true, int random_seed=1, bool debug_algorithm=false)
Constructor.
Definition: GreedySampler.cpp:90
CAROM::GreedySampler::printConvergenceAchieved
void printConvergenceAchieved()
Print that convergence was achieved.
Definition: GreedySampler.cpp:816
CAROM::GreedySampler::setPointRelativeError
void setPointRelativeError(double error)
Set the relative error of the specified point.
Definition: GreedySampler.cpp:848
CAROM::GreedySampler::generateConvergenceSubset
void generateConvergenceSubset()
Generate a new set of convergence points.
Definition: GreedySampler.cpp:1181
CAROM::GreedySampler::getNextPointRequiringRelativeError
GreedyErrorIndicatorPoint getNextPointRequiringRelativeError()
Returns the next parameter point for which a relative error is required.
Definition: GreedySampler.cpp:585
CAROM::GreedySampler::d_max_error
double d_max_error
The current max error of the parameter points of the current iteration.
Definition: GreedySampler.h:565
CAROM::GreedySampler::getNearestNonSampledPoint
int getNearestNonSampledPoint(const Vector &point)
Returns the index of the nearest non-sampled parameter point to the given point.
Definition: GreedySampler.cpp:1263
CAROM::GreedySampler::save
virtual void save(std::string base_file_name)
Save the object state to a file.
Definition: GreedySampler.cpp:1339
CAROM::GreedySampler::d_subset_created
bool d_subset_created
Whether the database has already created a random subset for this iteration.
Definition: GreedySampler.h:656
CAROM::GreedySampler::d_next_point_requiring_error_indicator
int d_next_point_requiring_error_indicator
The next parameter point requiring a error indicator.
Definition: GreedySampler.h:615
CAROM::GreedySampler::printSamplingType
void printSamplingType(std::string sampling_type)
Print the sampling type.
Definition: GreedySampler.cpp:805
CAROM::GreedySampler::setSubsetErrorIndicator
void setSubsetErrorIndicator(double proc_errors)
Set the error indicator for a subset point.
Definition: GreedySampler.cpp:1049
CAROM::GreedySampler::d_subset_counter
int d_subset_counter
An internal subset counter.
Definition: GreedySampler.h:674
CAROM::GreedySampler::d_max_clamp
double d_max_clamp
The max_clamp constant.
Definition: GreedySampler.h:580
CAROM::GreedySampler::d_use_centroid
bool d_use_centroid
Whether the use the centroid heuristic for obtaining the first parameter point.
Definition: GreedySampler.h:621
CAROM::GreedySampler::d_point_requiring_error_indicator_computed
bool d_point_requiring_error_indicator_computed
Whether the database has already computed a new paramter point requiring a error indicator.
Definition: GreedySampler.h:650
CAROM::GreedySampler::setConvergenceErrorIndicator
void setConvergenceErrorIndicator(double proc_errors)
Set the error indicator for a convergence point.
Definition: GreedySampler.cpp:1145
CAROM::GreedySampler::isComplete
bool isComplete()
Check if the greedy algorithm procedure is complete.
Definition: GreedySampler.cpp:1446
CAROM::GreedySampler::d_random_seed
int d_random_seed
Then random seed used to generate subsets.
Definition: GreedySampler.h:689
CAROM::GreedySampler::d_subset_size
int d_subset_size
The size of the subset of parameter points used per iteration.
Definition: GreedySampler.h:600
CAROM::GreedySampler::d_debug_algorithm
bool d_debug_algorithm
Turn off randomness for debugging purposes.
Definition: GreedySampler.h:662
CAROM::GreedySampler::getParameterPointDomain
std::vector< Vector > getParameterPointDomain()
Get the domain of the parameter points.
Definition: GreedySampler.cpp:1322
CAROM::GreedySampler::d_relative_error_tol
double d_relative_error_tol
The relative error tolerance used to terminate the greedy procedure.
Definition: GreedySampler.h:590
CAROM::GreedySampler::d_output_log_path
std::string d_output_log_path
Output log path.
Definition: GreedySampler.h:560
CAROM::GreedySampler::constructObject
void constructObject(bool check_local_rom, double relative_error_tolerance, double alpha, double max_clamp, int subset_size, int convergence_subset_size, std::string output_log_path, bool use_centroid, int random_seed, bool debug_algorithm)
Construct the GreedySampler object.
Definition: GreedySampler.cpp:442
CAROM::GreedySampler::d_check_local_rom
bool d_check_local_rom
Whether the check the last sampled local ROM's error indicator after each iteration.
Definition: GreedySampler.h:627
CAROM::GreedySampler::d_parameter_point_random_indices
std::vector< int > d_parameter_point_random_indices
The parameter point indices (used to generate the random subsets).
Definition: GreedySampler.h:545
CAROM::GreedySampler::getNextParameterPointAfterConvergenceFailure
virtual void getNextParameterPointAfterConvergenceFailure()=0
Get the next parameter point to sample after a convergence failure.
CAROM::GreedySampler::d_convergence_points
std::vector< Vector > d_convergence_points
The convergence parameter points to explore.
Definition: GreedySampler.h:525
CAROM::GreedySampler::getNextSubsetPointRequiringErrorIndicator
GreedyErrorIndicatorPoint getNextSubsetPointRequiringErrorIndicator()
Get the next subset point requiring an error indicator.
Definition: GreedySampler.cpp:637
CAROM::GreedySampler::d_parameter_points
std::vector< Vector > d_parameter_points
The parameter points to explore.
Definition: GreedySampler.h:520
CAROM::GreedySampler::d_rank
int d_rank
The rank of the given processor.
Definition: GreedySampler.h:684
CAROM::GreedySampler::d_curr_relative_error
double d_curr_relative_error
The current relative error of the current iteration.
Definition: GreedySampler.h:570
CAROM::GreedySampler::d_parameter_point_errors
std::vector< double > d_parameter_point_errors
The current errors of the parameter points.
Definition: GreedySampler.h:550
CAROM::GreedySampler::initializeParameterPoints
void initializeParameterPoints()
Initialize the list of parameter point candidates to sample.
Definition: GreedySampler.cpp:511
CAROM::GreedySampler::d_counter
int d_counter
An internal counter.
Definition: GreedySampler.h:668
CAROM::GreedySampler::d_next_parameter_point_computed
bool d_next_parameter_point_computed
Whether the database has already computed a new parameter point to sample.
Definition: GreedySampler.h:644
CAROM::GreedySampler::~GreedySampler
virtual ~GreedySampler()
Destructor.
Definition: GreedySampler.cpp:1462
CAROM::GreedySampler::setPointErrorIndicator
void setPointErrorIndicator(double error, int vec_size)
Set the error indicator error of the specified parameter point.
Definition: GreedySampler.cpp:965
CAROM::GreedySampler::rng
std::default_random_engine rng
Random engine used to generate subsets.
Definition: GreedySampler.h:694
CAROM::GreedySampler::getNextConvergencePointRequiringErrorIndicator
GreedyErrorIndicatorPoint getNextConvergencePointRequiringErrorIndicator()
Get the next convergence point requiring an error indicator.
Definition: GreedySampler.cpp:764
CAROM::GreedySampler::d_parameter_point_local_rom
std::vector< int > d_parameter_point_local_rom
The local ROMs used to obtain the current errors of the parameter points.
Definition: GreedySampler.h:555
CAROM::GreedySampler::getNextPointRequiringErrorIndicator
GreedyErrorIndicatorPoint getNextPointRequiringErrorIndicator()
Returns the next parameter point for which an error indicator is required.
Definition: GreedySampler.cpp:615
CAROM::GreedySampler::d_next_point_to_sample
int d_next_point_to_sample
The next parameter point to sample.
Definition: GreedySampler.h:610
CAROM::GreedySampler::agnosticPrint
void agnosticPrint(std::string str)
Print to output_log_file or cout.
Definition: GreedySampler.cpp:1021
CAROM::GreedySampler::d_num_parameter_points
int d_num_parameter_points
The maximum number of parameter points to sample.
Definition: GreedySampler.h:595
CAROM::GreedySampler::getNearestROM
std::shared_ptr< Vector > getNearestROM(const Vector &point)
Returns the nearest local ROM to the specified parameter point.
Definition: GreedySampler.cpp:1234
CAROM::GreedySampler::d_error_indicator_tol
double d_error_indicator_tol
The error indicator tolerance used to terminate the greedy procedure.
Definition: GreedySampler.h:585
CAROM::GreedySampler::load
virtual void load(std::string base_file_name)
Load the object state from a file.
Definition: GreedySampler.cpp:264
CAROM::GreedySampler::getNearestROMIndexToParameterPoint
int getNearestROMIndexToParameterPoint(int index, bool ignore_self)
Returns the index to the nearest local ROM to the specified parameter point index in the parameter po...
Definition: GreedySampler.cpp:1288
CAROM::GreedySampler::d_procedure_completed
bool d_procedure_completed
Whether the greedy procedure has completed.
Definition: GreedySampler.h:679
CAROM::GreedySampler::d_convergence_subset_size
int d_convergence_subset_size
The size of the subset of parameter points used to check convergence.
Definition: GreedySampler.h:605
CAROM::GreedySampler::d_convergence_started
bool d_convergence_started
Whether the database is in the convergence verifying phase.
Definition: GreedySampler.h:638
CAROM::GreedySampler::printErrorIndicatorToleranceNotMet
void printErrorIndicatorToleranceNotMet()
Print that the error indicator was not met.
Definition: GreedySampler.cpp:1037
CAROM::GreedySampler::d_alpha
double d_alpha
The alpha constant.
Definition: GreedySampler.h:575
CAROM::GreedySampler::getSampledParameterPoints
std::vector< Vector > getSampledParameterPoints()
Get the sampled parameter points.
Definition: GreedySampler.cpp:1328
CAROM::GreedySampler::d_max_param_point
Vector d_max_param_point
The maximum value of the parameter space.
Definition: GreedySampler.h:535
CAROM::GreedySampler::d_parameter_sampled_indices
std::set< int > d_parameter_sampled_indices
The parameter points that were already sampled.
Definition: GreedySampler.h:540
CAROM::GreedySampler::startConvergence
void startConvergence()
Switch to convergence mode.
Definition: GreedySampler.cpp:1188
CAROM::GreedySampler::d_min_param_point
Vector d_min_param_point
The minimum value of the parameter space.
Definition: GreedySampler.h:530
CAROM::GreedySampler::addDatabaseFromFile
void addDatabaseFromFile(std::string const &warm_start_file_name)
Do a warm start by adding the sampled parameters from the database in a file to the current database.
Definition: GreedySampler.cpp:226
CAROM::GreedySampler::getNextParameterPoint
std::shared_ptr< Vector > getNextParameterPoint()
Returns the next parameter point for which sampling is required.
Definition: GreedySampler.cpp:531
CAROM::GreedySampler::d_iteration_started
bool d_iteration_started
Whether the database has already computed a new parameter point for the current iteration.
Definition: GreedySampler.h:633
CAROM::GreedySampler::checkParameterPointInput
void checkParameterPointInput()
Construct the list of parameter point candidates to sample.
Definition: GreedySampler.cpp:393
CAROM::GreedySampler::generateRandomPoints
std::vector< Vector > generateRandomPoints(int num_points)
Generate a vector of random points.
Definition: GreedySampler.cpp:1210
CAROM::HDFDatabase::getDoubleArray
virtual void getDoubleArray(const std::string &key, double *data, int nelements, const bool distributed=false) override
Reads an array of doubles associated with the supplied key from the currently open HDF5 database file...
Definition: HDFDatabase.cpp:323
CAROM::HDFDatabase::putIntegerArray
virtual void putIntegerArray(const std::string &key, const int *const data, int nelements, const bool distributed=false) override
Writes an array of integers associated with the supplied key to the currently open HDF5 database file...
Definition: HDFDatabase.cpp:136
CAROM::HDFDatabase::getIntegerArray
virtual void getIntegerArray(const std::string &key, int *data, int nelements, const bool distributed=false) override
Reads an array of integers associated with the supplied key from the currently open HDF5 database fil...
Definition: HDFDatabase.cpp:252
CAROM::HDFDatabase::create
virtual bool create(const std::string &file_name, const MPI_Comm comm=MPI_COMM_NULL) override
Creates a new HDF5 database file with the supplied name.
Definition: HDFDatabase.cpp:45
CAROM::HDFDatabase::open
virtual bool open(const std::string &file_name, const std::string &type, const MPI_Comm comm=MPI_COMM_NULL) override
Opens an existing HDF5 database file with the supplied name.
Definition: HDFDatabase.cpp:76
CAROM::HDFDatabase::putDoubleArray
virtual void putDoubleArray(const std::string &key, const double *const data, int nelements, const bool distributed=false) override
Writes an array of doubles associated with the supplied key to the currently open HDF5 database file.
Definition: HDFDatabase.cpp:189
CAROM::HDFDatabase::close
virtual bool close() override
Closes the currently open HDF5 database file.
Definition: HDFDatabase.cpp:118
CAROM::Vector::norm
double norm() const
Form the norm of this.
Definition: Vector.cpp:216
CAROM::Vector::minus
std::unique_ptr< Vector > minus(const Vector &other) const
Subtracts other and this and returns the result.
Definition: Vector.h:402
CAROM::Vector::write
void write(const std::string &base_file_name)
write Vector into (a) HDF file(s).
Definition: Vector.cpp:325
CAROM::Vector::item
const double & item(int i) const
Const Vector member access.
Definition: Vector.h:468
CAROM::Vector::local_read
void local_read(const std::string &base_file_name, int rank)
read read a single rank of a distributed Vector from (a) HDF file(s).
Definition: Vector.cpp:412
CAROM::Vector::dim
int dim() const
Returns the dimension of the Vector on this processor.
Definition: Vector.h:251
CAROM::getNearestPointIndex
int getNearestPointIndex(const std::vector< Vector > &param_points, Vector &point)
Given a vector, find the nearest point in a domain.
Definition: GreedySampler.cpp:41
CAROM::createGreedyErrorIndicatorPoint
GreedyErrorIndicatorPoint createGreedyErrorIndicatorPoint(const std::shared_ptr< Vector > &point, const std::shared_ptr< Vector > &localROM)
Create a greedy error indicator point.
Definition: GreedySampler.cpp:25
CAROM::getCenterPoint
int getCenterPoint(const std::vector< const Vector * > &points, bool use_centroid)
Get center point of a group of points.
Definition: Vector.cpp:519
Generated by doxygen 1.9.1