Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
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   /* Create the evaluators for the children first so that I can refer to
00163    * them when I create my own evaluator */
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   /* count self */
00293   int count = EvaluatableExpr::countNodes();
00294 
00295   /* count children */
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   /* we'll dealt with zero order derivatives specially */
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   /* now do arbitrary order derivatives with the multivariable chain rule*/
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   /* we'll dealt with zero order derivatives specially */
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   /* now do arbitrary order derivatives with the multivariable chain rule*/
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   /* we'll deal with zero order derivatives specially */
00550   if (order==0) 
00551   {
00552     Tabs tab1;
00553     SUNDANCE_MSG2(verb, tab1 << "** CASE: order=0");  
00554     /* If I am an arbitrary nonlinear expression, I cannot be known to
00555      * be zero regardless of the state of my arguments. Return the
00556      * zeroth-order derivative. */
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     /* At this point, I've dealt with arbitrary nonlinear exprs so 
00569      * I know I'm either a product or a linear expr */
00570 
00571     SUNDANCE_MSG3(verb, tab1 << "getting Q");  
00572     const Set<MultiSet<int> >& Q = findQ_W(0, context);
00573 
00574     /* If there are no nonzero terms, a linear combination or a product
00575      * will be zero. Return the empty set. */
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     /* if I'm a linear combination and any term is nonzero, I am nonzero */
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     /* The only possibility remaining is that I'm a product.
00598      * If any term is zero, I am zero. If a term is 
00599      * known to be zero, it's index will not appear in Q, so that 
00600      * comparing the size of Q and the number of children tells me
00601      * if I'm zero.
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   /* now do arbitrary order derivatives with the multivariable chain rule*/
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     /* first derivatives of the sum wrt the arguments are 
00767      * always nonzero */
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     /* zeroth derivatives are nonzero if terms are nonzero */
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     //      int i = *(mi.begin());
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   /* convert to tuple representation of b */
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   /* count orders of each functional deriv in c */
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   /* We'll generate the subsets by traversing them in bitwise order.
01167    * For a multiset having N elements, there are up to 2^N subsets each
01168    * of which can be described by a N-bit number with the i-th
01169    * bit indicating whether the i-th element is in the subset. Note that
01170    * with a multiset, repetitions can occur so we need to record the
01171    * results in a Set object to eliminate duplicates. 
01172    */
01173 
01174   /* Make an indexable array of the elements. This will be convenient
01175    * because we'll need to access the i-th element after reading
01176    * the i-th bit. */
01177   Array<Deriv> elements = nu.elements();
01178 
01179   /* Compute the maximum number of subsets. This number will be reached
01180    * only in the case of no repetitions. */
01181   int n = elements.size();
01182   int maxNumSubsets = pow2(n);
01183     
01184   Set<MultipleDeriv> rtn;
01185     
01186     
01187   /* Loop over subsets in bitwise order. We start the count at 1 
01188      to avoid including the empty subset */
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
SundanceExprWithChildren.cpp
