MueLu  Version of the Day
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 // ***********************************************************************
59 /* Notional Parameterlist Structure
60  "avatar: filestem" "{'mystem1','mystem2'}"
61  "avatar: decision tree files" "{'mystem1.trees','mystem2.trees'}"
62  "avatar: names files" "{'mystem1.names','mystem2.names'}"
63  "avatar: good class" "1"
64  "avatar: heuristic" "1"
65  "avatar: bounds file" "{'bounds.data'}"
66  "avatar: muelu parameter mapping"
67  - "param0'
68  - "muelu parameter" "aggregation: threshold"
69  - "avatar parameter" "DROP_TOL"
70  - "muelu values" "{0,1e-4,1e-3,1e-2}"
71  - "avatar values" "{1,2,3,4}
72  - "param1'
73  - "muelu parameter" "smoother: sweeps"
74  - "avatar parameter" "SWEEPS"
75  - "muelu values" "{1,2,3}"
76  - "avatar values" "{1,2,3}"
77 
78 
79  Notional SetMueLuParameters "problemFeatures" Structure
80  "my feature #1" "246.01"
81  "my feature #2" "123.45"
82 
83  */
84 
85 
86 /*
87 TODO List:
88 Modify MueLu
89  Parameter name checking (make sure names match between Avatar’s names file and the parameter / feature names that MueLu sees).
90 */
91 
92 #ifdef HAVE_MUELU_AVATAR
93 
94 extern "C" {
95 #include "avatar_api.h"
96 }
97 
98 
99 
100 namespace MueLu {
101 
102 
103 // ***********************************************************************
104 RCP<const ParameterList> AvatarInterface::GetValidParameterList() const {
105  RCP<ParameterList> validParamList = rcp(new ParameterList());
106 
107  Teuchos::ParameterList pl_dummy;
109  int int_dummy;
110 
111  // Files from which to read Avatar trees
112  validParamList->set<Teuchos::Array<std::string> >("avatar: decision tree files",ar_dummy,"Names of Avatar decision tree files");
113 
114  // Files from which to read Avatar names
115  validParamList->set<Teuchos::Array<std::string> >("avatar: names files",ar_dummy,"Names of Avatar decision names files");
116 
117  // Filestem arg for Avatar
118  validParamList->set<Teuchos::Array<std::string> >("avatar: filestem",ar_dummy,"Filestem for the files Avatar requires");
119 
120  // This should be a MueLu parameter-to-Avatar parameter mapping (e.g. if Avatar doesn't like spaces)
121  validParamList->set<Teuchos::ParameterList>("avatar: muelu parameter mapping",pl_dummy,"Mapping of MueLu to Avatar Parameters");
122 
123  // "Good" Class ID for Avatar
124  validParamList->set<int>("avatar: good class",int_dummy,"Numeric code for class Avatar considers to be good");
125 
126  // Which drop tol choice heuristic to use
127  validParamList->set<int>("avatar: heuristic",int_dummy,"Numeric code for which heurisitc we want to use");
128 
129  // Bounds file for extrapolation risk
130  validParamList->set<Teuchos::Array<std::string> >("avatar: bounds file",ar_dummy,"Bounds file for Avatar extrapolation risk");
131 
132  return validParamList;
133 }
134 
135 
136 // ***********************************************************************
137 Teuchos::ArrayRCP<std::string> AvatarInterface::ReadFromFiles(const char * paramName) const {
138  // const Teuchos::Array<std::string> & tf = params_.get<const Teuchos::Array<std::string> >(paramName);
139  Teuchos::Array<std::string> & tf = params_.get<Teuchos::Array<std::string> >(paramName);
141  // Only Proc 0 will read the files and print the strings
142  if (comm_->getRank() == 0) {
143  treelist.resize(tf.size());
144  for(Teuchos_Ordinal i=0; i<tf.size(); i++) {
145  std::fstream file;
146  std::stringstream ss;
147  file.open(tf[i]);
148  ss << file.rdbuf();
149  treelist[i] = ss.str();
150  file.close();
151  }
152  }
153  return treelist;
154 }
155 
156 
157 
158 // ***********************************************************************
159 void AvatarInterface::Setup() {
160  // Sanity check
161  if(comm_.is_null()) throw std::runtime_error("MueLu::AvatarInterface::Setup(): Communicator cannot be null");
162 
163  // Get the avatar strings (NOTE: Only exist on proc 0)
164  avatarStrings_ = ReadFromFiles("avatar: decision tree files");
165  namesStrings_ = ReadFromFiles("avatar: names files");
166  boundsString_ = ReadFromFiles("avatar: bounds file");
167  filestem_ = params_.get<Teuchos::Array<std::string>>("avatar: filestem");
168 
169 
170  if(comm_->getRank() == 0) {
171  // Now actually set up avatar - Avatar's cleanup routine will free the pointer
172  // 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);
173  const int namesfile_is_a_string = 1;
174  const int treesfile_is_a_string = 1;
175  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);
176 
177  }
178 
179  // Which class does Avatar consider "good"
180  avatarGoodClass_ = params_.get<int>("avatar: good class");
181 
182  heuristicToUse_ = params_.get<int>("avatar: heuristic");
183 
184  // Unpack the MueLu Mapping into something actionable
185  UnpackMueLuMapping();
186 
187 }
188 
189 // ***********************************************************************
190 void AvatarInterface::Cleanup() {
191  avatar_cleanup(avatarHandle_);
192  avatarHandle_=0;
193 }
194 
195 
196 // ***********************************************************************
197 void AvatarInterface::GenerateFeatureString(const Teuchos::ParameterList & problemFeatures, std::string & featureString) const {
198  // NOTE: Assumes that the features are in the same order Avatar wants them.
199  std::stringstream ss;
200  for(Teuchos::ParameterList::ConstIterator i=problemFeatures.begin(); i != problemFeatures.end(); i++) {
201  // const std::string& name = problemFeatures.name(i);
202  const Teuchos::ParameterEntry& entry = problemFeatures.entry(i);
203  if(i!=problemFeatures.begin()) ss<<",";
204  entry.leftshift(ss,false); // Because ss<<entry prints out '[unused]' and we don't want that.
205  }
206  featureString = ss.str();
207 }
208 
209 // ***********************************************************************
210 void AvatarInterface::UnpackMueLuMapping() {
211  const Teuchos::ParameterList & mapping = params_.get<Teuchos::ParameterList>("avatar: muelu parameter mapping");
212  // Each MueLu/Avatar parameter pair gets its own sublist. These must be numerically ordered with no gap
213 
214  bool done=false;
215  int idx=0;
216  int numParams = mapping.numParams();
217 
218  mueluParameterName_.resize(numParams);
219  avatarParameterName_.resize(numParams);
220  mueluParameterValues_.resize(numParams);
221  avatarParameterValues_.resize(numParams);
222 
223  while(!done) {
224  std::stringstream ss;
225  ss << "param" << idx;
226  if(mapping.isSublist(ss.str())) {
227  const Teuchos::ParameterList & sublist = mapping.sublist(ss.str());
228 
229  // Get the names
230  mueluParameterName_[idx] = sublist.get<std::string>("muelu parameter");
231  avatarParameterName_[idx] = sublist.get<std::string>("avatar parameter");
232 
233  // Get the values
234  //FIXME: For now we assume that all of these guys are doubles and their Avatar analogues are doubles
235  mueluParameterValues_[idx] = sublist.get<Teuchos::Array<double> >("muelu values");
236  avatarParameterValues_[idx] = sublist.get<Teuchos::Array<double> >("avatar values");
237 
238  idx++;
239  }
240  else {
241  done=true;
242  }
243  }
244 
245  if(idx!=numParams)
246  throw std::runtime_error("MueLu::AvatarInterface::UnpackMueLuMapping(): 'avatar: muelu parameter mapping' has unknown fields");
247 }
248 // ***********************************************************************
249 std::string AvatarInterface::ParamsToString(const std::vector<int> & indices) const {
250  std::stringstream ss;
251  for(Teuchos_Ordinal i=0; i<avatarParameterValues_.size(); i++) {
252  ss << "," << avatarParameterValues_[i][indices[i]];
253  }
254  return ss.str();
255 }
256 
257 // ***********************************************************************
258 void AvatarInterface::SetIndices(int id,std::vector<int> & indices) const {
259  // The combo numbering here starts with the first guy
260  int numParams = (int)avatarParameterValues_.size();
261  int curr_id = id;
262  for(int i=0; i<numParams; i++) {
263  int div = avatarParameterValues_[i].size();
264  int mod = curr_id % div;
265  indices[i] = mod;
266  curr_id = (curr_id - mod)/div;
267  }
268 }
269 
270 
271 
272 // ***********************************************************************
273 void AvatarInterface::GenerateMueLuParametersFromIndex(int id,Teuchos::ParameterList & pl) const {
274  // The combo numbering here starts with the first guy
275  int numParams = (int)avatarParameterValues_.size();
276  int curr_id = id;
277  for(int i=0; i<numParams; i++) {
278  int div = avatarParameterValues_[i].size();
279  int mod = curr_id % div;
280  pl.set(mueluParameterName_[i],mueluParameterValues_[i][mod]);
281  curr_id = (curr_id - mod)/div;
282  }
283 }
284 
285 
286 // ***********************************************************************
287 void AvatarInterface::SetMueLuParameters(const Teuchos::ParameterList & problemFeatures, Teuchos::ParameterList & mueluParams, bool overwrite) const {
288  Teuchos::ParameterList avatarParams;
289  std::string paramString;
290 
291  if (comm_->getRank() == 0) {
292  // Only Rank 0 calls Avatar
293  if(!avatarHandle_) throw std::runtime_error("MueLu::AvatarInterface::SetMueLuParameters(): Setup has not been run");
294 
295  // Turn the problem features into a "trial" string for Avatar
296  std::string trialString;
297  GenerateFeatureString(problemFeatures,trialString);
298 
299  // Compute the number of things we need to test
300  int numParams = (int)avatarParameterValues_.size();
301  std::vector<int> indices(numParams);
302  std::vector<int> sizes(numParams);
303  int num_combos = 1;
304  for(int i=0; i<numParams; i++) {
305  sizes[i] = avatarParameterValues_[i].size();
306  num_combos *= avatarParameterValues_[i].size();
307  }
308  GetOStream(Runtime0)<< "MueLu::AvatarInterface: Testing "<< num_combos << " option combinations"<<std::endl;
309 
310  // For each input parameter to avatar we iterate over its allowable values and then compute the list of options which Avatar
311  // views as acceptable
312  // FIXME: Find alternative to hard coding malloc size (input deck?)
313  int* predictions = (int*)malloc(8 * sizeof(int));
314  float* probabilities = (float*)malloc(3 * 8 * sizeof(float));
315 
316  std::string testString;
317  for(int i=0; i<num_combos; i++) {
318  SetIndices(i,indices);
319  // Now we add the MueLu parameters into one, enormous Avatar trial string and run avatar once
320  testString += trialString + ParamsToString(indices) + ",0\n";
321  }
322 
323  std::cout<<"** Avatar TestString ***\n"<<testString<<std::endl;//DEBUG
324 
325  int bound_check = checkBounds(testString, boundsString_);
326 
327  // FIXME: Only send in first tree's string
328  //int* avatar_test(Avatar_handle* a, char* test_data_file, int test_data_is_a_string);
329  const int test_data_is_a_string = 1;
330  avatar_test(avatarHandle_,const_cast<char*>(testString.c_str()),test_data_is_a_string,predictions,probabilities);
331 
332  // Look at the list of acceptable combinations of options
333  std::vector<int> acceptableCombos; acceptableCombos.reserve(100);
334  for(int i=0; i<num_combos; i++) {
335  if(predictions[i] == avatarGoodClass_) acceptableCombos.push_back(i);
336  }
337  GetOStream(Runtime0)<< "MueLu::AvatarInterface: "<< acceptableCombos.size() << " acceptable option combinations found"<<std::endl;
338 
339  // Did we have any good combos at all?
340  int chosen_option_id = 0;
341  if(acceptableCombos.size() == 0) {
342  GetOStream(Runtime0) << "WARNING: MueLu::AvatarInterface: found *no* combinations of options which it believes will perform well on this problem" <<std::endl
343  << " An arbitrary set of options will be chosen instead"<<std::endl;
344  }
345  else {
346  // If there is only one acceptable combination, use it;
347  // otherwise, find the parameter choice with the highest
348  // probability of success
349  if(acceptableCombos.size() == 1){
350  chosen_option_id = acceptableCombos[0];
351  }
352  else {
353  switch (heuristicToUse_){
354  case 1:
355  chosen_option_id = hybrid(probabilities, acceptableCombos);
356  break;
357  case 2:
358  chosen_option_id = highProb(probabilities, acceptableCombos);
359  break;
360  case 3:
361  // Choose the first option in the list of acceptable
362  // combinations; the lowest drop tolerance among the
363  // acceptable combinations
364  chosen_option_id = acceptableCombos[0];
365  break;
366  case 4:
367  chosen_option_id = lowCrash(probabilities, acceptableCombos);
368  break;
369  case 5:
370  chosen_option_id = weighted(probabilities, acceptableCombos);
371  break;
372  }
373 
374  }
375  }
376 
377  // If mesh parameters are outside bounding box, set drop tolerance
378  // to 0, otherwise use avatar recommended drop tolerance
379  if (bound_check == 0){
380  GetOStream(Runtime0) << "WARNING: Extrapolation risk detected, setting drop tolerance to 0" <<std::endl;
381  GenerateMueLuParametersFromIndex(0,avatarParams);
382  } else {
383  GenerateMueLuParametersFromIndex(chosen_option_id,avatarParams);
384  }
385 
386  // Cleanup
387  free(predictions);
388  free(probabilities);
389  }
390 
391  Teuchos::updateParametersAndBroadcast(outArg(avatarParams),outArg(mueluParams),*comm_,0,overwrite);
392 
393 
394 }
395 
396 int AvatarInterface::checkBounds(std::string trialString, Teuchos::ArrayRCP<std::string> boundsString_) const {
397  std::stringstream ss(trialString);
398  std::vector<float> vect;
399 
400  int useNewFeatures = 0;
401 
402  float i;
403 
404  while (ss >> i)
405  {
406  vect.push_back(i);
407 
408  if (ss.peek() == ',')
409  ss.ignore();
410  }
411 
412  std::string bounds = const_cast<char*>(boundsString_[0].c_str());
413 
414  std::stringstream ssBounds(bounds);
415  std::vector<float> boundsVect;
416 
417  float b;
418 
419  while (ssBounds >> b)
420  {
421  boundsVect.push_back(b);
422 
423  if (ssBounds.peek() == ',')
424  ssBounds.ignore();
425  }
426 
427  if (vect.at(3) > boundsVect.at(0) || vect.at(3) < boundsVect.at(1))
428  return 0;
429 
430  if (vect.at(4) > boundsVect.at(2) || vect.at(4) < boundsVect.at(3))
431  return 0;
432 
433  if (vect.at(5) > boundsVect.at(4) || vect.at(5) < boundsVect.at(5))
434  return 0;
435 
436  if (vect.at(6) > boundsVect.at(6) || vect.at(6) < boundsVect.at(7))
437  return 0;
438 
439  if (useNewFeatures == 1){
440  if (vect.at(8) > boundsVect.at(8) || vect.at(8) < boundsVect.at(9))
441  return 0;
442 
443  if (vect.at(9) > boundsVect.at(10) || vect.at(9) < boundsVect.at(11))
444  return 0;
445  }
446 
447  return 1;
448 }
449 
450 int AvatarInterface::hybrid(float * probabilities, std::vector<int> acceptableCombos) const{
451  float low_crash = probabilities[0];
452  float best_prob = probabilities[2];
453  float diff;
454  int this_combo;
455  int chosen_option_id = acceptableCombos[0];
456  for(int x=0; x<acceptableCombos.size(); x++){
457  this_combo = acceptableCombos[x] * 3;
458  diff = probabilities[this_combo] - low_crash;
459  // If this parameter combination has a crash
460  // probability .2 lower than the current "best", we
461  // will use this drop tolerance
462  if(diff < -.2){
463  low_crash = probabilities[this_combo];
464  best_prob = probabilities[this_combo + 2];
465  chosen_option_id = acceptableCombos[x];
466  }
467  // If this parameter combination has the same
468  // or slightly lower crash probability than the
469  // current best, we compare their "GOOD" probabilities
470  else if(diff <= 0 && probabilities[this_combo + 2] > best_prob){
471  low_crash = probabilities[this_combo];
472  best_prob = probabilities[this_combo + 2];
473  chosen_option_id = acceptableCombos[x];
474  }
475  }
476  return chosen_option_id;
477 }
478 
479 int AvatarInterface::highProb(float * probabilities, std::vector<int> acceptableCombos) const{
480  float high_prob = probabilities[2];
481  int this_combo;
482  int chosen_option_id = acceptableCombos[0];
483  for(int x=0; x<acceptableCombos.size(); x++){
484  this_combo = acceptableCombos[x] * 3;
485  // If this parameter combination has a higher "GOOD"
486  // probability, use this combination
487  if(probabilities[this_combo + 2] > high_prob){
488  high_prob = probabilities[this_combo + 2];
489  chosen_option_id = acceptableCombos[x];
490  }
491  }
492  return chosen_option_id;
493 }
494 
495 int AvatarInterface::lowCrash(float * probabilities, std::vector<int> acceptableCombos) const{
496  float low_crash = probabilities[0];
497  int this_combo;
498  int chosen_option_id = acceptableCombos[0];
499  for(int x=0; x<acceptableCombos.size(); x++){
500  this_combo = acceptableCombos[x] * 3;
501  // If this parameter combination has a lower "CRASH"
502  // probability, use this combination
503  if(probabilities[this_combo] < low_crash){
504  low_crash = probabilities[this_combo];
505  chosen_option_id = acceptableCombos[x];
506  }
507  }
508  return chosen_option_id;
509 }
510 
511 int AvatarInterface::weighted(float * probabilities, std::vector<int> acceptableCombos) const{
512  float low_crash = probabilities[0];
513  float best_prob = probabilities[2];
514  float diff;
515  int this_combo;
516  int chosen_option_id = acceptableCombos[0];
517  for(int x=0; x<acceptableCombos.size(); x++){
518  this_combo = acceptableCombos[x] * 3;
519  diff = probabilities[this_combo] - low_crash;
520  // If this parameter combination has a crash
521  // probability .2 lower than the current "best", we
522  // will use this drop tolerance
523  if(diff < -.2){
524  low_crash = probabilities[this_combo];
525  best_prob = probabilities[this_combo + 2];
526  chosen_option_id = acceptableCombos[x];
527  }
528  // If this parameter combination is within .1
529  // or has a slightly lower crash probability than the
530  // current best, we compare their "GOOD" probabilities
531  else if(diff <= .1 && probabilities[this_combo + 2] > best_prob){
532  low_crash = probabilities[this_combo];
533  best_prob = probabilities[this_combo + 2];
534  chosen_option_id = acceptableCombos[x];
535  }
536  }
537  return chosen_option_id;
538 }
539 
540 
541 }// namespace MueLu
542 
543 #endif// HAVE_MUELU_AVATAR
544 
ConstIterator begin() 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)
const ParameterEntry & entry(ConstIterator i) const
bool isSublist(const std::string &name) const
Ordinal numParams() const
One-liner description of what is happening.
Namespace for MueLu classes and methods.
ConstIterator end() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_type size() const
std::ostream & leftshift(std::ostream &os, bool printFlags=true) const
params_t::ConstIterator ConstIterator
void resize(size_type new_size, const value_type &x=value_type())
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")