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 "SundanceSymbPreprocessor.hpp"
00032 #include "SundanceEvaluatorFactory.hpp"
00033 #include "SundanceEvaluator.hpp"
00034 #include "SundanceEvalManager.hpp"
00035 #include "SundanceExpr.hpp"
00036 #include "PlayaTabs.hpp"
00037 #include "SundanceZeroExpr.hpp"
00038 #include "SundanceDiscreteFuncElement.hpp"
00039 #include "SundanceDiscreteFunctionStub.hpp"
00040 #include "SundanceTestFuncElement.hpp"
00041 #include "SundanceUnknownFuncElement.hpp"
00042 #include "SundanceUnknownParameterElement.hpp"
00043 #include "SundanceUnknownFunctionStub.hpp"
00044 #include "SundanceTestFunctionStub.hpp"
00045
00046 #include "SundanceOut.hpp"
00047 #include "Teuchos_Utils.hpp"
00048
00049 using namespace Sundance;
00050 using namespace Sundance;
00051 using namespace Teuchos;
00052
00053
00054
00055
00056
00057 DerivSet SymbPreprocessor::setupFwdProblem(const Expr& expr,
00058 const Expr& tests,
00059 const Expr& unks,
00060 const Expr& unkEvalPts,
00061 const Expr& unkParams,
00062 const Expr& unkParamEvalPts,
00063 const Expr& fixedParams,
00064 const Expr& fixedParamEvalPts,
00065 const Expr& fixedFields,
00066 const Expr& fixedFieldEvalPts,
00067 const EvalContext& context,
00068 const ComputationType& compType)
00069 {
00070 Expr zero;
00071 Expr v = tests.flatten();
00072 Array<Expr> z(v.size());
00073 for (int i=0; i<v.size(); i++) z[i] = new ZeroExpr();
00074 zero = new ListExpr(z);
00075
00076 return setupVariations(expr,
00077 tests, zero,
00078 unks, unkEvalPts,
00079 unkParams, unkParamEvalPts,
00080 fixedFields, fixedFieldEvalPts,
00081 fixedParams, fixedParamEvalPts,
00082 context, compType);
00083 }
00084
00085
00086
00087
00088 DerivSet SymbPreprocessor::setupSensitivities(const Expr& expr,
00089 const Expr& tests,
00090 const Expr& unks,
00091 const Expr& unkEvalPts,
00092 const Expr& unkParams,
00093 const Expr& unkParamEvalPts,
00094 const Expr& fixedParams,
00095 const Expr& fixedParamEvalPts,
00096 const Expr& fixedFields,
00097 const Expr& fixedFieldEvalPts,
00098 const EvalContext& context,
00099 const ComputationType& compType)
00100 {
00101 Expr zero;
00102 Expr v = tests.flatten();
00103 Array<Expr> z(v.size());
00104 for (int i=0; i<v.size(); i++) z[i] = new ZeroExpr();
00105 zero = new ListExpr(z);
00106
00107 return setupVariations(expr,
00108 tests, zero,
00109 unks, unkEvalPts,
00110 unkParams, unkParamEvalPts,
00111 fixedFields, fixedFieldEvalPts,
00112 fixedParams, fixedParamEvalPts,
00113 context, compType);
00114 }
00115
00116
00117 DerivSet SymbPreprocessor::setupFunctional(const Expr& expr,
00118 const Expr& fixedParams,
00119 const Expr& fixedParamEvalPts,
00120 const Expr& fixedFields,
00121 const Expr& fixedFieldEvalPts,
00122 const EvalContext& context,
00123 const ComputationType& compType)
00124 {
00125 Expr vars;
00126 Expr varEvalPts;
00127 Expr unks;
00128 Expr unkEvalPts;
00129 Expr unkParams;
00130 Expr unkParamEvalPts;
00131
00132 return setupVariations(expr,
00133 vars, varEvalPts,
00134 unks, unkEvalPts,
00135 unkParams, unkParamEvalPts,
00136 fixedFields, fixedFieldEvalPts,
00137 fixedParams, fixedParamEvalPts,
00138 context,compType);
00139 }
00140
00141
00142
00143
00144
00145 DerivSet SymbPreprocessor::setupGradient(const Expr& expr,
00146 const Expr& vars,
00147 const Expr& varEvalPts,
00148 const Expr& fixedParams,
00149 const Expr& fixedParamEvalPts,
00150 const Expr& fixedFields,
00151 const Expr& fixedFieldEvalPts,
00152 const EvalContext& context,
00153 const ComputationType& compType)
00154 {
00155 Expr unks;
00156 Expr unkEvalPts;
00157 Expr unkParams;
00158 Expr unkParamEvalPts;
00159
00160 return setupVariations(expr,
00161 vars, varEvalPts,
00162 unks, unkEvalPts,
00163 unkParams, unkParamEvalPts,
00164 fixedFields, fixedFieldEvalPts,
00165 fixedParams, fixedParamEvalPts,
00166 context, compType);
00167 }
00168
00169
00170
00171
00172 DerivSet SymbPreprocessor::setupVariations(const Expr& expr,
00173 const Expr& vars,
00174 const Expr& varEvalPts,
00175 const Expr& unks,
00176 const Expr& unkEvalPts,
00177 const Expr& unkParams,
00178 const Expr& unkParamEvalPts,
00179 const Expr& fixedFields,
00180 const Expr& fixedFieldEvalPts,
00181 const Expr& fixedParams,
00182 const Expr& fixedParamEvalPts,
00183 const EvalContext& context,
00184 const ComputationType& compType)
00185 {
00186 TimeMonitor t(preprocTimer());
00187 Tabs tab;
00188
00189 const EvaluatableExpr* e
00190 = dynamic_cast<const EvaluatableExpr*>(expr.ptr().get());
00191
00192 Array<Set<MultiSet<int> > > funcDerivs(3);
00193 Array<Set<MultiIndex> > spatialDerivs(3);
00194
00195 int verb=context.setupVerbosity();
00196 SUNDANCE_BANNER1(verb, tab, "in setupVariations()");
00197 verbosity<EvaluatableExpr>() = verb;
00198 SUNDANCE_MSG1(verb,
00199 tab << "************ setting up variations of expr: "
00200 << expr
00201 << std::endl << tab << "context is " << context
00202 << std::endl << tab << "conp type is " << compType
00203 << std::endl << tab << "vars are " << vars
00204 << std::endl << tab << "unks are " << unks
00205 << std::endl << tab << "unk parameters " << unkParams
00206 << std::endl << tab << "fixed parameters " << fixedParams
00207 << std::endl << tab << "the eval points for the vars are "
00208 << varEvalPts
00209 << std::endl << tab << "the eval points for the unks are "
00210 << unkEvalPts
00211 << std::endl << tab
00212 << "the eval points for the unknown parameters are "
00213 << unkParamEvalPts
00214 << std::endl << tab
00215 << "the eval points for the fixed parameters are "
00216 << fixedParamEvalPts
00217 << tab << std::endl);
00218
00219 TEUCHOS_TEST_FOR_EXCEPTION(e==0, std::logic_error,
00220 "Non-evaluatable expr " << expr.toString()
00221 << " given to SymbPreprocessor::setupExpr()");
00222
00223
00224 Expr v = vars.flatten();
00225 Expr v0 = varEvalPts.flatten();
00226 Expr u = unks.flatten();
00227 Expr u0 = unkEvalPts.flatten();
00228 Expr alpha = unkParams.flatten();
00229 Expr alpha0 = unkParamEvalPts.flatten();
00230 Expr beta = fixedParams.flatten();
00231 Expr beta0 = fixedParamEvalPts.flatten();
00232 Expr f = fixedFields.flatten();
00233 Expr f0 = fixedFieldEvalPts.flatten();
00234
00235
00236
00237 Set<int> varID = processInputFuncs<SymbolicFuncElement>(v, v0);
00238
00239 Set<int> unkID = processInputFuncs<UnknownFuncElement>(u, u0);
00240
00241 Set<int> fixedID = processInputFuncs<UnknownFuncElement>(f, f0);
00242
00243 Set<int> unkParamID
00244 = processInputParams<UnknownParameterElement>(alpha, alpha0);
00245
00246 Set<int> fixedParamID
00247 = processInputParams<UnknownParameterElement>(beta, beta0);
00248
00249
00250
00251
00252
00253
00254
00255 SUNDANCE_MSG2(verb, tab << "forming active set");
00256 Array<Sundance::Set<MultiSet<int> > > activeFuncIDs(3);
00257 if (context.needsDerivOrder(0)) activeFuncIDs[0].put(MultiSet<int>());
00258 if (context.topLevelDiffOrder() >= 1)
00259 {
00260 for (Set<int>::const_iterator i=varID.begin(); i != varID.end(); i++)
00261 {
00262 if (context.needsDerivOrder(1)) activeFuncIDs[1].put(makeMultiSet<int>(*i));
00263 if (context.topLevelDiffOrder()==2)
00264 {
00265 for (Set<int>::const_iterator j=unkID.begin(); j != unkID.end(); j++)
00266 {
00267 activeFuncIDs[2].put(makeMultiSet<int>(*i, *j));
00268 }
00269 if (compType==MatrixAndVector)
00270 {
00271 for (Set<int>::const_iterator
00272 j=unkParamID.begin(); j != unkParamID.end(); j++)
00273 {
00274 activeFuncIDs[2].put(makeMultiSet<int>(*i, *j));
00275 }
00276 }
00277 else if (compType==Sensitivities)
00278 {
00279 for (Set<int>::const_iterator
00280 j=fixedParamID.begin(); j != fixedParamID.end(); j++)
00281 {
00282 activeFuncIDs[2].put(makeMultiSet<int>(*i, *j));
00283 }
00284 }
00285 }
00286 }
00287 }
00288 SUNDANCE_MSG1(verb, tab << std::endl << tab
00289 << " ************* Finding nonzeros for expr " << std::endl << tab);
00290 for (int i=0; i<=context.topLevelDiffOrder(); i++)
00291 {
00292 Tabs tab2;
00293 SUNDANCE_MSG4(verb, tab2 << "diff order=" << i << ", active funcs="
00294 << activeFuncIDs[i]);
00295 }
00296
00297 Set<MultiIndex> miSet;
00298 miSet.put(MultiIndex());
00299 e->registerSpatialDerivs(context, miSet);
00300
00301 SUNDANCE_MSG2(verb,
00302 tab << std::endl << tab
00303 << " ************* finding required functions" << std::endl << tab);
00304 SUNDANCE_MSG2(verb,
00305 tab << "activeFuncIDs are = " << activeFuncIDs);
00306 SUNDANCE_MSG2(verb,
00307 tab << "spatial derivs are = " << spatialDerivs);
00308 Array<Set<MultipleDeriv> > RInput
00309 = e->computeInputR(context, activeFuncIDs, spatialDerivs);
00310
00311 SUNDANCE_MSG3(verb,
00312 tab << std::endl << tab
00313 << " ************* Top-level required funcs are " << RInput << std::endl << tab);
00314
00315
00316 SUNDANCE_MSG2(verb,
00317 tab << std::endl << tab
00318 << " ************* Calling determineR()" << std::endl << tab);
00319
00320 e->determineR(context, RInput);
00321
00322
00323 SUNDANCE_MSG1(verb,
00324 tab << std::endl << tab
00325 << " ************* Finding sparsity structure " << std::endl << tab);
00326 DerivSet derivs = e->sparsitySuperset(context)->derivSet();
00327 SUNDANCE_MSG1(verb,
00328 tab << std::endl << tab
00329 << "Nonzero deriv set = " << derivs);
00330
00331 SUNDANCE_MSG1(verb,
00332 tab << std::endl << tab
00333 << " ************* Setting up evaluators for expr " << std::endl << tab);
00334
00335 int saveVerb = context.setupVerbosity();
00336 context.setSetupVerbosity(0);
00337 e->setupEval(context);
00338
00339 if (verb>1)
00340 {
00341 SUNDANCE_MSG1(verb,
00342 tab << std::endl << tab
00343 << " ************* Nonzeros are:");
00344 e->displayNonzeros(Out::os(), context);
00345
00346 }
00347
00348
00349 context.setSetupVerbosity(saveVerb);
00350 return derivs;
00351 }
00352
00353
00354 namespace Sundance
00355 {
00356 Expr makeZeros(const Expr& v)
00357 {
00358 Array<Expr> z(v.size());
00359 for (int i=0; i<v.size(); i++)
00360 {
00361 if (v[i].size()==1) z[i] = new ZeroExpr();
00362 else z[i] = makeZeros(v[i]);
00363 }
00364 return new ListExpr(z);
00365 }
00366 }
