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 "SundanceFunctionalPolynomial.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "SundanceExpr.hpp"
00035 #include "SundanceEvaluatorFactory.hpp"
00036 #include "SundanceEvaluator.hpp"
00037 #include "SundanceSymbolicFuncElement.hpp"
00038 #include "SundanceDerivOfSymbFunc.hpp"
00039 #include "SundanceUnaryExpr.hpp"
00040
00041
00042 using namespace Sundance;
00043 using namespace Sundance;
00044
00045 using namespace Sundance;
00046 using namespace Teuchos;
00047
00048
00049
00050 FunctionalPolynomial::FunctionalPolynomial(const RCP<ScalarExpr>& expr)
00051 : funcs_(),
00052 funcMultiIndices_(),
00053 coeffs_(),
00054 keys_()
00055 {
00056 TEUCHOS_TEST_FOR_EXCEPTION(!isConvertibleToPoly(expr.get()), std::logic_error,
00057 "FunctionalPolynomial ctor called with an argument that "
00058 "cannot be converted to a polynomial");
00059 int funcID;
00060 MultiIndex mi;
00061 Deriv deriv;
00062
00063 const SymbolicFuncElement* s
00064 = dynamic_cast<const SymbolicFuncElement*>(expr.get());
00065 if (s != 0)
00066 {
00067 mi = MultiIndex();
00068 deriv = funcDeriv(s);
00069 }
00070
00071 const DerivOfSymbFunc* d
00072 = dynamic_cast<const DerivOfSymbFunc*>(expr.get());
00073 if (d != 0)
00074 {
00075 deriv = d->representMeAsFunctionalDeriv();
00076 }
00077
00078 MultipleDeriv mu;
00079 mu.put(deriv);
00080
00081
00082 if (d != 0 || s != 0)
00083 {
00084 funcs_.put(funcID, expr);
00085 Set<MultiIndex> miSet;
00086 miSet.put(mi);
00087 funcMultiIndices_.put(funcID, miSet);
00088
00089 Expr coeff = 1.0;
00090 RCP<ScalarExpr> cPtr = rcp_dynamic_cast<ScalarExpr>(coeff.ptr());
00091
00092 coeffs_.resize(2);
00093 keys_.resize(2);
00094 coeffs_[1].put(mu, cPtr);
00095 keys_[1].put(mu);
00096 }
00097 else
00098 {
00099 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00100 "impossible case in FunctionalPolynomial ctor");
00101 }
00102 }
00103
00104
00105 FunctionalPolynomial::FunctionalPolynomial(const Map<int, RCP<ScalarExpr> >& funcs,
00106 const Map<int, Set<MultiIndex> >& funcMultiIndices,
00107 const Array<Map<MultipleDeriv, RCP<ScalarExpr> > > & coeffs)
00108 : funcs_(funcs),
00109 funcMultiIndices_(funcMultiIndices),
00110 coeffs_(coeffs),
00111 keys_(coeffs.size())
00112 {
00113 typedef Map<MultipleDeriv, RCP<ScalarExpr> > termMap;
00114
00115 for (int i=0; i < coeffs_.size(); i++)
00116 {
00117 for (termMap::const_iterator
00118 j = coeffs_[i].begin(); j != coeffs_[i].end(); j++)
00119 {
00120 const MultipleDeriv& key = j->first;
00121 keys_[i].put(key);
00122 }
00123 }
00124 }
00125
00126
00127 Set<MultipleDeriv>
00128 FunctionalPolynomial::internalFindW(int order, const EvalContext& context) const
00129 {
00130 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00131 "FunctionalPolynomial not implemented");
00132 return Set<MultipleDeriv> ();
00133 }
00134
00135
00136 RCP<FunctionalPolynomial> FunctionalPolynomial::
00137 addPoly(const FunctionalPolynomial* other, int sign) const
00138 {
00139 typedef Map<MultipleDeriv, RCP<ScalarExpr> > termMap;
00140 Map<int, RCP<ScalarExpr> > newFuncs = funcs_;
00141 Map<int, Set<MultiIndex> > newFuncMultiIndices = funcMultiIndices_;
00142 Array<Map<MultipleDeriv, RCP<ScalarExpr> > > newCoeffs = coeffs_;
00143
00144 if (other->coeffs_.size() > coeffs_.size())
00145 newCoeffs.resize(other->coeffs_.size());
00146
00147
00148 for (int i=0; i < other->coeffs_.size(); i++)
00149 {
00150
00151 for (termMap::const_iterator
00152 j = other->coeffs_[i].begin(); j != other->coeffs_[i].end(); j++)
00153 {
00154 const MultipleDeriv& key = j->first;
00155 Expr right = Expr::handle(j->second);
00156 Expr sum;
00157
00158 if (newCoeffs[i].containsKey(key))
00159 {
00160 Expr left = Expr::handle(newCoeffs[i].get(key));
00161
00162 if (sign > 0) sum = left + right;
00163 else sum = left - right;
00164 }
00165 else
00166 {
00167 if (sign > 0) sum = right;
00168 else sum = -right;
00169 }
00170
00171 const ScalarExpr* se = dynamic_cast<const ScalarExpr*>(sum.ptr().get());
00172 TEUCHOS_TEST_FOR_EXCEPTION(se==0, std::logic_error,
00173 "Sum could not be cast to scalar expression");
00174 newCoeffs[i].put(key, rcp_dynamic_cast<ScalarExpr>(sum.ptr()));
00175 }
00176 }
00177
00178 for (Map<int, RCP<ScalarExpr> >::const_iterator
00179 i = other->funcs_.begin(); i != other->funcs_.end(); i++)
00180 {
00181 newFuncs.put(i->first, i->second);
00182 }
00183
00184 for (Map<int, Set<MultiIndex> >::const_iterator
00185 i = other->funcMultiIndices_.begin();
00186 i != other->funcMultiIndices_.end(); i++)
00187 {
00188 newFuncMultiIndices.put(i->first, i->second);
00189 }
00190
00191
00192
00193 return rcp(new FunctionalPolynomial(newFuncs, newFuncMultiIndices, newCoeffs));
00194 }
00195
00196 RCP<FunctionalPolynomial> FunctionalPolynomial::
00197 multiplyPoly(const FunctionalPolynomial* other) const
00198 {
00199 typedef Map<MultipleDeriv, RCP<ScalarExpr> > termMap;
00200 Map<int, RCP<ScalarExpr> > newFuncs;
00201 Map<int, Set<MultiIndex> > newFuncMultiIndices;
00202 Array<Map<MultipleDeriv, RCP<ScalarExpr> > > newCoeffs;
00203
00204 newCoeffs.resize(coeffs_.size() + other->coeffs_.size() - 1);
00205
00206 for (int i=0; i < coeffs_.size(); i++)
00207 {
00208 for (int j = 0; j<other->coeffs_.size(); j++)
00209 {
00210 for (termMap::const_iterator
00211 me = coeffs_[i].begin(); me != coeffs_[i].end(); me++)
00212 {
00213 const MultipleDeriv& myKey = me->first;
00214 Expr myExpr = Expr::handle(me->second);
00215 for (termMap::const_iterator
00216 you = other->coeffs_[j].begin();
00217 you != other->coeffs_[j].end(); you++)
00218 {
00219 const MultipleDeriv& yourKey = you->first;
00220 Expr yourExpr = Expr::handle(you->second);
00221 MultipleDeriv newKey = myKey.product(yourKey);
00222 int order = i+j;
00223 Expr newTerm = myExpr*yourExpr;
00224 if (newCoeffs[order].containsKey(newKey))
00225 {
00226 Expr old = Expr::handle(newCoeffs[i].get(newKey));
00227 newTerm = newTerm + old;
00228 }
00229 const ScalarExpr* se
00230 = dynamic_cast<const ScalarExpr*>(newTerm.ptr().get());
00231 TEUCHOS_TEST_FOR_EXCEPTION(se==0, std::logic_error,
00232 "New coeff could not be cast to scalar expression");
00233 newCoeffs[order].put(newKey,
00234 rcp_dynamic_cast<ScalarExpr>(newTerm.ptr()));
00235 }
00236 }
00237 }
00238 }
00239
00240 for (Map<int, RCP<ScalarExpr> >::const_iterator
00241 i = funcs_.begin(); i != funcs_.end(); i++)
00242 {
00243 newFuncs.put(i->first, i->second);
00244 }
00245 for (Map<int, RCP<ScalarExpr> >::const_iterator
00246 i = other->funcs_.begin(); i != other->funcs_.end(); i++)
00247 {
00248 newFuncs.put(i->first, i->second);
00249 }
00250
00251 for (Map<int, Set<MultiIndex> >::const_iterator
00252 i = funcMultiIndices_.begin();
00253 i != funcMultiIndices_.end(); i++)
00254 {
00255 newFuncMultiIndices.put(i->first, i->second);
00256 }
00257 for (Map<int, Set<MultiIndex> >::const_iterator
00258 i = other->funcMultiIndices_.begin();
00259 i != other->funcMultiIndices_.end(); i++)
00260 {
00261 Set<MultiIndex> miSet = i->second;
00262 if (newFuncMultiIndices.containsKey(i->first))
00263 {
00264 miSet.merge(newFuncMultiIndices.get(i->first));
00265 }
00266 newFuncMultiIndices.put(i->first, miSet);
00267 }
00268
00269 return rcp(new FunctionalPolynomial(newFuncs, newFuncMultiIndices, newCoeffs));
00270 }
00271
00272 RCP<FunctionalPolynomial> FunctionalPolynomial::
00273 addFunction(const RCP<ScalarExpr>& u, int sign) const
00274 {
00275 RCP<FunctionalPolynomial> other = rcp(new FunctionalPolynomial(u));
00276 return addPoly(other.get(), sign);
00277 }
00278
00279 RCP<FunctionalPolynomial> FunctionalPolynomial::
00280 multiplyScalar(const RCP<ScalarExpr>& alpha) const
00281 {
00282 typedef Map<MultipleDeriv, RCP<ScalarExpr> > termMap;
00283
00284 Array<Map<MultipleDeriv, RCP<ScalarExpr> > > newCoeffs = coeffs_;
00285
00286 Expr alphaExpr = Expr::handle(alpha);
00287
00288 for (int i=0; i < coeffs_.size(); i++)
00289 {
00290 for (termMap::const_iterator
00291 j = coeffs_[i].begin(); j != coeffs_[i].end(); j++)
00292 {
00293 const MultipleDeriv& key = j->first;
00294 Expr oldCoeff = Expr::handle(j->second);
00295 Expr newCoeff = alphaExpr*oldCoeff;
00296
00297 const ScalarExpr* se
00298 = dynamic_cast<const ScalarExpr*>(newCoeff.ptr().get());
00299 TEUCHOS_TEST_FOR_EXCEPTION(se==0, std::logic_error,
00300 "Coefficient could not be cast to "
00301 "scalar expression");
00302
00303 newCoeffs[i].put(key, rcp_dynamic_cast<ScalarExpr>(newCoeff.ptr()));
00304 }
00305 }
00306
00307 return rcp(new FunctionalPolynomial(funcs_, funcMultiIndices_, newCoeffs));
00308 }
00309
00310
00311 Evaluator* FunctionalPolynomial::createEvaluator(const EvaluatableExpr* expr,
00312 const EvalContext& context) const
00313 {
00314 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "poly eval not ready");
00315 return (Evaluator*) 0;
00316 }
00317
00318
00319
00320 bool FunctionalPolynomial::isConvertibleToPoly(const ScalarExpr* expr)
00321 {
00322 const SymbolicFuncElement* s
00323 = dynamic_cast<const SymbolicFuncElement*>(expr);
00324
00325 const DerivOfSymbFunc* d
00326 = dynamic_cast<const DerivOfSymbFunc*>(expr);
00327
00328 const FunctionalPolynomial* p
00329 = dynamic_cast<const FunctionalPolynomial*>(expr);
00330
00331 return (s != 0 || d != 0 || p != 0);
00332 }
00333
00334
00335 RCP<FunctionalPolynomial> FunctionalPolynomial::toPoly(const RCP<ScalarExpr>& expr)
00336 {
00337 const FunctionalPolynomial* p
00338 = dynamic_cast<const FunctionalPolynomial*>(expr.get());
00339
00340 if (p != 0) return rcp_dynamic_cast<FunctionalPolynomial>(expr);
00341
00342 Expr rtn = new FunctionalPolynomial(expr);
00343 return rcp_dynamic_cast<FunctionalPolynomial>(rtn.ptr());
00344 }
00345
00346
00347 std::ostream& FunctionalPolynomial::toText(std::ostream& os, bool paren) const
00348 {
00349 os << evalString();
00350 return os;
00351 }
00352
00353
00354 XMLObject FunctionalPolynomial::toXML() const
00355 {
00356 XMLObject rtn("Polynomial");
00357 for (int order=0; order<coeffs_.size(); order++)
00358 {
00359 for (Map<MultipleDeriv, RCP<ScalarExpr> >::const_iterator
00360 i = coeffs_[order].begin(); i != coeffs_[order].end(); i++)
00361 {
00362 const MultipleDeriv& key = i->first;
00363 Expr e = Expr::handle(i->second);
00364 XMLObject term("Term");
00365 XMLObject coeff("Coeff");
00366 coeff.addChild(e.toXML());
00367 term.addChild(coeff);
00368 for (MultipleDeriv::const_iterator
00369 j = key.begin(); j != key.end(); j++)
00370 {
00371 XMLObject f("FunctionalDeriv");
00372 f.addAttribute("mu", j->toString());
00373 term.addChild(f);
00374 }
00375 rtn.addChild(term);
00376 }
00377 }
00378 return rtn;
00379 }
00380
00381 Set<Deriv> FunctionalPolynomial
00382 ::findFuncsForSummation(const Set<MultipleDeriv>& prevSet,
00383 const MultipleDeriv& currentDeriv) const
00384 {
00385 Set<Deriv> rtn;
00386
00387 for (Set<MultipleDeriv>::const_iterator
00388 i = prevSet.begin(); i != prevSet.end(); i++)
00389 {
00390 const MultipleDeriv& mdPrev = *i;
00391 TEUCHOS_TEST_FOR_EXCEPTION(currentDeriv.size()+1 != mdPrev.size(),
00392 std::logic_error,
00393 "deriv orders must differ by 1. Found "
00394 "currentDeriv.size()=" << currentDeriv.size()
00395 << " while mdPrev.size()=" << mdPrev.size());
00396
00397
00398
00399
00400
00401
00402 Set<Deriv> tmp;
00403 set_difference(mdPrev.begin(), mdPrev.end(),
00404 currentDeriv.begin(), currentDeriv.end(),
00405 inserter(tmp, tmp.begin()));
00406 if (tmp.size()==1)
00407 {
00408 const Deriv& d = *(tmp.begin());
00409 if (currentDeriv.upper_bound(d) == currentDeriv.end()) rtn.put(d);
00410 }
00411 }
00412 return rtn;
00413 }
00414
00415
00416 MultipleDeriv FunctionalPolynomial::successorTerm(const MultipleDeriv& md) const
00417 {
00418 MultipleDeriv rtn;
00419
00420 unsigned int k = 0;
00421 for (MultipleDeriv::const_iterator i=md.begin(); i!=md.end(); i++, k++)
00422 {
00423 if (k < md.size()-1) rtn.put(*i);
00424 }
00425 return rtn;
00426 }
00427
00428 void FunctionalPolynomial
00429 ::stepRecurrence(int level, const Map<MultipleDeriv, std::string>& sPrev,
00430 Map<MultipleDeriv, std::string>& sCurr) const
00431 {
00432 Set<MultipleDeriv> allKeys;
00433 Set<MultipleDeriv> inducedKeys;
00434 Set<MultipleDeriv> prevKeys;
00435 for (Map<MultipleDeriv, std::string>::const_iterator
00436 i = sPrev.begin(); i != sPrev.end(); i++)
00437 {
00438 inducedKeys.put(successorTerm(i->first));
00439 }
00440 for (Map<MultipleDeriv, std::string>::const_iterator
00441 i = sPrev.begin(); i != sPrev.end(); i++)
00442 {
00443 prevKeys.put(i->first);
00444 }
00445
00446 Out::os() << "keys from prev level are " << prevKeys << std::endl;
00447 Out::os() << "induced keys are " << inducedKeys << std::endl;
00448 Out::os() << "natural keys are " << keys_[level] << std::endl;
00449
00450 allKeys.merge(inducedKeys);
00451 allKeys.merge(keys_[level]);
00452
00453 Out::os() << "keys active at this level are " << allKeys << std::endl;
00454
00455 for (Set<MultipleDeriv>::const_iterator
00456 i = allKeys.begin(); i != allKeys.end(); i++)
00457 {
00458 const MultipleDeriv& key = *i;
00459 Out::os() << "working on key " << key << std::endl;
00460 std::string str;
00461 if (coeffs_[level].containsKey(key))
00462 {
00463 str = coeffs_[level].get(key)->toString();
00464 }
00465
00466 Set<Deriv> funcs = findFuncsForSummation(prevKeys, key);
00467 Out::os() << "funcs to sum are " << funcs << std::endl;
00468 for (Set<Deriv>::const_iterator
00469 j = funcs.begin(); j != funcs.end(); j++)
00470 {
00471 MultipleDeriv prevKey = key;
00472 Out::os() << "prev key = " << prevKey << std::endl;
00473 Out::os() << "func = " << *j << std::endl;
00474
00475
00476
00477
00478
00479 prevKey.put(*j);
00480 TEUCHOS_TEST_FOR_EXCEPTION(!sPrev.containsKey(prevKey), std::logic_error,
00481 "inconsisent key lookup");
00482 const std::string& prevStr = sPrev.get(prevKey);
00483 std::string funcStr = (*j).toString();
00484 if (str.length() > 0) str += " + ";
00485 str += funcStr + "*(" + prevStr + ")";
00486 }
00487 if (str.length() > 0) sCurr.put(key, str);
00488 }
00489 }
00490
00491 string FunctionalPolynomial::evalString() const
00492 {
00493 int maxLevel = coeffs_.size()-1;
00494
00495 Map<MultipleDeriv, std::string> sPrev;
00496 Map<MultipleDeriv, std::string> sCurr;
00497
00498 for (int level=maxLevel; level>=0; level--)
00499 {
00500 Out::os() << "********* Recurrence level = " << level << " ***************"
00501 << std::endl;
00502 sCurr = Map<MultipleDeriv, std::string>();
00503 stepRecurrence(level, sPrev, sCurr);
00504 sPrev = sCurr;
00505
00506 for (Map<MultipleDeriv, std::string>::const_iterator
00507 j = sCurr.begin(); j != sCurr.end(); j++)
00508 {
00509 Out::os() << "key=" << j->first << std::endl;
00510 Out::os() << "value=" << j->second << std::endl;
00511 }
00512 }
00513
00514
00515
00516
00517 return sCurr.begin()->second;
00518 }
00519
00520
00521
00522 bool FunctionalPolynomial::lessThan(const ScalarExpr* other) const
00523 {
00524 const FunctionalPolynomial* f = dynamic_cast<const FunctionalPolynomial*>(other);
00525 TEUCHOS_TEST_FOR_EXCEPTION(f==0, std::logic_error, "cast should never fail at this point");
00526
00527 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "FunctionalPolynomial::lessThan() not "
00528 "implemented");
00529 }
00530
00531