queso-0.53.0
GcmJointTildeInfo.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2015 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/GcmJointTildeInfo.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 template <class S_V,class S_M,class D_V,class D_M,class P_V,class P_M,class Q_V,class Q_M>
33  const GpmsaComputerModelOptions& gcmOptionsObj,
36  :
37  m_env (jj.m_env),
38  m_Bmat_tilde (m_env,e.m_y_space.map(),jj.m_Bmat_rank),
39  m_Bmat_tilde_rank (0),
40  m_vu_tilde_space (m_env, "vu_tilde_", jj.m_Bmat_rank, NULL), // rr0 check
41  m_Lbmat (m_env,m_vu_tilde_space.map(),jj.m_Bmat_with_permut->numCols()), // rr0 check
42  m_Btildet_Wy_Btilde (m_vu_tilde_space.zeroVector()),
43  m_Btildet_Wy_Btilde_inv(m_vu_tilde_space.zeroVector()),
44  m_Zvec_tilde_hat_vu (m_vu_tilde_space.zeroVector()),
45  m_a_y_modifier_tilde (0),
46  m_b_y_modifier_tilde (0)
47 {
48  std::set<unsigned int> tmpSet;
49  tmpSet.insert(m_env.subId());
50  //******************************************************************************
51  // Tilde situation: form 'm_Bmat_tilde'
52  // Tilde situation: form 'm_vu_tilde_space'
53  // Tilde situation: form 'm_Lbmat'
54  //******************************************************************************
55  if (jj.m_Bmat_with_permut->numRowsGlobal() >= jj.m_Bmat_with_permut->numCols()) {
56  D_M matbU(m_env,e.m_y_space.map(),jj.m_vu_size); // same as m_Bmat
57  matbU = jj.m_Bmat_with_permut->svdMatU();
58  unsigned int buMatRank = matbU.rank(0.,1.e-8 ); // todo: should be an option
59  unsigned int buMatRank14 = matbU.rank(0.,1.e-14);
60  if (m_env.subDisplayFile()) {
61  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
62  << ": matbU.numRowsLocal() = " << matbU.numRowsLocal()
63  << ", matbU.numCols() = " << matbU.numCols()
64  << ", matbU.rank(0.,1.e-8) = " << buMatRank
65  << ", matbU.rank(0.,1.e-14) = " << buMatRank14
66  << std::endl;
67  }
68 
69  if (m_env.checkingLevel() >= 1) {
70  D_M matbUcheck(jj.m_vu_space.zeroVector());
71  D_V vecI(e.m_y_space.zeroVector());
72  D_V vecJ(e.m_y_space.zeroVector());
73  for (unsigned int i = 0; i < matbU.numCols(); ++i) {
74  matbU.getColumn(i,vecI);
75  for (unsigned int j = i; j < matbU.numCols(); ++j) {
76  matbU.getColumn(j,vecJ);
77  matbUcheck(i,j) = scalarProduct(vecI,vecJ);
78  }
79  }
80  matbUcheck.setPrintHorizontally(false);
81  if (m_env.subDisplayFile()) {
82  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
83  << ": m_Bmat_with_permut->numRowsLocal() = " << jj.m_Bmat_with_permut->numRowsLocal()
84  << ", m_Bmat_with_permut->numCols() = " << jj.m_Bmat_with_permut->numCols()
85  << ", m_Bmat_rank = " << jj.m_Bmat_rank
86  << ", matbU.numRowsLocal() = " << matbU.numRowsLocal()
87  << ", matbU.numCols() = " << matbU.numCols()
88  << ", buMatrank(0.,1.e-8) = " << buMatRank
89  << ", buMatrank(0.,1.e-14) = " << buMatRank14
90  << ", matbUcheck.numRowsLocal() = " << matbUcheck.numRowsLocal()
91  << ", matbUcheck.numCols() = " << matbUcheck.numCols()
92  << ", matbUcheck =\n" << matbUcheck
93  << std::endl;
94  }
95  }
96 
97  D_V vecbJ(e.m_y_space.zeroVector());
98  for (unsigned int j = 0; j < jj.m_Bmat_rank; ++j) {
99  matbU.getColumn(j,vecbJ);
100  m_Bmat_tilde.setColumn(j,vecbJ);
101  }
102  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
103  m_Bmat_tilde.subWriteContents("Btilde",
104  "Btilde1",
105  "m",
106  tmpSet);
107  }
108  if (m_env.subDisplayFile()) {
109  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
110  << ": m_Bmat_tilde computed (1)"
111  << std::endl;
112  }
113  }
114  else {
115  D_M Bmat_t(m_env,jj.m_vu_space.map(),e.m_paper_n_y); // same as m_Bmat^T
116  Bmat_t.fillWithTranspose(0,0,*jj.m_Bmat_with_permut,true,true);
117 
118  D_M matbV(e.m_y_space.zeroVector());
119  matbV = Bmat_t.svdMatV();
120  unsigned int bvMatRank = matbV.rank(0.,1.e-8); // todo: should be an option
121  unsigned int bvMatRank14 = matbV.rank(0.,1.e-14);
122  if (m_env.subDisplayFile()) {
123  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
124  << ": matbV.numRowsLocal() = " << matbV.numRowsLocal()
125  << ", matbV.numCols() = " << matbV.numCols()
126  << ", matbV.rank(0.,1.e-8) = " << bvMatRank
127  << ", matbV.rank(0.,1.e-14) = " << bvMatRank14
128  << std::endl;
129  }
130 
131  if (m_env.checkingLevel() >= 1) {
132  D_M matbVcheck(e.m_y_space.zeroVector());
133  D_V vecI(e.m_y_space.zeroVector());
134  D_V vecJ(e.m_y_space.zeroVector());
135  for (unsigned int i = 0; i < matbV.numCols(); ++i) {
136  matbV.getColumn(i,vecI);
137  for (unsigned int j = i; j < matbV.numCols(); ++j) {
138  matbV.getColumn(j,vecJ);
139  matbVcheck(i,j) = scalarProduct(vecI,vecJ);
140  }
141  }
142  matbVcheck.setPrintHorizontally(false);
143  if (m_env.subDisplayFile()) {
144  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
145  << ": m_Bmat_with_permut->numRowsLocal() = " << jj.m_Bmat_with_permut->numRowsLocal()
146  << ", m_Bmat_with_permut->numCols() = " << jj.m_Bmat_with_permut->numCols()
147  << ", m_Bmat_rank = " << jj.m_Bmat_rank
148  << ", matbV.numRowsLocal() = " << matbV.numRowsLocal()
149  << ", matbV.numCols() = " << matbV.numCols()
150  << ", bvMatrank(0.,1.e-8) = " << bvMatRank
151  << ", bvMatrank(0.,1.e-14) = " << bvMatRank14
152  << ", matbVcheck.numRowsLocal() = " << matbVcheck.numRowsLocal()
153  << ", matbVcheck.numCols() = " << matbVcheck.numCols()
154  << ", matbVcheck =\n" << matbVcheck
155  << std::endl;
156  }
157  }
158 
159  D_V vecbJ(e.m_y_space.zeroVector());
160  for (unsigned int j = 0; j < jj.m_Bmat_rank; ++j) {
161  matbV.getColumn(j,vecbJ);
162  m_Bmat_tilde.setColumn(j,vecbJ);
163  }
164  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
165  m_Bmat_tilde.subWriteContents("Btilde",
166  "Btilde2",
167  "m",
168  tmpSet);
169  }
170  if (m_env.subDisplayFile()) {
171  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
172  << ": m_Bmat_tilde computed (2)"
173  << std::endl;
174  }
175  }
176 
177  if (m_env.subDisplayFile()) {
178  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
179  << ": finished forming 'm_Bmat_tilde'"
180  << ", m_Bmat_tilde.numRowsLocal() = " << m_Bmat_tilde.numRowsLocal()
181  << ", m_Bmat_tilde.numCols() = " << m_Bmat_tilde.numCols()
182  << std::endl;
183  }
184 
186  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
187  m_Lbmat.subWriteContents("Lbmat",
188  "Lbmat",
189  "m",
190  tmpSet);
191  }
192  if (m_env.subDisplayFile()) {
193  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
194  << ": m_Lbmat_tilde computed"
195  << std::endl;
196  }
197 
198  //********************************************************************************
199  // Form 'Btilde^T' matrix
200  //********************************************************************************
201  D_M Btildet(m_env,m_vu_tilde_space.map(),e.m_paper_n_y);
202  Btildet.fillWithTranspose(0,0,m_Bmat_tilde,true,true);
203 
204  if (m_env.checkingLevel() >= 1) {
205  // Check transpose operation
206  D_M Btildett(m_env,e.m_y_space.map(),m_vu_tilde_space.dimGlobal());
207  Btildett.fillWithTranspose(0,0,Btildet,true,true);
208  Btildett -= m_Bmat_tilde;
209  double btDiffNorm = Btildett.normFrob();
210  if (m_env.subDisplayFile()) {
211  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
212  << ": ||Btildett - Btilde||_2 = " << btDiffNorm
213  << std::endl;
214  }
215  }
216 
217  //********************************************************************************
218  // Compute 'Btilde' rank
219  //********************************************************************************
220  double bTildeRank14 = 0.;
221  if (m_Bmat_tilde.numRowsGlobal() >= m_Bmat_tilde.numCols()) {
222  m_Bmat_tilde_rank = m_Bmat_tilde.rank(0.,1.e-8 ); // todo: should be an option
223  bTildeRank14 = m_Bmat_tilde.rank(0.,1.e-14);
224  }
225  else {
226  m_Bmat_tilde_rank = Btildet.rank(0.,1.e-8 ); // todo: should be an option
227  bTildeRank14 = Btildet.rank(0.,1.e-14);
228  }
229 
230  if (m_env.subDisplayFile()) {
231  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
232  << ": m_Bmat_tilde.numRowsLocal() = " << m_Bmat_tilde.numRowsLocal()
233  << ", m_Bmat_tilde.numCols() = " << m_Bmat_tilde.numCols()
234  << ", m_Bmat_tilde.rank(0.,1.e-8) = " << m_Bmat_tilde_rank
235  << ", m_Bmat_tilde.rank(0.,1.e-14) = " << bTildeRank14
236  << std::endl;
237  }
238  queso_require_equal_to_msg(m_Bmat_tilde_rank, std::min(m_Bmat_tilde.numRowsGlobal(),m_Bmat_tilde.numCols()), "'m_Bmat_tilde' does not have a proper rank");
239 
240  //******************************************************************************
241  // Tilde situation: compute 'Btilde^T W_y Btilde' matrix, and its inverse
242  //******************************************************************************
243  m_Btildet_Wy_Btilde = Btildet * (*e.m_Wy * m_Bmat_tilde); // todo: add 1.e-4 to diagonal
244  if (m_env.subDisplayFile()) {
245  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
246  << ": finished computing 'm_Btildet_Wy_Btilde'"
247  << std::endl;
248  }
249  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
250  m_Btildet_Wy_Btilde.subWriteContents("Btildet_Wy_Btilde",
251  "Btildet_Wy_Btilde",
252  "m",
253  tmpSet);
254  }
255  if (m_env.subDisplayFile()) {
256  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
257  << ": m_Btildet_Wy_Btilde computed"
258  << std::endl;
259  }
260 
261  double btildetWyBtildeLnDeterminant = m_Btildet_Wy_Btilde.lnDeterminant();
262  unsigned int btildetWyBtildeRank = m_Btildet_Wy_Btilde.rank(0.,1.e-8 ); // todo: should be an option
263  unsigned int btildetWyBtildeRank14 = m_Btildet_Wy_Btilde.rank(0.,1.e-14);
264  if (m_env.subDisplayFile()) {
265  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
266  << ": m_Btildet_Wy_Btilde.numRowsLocal() = " << m_Btildet_Wy_Btilde.numRowsLocal()
267  << ", m_Btildet_Wy_Btilde.numCols() = " << m_Btildet_Wy_Btilde.numCols()
268  << ", m_Btildet_Wy_Btilde.lnDeterminant() = " << btildetWyBtildeLnDeterminant
269  << ": m_Btildet_Wy_Btilde.rank(0.,1.e-8) = " << btildetWyBtildeRank
270  << ": m_Btildet_Wy_Btilde.rank(0.,1.e-14) = " << btildetWyBtildeRank14
271  << std::endl;
272  }
273 
274  m_Btildet_Wy_Btilde_inv = m_Btildet_Wy_Btilde.inverse(); // todo: add 1.e-6 to diagonal
275  if (m_env.subDisplayFile()) {
276  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
277  << ": finished computing 'm_Btildet_Wy_Btilde_inv'"
278  << ", m_Btildet_Wy_Btilde_inv.lnDeterminant() = " << m_Btildet_Wy_Btilde_inv.lnDeterminant()
279  << std::endl;
280  }
281  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
282  m_Btildet_Wy_Btilde_inv.subWriteContents("Btildet_Wy_Btilde_inv",
283  "Btildet_Wy_Btilde_inv",
284  "m",
285  tmpSet);
286  }
287  if (m_env.subDisplayFile()) {
288  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
289  << ": m_Btildet_Wy_Btilde_inv computed"
290  << std::endl;
291  }
292 
293  double btildetWyBtildeInvLnDeterminant = m_Btildet_Wy_Btilde_inv.lnDeterminant();
294  unsigned int btildetWyBtildeInvRank = m_Btildet_Wy_Btilde_inv.rank(0.,1.e-8 ); // todo: should be an option
295  unsigned int btildetWyBtildeInvRank14 = m_Btildet_Wy_Btilde_inv.rank(0.,1.e-14);
296  if (m_env.subDisplayFile()) {
297  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
298  << ": m_Btildet_Wy_Btilde_inv.numRowsLocal() = " << m_Btildet_Wy_Btilde_inv.numRowsLocal()
299  << ", m_Btildet_Wy_Btilde_inv.numCols() = " << m_Btildet_Wy_Btilde_inv.numCols()
300  << ": m_Btildet_Wy_Btilde_inv.lnDeterminant() = " << btildetWyBtildeInvLnDeterminant
301  << ": m_Btildet_Wy_Btilde_inv.rank(0.,1.e-8) = " << btildetWyBtildeInvRank
302  << ": m_Btildet_Wy_Btilde_inv.rank(0.,1.e-14) = " << btildetWyBtildeInvRank14
303  << std::endl;
304  }
305 
306  //********************************************************************************
307  // Compute 'tilde' exponent modifiers
308  //********************************************************************************
309  m_a_y_modifier_tilde = ((double) (e.m_paper_n_y - m_Bmat_tilde_rank)) / 2.;
310  if (m_env.subDisplayFile()) {
311  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
312  << ": m_a_y_modifier_tilde = " << m_a_y_modifier_tilde
313  << std::endl;
314  }
315 
316  D_V yVec_transformed(e.m_experimentStorage.yVec_transformed());
317  if (m_env.subDisplayFile()) {
318  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
319  << ": m_Zvec_tilde_hat_vu.sizeLocal() = " << m_Zvec_tilde_hat_vu.sizeLocal()
320  << ", m_Btildet_Wy_Btilde_inv.numRowsLocal() = " << m_Btildet_Wy_Btilde_inv.numRowsLocal()
321  << ", m_Btildet_Wy_Btilde_inv.numCols() = " << m_Btildet_Wy_Btilde_inv.numCols()
322  << ", Btildet.numRowsLocal() = " << Btildet.numRowsLocal()
323  << ", Btildet.numCols() = " << Btildet.numCols()
324  << ", m_Wy->numRowsLocal() = " << e.m_Wy->numRowsLocal()
325  << ", m_Wy->numCols() = " << e.m_Wy->numCols()
326  << ", yVec_transformed.sizeLocal() = " << yVec_transformed.sizeLocal()
327  << std::endl;
328  }
329  m_Zvec_tilde_hat_vu = m_Btildet_Wy_Btilde_inv * (Btildet * (*e.m_Wy * yVec_transformed));
330  D_V tmpVec2(yVec_transformed - (m_Bmat_tilde * m_Zvec_tilde_hat_vu));
331  tmpVec2 = *e.m_Wy * tmpVec2;
332  m_b_y_modifier_tilde = scalarProduct(yVec_transformed,tmpVec2) / 2.;
333 
334  if (m_env.subDisplayFile()) {
335  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
336  << ": m_b_y_modifier_tilde = " << m_b_y_modifier_tilde
337  << std::endl;
338  }
339 
340  if (m_env.subDisplayFile()) {
341  *m_env.subDisplayFile() << "In GcmJointTildeInfo<S_V,S_M,D_V,D_M,P_V,P_M,Q_V,Q_M>::constructor()"
342  << ": finished computing 'tilde' exponent modifiers"
343  << std::endl;
344  }
345 }
346 
347 template <class S_V,class S_M,class D_V,class D_M,class P_V,class P_M,class Q_V,class Q_M>
349 {
350 }
351 
352 } // End namespace QUESO
353 
const BaseEnvironment & m_env
const Map & map() const
Map.
Definition: VectorSpace.C:157
VectorSpace< D_V, D_M > m_vu_tilde_space
std::set< unsigned int > m_dataOutputAllowedSet
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:189
unsigned int dimGlobal() const
Definition: VectorSpace.C:176
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
const ExperimentStorage< S_V, S_M, D_V, D_M > & m_experimentStorage
double scalarProduct(const GslVector &x, const GslVector &y)
Definition: GslVector.C:1150
VectorSpace< D_V, D_M > m_vu_space
Definition: GcmJointInfo.h:72
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:274
GcmJointTildeInfo(const GpmsaComputerModelOptions &gcmOptionsObj, const GcmExperimentInfo< S_V, S_M, D_V, D_M, P_V, P_M > &e, const GcmJointInfo< S_V, S_M, D_V, D_M, P_V, P_M, Q_V, Q_M > &jj)
unsigned int m_Bmat_rank
Definition: GcmJointInfo.h:111
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:295
VectorSpace< D_V, D_M > m_y_space
unsigned int checkingLevel() const
Access function to private attribute m_checkingLevel.
Definition: Environment.C:410
unsigned int m_vu_size
Definition: GcmJointInfo.h:71

Generated on Thu Jun 11 2015 13:52:32 for queso-0.53.0 by  doxygen 1.8.5