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