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 "SundanceNonlinearUnaryOpEvaluator.hpp"
00032 #include "SundanceEvalManager.hpp"
00033
00034 #include "PlayaTabs.hpp"
00035 #include "SundanceOut.hpp"
00036 #include "SundanceNonlinearUnaryOp.hpp"
00037
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042
00043
00044
00045
00046
00047 NonlinearUnaryOpEvaluator
00048 ::NonlinearUnaryOpEvaluator(const NonlinearUnaryOp* expr,
00049 const EvalContext& context)
00050 : ChainRuleEvaluator(expr, context),
00051 op_(expr->op()),
00052 maxOrder_(-1),
00053 argIsConstant_(true),
00054 argValueIndex_(-1)
00055 {
00056 Tabs tabs;
00057 SUNDANCE_VERB_LOW(tabs << "NonlinearUnaryOp evaluator ctor for "
00058 << expr->toString());
00059
00060
00061 const SparsitySuperset* sArg = childSparsity(0);
00062 SUNDANCE_VERB_HIGH(tabs << "arg sparsity " << *sArg);
00063 maxOrder_ = expr->sparsitySuperset(context)->maxOrder();
00064
00065
00066 for (int i=0; i<sArg->numDerivs(); i++)
00067 {
00068 if (sArg->state(i) == VectorDeriv)
00069 {
00070 argIsConstant_ = false;
00071 break;
00072 }
00073 }
00074
00075
00076
00077 for (int order=0; order<=maxOrder_; order++)
00078 {
00079 MultiSet<int> df;
00080
00081
00082 for (int i=0; i<order; i++) df.put(0);
00083
00084 if (argIsConstant_) addConstArgDeriv(df, order);
00085 else addVarArgDeriv(df, order);
00086 }
00087
00088
00089
00090
00091 MultipleDeriv d0;
00092 TEUCHOS_TEST_FOR_EXCEPTION(!sArg->containsDeriv(d0), std::logic_error,
00093 "NonlinearUnaryOpEvaluator::ctor did not find zeroth-order "
00094 "derivative of argument");
00095 int d0Index = sArg->getIndex(d0);
00096 const Evaluator* argEv = childEvaluator(0);
00097 if (argIsConstant_) argValueIndex_ = argEv->constantIndexMap().get(d0Index);
00098 else argValueIndex_ = argEv->vectorIndexMap().get(d0Index);
00099
00100
00101 SUNDANCE_VERB_HIGH(tabs << "arg is constant: " << argIsConstant_);
00102 SUNDANCE_VERB_HIGH(tabs << "arg value index: " << argValueIndex_);
00103
00104 init(expr, context);
00105 }
00106
00107
00108 void NonlinearUnaryOpEvaluator
00109 ::evalArgDerivs(const EvalManager& mgr,
00110 const Array<RCP<Array<double> > >& constArgRes,
00111 const Array<RCP<Array<RCP<EvalVector> > > >& vArgResults,
00112 Array<double>& constArgDerivs,
00113 Array<RCP<EvalVector> >& varArgDerivs) const
00114 {
00115 Tabs tabs;
00116 SUNDANCE_MSG1(mgr.verb(), tabs
00117 << "NonlinearUnaryOpEvaluator::evalArgDerivs() for "
00118 << expr()->toString());
00119
00120 if (argIsConstant_)
00121 {
00122 double argValue = (*(constArgRes[0]))[argValueIndex_];
00123 constArgDerivs.resize(maxOrder_+1);
00124 switch(maxOrder_)
00125 {
00126 case 0:
00127 op_->eval0(&argValue, 1, &(constArgDerivs[0]));
00128 break;
00129 case 1:
00130 op_->eval1(&argValue, 1, &(constArgDerivs[0]), &(constArgDerivs[1]));
00131 break;
00132 case 2:
00133 op_->eval2(&argValue, 1, &(constArgDerivs[0]), &(constArgDerivs[1]),
00134 &(constArgDerivs[2]));
00135 break;
00136 case 3:
00137 op_->eval3(&argValue, 1, &(constArgDerivs[0]), &(constArgDerivs[1]),
00138 &(constArgDerivs[2]), &(constArgDerivs[3]));
00139 break;
00140 default:
00141 TEUCHOS_TEST_FOR_EXCEPT(true);
00142 }
00143 }
00144 else
00145 {
00146 Tabs tab1;
00147 SUNDANCE_MSG2(mgr.verb(),
00148 tab1 << "argument is a vector: argValIndex=" << argValueIndex_);
00149 SUNDANCE_MSG2(mgr.verb(),
00150 tab1 << "num vector results =" << vArgResults.size());
00151
00152 const Array<RCP<EvalVector> >& av = *(vArgResults[0]);
00153
00154 SUNDANCE_MSG2(mgr.verb(),
00155 tab1 << "av size=" << av.size());
00156
00157
00158 const double* argValue = av[argValueIndex_]->start();
00159 varArgDerivs.resize(maxOrder_+1);
00160 varArgDerivs[0] = mgr.stack().popVector();
00161 int nx = varArgDerivs[0]->length();
00162
00163 const std::string& argStr = (*(vArgResults[0]))[argValueIndex_]->str();
00164 if (EvalVector::shadowOps())
00165 {
00166 varArgDerivs[0]->setString(op_->name() + "(" + argStr + ")");
00167 }
00168
00169 double* f = varArgDerivs[0]->start();
00170 double* df_dArg = 0 ;
00171 if (maxOrder_ >= 1)
00172 {
00173 varArgDerivs[1] = mgr.stack().popVector();
00174 if (EvalVector::shadowOps())
00175 {
00176 varArgDerivs[1]->setString(op_->name() + "'(" + argStr + ")");
00177 }
00178 df_dArg = varArgDerivs[1]->start();
00179 }
00180 double* d2f_dArg2 = 0 ;
00181 if (maxOrder_ >= 2)
00182 {
00183 varArgDerivs[2] =mgr.stack().popVector();
00184 if (EvalVector::shadowOps())
00185 {
00186 varArgDerivs[2]->setString(op_->name() + "''(" + argStr + ")");
00187 }
00188 d2f_dArg2 = varArgDerivs[2]->start();
00189 }
00190 double* d3f_dArg3 = 0 ;
00191 if (maxOrder_ >= 3)
00192 {
00193 UnaryFunctor::fdStep()=1.0e-3;
00194 varArgDerivs[3] = mgr.stack().popVector();
00195 if (EvalVector::shadowOps())
00196 {
00197 varArgDerivs[3]->setString(op_->name() + "'''(" + argStr + ")");
00198 }
00199 d3f_dArg3 = varArgDerivs[3]->start();
00200 }
00201
00202 switch(maxOrder_)
00203 {
00204 case 0:
00205 op_->eval0(argValue, nx, f);
00206 break;
00207 case 1:
00208 TEUCHOS_TEST_FOR_EXCEPT(df_dArg==0);
00209 op_->eval1(argValue, nx, f, df_dArg);
00210 break;
00211 case 2:
00212 TEUCHOS_TEST_FOR_EXCEPT(df_dArg==0);
00213 TEUCHOS_TEST_FOR_EXCEPT(d2f_dArg2==0);
00214 op_->eval2(argValue, nx, f, df_dArg, d2f_dArg2);
00215 break;
00216 case 3:
00217 TEUCHOS_TEST_FOR_EXCEPT(df_dArg==0);
00218 TEUCHOS_TEST_FOR_EXCEPT(d2f_dArg2==0);
00219 TEUCHOS_TEST_FOR_EXCEPT(d3f_dArg3==0);
00220 op_->eval3(argValue, nx, f, df_dArg, d2f_dArg2, d3f_dArg3);
00221 break;
00222 default:
00223 TEUCHOS_TEST_FOR_EXCEPT(true);
00224 }
00225 }
00226
00227 SUNDANCE_MSG1(mgr.verb(),
00228 tabs << "done NonlinearUnaryOpEvaluator::evalArgDerivs() for "
00229 << expr()->toString());
00230 }
00231
00232
00233