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 "SundanceExprWithChildren.hpp"
00032 #include "SundanceCombinatorialUtils.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceOut.hpp"
00035 #include "SundanceExpr.hpp"
00036 #include "SundanceEvaluatorFactory.hpp"
00037 #include "SundanceEvaluator.hpp"
00038 #include "SundanceNullEvaluator.hpp"
00039 #include "SundanceUnknownFuncElement.hpp"
00040 #include "SundanceTestFuncElement.hpp"
00041 #include "SundanceUnaryExpr.hpp"
00042
00043 using namespace Sundance;
00044 using namespace Sundance;
00045
00046 using namespace Sundance;
00047 using namespace Teuchos;
00048
00049
00050
00051
00052 ExprWithChildren::ExprWithChildren(const Array<RCP<ScalarExpr> >& children)
00053 : EvaluatableExpr(),
00054 children_(children),
00055 contextToQWMap_(4),
00056 contextToQVMap_(4),
00057 contextToQCMap_(4)
00058 {}
00059
00060 bool ExprWithChildren::lessThan(const ScalarExpr* other) const
00061 {
00062 const ExprWithChildren* c = dynamic_cast<const ExprWithChildren*>(other);
00063 TEUCHOS_TEST_FOR_EXCEPTION(c==0, std::logic_error, "cast should never fail at this point");
00064
00065 if (children_.size() < c->children_.size()) return true;
00066 if (children_.size() > c->children_.size()) return false;
00067
00068 for (int i=0; i<children_.size(); i++)
00069 {
00070 Expr me = Expr::handle(children_[i]);
00071 Expr you = Expr::handle(c->children_[i]);
00072 if (me.lessThan(you)) return true;
00073 if (you.lessThan(me)) return false;
00074 }
00075 return false;
00076 }
00077
00078
00079
00080 bool ExprWithChildren::isConstant() const
00081 {
00082 for (int i=0; i<children_.size(); i++)
00083 {
00084 if (!children_[i]->isConstant()) return false;
00085 }
00086 return true;
00087 }
00088
00089 bool ExprWithChildren::isIndependentOf(const Expr& u) const
00090 {
00091 for (int i=0; i<children_.size(); i++)
00092 {
00093 if (!children_[i]->isIndependentOf(u)) return false;
00094 }
00095 return true;
00096 }
00097
00098 const EvaluatableExpr* ExprWithChildren::evaluatableChild(int i) const
00099 {
00100 const EvaluatableExpr* e
00101 = dynamic_cast<const EvaluatableExpr*>(children_[i].get());
00102
00103 TEUCHOS_TEST_FOR_EXCEPTION(e==0, std::logic_error,
00104 "ExprWithChildren: cast of child ["
00105 << children_[i]->toString()
00106 << " to evaluatable expr failed");
00107
00108 return e;
00109 }
00110
00111 int ExprWithChildren::maxDiffOrderOnDiscreteFunctions() const
00112 {
00113 int biggest = -1;
00114 for (int i=0; i<numChildren(); i++)
00115 {
00116 int x = evaluatableChild(i)->maxDiffOrderOnDiscreteFunctions();
00117 if (x > biggest) biggest = x;
00118 }
00119 return biggest;
00120 }
00121
00122 bool ExprWithChildren::hasDiscreteFunctions() const
00123 {
00124 for (int i=0; i<numChildren(); i++)
00125 {
00126 if (evaluatableChild(i)->hasDiscreteFunctions()) return true;
00127 }
00128 return false;
00129 }
00130
00131
00132 void ExprWithChildren::accumulateFuncSet(Set<int>& funcIDs,
00133 const Set<int>& activeSet) const
00134 {
00135 for (int i=0; i<children_.size(); i++)
00136 {
00137 children_[i]->accumulateFuncSet(funcIDs, activeSet);
00138 }
00139 }
00140
00141 void ExprWithChildren::registerSpatialDerivs(const EvalContext& context,
00142 const Set<MultiIndex>& miSet) const
00143 {
00144 for (int i=0; i<children_.size(); i++)
00145 {
00146 evaluatableChild(i)->registerSpatialDerivs(context, miSet);
00147 }
00148 EvaluatableExpr::registerSpatialDerivs(context, miSet);
00149 }
00150
00151 void ExprWithChildren::setupEval(const EvalContext& context) const
00152 {
00153 Tabs tabs;
00154 int verb = context.evalSetupVerbosity();
00155 SUNDANCE_MSG1(verb, tabs << "expr " + toString()
00156 << ": creating evaluators for children");
00157 RCP<SparsitySuperset> ss = sparsitySuperset(context);
00158 SUNDANCE_MSG2(verb, tabs << "my sparsity superset = "
00159 << std::endl << *ss);
00160
00161
00162
00163
00164 if (ss->numDerivs() > 0)
00165 {
00166 for (int i=0; i<children_.size(); i++)
00167 {
00168 Tabs tabs1;
00169 SUNDANCE_MSG1(verb, tabs1 << "creating evaluator for child "
00170 << evaluatableChild(i)->toString());
00171 evaluatableChild(i)->setupEval(context);
00172 }
00173 }
00174
00175 if (!evaluators().containsKey(context))
00176 {
00177 RCP<Evaluator> eval;
00178 if (ss->numDerivs()>0)
00179 {
00180 eval = rcp(createEvaluator(this, context));
00181 }
00182 else
00183 {
00184 SUNDANCE_MSG2(verb,
00185 tabs << "EWC: no results needed... creating null evaluator");
00186 eval = rcp(new NullEvaluator());
00187 }
00188 evaluators().put(context, eval);
00189 }
00190 }
00191
00192 void ExprWithChildren::showSparsity(std::ostream& os,
00193 const EvalContext& context) const
00194 {
00195 Tabs tab0;
00196 os << tab0 << "Node: " << toString() << std::endl;
00197 sparsitySuperset(context)->print(os);
00198 for (int i=0; i<children_.size(); i++)
00199 {
00200 Tabs tab1;
00201 os << tab1 << "Child " << i << std::endl;
00202 evaluatableChild(i)->showSparsity(os, context);
00203 }
00204 }
00205
00206
00207 bool ExprWithChildren::everyTermHasTestFunctions() const
00208 {
00209 for (int i=0; i<children_.size(); i++)
00210 {
00211 if (evaluatableChild(i)->everyTermHasTestFunctions()) return true;
00212 }
00213 return false;
00214 }
00215
00216
00217 bool ExprWithChildren::hasTestFunctions() const
00218 {
00219 for (int i=0; i<children_.size(); i++)
00220 {
00221 if (evaluatableChild(i)->hasTestFunctions()) return true;
00222 }
00223 return false;
00224 }
00225
00226
00227 bool ExprWithChildren::hasUnkFunctions() const
00228 {
00229 for (int i=0; i<children_.size(); i++)
00230 {
00231 if (evaluatableChild(i)->hasUnkFunctions()) return true;
00232 }
00233 return false;
00234 }
00235
00236
00237 void ExprWithChildren::getUnknowns(Set<int>& unkID, Array<Expr>& unks) const
00238 {
00239 for (int i=0; i<children_.size(); i++)
00240 {
00241 const RCP<ExprBase>& e = children_[i];
00242 const UnknownFuncElement* u
00243 = dynamic_cast<const UnknownFuncElement*>(e.get());
00244 if (u != 0)
00245 {
00246 Expr expr(e);
00247 if (!unkID.contains(u->fid().dofID()))
00248 {
00249 unks.append(expr);
00250 unkID.put(u->fid().dofID());
00251 }
00252 }
00253 else
00254 {
00255 evaluatableChild(i)->getUnknowns(unkID, unks);
00256 }
00257 }
00258 }
00259
00260 void ExprWithChildren::getTests(Set<int>& varID, Array<Expr>& vars) const
00261 {
00262 for (int i=0; i<children_.size(); i++)
00263 {
00264 const RCP<ExprBase>& e = children_[i];
00265 const TestFuncElement* u
00266 = dynamic_cast<const TestFuncElement*>(e.get());
00267 if (u != 0)
00268 {
00269 Expr expr(e);
00270 if (!varID.contains(u->fid().dofID()))
00271 {
00272 vars.append(expr);
00273 varID.put(u->fid().dofID());
00274 }
00275
00276 }
00277 else
00278 {
00279 evaluatableChild(i)->getTests(varID, vars);
00280 }
00281 }
00282 }
00283
00284
00285 int ExprWithChildren::countNodes() const
00286 {
00287 if (nodesHaveBeenCounted())
00288 {
00289 return 0;
00290 }
00291
00292
00293 int count = EvaluatableExpr::countNodes();
00294
00295
00296 for (int i=0; i<children_.size(); i++)
00297 {
00298 if (!evaluatableChild(i)->nodesHaveBeenCounted())
00299 {
00300 count += evaluatableChild(i)->countNodes();
00301 }
00302 }
00303 return count;
00304 }
00305
00306 Set<MultipleDeriv>
00307 ExprWithChildren::product(const Array<int>& J, const Array<int>& K,
00308 DerivSubsetSpecifier dss,
00309 const EvalContext& context) const
00310 {
00311 TEUCHOS_TEST_FOR_EXCEPTION(J.size() != K.size(), std::logic_error,
00312 "mismatched index set sizes");
00313 TEUCHOS_TEST_FOR_EXCEPTION(J.size() == 0, std::logic_error,
00314 "empty index set");
00315
00316 Set<MultipleDeriv> rtn
00317 = evaluatableChild(J[0])->findDerivSubset(K[0], dss, context);
00318
00319 for (int i=1; i<J.size(); i++)
00320 {
00321 const Set<MultipleDeriv>& S
00322 = evaluatableChild(J[i])->findDerivSubset(K[i], dss, context);
00323 rtn = setProduct(rtn, S);
00324 }
00325
00326 return rtn;
00327 }
00328
00329 Set<MultipleDeriv>
00330 ExprWithChildren::internalFindV(int order, const EvalContext& context) const
00331 {
00332 Tabs tab0;
00333 int verb = context.setupVerbosity();
00334 SUNDANCE_MSG2(verb, tab0 << "EWC::internalFindV(" << order
00335 << ") for " << toString());
00336
00337 Set<MultipleDeriv> rtn;
00338
00339
00340
00341
00342 if (order==0)
00343 {
00344 Tabs tab1;
00345 SUNDANCE_MSG3(verb, tab1 << "case: order=0");
00346 for (int i=0; i<numChildren(); i++)
00347 {
00348 if (!childIsRequired(i, order, context)) continue;
00349 const Set<MultipleDeriv>& childV0
00350 = evaluatableChild(i)->findV(0, context);
00351 rtn.merge(childV0);
00352 }
00353 const Set<MultipleDeriv>& R0 = findR(0, context);
00354 rtn = rtn.intersection(R0);
00355
00356 SUNDANCE_MSG3(verb, tab0 << "V[" << order << "]=" << rtn);
00357 SUNDANCE_MSG3(verb, tab0 << "done with EWC::internalFindV(" << order
00358 << ") for " << toString());
00359 return rtn;
00360 }
00361
00362
00363 SUNDANCE_MSG3(verb, tab0 << "getting required set");
00364 const Set<MultipleDeriv>& RM = findR(order, context);
00365 SUNDANCE_MSG3(verb, tab0 << "RM=" << RM);
00366 Array<Array<Array<int> > > comp = compositions(order);
00367 for (int i=1; i<=order; i++)
00368 {
00369 Tabs tab1;
00370 SUNDANCE_MSG3(verb, tab1 << "i=" << i << " out of order=" << order);
00371 const Set<MultiSet<int> >& QC = findQ_C(i, context);
00372 SUNDANCE_MSG3(verb, tab1 << "QC[" << i << "] = " << QC);
00373 for (Set<MultiSet<int> >::const_iterator j=QC.begin(); j!=QC.end(); j++)
00374 {
00375 Tabs tab2;
00376 Array<int> J = j->elements();
00377 const Array<Array<int> >& K = comp[J.size()-1];
00378 SUNDANCE_MSG3(verb, tab2 << "J=" << J);
00379 SUNDANCE_MSG3(verb, tab2 << "K=" << K);
00380
00381 for (int k=0; k<K.size(); k++)
00382 {
00383 Tabs tab3;
00384 Set<MultipleDeriv> RProd = product(J, K[k],
00385 RequiredNonzeros, context);
00386 Set<MultipleDeriv> CProd = product(J, K[k],
00387 ConstantNonzeros, context);
00388 SUNDANCE_MSG3(verb, tab3 << "CProd = " << CProd);
00389 SUNDANCE_MSG3(verb, tab3 << "RProd = " << RProd);
00390 rtn.merge(RProd.setDifference(CProd));
00391 }
00392 }
00393
00394 const Set<MultiSet<int> >& QV = findQ_V(i, context);
00395 SUNDANCE_MSG3(verb, tab1 << "QV[" << i << "] = " << QV);
00396 for (Set<MultiSet<int> >::const_iterator j=QV.begin(); j!=QV.end(); j++)
00397 {
00398 Tabs tab2;
00399 Array<int> J = j->elements();
00400 const Array<Array<int> >& K = comp[J.size()-1];
00401 SUNDANCE_MSG3(verb, tab2 << "J=" << J);
00402 SUNDANCE_MSG3(verb, tab2 << "K=" << K);
00403
00404 for (int k=0; k<K.size(); k++)
00405 {
00406 Tabs tab3;
00407 Set<MultipleDeriv> RProd = product(J, K[k],
00408 RequiredNonzeros, context);
00409 SUNDANCE_MSG3(verb, tab3 << "RProd = " << RProd);
00410 rtn.merge(RProd);
00411 }
00412 }
00413 }
00414
00415 SUNDANCE_MSG3(verb, tab0 << "rtn before intersection = " << rtn);
00416 rtn = rtn.intersection(RM);
00417 SUNDANCE_MSG2(verb, tab0 << "V[" << order << "]=" << rtn);
00418 SUNDANCE_MSG2(verb, tab0 << "done with EWC::internalFindV(" << order
00419 << ") for " << toString());
00420
00421 return rtn;
00422 }
00423
00424
00425 Set<MultipleDeriv>
00426 ExprWithChildren::internalFindC(int order, const EvalContext& context) const
00427 {
00428 Tabs tab0;
00429 Set<MultipleDeriv> rtn;
00430 bool started = false;
00431 int verb = context.setupVerbosity();
00432
00433 SUNDANCE_MSG2(verb, tab0 << "EWC::internalFindC(" << order
00434 << ") for " << toString());
00435
00436
00437 if (order==0)
00438 {
00439 for (int i=0; i<numChildren(); i++)
00440 {
00441 if (!childIsRequired(i, 0, context)) continue;
00442 const Set<MultipleDeriv>& childV0
00443 = evaluatableChild(i)->findV(0, context);
00444 rtn.merge(childV0);
00445 }
00446 const Set<MultipleDeriv>& R0 = findR(0, context);
00447 rtn = R0.setDifference(rtn);
00448 SUNDANCE_MSG3(verb, tab0 << "C[" << order << "]=" << rtn);
00449 SUNDANCE_MSG3(verb, tab0 << "done with EWC::internalFindC(" << order
00450 << ") for " << toString());
00451 return rtn;
00452 }
00453
00454
00455
00456 Array<Array<Array<int> > > comp = compositions(order);
00457 const Set<MultipleDeriv>& RM = findR(order, context);
00458 for (int i=1; i<=order; i++)
00459 {
00460 Tabs tab1;
00461 SUNDANCE_MSG3(verb, tab1 << "i=" << i << " up to order=" << order);
00462
00463
00464 const Set<MultiSet<int> >& QC = findQ_C(i, context);
00465 SUNDANCE_MSG5(verb, tab1 << "finding CProd union (R\\RProd) over QC");
00466 SUNDANCE_MSG5(verb, tab1 << "QC = " << QC);
00467
00468
00469 for (Set<MultiSet<int> >::const_iterator j=QC.begin(); j!=QC.end(); j++)
00470 {
00471 Tabs tab2;
00472 Array<int> J = j->elements();
00473 const Array<Array<int> >& K = comp[J.size()-1];
00474 SUNDANCE_MSG5(verb, tab2 << "J=" << J);
00475 SUNDANCE_MSG5(verb, tab2 << "K=" << K);
00476
00477 for (int k=0; k<K.size(); k++)
00478 {
00479 Tabs tab3;
00480 Set<MultipleDeriv> RProd = product(J, K[k],
00481 RequiredNonzeros, context);
00482 Set<MultipleDeriv> CProd = product(J, K[k],
00483 ConstantNonzeros, context);
00484 SUNDANCE_MSG5(verb, tab3 << "CProd = " << CProd);
00485 SUNDANCE_MSG5(verb, tab3 << "RProd = " << RProd);
00486 Set<MultipleDeriv> X = CProd.setUnion(RM.setDifference(RProd));
00487 if (!started)
00488 {
00489 rtn = X;
00490 started = true;
00491 }
00492 else
00493 {
00494 rtn = rtn.intersection(X);
00495 }
00496 }
00497 }
00498
00499 const Set<MultiSet<int> >& QV = findQ_V(i, context);
00500 SUNDANCE_MSG5(verb, tab1 << "finding R\\RProd over QV");
00501 SUNDANCE_MSG5(verb, tab1 << "QV = " << QV);
00502
00503 for (Set<MultiSet<int> >::const_iterator j=QV.begin(); j!=QV.end(); j++)
00504 {
00505 Tabs tab2;
00506 Array<int> J = j->elements();
00507 const Array<Array<int> >& K = comp[J.size()-1];
00508 SUNDANCE_MSG5(verb, tab2 << "J=" << J);
00509 SUNDANCE_MSG5(verb, tab2 << "K=" << K);
00510
00511 for (int k=0; k<K.size(); k++)
00512 {
00513 Tabs tab3;
00514 Set<MultipleDeriv> RProd = product(J, K[k],
00515 RequiredNonzeros, context);
00516 Set<MultipleDeriv> X = RM.setDifference(RProd);
00517 if (!started)
00518 {
00519 rtn = X;
00520 started = true;
00521 }
00522 else
00523 {
00524 rtn = rtn.intersection(X);
00525 }
00526 SUNDANCE_MSG5(verb, tab3 << "RProd = " << RProd);
00527 }
00528 }
00529 }
00530
00531 rtn = rtn.intersection(RM);
00532 SUNDANCE_MSG2(verb, tab0 << "C[" << order << "]=" << rtn);
00533 SUNDANCE_MSG2(verb, tab0 << "done with EWC::internalFindC(" << order
00534 << ") for " << toString());
00535 return rtn;
00536 }
00537
00538
00539
00540 Set<MultipleDeriv>
00541 ExprWithChildren::internalFindW(int order, const EvalContext& context) const
00542 {
00543 Tabs tab0;
00544 int verb = context.setupVerbosity();
00545 Set<MultipleDeriv> rtn;
00546 SUNDANCE_MSG2(verb, tab0 << "EWC::internalFindW(order="
00547 << order << ") for " << toString());
00548 Tabs tabz;
00549
00550 if (order==0)
00551 {
00552 Tabs tab1;
00553 SUNDANCE_MSG2(verb, tab1 << "** CASE: order=0");
00554
00555
00556
00557 if (!(isLinear() || isProduct()))
00558 {
00559 Tabs tab2;
00560 SUNDANCE_MSG3(verb, tab2 << "I am neither product nor linear");
00561 rtn.put(MultipleDeriv());
00562 SUNDANCE_MSG3(verb, tab2 << "W[" << order << "]=" << rtn);
00563 SUNDANCE_MSG3(verb, tab2 << "done with EWC::internalFindW for "
00564 << toString());
00565 return rtn;
00566 }
00567
00568
00569
00570
00571 SUNDANCE_MSG3(verb, tab1 << "getting Q");
00572 const Set<MultiSet<int> >& Q = findQ_W(0, context);
00573
00574
00575
00576 if (Q.size()==0)
00577 {
00578 Tabs tab2;
00579 SUNDANCE_MSG3(verb, tab2 << "W[" << order << "]=" << rtn);
00580 SUNDANCE_MSG3(verb, tab2 << "done with EWC::internalFindW for "
00581 << toString());
00582 return rtn;
00583 }
00584
00585
00586 if (isLinear())
00587 {
00588 Tabs tab2;
00589 SUNDANCE_MSG3(verb, tab2 << "I am a linear combination");
00590 rtn.put(MultipleDeriv());
00591 SUNDANCE_MSG3(verb, tab2 << "W[" << order << "]=" << rtn);
00592 SUNDANCE_MSG3(verb, tab2 << "done with EWC::internalFindW for "
00593 << toString());
00594 return rtn;
00595 }
00596
00597
00598
00599
00600
00601
00602
00603 if (((int) Q.size()) == numChildren())
00604 {
00605 rtn.put(MultipleDeriv());
00606 }
00607 SUNDANCE_MSG2(verb, tab1 << "I am a product");
00608 SUNDANCE_MSG2(verb, tab1 << "W[" << order << "]=" << rtn);
00609 SUNDANCE_MSG2(verb, tab0 << "done with EWC::internalFindW for "
00610 << toString());
00611 return rtn;
00612 }
00613
00614
00615
00616 Array<Array<Array<int> > > comp = compositions(order);
00617 for (int i=1; i<=order; i++)
00618 {
00619 Tabs tab1;
00620 SUNDANCE_MSG2(verb, tab1 << "** CASE order=" << i);
00621
00622 const Set<MultiSet<int> >& QW = findQ_W(i, context);
00623
00624 SUNDANCE_MSG3(verb, tab1 << "QW=" << QW);
00625
00626 for (Set<MultiSet<int> >::const_iterator j=QW.begin(); j!=QW.end(); j++)
00627 {
00628 Array<int> J = j->elements();
00629 const Array<Array<int> >& K = comp[J.size()-1];
00630
00631 for (int k=0; k<K.size(); k++)
00632 {
00633 Set<MultipleDeriv> WProd = product(J, K[k], AllNonzeros, context);
00634 rtn.merge(WProd);
00635 }
00636 }
00637 }
00638
00639 SUNDANCE_MSG2(verb, tabz << "W[" << order << "]=" << rtn);
00640 SUNDANCE_MSG2(verb, tab0 << "done with EWC::internalFindW for "
00641 << toString());
00642 return rtn;
00643 }
00644
00645 const Set<MultiSet<int> >&
00646 ExprWithChildren::findQ_W(int order,
00647 const EvalContext& context) const
00648 {
00649 Tabs tab1;
00650 int verb = context.setupVerbosity();
00651 SUNDANCE_MSG2(verb, tab1 << "finding Q_W");
00652 if (!contextToQWMap_[order].containsKey(context))
00653 {
00654 contextToQWMap_[order].put(context, internalFindQ_W(order, context));
00655 }
00656 else
00657 {
00658 SUNDANCE_MSG5(verb, tab1 << "using previously computed Q_W");
00659 }
00660 return contextToQWMap_[order].get(context);
00661
00662 }
00663
00664 const Set<MultiSet<int> >&
00665 ExprWithChildren::findQ_V(int order,
00666 const EvalContext& context) const
00667 {
00668 if (!contextToQVMap_[order].containsKey(context))
00669 {
00670 contextToQVMap_[order].put(context, internalFindQ_V(order, context));
00671 }
00672 return contextToQVMap_[order].get(context);
00673 }
00674
00675 const Set<MultiSet<int> >&
00676 ExprWithChildren::findQ_C(int order,
00677 const EvalContext& context) const
00678 {
00679 if (!contextToQCMap_[order].containsKey(context))
00680 {
00681 contextToQCMap_[order].put(context, internalFindQ_C(order, context));
00682 }
00683 return contextToQCMap_[order].get(context);
00684 }
00685
00686
00687 const Set<MultiSet<int> >& ExprWithChildren::getI_N() const
00688 {
00689 if (!cachedI_N().containsKey(numChildren()))
00690 {
00691 Set<MultiSet<int> > x;
00692 for (int i=0; i<numChildren(); i++)
00693 {
00694 x.put(makeMultiSet<int>(i));
00695 }
00696 cachedI_N().put(numChildren(), x);
00697 }
00698 return cachedI_N().get(numChildren());
00699 }
00700
00701 Set<MultiSet<int> > ExprWithChildren::indexSetProduct(const Set<MultiSet<int> >& a,
00702 const Set<MultiSet<int> >& b) const
00703 {
00704 Set<MultiSet<int> > rtn;
00705 for (Set<MultiSet<int> >::const_iterator i=a.begin(); i!=a.end(); i++)
00706 {
00707 for (Set<MultiSet<int> >::const_iterator j=b.begin(); j!=b.end(); j++)
00708 {
00709 MultiSet<int> ab = (*i).merge(*j);
00710 rtn.put(ab);
00711 }
00712 }
00713 return rtn;
00714 }
00715
00716
00717 Set<MultiSet<int> > ExprWithChildren::internalFindQ_V(int order,
00718 const EvalContext& context) const
00719 {
00720 Tabs tab0;
00721 int verb = context.setupVerbosity();
00722 SUNDANCE_MSG2(verb, tab0 << "EWC::internalFindQ_V(order=" << order <<")");
00723 Set<MultiSet<int> > rtn;
00724
00725 if (!isLinear())
00726 {
00727 bool isVar = false;
00728 for (int i=0; i<numChildren(); i++)
00729 {
00730 if (childIsRequired(i,order,context) && evaluatableChild(i)->findV(0, context).size() != 0)
00731 {
00732 isVar=true;
00733 break;
00734 }
00735 }
00736 if (isVar) rtn = findQ_W(order, context);
00737 }
00738
00739 SUNDANCE_MSG2(verb, tab0 << "Q_V = " << rtn);
00740 return rtn;
00741 }
00742
00743 Set<MultiSet<int> > ExprWithChildren::internalFindQ_C(int order,
00744 const EvalContext& context) const
00745 {
00746 if (isLinear()) return findQ_W(order,context);
00747
00748 return findQ_W(order,context).setDifference(findQ_V(order, context));
00749 }
00750
00751
00752 Set<MultiSet<int> > ExprWithChildren
00753 ::internalFindQ_W(int order,
00754 const EvalContext& context) const
00755 {
00756 Tabs tab0;
00757 int verb = context.setupVerbosity();
00758 SUNDANCE_MSG2(verb, tab0 << "in internalFindQ_W");
00759 Set<MultiSet<int> > rtn;
00760 const Set<MultiSet<int> >& I_N = getI_N();
00761 SUNDANCE_MSG2(verb, tab0 << "I_N=" << I_N);
00762
00763 if (isLinear())
00764 {
00765 Tabs tab1;
00766
00767
00768 if (order==1)
00769 {
00770 SUNDANCE_MSG3(verb, tab1 << "is linear, order=1, so Q_W = I_N");
00771 return I_N;
00772 }
00773
00774 if (order==0)
00775 {
00776 SUNDANCE_MSG3(verb, tab1 << "is linear, order=0");
00777 for (int i=0; i<numChildren(); i++)
00778 {
00779 const Set<MultipleDeriv>& W_i = evaluatableChild(i)->findW(0,context);
00780 SUNDANCE_MSG3(verb, tab1 << "is linear, order=0");
00781 if (W_i.size() > 0)
00782 {
00783 Tabs tab2;
00784 SUNDANCE_MSG3(verb, tab2 << "child " << i << " is nonzero");
00785 rtn.put(makeMultiSet(i));
00786 }
00787 else
00788 {
00789 Tabs tab2;
00790 SUNDANCE_MSG3(verb, tab2 << "child " << i << " is zero");
00791 }
00792 }
00793 }
00794 }
00795 else
00796 {
00797 Tabs tab1;
00798 SUNDANCE_MSG3(verb, tab1 << "is nonlinear, using all index set products");
00799 rtn = I_N;
00800
00801 for (int i=1; i<order; i++)
00802 {
00803 rtn = indexSetProduct(rtn, I_N);
00804 }
00805 }
00806 SUNDANCE_MSG2(verb, tab0 << "rtn=" << rtn);
00807 return rtn;
00808 }
00809
00810 bool ExprWithChildren::childIsRequired(int index, int diffOrder,
00811 const EvalContext& context) const
00812 {
00813 const Set<MultiSet<int> >& Q = findQ_W(diffOrder, context);
00814 for (Set<MultiSet<int> >::const_iterator it=Q.begin(); it != Q.end(); it++)
00815 {
00816 if (it->contains(index)) return true;
00817 }
00818 return true;
00819 }
00820
00821 RCP<Array<Set<MultipleDeriv> > > ExprWithChildren
00822 ::internalDetermineR(const EvalContext& context,
00823 const Array<Set<MultipleDeriv> >& RInput) const
00824 {
00825 Tabs tab0;
00826 int verb = context.setupVerbosity();
00827 RCP<Array<Set<MultipleDeriv> > > rtn
00828 = rcp(new Array<Set<MultipleDeriv> >(RInput.size()));
00829
00830 SUNDANCE_MSG2(verb, tab0 << "in internalDetermineR() for " << toString());
00831 SUNDANCE_MSG2(verb, tab0 << "RInput = " << RInput);
00832
00833
00834
00835 for (int i=0; i<RInput.size(); i++)
00836 {
00837 if (RInput[i].size()==0) continue;
00838 const Set<MultipleDeriv>& Wi = findW(i, context);
00839 (*rtn)[i] = RInput[i].intersection(Wi);
00840 }
00841
00842 int maxOrder = rtn->size()-1;
00843
00844 const Set<MultiSet<int> >& Q1 = findQ_W(1, context);
00845 const Set<MultiSet<int> >& Q2 = findQ_W(2, context);
00846 const Set<MultiSet<int> >& Q3 = findQ_W(3, context);
00847
00848 SUNDANCE_MSG3(verb, tab0 << "Q1 = " << Q1);
00849 SUNDANCE_MSG3(verb, tab0 << "Q2 = " << Q2);
00850 SUNDANCE_MSG3(verb, tab0 << "Q3 = " << Q3);
00851
00852 for (int i=0; i<numChildren(); i++)
00853 {
00854 Tabs tab1;
00855 MultiSet<int> mi = makeMultiSet(i);
00856 SUNDANCE_MSG4(verb, tab1 << "i=" << i << ", Q1_i = " << mi );
00857 TEUCHOS_TEST_FOR_EXCEPTION(mi.size() != 1, std::logic_error, "unexpected multiset size");
00858
00859 Set<MultipleDeriv> R11;
00860 Set<MultipleDeriv> R12;
00861 Set<MultipleDeriv> R13;
00862 Set<MultipleDeriv> R22;
00863 Set<MultipleDeriv> R23;
00864 Set<MultipleDeriv> R33;
00865 if (maxOrder >= 1)
00866 {
00867 R11 = (*rtn)[1];
00868 if (maxOrder >=2) R22 = (*rtn)[2];
00869 if (maxOrder >=3) R33 = (*rtn)[3];
00870 }
00871 if (maxOrder >= 2)
00872 {
00873 Tabs tab2;
00874 Set<MultiSet<int> > jSet = setDivision(Q2, i);
00875 SUNDANCE_MSG5(verb, tab2 << "Q2/i = " << jSet);
00876 for (Set<MultiSet<int> >::const_iterator
00877 j=jSet.begin(); j!=jSet.end(); j++)
00878 {
00879 Tabs tab3;
00880 TEUCHOS_TEST_FOR_EXCEPTION(j->size()!=1, std::logic_error,
00881 "unexpected set size");
00882 int jIndex = *(j->begin());
00883 SUNDANCE_MSG5(verb, tab3 << "j=" << jIndex );
00884 const Set<MultipleDeriv>& W1j = evaluatableChild(jIndex)->findW(1, context);
00885 Set<MultipleDeriv> ROverW = setDivision((*rtn)[2], W1j);
00886 R12.merge(ROverW);
00887 SUNDANCE_MSG5(verb, tab3 << "R2=" << (*rtn)[2] );
00888 SUNDANCE_MSG5(verb, tab3 << "W1(j)=" << W1j );
00889 SUNDANCE_MSG5(verb, tab3 << "R2/W1(j)=" << ROverW );
00890
00891 if (maxOrder>=3)
00892 {
00893 const Set<MultipleDeriv>& W2j
00894 = evaluatableChild(jIndex)->findW(2, context);
00895 R13.merge(setDivision((*rtn)[3], W2j));
00896 R23.merge(setDivision((*rtn)[3], W1j));
00897 }
00898 }
00899 }
00900 if (maxOrder >= 3)
00901 {
00902 Set<MultiSet<int> > jkSet = setDivision(Q3, i);
00903 for (Set<MultiSet<int> >::const_iterator
00904 jk=jkSet.begin(); jk!=jkSet.end(); jk++)
00905 {
00906 TEUCHOS_TEST_FOR_EXCEPTION(jk->size()!=2, std::logic_error,
00907 "unexpected set size");
00908 Array<int> jka = jk->elements();
00909 int j = jka[0];
00910 int k = jka[1];
00911 const Set<MultipleDeriv>& W1j = evaluatableChild(j)->findW(1, context);
00912 const Set<MultipleDeriv>& W1k = evaluatableChild(k)->findW(1, context);
00913 R13.merge(setDivision((*rtn)[3], setProduct(W1j, W1k)));
00914 }
00915 }
00916 SUNDANCE_MSG5(verb, tab1 << "R11 = " << R11 );
00917 SUNDANCE_MSG5(verb, tab1 << "R12 = " << R12 );
00918 SUNDANCE_MSG5(verb, tab1 << "R13 = " << R13 );
00919 SUNDANCE_MSG5(verb, tab1 << "R22 = " << R22 );
00920 SUNDANCE_MSG5(verb, tab1 << "R23 = " << R23 );
00921 SUNDANCE_MSG5(verb, tab1 << "R33 = " << R33 );
00922 Set<MultipleDeriv> R1 = R11;
00923 R1.merge(R12);
00924 R1.merge(R13);
00925 Set<MultipleDeriv> R2 = R22;
00926 R2.merge(R23);
00927 Set<MultipleDeriv> R3 = R33;
00928
00929 Set<MultipleDeriv> R0;
00930 bool childIsNeeded = (*rtn)[0].size() > 0;
00931 if (!childIsNeeded && R1.size() > 0) childIsNeeded = childIsRequired(i, 2, context);
00932 if (!childIsNeeded && R2.size() > 0) childIsNeeded = childIsRequired(i, 3, context);
00933 if (!childIsNeeded && R3.size() > 0) childIsNeeded = childIsRequired(i, 4, context);
00934 if (childIsNeeded) R0.put(MultipleDeriv());
00935
00936 Array<Set<MultipleDeriv> > RChild;
00937
00938 RChild.append(R0);
00939 if (maxOrder >= 1) RChild.append(R1);
00940 if (maxOrder >= 2) RChild.append(R2);
00941 if (maxOrder >= 3) RChild.append(R3);
00942 SUNDANCE_MSG4(verb, tab1 << "RChild = " << RChild );
00943 evaluatableChild(i)->determineR(context, RChild);
00944 }
00945
00946 printR(verb, rtn);
00947 SUNDANCE_MSG2(verb, tab0 << "done with EWC::internalDetermineR for "
00948 << toString());
00949
00950 return rtn;
00951 }
00952
00953
00954
00955
00956 void ExprWithChildren::displayNonzeros(std::ostream& os, const EvalContext& context) const
00957 {
00958 Tabs tabs0(0);
00959 os << tabs0 << "Nonzeros of " << toString() << std::endl;
00960 if (context.setupVerbosity() > 4)
00961 {
00962 os << tabs0 << "Diving into children " << std::endl;
00963 }
00964
00965 for (int i=0; i<numChildren(); i++)
00966 {
00967 Tabs tab1;
00968 if (context.setupVerbosity() > 4) os << tab1 << "Child " << i << std::endl;
00969 evaluatableChild(i)->displayNonzeros(os, context);
00970 }
00971
00972 if (context.setupVerbosity() > 4)
00973 os << tabs0 << "Printing nonzeros for parent " << toString() << std::endl;
00974 const Set<MultipleDeriv>& W = findW(context);
00975 const Set<MultipleDeriv>& R = findR(context);
00976 const Set<MultipleDeriv>& C = findC(context);
00977
00978
00979 for (Set<MultipleDeriv>::const_iterator i=W.begin(); i != W.end(); i++)
00980 {
00981 Tabs tab1;
00982 std::string state = "Variable";
00983 if (C.contains(*i)) state = "Constant";
00984 if (!R.contains(*i)) state = "Not Required";
00985 os << tab1 << std::setw(25) << std::left << i->toString() << ": " << state << std::endl;
00986 }
00987 }
00988
00989
00990 namespace Sundance {
00991
00992 Array<Array<std::pair<int, Array<MultipleDeriv> > > >
00993 chainRuleDerivsOfArgs(int nArgs,
00994 const MultiSet<int>& bSet,
00995 const MultipleDeriv& c)
00996 {
00997 Array<Array<std::pair<int, Array<MultipleDeriv> > > > rtn;
00998
00999
01000 Array<int> b(nArgs, 0);
01001 int J = 0;
01002 for (MultiSet<int>::const_iterator i=bSet.begin(); i!=bSet.end(); i++)
01003 {
01004 b[*i]++;
01005 J++;
01006 }
01007
01008
01009 Sundance::Map<Deriv, int> counts;
01010 Array<Deriv> d;
01011 typedef Sundance::Map<Deriv, int>::const_iterator iter;
01012 for (MultipleDeriv::const_iterator i=c.begin(); i!=c.end(); i++)
01013 {
01014 if (!counts.containsKey(*i)) counts.put(*i, 1);
01015 else counts[*i]++;
01016 d.append(*i);
01017 }
01018
01019 Array<Array<Array<Array<int> > > > a;
01020 Array<int> s;
01021 for (iter i=counts.begin(); i!=counts.end(); i++)
01022 {
01023 Array<Array<int> > tmp = nonNegCompositions(i->second, J);
01024 Array<Array<Array<int> > > ai = bStructure(b, tmp);
01025 a.append(ai);
01026 s.append(ai.size());
01027 }
01028
01029 Array<Array<int> > ic = indexCombinations(s);
01030
01031
01032
01033 Array<Array<Array<Array<int> > > > all;
01034 for (int i=0; i<ic.size(); i++)
01035 {
01036 bool good = true;
01037 int numFuncs = ic[i].size();
01038 Array<Array<Array<int> > > tmp;
01039
01040 Array<Array<int> > aTot(nArgs);
01041 for (int j=0; j<nArgs; j++)
01042 {
01043 if (b[j] > 0) aTot[j].resize(b[j]);
01044 for (int k=0; k<b[j]; k++) aTot[j][k] = 0;
01045 }
01046 for (int f=0; f<numFuncs; f++)
01047 {
01048 bool skip = false;
01049 const Array<Array<int> >& e = a[f][ic[i][f]];
01050 for (int j=0; j<nArgs; j++)
01051 {
01052 for (int k=0; k<b[j]; k++)
01053 {
01054 aTot[j][k] += e[j][k];
01055 }
01056 }
01057 if (!skip) tmp.append(e);
01058 }
01059
01060 for (int j=0; j<nArgs; j++)
01061 {
01062 for (int k=0; k<b[j]; k++)
01063 {
01064 if (aTot[j][k] == 0) good = false;
01065 }
01066 }
01067 if (good) all.append(tmp);
01068 }
01069
01070
01071 for (int p=0; p<all.size(); p++)
01072 {
01073 Array<std::pair<int, Array<MultipleDeriv> > > terms;
01074 for (int j=0; j<nArgs; j++)
01075 {
01076 std::pair<int, Array<MultipleDeriv> > factors;
01077 factors.first = j;
01078 for (int k=0; k<b[j]; k++)
01079 {
01080 MultipleDeriv md;
01081 for (int i=0; i<all[p].size(); i++)
01082 {
01083 int order = all[p][i][j][k];
01084 if (order > 0)
01085 {
01086 md.put(d[i]);
01087 }
01088 }
01089 if (md.order() > 0) factors.second.append(md);
01090 }
01091 if (factors.second.size() > 0) terms.append(factors);
01092 }
01093 rtn.append(terms);
01094 }
01095
01096 return rtn;
01097 }
01098
01099
01100 Array<Array<Array<int> > > bStructure(const Array<int>& b,
01101 const Array<Array<int> >& tmp)
01102 {
01103 Array<Array<Array<int> > > rtn(tmp.size());
01104
01105 for (int p=0; p<tmp.size(); p++)
01106 {
01107 int count=0;
01108 rtn[p].resize(b.size());
01109 for (int j=0; j<b.size(); j++)
01110 {
01111 rtn[p][j].resize(b[j]);
01112 for (int k=0; k<b[j]; k++, count++)
01113 {
01114 rtn[p][j][k] = tmp[p][count];
01115 }
01116 }
01117 }
01118 return rtn;
01119 }
01120
01121 Array<OrderedPair<Array<MultiSet<int> >, Array<MultipleDeriv> > >
01122 chainRuleTerms(int s,
01123 const MultiSet<int>& lambda,
01124 const MultipleDeriv& nu)
01125 {
01126 Array<Array<MultiSet<int> > > allK = multisetCompositions(s, lambda);
01127
01128 Array<MultipleDeriv> allL = multisetSubsets(nu).elements();
01129
01130 Array<Array<int> > indexTuples = distinctIndexTuples(s, allL.size());
01131
01132 Array<Array<MultipleDeriv> > allLTuples;
01133 for (int i=0; i<indexTuples.size(); i++)
01134 {
01135 Array<MultipleDeriv> t(s);
01136 for (int p=0; p<s; p++) t[p] = allL[indexTuples[i][p]];
01137 allLTuples.append(t);
01138 }
01139
01140 Array<OrderedPair<Array<MultiSet<int> >, Array<MultipleDeriv> > > rtn;
01141 for (int i=0; i<allLTuples.size(); i++)
01142 {
01143 for (int j=0; j<allK.size(); j++)
01144 {
01145 MultipleDeriv result;
01146 for (int p=0; p<s; p++)
01147 {
01148 for (unsigned int q=0; q<allK[j][p].size(); q++)
01149 {
01150 result = result.product(allLTuples[i][p]);
01151 }
01152 }
01153 if (result==nu)
01154 {
01155 OrderedPair<Array<MultiSet<int> >, Array<MultipleDeriv> >
01156 kl(allK[j], allLTuples[i]);
01157 rtn.append(kl);
01158 }
01159 }
01160 }
01161 return rtn;
01162 }
01163
01164 Set<MultipleDeriv> multisetSubsets(const MultipleDeriv& nu)
01165 {
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177 Array<Deriv> elements = nu.elements();
01178
01179
01180
01181 int n = elements.size();
01182 int maxNumSubsets = pow2(n);
01183
01184 Set<MultipleDeriv> rtn;
01185
01186
01187
01188
01189 for (int i=1; i<maxNumSubsets; i++)
01190 {
01191 Array<int> bits = bitsOfAnInteger(i, n);
01192 MultipleDeriv md;
01193 for (int j=0; j<n; j++)
01194 {
01195 if (bits[j] == 1) md.put(elements[j]);
01196 }
01197 rtn.put(md);
01198 }
01199 return rtn;
01200 }
01201
01202
01203
01204 int chainRuleMultiplicity(const MultipleDeriv& nu,
01205 const Array<MultiSet<int> >& K,
01206 const Array<MultipleDeriv>& L)
01207 {
01208 int rtn = factorial(nu);
01209 for (int i=0; i<K.size(); i++)
01210 {
01211 rtn = rtn/factorial(K[i]);
01212 int lFact = factorial(L[i]);
01213 int lPow = 1;
01214 int kNorm = K[i].size();
01215 for (int j=0; j<kNorm; j++) lPow *= lFact;
01216 TEUCHOS_TEST_FOR_EXCEPT(rtn % lPow != 0);
01217 rtn = rtn/lPow;
01218 }
01219 return rtn;
01220 }
01221
01222
01223 }
01224
01225
01226