00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceEvaluatableExpr.hpp"
00032 #include "SundanceEvaluatorFactory.hpp"
00033 #include "SundanceEvaluator.hpp"
00034 #include "SundanceNullEvaluator.hpp"
00035 #include "SundanceEvalManager.hpp"
00036 #include "SundanceSymbolicFuncElement.hpp"
00037 #include "SundanceExpr.hpp"
00038 #include "PlayaTabs.hpp"
00039 #include "SundanceOut.hpp"
00040 #include "Teuchos_Utils.hpp"
00041
00042 using namespace Sundance;
00043 using namespace Teuchos;
00044 using std::endl;
00045
00046 TEUCHOS_TIMER(evalTimer, "Symbolic Evaluation")
00047
00048 EvaluatableExpr::EvaluatableExpr()
00049 : ScalarExpr(),
00050 evaluators_(),
00051 sparsity_(),
00052 nodesHaveBeenCounted_(false),
00053 contextToDSSMap_(numDerivSubsetTypes(), contextToDSSMap_ele_t(maxFuncDiffOrder()+1)),
00054 rIsReady_(false),
00055 allOrdersMap_(numDerivSubsetTypes())
00056 {}
00057
00058
00059 Time& EvaluatableExpr::evaluationTimer()
00060 {
00061 return evalTimer();
00062 }
00063
00064 string EvaluatableExpr::derivType(const DerivSubsetSpecifier& dss) const
00065 {
00066 switch(dss)
00067 {
00068 case RequiredNonzeros:
00069 return "R";
00070 case ConstantNonzeros:
00071 return "C";
00072 case VariableNonzeros:
00073 return "V";
00074 default:
00075 return "W";
00076 }
00077 }
00078
00079 void EvaluatableExpr::registerSpatialDerivs(const EvalContext& context,
00080 const Set<MultiIndex>& miSet) const
00081 {
00082 Tabs tab;
00083 if (activeSpatialDerivMap_.containsKey(context))
00084 {
00085 Set<MultiIndex> s = activeSpatialDerivMap_[context];
00086 s.merge(miSet);
00087 activeSpatialDerivMap_.put(context, s);
00088 }
00089 else
00090 {
00091 activeSpatialDerivMap_.put(context, miSet);
00092 }
00093 }
00094
00095 const Set<MultiIndex>& EvaluatableExpr
00096 ::activeSpatialDerivs(const EvalContext& context) const
00097 {
00098 TEUCHOS_TEST_FOR_EXCEPTION(!activeSpatialDerivMap_.containsKey(context),
00099 std::logic_error,
00100 "Unknown context " << context);
00101 return activeSpatialDerivMap_[context];
00102 }
00103
00104 RCP<SparsitySuperset>
00105 EvaluatableExpr::sparsitySuperset(const EvalContext& context) const
00106 {
00107 Tabs tab;
00108
00109 SUNDANCE_MSG2(context.setupVerbosity(),
00110 tab << "getting sparsity superset for " << toString());
00111
00112 RCP<SparsitySuperset> rtn;
00113
00114 if (sparsity_.containsKey(context))
00115 {
00116 Tabs tab1;
00117 SUNDANCE_MSG2(context.setupVerbosity(),
00118 tab1 << "reusing previously computed data...");
00119 rtn = sparsity_.get(context);
00120 }
00121 else
00122 {
00123 Tabs tab1;
00124 SUNDANCE_MSG2(context.setupVerbosity(),
00125 tab1 << "computing from scratch...");
00126 const Set<MultipleDeriv>& R = findR(context);
00127 const Set<MultipleDeriv>& C = findC(context);
00128 const Set<MultipleDeriv>& V = findV(context);
00129 if (context.setupVerbosity() > 4)
00130 {
00131 Out::os() << tab1 << "R=" << R << endl;
00132 Out::os() << tab1 << "C=" << C << endl;
00133 Out::os() << tab1 << "V=" << V << endl;
00134 }
00135 rtn = rcp(new SparsitySuperset(C.intersection(R), V.intersection(R)));
00136 sparsity_.put(context, rtn);
00137 }
00138 return rtn;
00139 }
00140
00141
00142
00143 const RCP<Evaluator>&
00144 EvaluatableExpr::evaluator(const EvalContext& context) const
00145 {
00146 TEUCHOS_TEST_FOR_EXCEPTION(!evaluators_.containsKey(context), std::runtime_error,
00147 "Evaluator not found for context " << context);
00148 return evaluators_.get(context);
00149 }
00150
00151 void EvaluatableExpr::evaluate(const EvalManager& mgr,
00152 Array<double>& constantResults,
00153 Array<RCP<EvalVector> >& vectorResults) const
00154 {
00155 TimeMonitor timer(evalTimer());
00156 evaluator(mgr.getRegion())->eval(mgr, constantResults, vectorResults);
00157 }
00158
00159 void EvaluatableExpr::setupEval(const EvalContext& context) const
00160 {
00161 Tabs tabs0(0);
00162 int verb = context.evalSetupVerbosity();
00163 SUNDANCE_MSG1(verb, tabs0 << "setupEval() for " << this->toString());
00164 if (!evaluators_.containsKey(context))
00165 {
00166 Tabs tabs;
00167 SUNDANCE_MSG2(verb, tabs << "creating new evaluator...");
00168 SUNDANCE_MSG2(verb, tabs << "my sparsity superset = " << std::endl
00169 << *sparsitySuperset(context));
00170 RCP<Evaluator> eval;
00171 if (sparsitySuperset(context)->numDerivs()>0)
00172 {
00173 SUNDANCE_MSG2(verb, tabs << "calling createEvaluator()");
00174 eval = rcp(createEvaluator(this, context));
00175 }
00176 else
00177 {
00178 SUNDANCE_MSG2(verb,
00179 tabs << "EE: no results needed... creating null evaluator");
00180 eval = rcp(new NullEvaluator());
00181 }
00182 evaluators_.put(context, eval);
00183 }
00184 else
00185 {
00186 Tabs tabs;
00187 SUNDANCE_MSG2(verb, tabs << "reusing existing evaluator...");
00188 }
00189 }
00190
00191 void EvaluatableExpr::showSparsity(std::ostream& os,
00192 const EvalContext& context) const
00193 {
00194 Tabs tab0;
00195 os << tab0 << "Node: " << toString() << std::endl;
00196 sparsitySuperset(context)->print(os);
00197 }
00198
00199
00200
00201 int EvaluatableExpr::maxOrder(const Set<MultiIndex>& m) const
00202 {
00203 int rtn = 0;
00204 for (Set<MultiIndex>::const_iterator i=m.begin(); i != m.end(); i++)
00205 {
00206 rtn = std::max(i->order(), rtn);
00207 }
00208 return rtn;
00209 }
00210
00211 const EvaluatableExpr* EvaluatableExpr::getEvalExpr(const Expr& expr)
00212 {
00213 const EvaluatableExpr* rtn
00214 = dynamic_cast<const EvaluatableExpr*>(expr[0].ptr().get());
00215 TEUCHOS_TEST_FOR_EXCEPTION(rtn==0, std::logic_error,
00216 "cast of " << expr
00217 << " failed in EvaluatableExpr::getEvalExpr()");
00218 TEUCHOS_TEST_FOR_EXCEPTION(expr.size() != 1, std::logic_error,
00219 "non-scalar expression " << expr
00220 << " in EvaluatableExpr::getEvalExpr()");
00221
00222 return rtn;
00223 }
00224
00225
00226 bool EvaluatableExpr::isEvaluatable(const ExprBase* expr)
00227 {
00228 return dynamic_cast<const EvaluatableExpr*>(expr) != 0;
00229 }
00230
00231
00232 int EvaluatableExpr::countNodes() const
00233 {
00234 nodesHaveBeenCounted_ = true;
00235 return 1;
00236 }
00237
00238
00239
00240 const Set<MultipleDeriv>&
00241 EvaluatableExpr::findW(int order,
00242 const EvalContext& context) const
00243 {
00244 return findDerivSubset(order, AllNonzeros, context);
00245 }
00246
00247 const Set<MultipleDeriv>&
00248 EvaluatableExpr::findV(int order,
00249 const EvalContext& context) const
00250 {
00251 return findDerivSubset(order, VariableNonzeros, context);
00252 }
00253
00254 const Set<MultipleDeriv>&
00255 EvaluatableExpr::findC(int order,
00256 const EvalContext& context) const
00257 {
00258 return findDerivSubset(order, ConstantNonzeros, context);
00259 }
00260
00261 const Set<MultipleDeriv>&
00262 EvaluatableExpr::findR(int order,
00263 const EvalContext& context) const
00264 {
00265 TEUCHOS_TEST_FOR_EXCEPTION(!rIsReady_, std::logic_error,
00266 "findR() cannot be used for initial computation of the "
00267 "R subset. Calling object is " << toString());
00268 return findDerivSubset(order, RequiredNonzeros, context);
00269 }
00270
00271
00272 const Set<MultipleDeriv>&
00273 EvaluatableExpr::findDerivSubset(int order,
00274 const DerivSubsetSpecifier& dss,
00275 const EvalContext& context) const
00276 {
00277 Tabs tabs;
00278 int verb = context.setupVerbosity();
00279 SUNDANCE_MSG2(verb, tabs << "finding " << derivType(dss)
00280 << "[" << order << "] for " << toString());
00281 if (contextToDSSMap_[dss][order].containsKey(context))
00282 {
00283 Tabs tab1;
00284 SUNDANCE_MSG3(verb, tab1 << "...reusing previously computed data");
00285 }
00286 else
00287 {
00288 Tabs tab1;
00289 SUNDANCE_MSG3(verb, tab1 << "...need to call find<Set>");
00290 switch(dss)
00291 {
00292 case AllNonzeros:
00293 contextToDSSMap_[dss][order].put(context, internalFindW(order, context));
00294 break;
00295 case RequiredNonzeros:
00296 break;
00297 case VariableNonzeros:
00298 contextToDSSMap_[dss][order].put(context, internalFindV(order, context));
00299 break;
00300 case ConstantNonzeros:
00301 contextToDSSMap_[dss][order].put(context, internalFindC(order, context));
00302 break;
00303 default:
00304 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "this should never happen");
00305 }
00306 }
00307
00308
00309 if (!contextToDSSMap_[dss][order].containsKey(context))
00310 {
00311 contextToDSSMap_[dss][order].put(context, Set<MultipleDeriv>());
00312 }
00313 const Set<MultipleDeriv>& rtn = contextToDSSMap_[dss][order].get(context);
00314 SUNDANCE_MSG2(verb, tabs << "found " << rtn);
00315 return rtn;
00316 }
00317
00318
00319 Set<MultipleDeriv>
00320 EvaluatableExpr::setProduct(const Set<MultipleDeriv>& a,
00321 const Set<MultipleDeriv>& b) const
00322 {
00323 Set<MultipleDeriv> rtn;
00324 for (Set<MultipleDeriv>::const_iterator i=a.begin(); i!=a.end(); i++)
00325 {
00326 for (Set<MultipleDeriv>::const_iterator j=b.begin(); j!=b.end(); j++)
00327 {
00328 rtn.put(i->product(*j));
00329 }
00330 }
00331 return rtn;
00332 }
00333
00334 Set<MultiSet<int> > EvaluatableExpr::setDivision(const Set<MultiSet<int> >& s,
00335 int index) const
00336 {
00337 Set<MultiSet<int> > rtn;
00338 for (Set<MultiSet<int> >::const_iterator
00339 i=s.begin(); i!=s.end(); i++)
00340 {
00341 if (i->contains(index))
00342 {
00343 MultiSet<int> e = *i;
00344 e.erase(e.find(index));
00345 rtn.put(e);
00346 }
00347 }
00348 return rtn;
00349 }
00350
00351
00352 Set<MultipleDeriv>
00353 EvaluatableExpr::setDivision(const Set<MultipleDeriv>& a,
00354 const Set<MultipleDeriv>& b) const
00355 {
00356 Set<MultipleDeriv> rtn;
00357 for (Set<MultipleDeriv>::const_iterator i=a.begin(); i!=a.end(); i++)
00358 {
00359 for (Set<MultipleDeriv>::const_iterator j=b.begin(); j!=b.end(); j++)
00360 {
00361 MultipleDeriv c = i->factorOutDeriv(*j);
00362 if (c.size() != 0) rtn.put(c);
00363 if (*i == *j) rtn.put(MultipleDeriv());
00364 }
00365 }
00366 return rtn;
00367 }
00368
00369 void EvaluatableExpr::determineR(const EvalContext& context,
00370 const Array<Set<MultipleDeriv> >& RInput) const
00371 {
00372 Tabs tabs;
00373 RCP<Array<Set<MultipleDeriv> > > additionToR
00374 = internalDetermineR(context, RInput);
00375
00376 for (int i=0; i<RInput.size(); i++)
00377 {
00378 if ((*additionToR)[i].size()==0) continue;
00379 if (!contextToDSSMap_[RequiredNonzeros][i].containsKey(context))
00380 {
00381 contextToDSSMap_[RequiredNonzeros][i].put(context, (*additionToR)[i]);
00382 }
00383 else
00384 {
00385 contextToDSSMap_[RequiredNonzeros][i][context].merge((*additionToR)[i]);
00386 }
00387 }
00388 rIsReady_ = true;
00389
00390 }
00391
00392 RCP<Array<Set<MultipleDeriv> > > EvaluatableExpr
00393 ::internalDetermineR(const EvalContext& context,
00394 const Array<Set<MultipleDeriv> >& RInput) const
00395 {
00396 Tabs tab0(0);
00397 int verb = context.setupVerbosity();
00398 RCP<Array<Set<MultipleDeriv> > > rtn
00399 = rcp(new Array<Set<MultipleDeriv> >(RInput.size()));
00400
00401 SUNDANCE_MSG2(verb, tab0 << "EE::internalDetermineR() for " << toString() );
00402 SUNDANCE_MSG2(verb, tab0 << "RInput = " << RInput );
00403
00404 for (int i=0; i<RInput.size(); i++)
00405 {
00406 Tabs tab1;
00407 if (RInput[i].size()==0) continue;
00408 const Set<MultipleDeriv>& Wi = findW(i, context);
00409 SUNDANCE_MSG3(verb, tab1 << "W[" << i << "] = " << Wi );
00410 (*rtn)[i] = RInput[i].intersection(Wi);
00411 }
00412
00413 printR(verb, rtn);
00414 return rtn;
00415 }
00416
00417
00418 Array<Set<MultipleDeriv> > EvaluatableExpr
00419 ::computeInputR(const EvalContext& context,
00420 const Array<Set<MultiSet<int> > >& funcIDCombinations,
00421 const Array<Set<MultiIndex> >& spatialDerivs) const
00422 {
00423 int verb = context.setupVerbosity();
00424 Tabs tabs(0);
00425 SUNDANCE_MSG2(verb, tabs << "in EE::computeInputR()");
00426 Array<Set<MultipleDeriv> > rtn(funcIDCombinations.size());
00427 for (int order=0; order<funcIDCombinations.size(); order++)
00428 {
00429 Tabs tabs1;
00430 SUNDANCE_MSG2(verb, tabs << "order=" << order);
00431 SUNDANCE_MSG2(verb, tabs1 << "funcID combs="
00432 << funcIDCombinations[order]);
00433 const Set<MultipleDeriv>& W = findW(order, context);
00434 SUNDANCE_MSG2(verb, tabs1 << "W[order=" << order << "]=" << W);
00435
00436 for (Set<MultipleDeriv>::const_iterator i=W.begin(); i!=W.end(); i++)
00437 {
00438 Tabs tab2;
00439 const MultipleDeriv& nonzeroDeriv = *i;
00440 if (nonzeroDeriv.isInRequiredSet(funcIDCombinations[order],
00441 spatialDerivs[order]))
00442 {
00443 SUNDANCE_MSG2(verb, tab2 << "md=" << nonzeroDeriv
00444 << " added to inputR");
00445 rtn[order].put(nonzeroDeriv);
00446 }
00447 else
00448 {
00449 SUNDANCE_MSG2(verb, tab2 << "md=" << nonzeroDeriv << " not needed");
00450 }
00451 }
00452 }
00453 SUNDANCE_MSG2(verb, tabs << "inputR=" << rtn);
00454 return rtn;
00455 }
00456
00457
00458 const Set<MultipleDeriv>& EvaluatableExpr::findW(const EvalContext& context) const
00459 {
00460 return findDerivSubset(AllNonzeros, context);
00461 }
00462
00463 const Set<MultipleDeriv>& EvaluatableExpr::findR(const EvalContext& context) const
00464 {
00465 return findDerivSubset(RequiredNonzeros, context);
00466 }
00467
00468 const Set<MultipleDeriv>& EvaluatableExpr::findC(const EvalContext& context) const
00469 {
00470 return findDerivSubset(ConstantNonzeros, context);
00471 }
00472
00473 const Set<MultipleDeriv>& EvaluatableExpr::findV(const EvalContext& context) const
00474 {
00475 return findDerivSubset(VariableNonzeros, context);
00476 }
00477
00478 const Set<MultipleDeriv>&
00479 EvaluatableExpr::findDerivSubset(const DerivSubsetSpecifier& dss,
00480 const EvalContext& context) const
00481 {
00482 if (!allOrdersMap_[dss].containsKey(context))
00483 {
00484 Set<MultipleDeriv> tmp;
00485 for (int i=0; i<=3; i++)
00486 {
00487 tmp.merge(findDerivSubset(i, dss, context));
00488 }
00489 allOrdersMap_[dss].put(context, tmp);
00490 }
00491
00492 return allOrdersMap_[dss].get(context);
00493 }
00494
00495 void EvaluatableExpr::displayNonzeros(std::ostream& os, const EvalContext& context) const
00496 {
00497 Tabs tabs0;
00498 os << tabs0 << "Nonzeros of " << toString() << std::endl;
00499
00500 const Set<MultipleDeriv>& W = findW(context);
00501 const Set<MultipleDeriv>& R = findR(context);
00502 const Set<MultipleDeriv>& C = findC(context);
00503 const Set<MultipleDeriv>& V = findV(context);
00504
00505 TEUCHOS_TEST_FOR_EXCEPT(C.intersection(V).size() != 0);
00506
00507 for (Set<MultipleDeriv>::const_iterator i=W.begin(); i != W.end(); i++)
00508 {
00509 Tabs tab1;
00510 std::string state = "Variable";
00511 if (C.contains(*i)) state = "Constant";
00512 if (!R.contains(*i)) state = "Not Required";
00513 os << tab1 << std::setw(25) << std::left << i->toString() << ": " << state << std::endl;
00514 }
00515 }
00516
00517
00518 void EvaluatableExpr::printR(int verb,
00519 const RCP<Array<Set<MultipleDeriv> > >& R) const
00520 {
00521 Tabs tab(0);
00522 const Array<Set<MultipleDeriv> >& r = *R;
00523 for (int order=0; order<r.size(); order++)
00524 {
00525 Tabs tab1;
00526 SUNDANCE_MSG2(verb, tab << "order=" << order);
00527 SUNDANCE_MSG2(verb, tab1 << "R[" << order << "]=" << r[order]);
00528 }
00529 }
