Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ResponseStatisticModelEvaluator.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
46#include "Teuchos_Assert.hpp"
47
49 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
50 const Teuchos::Array< Teuchos::RCP<const Epetra_Map> >& base_g_maps_,
51 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
52 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
53 const Teuchos::RCP<const Epetra_BlockMap>& block_map_) :
54 me(me_),
55 base_g_maps(base_g_maps_),
56 sg_basis(sg_basis_),
57 sg_comm(sg_comm_),
58 block_map(block_map_),
59 num_p(0),
60 num_g(0)
61{
62 InArgs me_inargs = me->createInArgs();
63 OutArgs me_outargs = me->createOutArgs();
64 num_p = me_inargs.Np();
65 num_g = me_outargs.Ng();
66}
67
68// Overridden from EpetraExt::ModelEvaluator
69
70Teuchos::RCP<const Epetra_Map>
72get_x_map() const
73{
74 return Teuchos::null;
75}
76
77Teuchos::RCP<const Epetra_Map>
79get_f_map() const
80{
81 return Teuchos::null;
82}
83
84Teuchos::RCP<const Epetra_Map>
86get_p_map(int l) const
87{
88 TEUCHOS_TEST_FOR_EXCEPTION(
89 l >= num_p || l < 0, std::logic_error,
90 std::endl <<
91 "Error! Stokhos::ResponseStatisticModelEvaluator::get_p_map(): " <<
92 "Invalid parameter index l = " << l << std::endl);
93
94 return me->get_p_map(l);
95}
96
97Teuchos::RCP<const Epetra_Map>
99get_g_map(int l) const
100{
101 TEUCHOS_TEST_FOR_EXCEPTION(
102 l >= 3*num_g || l < 0, std::logic_error,
103 std::endl <<
104 "Error! Stokhos::ResponseStatisticModelEvaluator::get_g_map(): " <<
105 "Invalid response index l = " << l << std::endl);
106
107 if (l < num_g)
108 return me->get_g_map(l);
109 else if (l < 2*num_g)
110 return base_g_maps[l-num_g];
111 else
112 return base_g_maps[l-2*num_g];
113}
114
115Teuchos::RCP<const Teuchos::Array<std::string> >
117get_p_names(int l) const
118{
119 TEUCHOS_TEST_FOR_EXCEPTION(
120 l >= num_p || l < 0, std::logic_error,
121 std::endl <<
122 "Error! Stokhos::ResponseStatisticModelEvaluator::get_p_names(): " <<
123 "Invalid parameter index l = " << l << std::endl);
124
125 return me->get_p_names(l);
126}
127
128Teuchos::RCP<const Epetra_Vector>
130get_p_init(int l) const
131{
132 TEUCHOS_TEST_FOR_EXCEPTION(
133 l >= num_p || l < 0, std::logic_error,
134 std::endl <<
135 "Error! Stokhos::ResponseStatisticModelEvaluator::get_p_init(): " <<
136 "Invalid parameter index l = " << l << std::endl);
137
138 return me->get_p_init(l);
139}
140
141EpetraExt::ModelEvaluator::InArgs
143{
144 InArgsSetup inArgs;
145 InArgs me_inargs = me->createInArgs();
146
147 inArgs.setModelEvalDescription(this->description());
148 inArgs.set_Np(num_p);
149
150 // Support pass-through of sg data
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));
156
157 return inArgs;
158}
159
160EpetraExt::ModelEvaluator::OutArgs
162{
163 OutArgsSetup outArgs;
164 OutArgs me_outargs = me->createOutArgs();
165
166 outArgs.setModelEvalDescription(this->description());
167
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));
173
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);
184 }
185 }
186
187 return outArgs;
188}
189
190void
192 const OutArgs& outArgs) const
193{
194 // Create underlying inargs
195 InArgs me_inargs = me->createInArgs();
196
197 // Pass parameters
198 for (int i=0; i<num_p; i++)
199 me_inargs.set_p(i, inArgs.get_p(i));
200
201 // Pass SG data
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());
208
209 // Create underlying outargs
210 OutArgs me_outargs = me->createOutArgs();
211
212 Teuchos::Array<bool> need_g_var(num_g);
213
214 // Responses
215 for (int i=0; i<num_g; i++) {
216
217 // See if we need g_var for d/dp(g_var) below
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;
223
224 // g
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) {
230 g = Teuchos::rcp(new Epetra_Vector(*(me->get_g_map(i))));
231 }
232 me_outargs.set_g(i, g);
233
234 // dg/dp
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) &&
243 dgdp.isEmpty()) {
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();
250 else
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)
256 dgdp =
257 Derivative(Teuchos::rcp(new Epetra_MultiVector(
258 *p_map,
259 g_map->NumMyElements())),
260 DERIV_TRANS_MV_BY_ROW);
261 else
262 dgdp =
263 Derivative(Teuchos::rcp(new Epetra_MultiVector(
264 *g_map,
265 p_map->NumMyElements())),
266 DERIV_MV_BY_COL);
267
268 }
269 me_outargs.set_DgDp(i, j, dgdp);
270
271 // We need g to compute d/dp(g_var)
272 if (dgdp_var != Teuchos::null && g == Teuchos::null) {
273 g = Teuchos::rcp(new Epetra_Vector(*(me->get_g_map(i))));
274 me_outargs.set_g(i, g);
275 }
276 }
277 }
278
279 }
280
281 // Compute the functions
282 me->evalModel(me_inargs, me_outargs);
283
284 // Compute statistics
285 for (int i=0; i<num_g; i++) {
286
287 // g
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)
292 g_var = Teuchos::rcp(new Epetra_Vector(*(base_g_maps[i])));
293 Teuchos::RCP<EpetraVectorOrthogPoly> g_sg;
294 if (g_mean != Teuchos::null || g_var != Teuchos::null) {
295 g_sg = Teuchos::rcp(new EpetraVectorOrthogPoly(
296 sg_basis, block_map, base_g_maps[i], me->get_g_map(i),
297 sg_comm, View, *g));
298 }
299
300 if (g_mean != Teuchos::null)
301 g_sg->computeMean(*g_mean);
302 if (g_var != Teuchos::null)
303 g_sg->computeVariance(*g_var);
304
305 // dg/dp
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();
317 EpetraMultiVectorOrthogPoly X_sg(sg_basis, block_map,
318 base_g_maps[i], sg_comm, n);
319 dgdp_op->SetUseTranspose(true);
320 if (dgdp_mean != Teuchos::null) {
321 X_sg.init(0.0);
322 for (int l=0; l<n; l++)
323 X_sg[0][l][l] = 1.0;
324 TEUCHOS_TEST_FOR_EXCEPTION(
325 outArgs.get_DgDp(i+num_g,j).getMultiVectorOrientation() == DERIV_MV_BY_COL,
326 std::logic_error,
327 "Error! ResponseStatisticModelEvaluator does not support " <<
328 " dg/dp DERIV_MV_BY_COL for underlying operator dg/dp");
329 dgdp_op->Apply(*(X_sg.getBlockMultiVector()), *dgdp_mean);
330 }
331 if (dgdp_var != Teuchos::null) {
332 X_sg.init(0.0);
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,
338 std::logic_error,
339 "Error! ResponseStatisticModelEvaluator does not support " <<
340 " dg/dp DERIV_MV_BY_COL for underlying operator dg/dp");
341 dgdp_op->Apply(*(X_sg.getBlockMultiVector()), *dgdp_var);
342 }
343
344 }
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);
351 else
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)
357 dgdp_sg.computeMean(*dgdp_mean);
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);
365 }
366 }
367 }
368 }
369 }
370 }
371}
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.
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.
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.
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.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)