46#include "Teuchos_Assert.hpp"
49 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
50 const Teuchos::Array< Teuchos::RCP<const Epetra_Map> >& base_g_maps_,
52 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
53 const Teuchos::RCP<const Epetra_BlockMap>& block_map_) :
55 base_g_maps(base_g_maps_),
58 block_map(block_map_),
62 InArgs me_inargs =
me->createInArgs();
63 OutArgs me_outargs =
me->createOutArgs();
64 num_p = me_inargs.Np();
65 num_g = me_outargs.Ng();
70Teuchos::RCP<const Epetra_Map>
77Teuchos::RCP<const Epetra_Map>
84Teuchos::RCP<const Epetra_Map>
88 TEUCHOS_TEST_FOR_EXCEPTION(
89 l >= num_p || l < 0, std::logic_error,
91 "Error! Stokhos::ResponseStatisticModelEvaluator::get_p_map(): " <<
92 "Invalid parameter index l = " << l << std::endl);
94 return me->get_p_map(l);
97Teuchos::RCP<const Epetra_Map>
101 TEUCHOS_TEST_FOR_EXCEPTION(
102 l >= 3*num_g || l < 0, std::logic_error,
104 "Error! Stokhos::ResponseStatisticModelEvaluator::get_g_map(): " <<
105 "Invalid response index l = " << l << std::endl);
108 return me->get_g_map(l);
109 else if (l < 2*num_g)
110 return base_g_maps[l-num_g];
112 return base_g_maps[l-2*num_g];
115Teuchos::RCP<const Teuchos::Array<std::string> >
119 TEUCHOS_TEST_FOR_EXCEPTION(
120 l >= num_p || l < 0, std::logic_error,
122 "Error! Stokhos::ResponseStatisticModelEvaluator::get_p_names(): " <<
123 "Invalid parameter index l = " << l << std::endl);
125 return me->get_p_names(l);
128Teuchos::RCP<const Epetra_Vector>
132 TEUCHOS_TEST_FOR_EXCEPTION(
133 l >= num_p || l < 0, std::logic_error,
135 "Error! Stokhos::ResponseStatisticModelEvaluator::get_p_init(): " <<
136 "Invalid parameter index l = " << l << std::endl);
138 return me->get_p_init(l);
141EpetraExt::ModelEvaluator::InArgs
145 InArgs me_inargs = me->createInArgs();
147 inArgs.setModelEvalDescription(this->description());
148 inArgs.set_Np(num_p);
151 inArgs.setSupports(IN_ARG_sg_basis, me_inargs.supports(IN_ARG_sg_basis));
152 inArgs.setSupports(IN_ARG_sg_quadrature,
153 me_inargs.supports(IN_ARG_sg_quadrature));
154 inArgs.setSupports(IN_ARG_sg_expansion,
155 me_inargs.supports(IN_ARG_sg_expansion));
160EpetraExt::ModelEvaluator::OutArgs
163 OutArgsSetup outArgs;
164 OutArgs me_outargs = me->createOutArgs();
166 outArgs.setModelEvalDescription(this->description());
168 outArgs.set_Np_Ng(num_p, 3*num_g);
169 for (
int i=0; i<num_g; i++) {
170 for (
int j=0;
j<num_p;
j++) {
171 outArgs.setSupports(OUT_ARG_DgDp, i,
j,
172 me_outargs.supports(OUT_ARG_DgDp, i,
j));
174 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDp, i,
j);
175 DerivativeSupport ds_stat;
176 if (ds.supports(DERIV_MV_BY_COL))
177 ds_stat.plus(DERIV_MV_BY_COL);
178 if (ds.supports(DERIV_TRANS_MV_BY_ROW))
179 ds_stat.plus(DERIV_TRANS_MV_BY_ROW);
180 if (ds.supports(DERIV_LINEAR_OP))
181 ds_stat.plus(DERIV_TRANS_MV_BY_ROW);
182 outArgs.setSupports(OUT_ARG_DgDp, i+num_g,
j, ds_stat);
183 outArgs.setSupports(OUT_ARG_DgDp, i+2*num_g,
j, ds_stat);
192 const OutArgs& outArgs)
const
195 InArgs me_inargs = me->createInArgs();
198 for (
int i=0; i<num_p; i++)
199 me_inargs.set_p(i, inArgs.get_p(i));
202 if (me_inargs.supports(IN_ARG_sg_basis))
203 me_inargs.set_sg_basis(inArgs.get_sg_basis());
204 if (me_inargs.supports(IN_ARG_sg_quadrature))
205 me_inargs.set_sg_quadrature(inArgs.get_sg_quadrature());
206 if (me_inargs.supports(IN_ARG_sg_expansion))
207 me_inargs.set_sg_expansion(inArgs.get_sg_expansion());
210 OutArgs me_outargs = me->createOutArgs();
212 Teuchos::Array<bool> need_g_var(num_g);
215 for (
int i=0; i<num_g; i++) {
218 need_g_var[i] =
false;
219 for (
int j=0;
j<num_p;
j++)
220 if (!outArgs.supports(OUT_ARG_DgDp, i,
j).none())
221 if (!outArgs.get_DgDp(i+2*num_g,
j).isEmpty())
222 need_g_var[i] =
true;
225 Teuchos::RCP<Epetra_Vector>
g = outArgs.get_g(i);
226 Teuchos::RCP<Epetra_Vector> g_mean = outArgs.get_g(i+num_g);
227 Teuchos::RCP<Epetra_Vector> g_var = outArgs.get_g(i+2*num_g);
228 if ((g_mean != Teuchos::null || g_var != Teuchos::null || need_g_var[i]) &&
229 g == Teuchos::null) {
232 me_outargs.set_g(i,
g);
235 for (
int j=0;
j<num_p;
j++) {
236 if (!outArgs.supports(OUT_ARG_DgDp, i,
j).none()) {
237 Derivative dgdp = outArgs.get_DgDp(i,
j);
238 Teuchos::RCP<Epetra_MultiVector> dgdp_mean =
239 outArgs.get_DgDp(i+num_g,
j).getMultiVector();
240 Teuchos::RCP<Epetra_MultiVector> dgdp_var =
241 outArgs.get_DgDp(i+2*num_g,
j).getMultiVector();
242 if ((dgdp_mean != Teuchos::null || dgdp_var != Teuchos::null) &&
244 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(i);
245 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(
j);
246 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDp,i,
j);
247 EDerivativeMultiVectorOrientation mvOrientation;
248 if (dgdp_mean != Teuchos::null)
249 mvOrientation = outArgs.get_DgDp(i+num_g,
j).getMultiVectorOrientation();
251 mvOrientation = outArgs.get_DgDp(i+2*num_g,
j).getMultiVectorOrientation();
252 if (mvOrientation == DERIV_TRANS_MV_BY_ROW &&
253 ds.supports(DERIV_LINEAR_OP))
254 dgdp = Derivative(me->create_DgDp_op(i,
j));
255 else if (mvOrientation == DERIV_TRANS_MV_BY_ROW)
259 g_map->NumMyElements())),
260 DERIV_TRANS_MV_BY_ROW);
265 p_map->NumMyElements())),
269 me_outargs.set_DgDp(i,
j, dgdp);
272 if (dgdp_var != Teuchos::null &&
g == Teuchos::null) {
274 me_outargs.set_g(i,
g);
282 me->evalModel(me_inargs, me_outargs);
285 for (
int i=0; i<num_g; i++) {
288 Teuchos::RCP<Epetra_Vector>
g = me_outargs.get_g(i);
289 Teuchos::RCP<Epetra_Vector> g_mean = outArgs.get_g(i+num_g);
290 Teuchos::RCP<Epetra_Vector> g_var = outArgs.get_g(i+2*num_g);
291 if (need_g_var[i] && g_var == Teuchos::null)
293 Teuchos::RCP<EpetraVectorOrthogPoly> g_sg;
294 if (g_mean != Teuchos::null || g_var != Teuchos::null) {
296 sg_basis, block_map, base_g_maps[i], me->get_g_map(i),
300 if (g_mean != Teuchos::null)
301 g_sg->computeMean(*g_mean);
302 if (g_var != Teuchos::null)
303 g_sg->computeVariance(*g_var);
306 for (
int j=0;
j<num_p;
j++) {
307 if (!outArgs.supports(OUT_ARG_DgDp, i,
j).none()) {
308 Derivative dgdp = me_outargs.get_DgDp(i,
j);
309 Teuchos::RCP<Epetra_MultiVector> dgdp_mean =
310 outArgs.get_DgDp(i+num_g,
j).getMultiVector();
311 Teuchos::RCP<Epetra_MultiVector> dgdp_var =
312 outArgs.get_DgDp(i+2*num_g,
j).getMultiVector();
313 if (dgdp_mean != Teuchos::null || dgdp_var != Teuchos::null) {
314 if (dgdp.getLinearOp() != Teuchos::null) {
315 Teuchos::RCP<Epetra_Operator> dgdp_op = dgdp.getLinearOp();
316 int n = base_g_maps[i]->NumMyElements();
318 base_g_maps[i], sg_comm, n);
319 dgdp_op->SetUseTranspose(
true);
320 if (dgdp_mean != Teuchos::null) {
322 for (
int l=0; l<n; l++)
324 TEUCHOS_TEST_FOR_EXCEPTION(
325 outArgs.get_DgDp(i+num_g,
j).getMultiVectorOrientation() == DERIV_MV_BY_COL,
327 "Error! ResponseStatisticModelEvaluator does not support " <<
328 " dg/dp DERIV_MV_BY_COL for underlying operator dg/dp");
331 if (dgdp_var != Teuchos::null) {
333 for (
int k=1; k<sg_basis->size(); k++)
334 for (
int l=0; l<n; l++)
335 X_sg[k][l][l] = 2.0*(*g_sg)[k][l]*sg_basis->norm_squared(k);
336 TEUCHOS_TEST_FOR_EXCEPTION(
337 outArgs.get_DgDp(i+2*num_g,
j).getMultiVectorOrientation() == DERIV_MV_BY_COL,
339 "Error! ResponseStatisticModelEvaluator does not support " <<
340 " dg/dp DERIV_MV_BY_COL for underlying operator dg/dp");
345 else if (dgdp.getMultiVector() != Teuchos::null) {
346 Teuchos::RCP<const Epetra_MultiVector> dgdp_mv =
347 dgdp.getMultiVector();
348 Teuchos::RCP<const Epetra_Map> row_map;
349 if (dgdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
350 row_map = me->get_g_map(i);
352 row_map = me->get_p_map(
j);
354 sg_basis, block_map, row_map, Teuchos::rcpFromRef(dgdp_mv->Map()),
355 sg_comm,
View, *dgdp_mv);
356 if (dgdp_mean != Teuchos::null)
358 if (dgdp_var != Teuchos::null) {
359 dgdp_var->PutScalar(0.0);
360 for (
int k=1; k<sg_basis->size(); k++)
361 for (
int m=0; m<dgdp_var->NumVectors(); m++)
362 for (
int l=0; l<row_map->NumMyElements(); l++)
363 (*dgdp_var)[m][l] += 2.0*(*g_sg)[k][l]*dgdp_sg[k][m][l]*
364 sg_basis->norm_squared(k);
A container class storing an orthogonal polynomial whose coefficients are vectors,...
void computeMean(Epetra_MultiVector &v) const
Compute mean.
A container class storing an orthogonal polynomial whose coefficients are vectors,...
Abstract base class for multivariate orthogonal polynomials.
void init(const value_type &val)
Initialize coefficients.
Teuchos::RCP< EpetraExt::BlockMultiVector > getBlockMultiVector()
Get block vector.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
int num_p
Number of parameters.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
InArgs createInArgs() const
Create InArgs.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
int num_g
Number of responses.
ResponseStatisticModelEvaluator(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me, const Teuchos::Array< Teuchos::RCP< const Epetra_Map > > &base_g_maps, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const EpetraExt::MultiComm > &sg_comm, const Teuchos::RCP< const Epetra_BlockMap > &block_map)
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
OutArgs createOutArgs() const
Create OutArgs.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)