MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_AvatarInterface.cpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
47 
48 #include <string>
49 #include <fstream>
50 #include <sstream>
51 #include <vector>
52 #include "Teuchos_Array.hpp"
53 #include "Teuchos_CommHelpers.hpp"
54 #include "MueLu_BaseClass.hpp"
55 #include "Teuchos_RawParameterListHelpers.hpp"
56 
57 // ***********************************************************************
58 /* Notional Parameterlist Structure
59  "avatar: filestem" "{'mystem1','mystem2'}"
60  "avatar: decision tree files" "{'mystem1.trees','mystem2.trees'}"
61  "avatar: names files" "{'mystem1.names','mystem2.names'}"
62  "avatar: good class" "1"
63  "avatar: heuristic" "1"
64  "avatar: bounds file" "{'bounds.data'}"
65  "avatar: muelu parameter mapping"
66  - "param0'
67  - "muelu parameter" "aggregation: threshold"
68  - "avatar parameter" "DROP_TOL"
69  - "muelu values" "{0,1e-4,1e-3,1e-2}"
70  - "avatar values" "{1,2,3,4}
71  - "param1'
72  - "muelu parameter" "smoother: sweeps"
73  - "avatar parameter" "SWEEPS"
74  - "muelu values" "{1,2,3}"
75  - "avatar values" "{1,2,3}"
76 
77 
78  Notional SetMueLuParameters "problemFeatures" Structure
79  "my feature #1" "246.01"
80  "my feature #2" "123.45"
81 
82  */
83 
84 /*
85 TODO List:
86 Modify MueLu
87  Parameter name checking (make sure names match between Avatar’s names file and the parameter / feature names that MueLu sees).
88 */
89 
90 #ifdef HAVE_MUELU_AVATAR
91 
92 extern "C" {
93 #include "avatar_api.h"
94 }
95 
96 namespace MueLu {
97 
98 // ***********************************************************************
99 RCP<const ParameterList> AvatarInterface::GetValidParameterList() const {
100  RCP<ParameterList> validParamList = rcp(new ParameterList());
101 
102  Teuchos::ParameterList pl_dummy;
104  int int_dummy;
105 
106  // Files from which to read Avatar trees
107  validParamList->set<Teuchos::Array<std::string>>("avatar: decision tree files", ar_dummy, "Names of Avatar decision tree files");
108 
109  // Strings from which to read Avatar trees
110  validParamList->set<Teuchos::Array<std::string>>("avatar: decision tree strings", ar_dummy, "Avatar decision tree strings");
111 
112  // Files from which to read Avatar names
113  validParamList->set<Teuchos::Array<std::string>>("avatar: names files", ar_dummy, "Names of Avatar decision names files");
114 
115  // Strings from which to read Avatar names
116  validParamList->set<Teuchos::Array<std::string>>("avatar: names strings", ar_dummy, "Avatar decision names strings");
117 
118  // Filestem arg for Avatar
119  validParamList->set<Teuchos::Array<std::string>>("avatar: filestem", ar_dummy, "Filestem for the files Avatar requires");
120 
121  // This should be a MueLu parameter-to-Avatar parameter mapping (e.g. if Avatar doesn't like spaces)
122  validParamList->set<Teuchos::ParameterList>("avatar: muelu parameter mapping", pl_dummy, "Mapping of MueLu to Avatar Parameters");
123 
124  // "Good" Class ID for Avatar
125  validParamList->set<int>("avatar: good class", int_dummy, "Numeric code for class Avatar considers to be good");
126 
127  // Which drop tol choice heuristic to use
128  validParamList->set<int>("avatar: heuristic", int_dummy, "Numeric code for which heuristic we want to use");
129 
130  // Bounds file for extrapolation risk
131  validParamList->set<Teuchos::Array<std::string>>("avatar: bounds file", ar_dummy, "Bounds file for Avatar extrapolation risk");
132 
133  // Add dummy variables at the start
134  validParamList->set<int>("avatar: initial dummy variables", int_dummy, "Number of dummy variables to add at the start");
135 
136  // Add dummy variables before the class
137  validParamList->set<int>("avatar: pre-class dummy variables", int_dummy, "Number of dummy variables to add at the before the class");
138 
139  // Value of the dummy variables
140  validParamList->set<int>("avatar: dummy value", int_dummy, "Value of the dummy variables to add at the before/after the class");
141 
142  return validParamList;
143 }
144 
145 // ***********************************************************************
146 Teuchos::ArrayRCP<std::string> AvatarInterface::ReadFromFiles(const char* paramName) const {
147  // const Teuchos::Array<std::string> & tf = params_.get<const Teuchos::Array<std::string> >(paramName);
148  Teuchos::Array<std::string>& tf = params_.get<Teuchos::Array<std::string>>(paramName);
150  // Only Proc 0 will read the files and print the strings
151  if (comm_->getRank() == 0) {
152  treelist.resize(tf.size());
153  for (Teuchos_Ordinal i = 0; i < tf.size(); i++) {
154  std::fstream file;
155  std::stringstream ss;
156  file.open(tf[i]);
157  ss << file.rdbuf();
158  treelist[i] = ss.str();
159  file.close();
160  }
161  }
162  return treelist;
163 }
164 
165 // ***********************************************************************
166 void AvatarInterface::Setup() {
167  // Sanity check
168  if (comm_.is_null()) throw std::runtime_error("MueLu::AvatarInterface::Setup(): Communicator cannot be null");
169 
170  // Get the avatar strings (NOTE: Only exist on proc 0)
171  if (params_.isParameter("avatar: decision tree strings"))
172  avatarStrings_ = Teuchos::arcpFromArray(params_.get<Teuchos::Array<std::string>>("avatar: decision tree strings"));
173  else
174  avatarStrings_ = ReadFromFiles("avatar: decision tree files");
175 
176  if (params_.isParameter("avatar: names strings"))
177  namesStrings_ = Teuchos::arcpFromArray(params_.get<Teuchos::Array<std::string>>("avatar: names strings"));
178  else
179  namesStrings_ = ReadFromFiles("avatar: names files");
180 
181  if (params_.isParameter("avatar: bounds file"))
182  boundsString_ = ReadFromFiles("avatar: bounds file");
183 
184  filestem_ = params_.get<Teuchos::Array<std::string>>("avatar: filestem");
185 
186  if (comm_->getRank() == 0) {
187  // Now actually set up avatar - Avatar's cleanup routine will free the pointer
188  // Avatar_handle* avatar_train(int argc, char **argv, char* names_file, int names_file_is_a_string, char* train_file, int train_file_is_a_string);
189  const int namesfile_is_a_string = 1;
190  const int treesfile_is_a_string = 1;
191  avatarHandle_ = avatar_load(const_cast<char*>(filestem_[0].c_str()), const_cast<char*>(namesStrings_[0].c_str()), namesfile_is_a_string, const_cast<char*>(avatarStrings_[0].c_str()), treesfile_is_a_string);
192  }
193 
194  // Which class does Avatar consider "good"
195  avatarGoodClass_ = params_.get<int>("avatar: good class");
196 
197  heuristicToUse_ = params_.get<int>("avatar: heuristic");
198 
199  // Unpack the MueLu Mapping into something actionable
200  UnpackMueLuMapping();
201 }
202 
203 // ***********************************************************************
204 void AvatarInterface::Cleanup() {
205  avatar_cleanup(avatarHandle_);
206  avatarHandle_ = 0;
207 }
208 
209 // ***********************************************************************
210 void AvatarInterface::GenerateFeatureString(const Teuchos::ParameterList& problemFeatures, std::string& featureString) const {
211  // NOTE: Assumes that the features are in the same order Avatar wants them.
212  std::stringstream ss;
213 
214  // Initial Dummy Variables
215  if (params_.isParameter("avatar: initial dummy variables")) {
216  int num_dummy = params_.get<int>("avatar: initial dummy variables");
217  int dummy_value = params_.get("avatar: dummy value", 666);
218 
219  for (int i = 0; i < num_dummy; i++)
220  ss << dummy_value << ",";
221  }
222 
223  for (Teuchos::ParameterList::ConstIterator i = problemFeatures.begin(); i != problemFeatures.end(); i++) {
224  // const std::string& name = problemFeatures.name(i);
225  const Teuchos::ParameterEntry& entry = problemFeatures.entry(i);
226  if (i != problemFeatures.begin()) ss << ",";
227  entry.leftshift(ss, false); // Because ss<<entry prints out '[unused]' and we don't want that.
228  }
229  featureString = ss.str();
230 }
231 
232 // ***********************************************************************
233 void AvatarInterface::UnpackMueLuMapping() {
234  const Teuchos::ParameterList& mapping = params_.get<Teuchos::ParameterList>("avatar: muelu parameter mapping");
235  // Each MueLu/Avatar parameter pair gets its own sublist. These must be numerically ordered with no gap
236 
237  bool done = false;
238  int idx = 0;
239  int numParams = mapping.numParams();
240 
241  mueluParameterName_.resize(numParams);
242  avatarParameterName_.resize(numParams);
243  mueluParameterValues_.resize(numParams);
244  avatarParameterValues_.resize(numParams);
245 
246  while (!done) {
247  std::stringstream ss;
248  ss << "param" << idx;
249  if (mapping.isSublist(ss.str())) {
250  const Teuchos::ParameterList& sublist = mapping.sublist(ss.str());
251 
252  // Get the names
253  mueluParameterName_[idx] = sublist.get<std::string>("muelu parameter");
254  avatarParameterName_[idx] = sublist.get<std::string>("avatar parameter");
255 
256  // Get the values
257  // FIXME: For now we assume that all of these guys are doubles and their Avatar analogues are doubles
258  mueluParameterValues_[idx] = sublist.get<Teuchos::Array<double>>("muelu values");
259  avatarParameterValues_[idx] = sublist.get<Teuchos::Array<double>>("avatar values");
260 
261  idx++;
262  } else {
263  done = true;
264  }
265  }
266 
267  if (idx != numParams)
268  throw std::runtime_error("MueLu::AvatarInterface::UnpackMueLuMapping(): 'avatar: muelu parameter mapping' has unknown fields");
269 }
270 // ***********************************************************************
271 std::string AvatarInterface::ParamsToString(const std::vector<int>& indices) const {
272  std::stringstream ss;
273  for (Teuchos_Ordinal i = 0; i < avatarParameterValues_.size(); i++) {
274  ss << "," << avatarParameterValues_[i][indices[i]];
275  }
276 
277  // Pre-Class dummy variables
278  if (params_.isParameter("avatar: pre-class dummy variables")) {
279  int num_dummy = params_.get<int>("avatar: pre-class dummy variables");
280  int dummy_value = params_.get("avatar: dummy value", 666);
281  for (int i = 0; i < num_dummy; i++)
282  ss << "," << dummy_value;
283  }
284 
285  return ss.str();
286 }
287 
288 // ***********************************************************************
289 void AvatarInterface::SetIndices(int id, std::vector<int>& indices) const {
290  // The combo numbering here starts with the first guy
291  int numParams = (int)avatarParameterValues_.size();
292  int curr_id = id;
293  for (int i = 0; i < numParams; i++) {
294  int div = avatarParameterValues_[i].size();
295  int mod = curr_id % div;
296  indices[i] = mod;
297  curr_id = (curr_id - mod) / div;
298  }
299 }
300 
301 // ***********************************************************************
302 void AvatarInterface::GenerateMueLuParametersFromIndex(int id, Teuchos::ParameterList& pl) const {
303  // The combo numbering here starts with the first guy
304  int numParams = (int)avatarParameterValues_.size();
305  int curr_id = id;
306  for (int i = 0; i < numParams; i++) {
307  int div = avatarParameterValues_[i].size();
308  int mod = curr_id % div;
309  pl.set(mueluParameterName_[i], mueluParameterValues_[i][mod]);
310  curr_id = (curr_id - mod) / div;
311  }
312 }
313 
314 // ***********************************************************************
315 void AvatarInterface::SetMueLuParameters(const Teuchos::ParameterList& problemFeatures, Teuchos::ParameterList& mueluParams, bool overwrite) const {
316  Teuchos::ParameterList avatarParams;
317  std::string paramString;
318 
319  if (comm_->getRank() == 0) {
320  // Only Rank 0 calls Avatar
321  if (!avatarHandle_) throw std::runtime_error("MueLu::AvatarInterface::SetMueLuParameters(): Setup has not been run");
322 
323  // Turn the problem features into a "trial" string for Avatar
324  std::string trialString;
325  GenerateFeatureString(problemFeatures, trialString);
326 
327  // Compute the number of things we need to test
328  int numParams = (int)avatarParameterValues_.size();
329  std::vector<int> indices(numParams);
330  std::vector<int> sizes(numParams);
331  int num_combos = 1;
332  for (int i = 0; i < numParams; i++) {
333  sizes[i] = avatarParameterValues_[i].size();
334  num_combos *= avatarParameterValues_[i].size();
335  }
336  GetOStream(Runtime0) << "MueLu::AvatarInterface: Testing " << num_combos << " option combinations" << std::endl;
337 
338  // For each input parameter to avatar we iterate over its allowable values and then compute the list of options which Avatar
339  // views as acceptable
340  // FIXME: Find alternative to hard coding malloc size (input deck?)
341  int num_classes = avatar_num_classes(avatarHandle_);
342  std::vector<int> predictions(num_combos, 0);
343  std::vector<float> probabilities(num_classes * num_combos, 0);
344 
345  std::string testString;
346  for (int i = 0; i < num_combos; i++) {
347  SetIndices(i, indices);
348  // Now we add the MueLu parameters into one, enormous Avatar trial string and run avatar once
349  testString += trialString + ParamsToString(indices) + ",0\n";
350  }
351 
352  std::cout << "** Avatar TestString ***\n"
353  << testString << std::endl; // DEBUG
354 
355  int bound_check = true;
356  if (params_.isParameter("avatar: bounds file"))
357  bound_check = checkBounds(testString, boundsString_);
358 
359  // FIXME: Only send in first tree's string
360  // int* avatar_test(Avatar_handle* a, char* test_data_file, int test_data_is_a_string);
361  const int test_data_is_a_string = 1;
362  avatar_test(avatarHandle_, const_cast<char*>(testString.c_str()), test_data_is_a_string, predictions.data(), probabilities.data());
363 
364  // Look at the list of acceptable combinations of options
365  std::vector<int> acceptableCombos;
366  acceptableCombos.reserve(100);
367  for (int i = 0; i < num_combos; i++) {
368  if (predictions[i] == avatarGoodClass_) acceptableCombos.push_back(i);
369  }
370  GetOStream(Runtime0) << "MueLu::AvatarInterface: " << acceptableCombos.size() << " acceptable option combinations found" << std::endl;
371 
372  // Did we have any good combos at all?
373  int chosen_option_id = 0;
374  if (acceptableCombos.size() == 0) {
375  GetOStream(Runtime0) << "WARNING: MueLu::AvatarInterface: found *no* combinations of options which it believes will perform well on this problem" << std::endl
376  << " An arbitrary set of options will be chosen instead" << std::endl;
377  } else {
378  // If there is only one acceptable combination, use it;
379  // otherwise, find the parameter choice with the highest
380  // probability of success
381  if (acceptableCombos.size() == 1) {
382  chosen_option_id = acceptableCombos[0];
383  } else {
384  switch (heuristicToUse_) {
385  case 1:
386  chosen_option_id = hybrid(probabilities.data(), acceptableCombos);
387  break;
388  case 2:
389  chosen_option_id = highProb(probabilities.data(), acceptableCombos);
390  break;
391  case 3:
392  // Choose the first option in the list of acceptable
393  // combinations; the lowest drop tolerance among the
394  // acceptable combinations
395  chosen_option_id = acceptableCombos[0];
396  break;
397  case 4:
398  chosen_option_id = lowCrash(probabilities.data(), acceptableCombos);
399  break;
400  case 5:
401  chosen_option_id = weighted(probabilities.data(), acceptableCombos);
402  break;
403  }
404  }
405  }
406 
407  // If mesh parameters are outside bounding box, set drop tolerance
408  // to 0, otherwise use avatar recommended drop tolerance
409  if (bound_check == 0) {
410  GetOStream(Runtime0) << "WARNING: Extrapolation risk detected, setting drop tolerance to 0" << std::endl;
411  GenerateMueLuParametersFromIndex(0, avatarParams);
412  } else {
413  GenerateMueLuParametersFromIndex(chosen_option_id, avatarParams);
414  }
415  }
416 
417  Teuchos::updateParametersAndBroadcast(outArg(avatarParams), outArg(mueluParams), *comm_, 0, overwrite);
418 }
419 
420 int AvatarInterface::checkBounds(std::string trialString, Teuchos::ArrayRCP<std::string> boundsString) const {
421  std::stringstream ss(trialString);
422  std::vector<double> vect;
423 
424  double b;
425  while (ss >> b) {
426  vect.push_back(b);
427  if (ss.peek() == ',') ss.ignore();
428  }
429 
430  std::stringstream ssBounds(boundsString[0]);
431  std::vector<double> boundsVect;
432 
433  while (ssBounds >> b) {
434  boundsVect.push_back(b);
435  if (ssBounds.peek() == ',') ssBounds.ignore();
436  }
437 
438  int min_idx = (int)std::min(vect.size(), boundsVect.size() / 2);
439 
440  bool inbounds = true;
441  for (int i = 0; inbounds && i < min_idx; i++)
442  inbounds = boundsVect[2 * i] <= vect[i] && vect[i] <= boundsVect[2 * i + 1];
443 
444  return (int)inbounds;
445 }
446 
447 int AvatarInterface::hybrid(float* probabilities, std::vector<int> acceptableCombos) const {
448  float low_crash = probabilities[0];
449  float best_prob = probabilities[2];
450  float diff;
451  int this_combo;
452  int chosen_option_id = acceptableCombos[0];
453  for (int x = 0; x < (int)acceptableCombos.size(); x++) {
454  this_combo = acceptableCombos[x] * 3;
455  diff = probabilities[this_combo] - low_crash;
456  // If this parameter combination has a crash
457  // probability .2 lower than the current "best", we
458  // will use this drop tolerance
459  if (diff < -.2) {
460  low_crash = probabilities[this_combo];
461  best_prob = probabilities[this_combo + 2];
462  chosen_option_id = acceptableCombos[x];
463  }
464  // If this parameter combination has the same
465  // or slightly lower crash probability than the
466  // current best, we compare their "GOOD" probabilities
467  else if (diff <= 0 && probabilities[this_combo + 2] > best_prob) {
468  low_crash = probabilities[this_combo];
469  best_prob = probabilities[this_combo + 2];
470  chosen_option_id = acceptableCombos[x];
471  }
472  }
473  return chosen_option_id;
474 }
475 
476 int AvatarInterface::highProb(float* probabilities, std::vector<int> acceptableCombos) const {
477  float high_prob = probabilities[2];
478  int this_combo;
479  int chosen_option_id = acceptableCombos[0];
480  for (int x = 0; x < (int)acceptableCombos.size(); x++) {
481  this_combo = acceptableCombos[x] * 3;
482  // If this parameter combination has a higher "GOOD"
483  // probability, use this combination
484  if (probabilities[this_combo + 2] > high_prob) {
485  high_prob = probabilities[this_combo + 2];
486  chosen_option_id = acceptableCombos[x];
487  }
488  }
489  return chosen_option_id;
490 }
491 
492 int AvatarInterface::lowCrash(float* probabilities, std::vector<int> acceptableCombos) const {
493  float low_crash = probabilities[0];
494  int this_combo;
495  int chosen_option_id = acceptableCombos[0];
496  for (int x = 0; x < (int)acceptableCombos.size(); x++) {
497  this_combo = acceptableCombos[x] * 3;
498  // If this parameter combination has a lower "CRASH"
499  // probability, use this combination
500  if (probabilities[this_combo] < low_crash) {
501  low_crash = probabilities[this_combo];
502  chosen_option_id = acceptableCombos[x];
503  }
504  }
505  return chosen_option_id;
506 }
507 
508 int AvatarInterface::weighted(float* probabilities, std::vector<int> acceptableCombos) const {
509  float low_crash = probabilities[0];
510  float best_prob = probabilities[2];
511  float diff;
512  int this_combo;
513  int chosen_option_id = acceptableCombos[0];
514  for (int x = 0; x < (int)acceptableCombos.size(); x++) {
515  this_combo = acceptableCombos[x] * 3;
516  diff = probabilities[this_combo] - low_crash;
517  // If this parameter combination has a crash
518  // probability .2 lower than the current "best", we
519  // will use this drop tolerance
520  if (diff < -.2) {
521  low_crash = probabilities[this_combo];
522  best_prob = probabilities[this_combo + 2];
523  chosen_option_id = acceptableCombos[x];
524  }
525  // If this parameter combination is within .1
526  // or has a slightly lower crash probability than the
527  // current best, we compare their "GOOD" probabilities
528  else if (diff <= .1 && probabilities[this_combo + 2] > best_prob) {
529  low_crash = probabilities[this_combo];
530  best_prob = probabilities[this_combo + 2];
531  chosen_option_id = acceptableCombos[x];
532  }
533  }
534  return chosen_option_id;
535 }
536 
537 } // namespace MueLu
538 
539 #endif // HAVE_MUELU_AVATAR
ConstIterator end() const
std::ostream & leftshift(std::ostream &os, bool printFlags=true) const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
One-liner description of what is happening.
Ordinal numParams() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isSublist(const std::string &name) const
params_t::ConstIterator ConstIterator
ConstIterator begin() const
void resize(size_type new_size, const value_type &x=value_type())
const ParameterEntry & entry(ConstIterator i) const
size_type size() const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")