MueLu Version of the Day
Loading...
Searching...
No Matches
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/*
87TODO List:
88Modify 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
94extern "C" {
95#include "avatar_api.h"
96}
97
98
99
100namespace MueLu {
101
102
103// ***********************************************************************
104RCP<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// ***********************************************************************
152Teuchos::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// ***********************************************************************
174void 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// ***********************************************************************
215void AvatarInterface::Cleanup() {
216 avatar_cleanup(avatarHandle_);
217 avatarHandle_=0;
218}
219
220
221// ***********************************************************************
222void 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// ***********************************************************************
245void 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// ***********************************************************************
284std::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// ***********************************************************************
302void 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// ***********************************************************************
317void 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// ***********************************************************************
331void 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
439int 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
466int 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
495int 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
511int 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
527int 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
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.