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