queso-0.51.1
VectorSequence.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,2009,2010,2011,2012 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 
26 #include <queso/VectorSequence.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
29 
30 namespace QUESO {
31 
32 // Default constructor -----------------------------
33 template <class V, class M>
34 BaseVectorSequence<V,M>::BaseVectorSequence(
35  const VectorSpace<V,M>& vectorSpace,
36  unsigned int subSequenceSize,
37  const std::string& name)
38  :
39  m_env (vectorSpace.env()),
40  m_vectorSpace (vectorSpace),
41  m_name (name),
42  m_fftObj (new Fft<double>(m_env)),
43  m_subMinPlain (NULL),
44  m_unifiedMinPlain (NULL),
45  m_subMaxPlain (NULL),
46  m_unifiedMaxPlain (NULL),
47  m_subMeanPlain (NULL),
48  m_unifiedMeanPlain (NULL),
49  m_subMedianPlain (NULL),
50  m_unifiedMedianPlain (NULL),
51  m_subSampleVariancePlain (NULL),
52  m_unifiedSampleVariancePlain(NULL),
53  m_subBoxPlain (NULL),
54  m_unifiedBoxPlain (NULL)
55 {
56  if (subSequenceSize) {}; // just to avoid compiler warning
57 }
58 // Destructor ---------------------------------------
59 template <class V, class M>
60 BaseVectorSequence<V,M>::~BaseVectorSequence()
61 {
62  //clear();
63  this->deleteStoredVectors();
64  if (m_fftObj != NULL) delete m_fftObj;
65 }
66 
67 // Sequence methods----------------------------------
68 template <class V, class M>
69 unsigned int
70 BaseVectorSequence<V,M>::unifiedSequenceSize() const
71 {
72  unsigned int unifiedNumSamples = 0;
73 
74  bool useOnlyInter0Comm = (m_vectorSpace.numOfProcsForStorage() == 1);
75 
76  if (useOnlyInter0Comm) {
77  if (m_env.inter0Rank() >= 0) {
78  unsigned int subNumSamples = this->subSequenceSize();
79  m_env.inter0Comm().Allreduce((void *) &subNumSamples, (void *) &unifiedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
80  "BaseVectorSequence<V,M>::unifiedSequenceSize()",
81  "failed MPI.Allreduce() for unifiedSequenceSize()");
82  }
83  else {
84  // Node not in the 'inter0' communicator
85  unifiedNumSamples = this->subSequenceSize();
86  }
87  }
88  else {
89  UQ_FATAL_TEST_MACRO((useOnlyInter0Comm == false),
90  m_env.worldRank(),
91  "BaseVectorSequence<V,M>::unifiedSequenceSize()",
92  "parallel vectors not supported yet");
93  }
94 
95  return unifiedNumSamples;
96 }
97 // --------------------------------------------------
98 template <class V, class M>
99 unsigned int
100 BaseVectorSequence<V,M>::vectorSizeLocal() const
101 {
102  return m_vectorSpace.dimLocal();
103 }
104 // --------------------------------------------------
105 template <class V, class M>
106 unsigned int
107 BaseVectorSequence<V,M>::vectorSizeGlobal() const
108 {
109  return m_vectorSpace.dimGlobal();
110 }
111 // --------------------------------------------------
112 template <class V, class M>
113 const VectorSpace<V,M>&
114 BaseVectorSequence<V,M>::vectorSpace() const
115 {
116  return m_vectorSpace;
117 }
118 // --------------------------------------------------
119 template <class V, class M>
120 const std::string&
121 BaseVectorSequence<V,M>::name() const
122 {
123  return m_name;
124 }
125 // --------------------------------------------------
126 template <class V, class M>
127 void
128 BaseVectorSequence<V,M>::setName(const std::string& newName)
129 {
130  m_name = newName;
131  return;
132 }
133 // --------------------------------------------------
134 template <class V, class M>
135 void
136 BaseVectorSequence<V,M>::clear()
137 {
138  unsigned int numPos = this->subSequenceSize();
139  if (numPos) {
140  this->resetValues(0,numPos);
141  this->resizeSequence(0);
142  }
143 
144  return;
145 }
146 // --------------------------------------------------
147 template <class V, class M>
148 const V&
149 BaseVectorSequence<V,M>::subMinPlain() const
150 {
151  if (m_subMinPlain == NULL) {
152  m_subMinPlain = m_vectorSpace.newVector();
153  if (m_subMaxPlain == NULL) m_subMaxPlain = m_vectorSpace.newVector();
154  subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
155  }
156 
157  return *m_subMinPlain;
158 }
159 // --------------------------------------------------
160 template <class V, class M>
161 const V&
162 BaseVectorSequence<V,M>::unifiedMinPlain() const
163 {
164  if (m_unifiedMinPlain == NULL) {
165  m_unifiedMinPlain = m_vectorSpace.newVector();
166  if (m_unifiedMaxPlain == NULL) m_unifiedMaxPlain = m_vectorSpace.newVector();
167  unifiedMinMaxExtra(0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
168  }
169 
170  return *m_unifiedMinPlain;
171 }
172 // --------------------------------------------------
173 template <class V, class M>
174 const V&
175 BaseVectorSequence<V,M>::subMaxPlain() const
176 {
177  if (m_subMaxPlain == NULL) {
178  if (m_subMinPlain == NULL) m_subMinPlain = m_vectorSpace.newVector();
179  m_subMaxPlain = m_vectorSpace.newVector();
180  subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
181  }
182 
183  return *m_subMaxPlain;
184 }
185 // --------------------------------------------------
186 template <class V, class M>
187 const V&
188 BaseVectorSequence<V,M>::unifiedMaxPlain() const
189 {
190  if (m_unifiedMaxPlain == NULL) {
191  if (m_unifiedMinPlain == NULL) m_unifiedMinPlain = m_vectorSpace.newVector();
192  m_unifiedMaxPlain = m_vectorSpace.newVector();
193  unifiedMinMaxExtra(0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
194  }
195 
196  return *m_unifiedMaxPlain;
197 }
198 // --------------------------------------------------
199 template <class V, class M>
200 const V&
201 BaseVectorSequence<V,M>::subMeanPlain() const
202 {
203  if (m_subMeanPlain == NULL) {
204  m_subMeanPlain = m_vectorSpace.newVector();
205  subMeanExtra(0,subSequenceSize(),*m_subMeanPlain);
206  }
207 
208  return *m_subMeanPlain;
209 }
210 // --------------------------------------------------
211 template <class V, class M>
212 const V&
213 BaseVectorSequence<V,M>::unifiedMeanPlain() const
214 {
215  if (m_unifiedMeanPlain == NULL) {
216  m_unifiedMeanPlain = m_vectorSpace.newVector();
217  unifiedMeanExtra(0,subSequenceSize(),*m_unifiedMeanPlain);
218  }
219 
220  return *m_unifiedMeanPlain;
221 }
222 // --------------------------------------------------
223 template <class V, class M>
224 const V&
225 BaseVectorSequence<V,M>::subMedianPlain() const
226 {
227  if (m_subMedianPlain == NULL) {
228  m_subMedianPlain = m_vectorSpace.newVector();
229  subMedianExtra(0, subSequenceSize(), *m_subMedianPlain);
230  }
231 
232  return *m_subMedianPlain;
233 }
234 // --------------------------------------------------
235 template <class V, class M>
236 const V&
237 BaseVectorSequence<V,M>::unifiedMedianPlain() const
238 {
239  if (m_unifiedMedianPlain == NULL) {
240  m_unifiedMedianPlain = m_vectorSpace.newVector();
241  unifiedMedianExtra(0, subSequenceSize(), *m_unifiedMedianPlain);
242  }
243 
244  return *m_unifiedMedianPlain;
245 }
246 // --------------------------------------------------
247 template <class V, class M>
248 const V&
249 BaseVectorSequence<V,M>::subSampleVariancePlain() const
250 {
251  if (m_subSampleVariancePlain == NULL) {
252  m_subSampleVariancePlain = m_vectorSpace.newVector();
253  subSampleVarianceExtra(0,subSequenceSize(),subMeanPlain(),*m_subSampleVariancePlain);
254  }
255 
256  return *m_subSampleVariancePlain;
257 }
258 // --------------------------------------------------
259 template <class V, class M>
260 const V&
261 BaseVectorSequence<V,M>::unifiedSampleVariancePlain() const
262 {
263  if (m_unifiedSampleVariancePlain == NULL) {
264  m_unifiedSampleVariancePlain = m_vectorSpace.newVector();
265  unifiedSampleVarianceExtra(0,subSequenceSize(),unifiedMeanPlain(),*m_unifiedSampleVariancePlain);
266  }
267 
268  return *m_unifiedSampleVariancePlain;
269 }
270 // --------------------------------------------------
271 template <class V, class M>
272 const BoxSubset<V,M>&
273 BaseVectorSequence<V,M>::subBoxPlain() const
274 {
275  if (m_subBoxPlain == NULL) {
276  m_subBoxPlain = new BoxSubset<V,M>(m_name.c_str(),
277  m_vectorSpace,
278  this->subMinPlain(),
279  this->subMaxPlain());
280  }
281 
282  return *m_subBoxPlain;
283 }
284 // --------------------------------------------------
285 template <class V, class M>
286 const BoxSubset<V,M>&
287 BaseVectorSequence<V,M>::unifiedBoxPlain() const
288 {
289  if (m_unifiedBoxPlain == NULL) {
290  m_unifiedBoxPlain = new BoxSubset<V,M>(m_name.c_str(),
291  m_vectorSpace,
292  this->unifiedMinPlain(),
293  this->unifiedMaxPlain());
294  }
295 
296  return *m_unifiedBoxPlain;
297 }
298 // --------------------------------------------------
299 template <class V, class M>
300 void
301 BaseVectorSequence<V,M>::deleteStoredVectors()
302 {
303  if (m_subMinPlain) {
304  delete m_subMinPlain;
305  m_subMinPlain = NULL;
306  }
307  if (m_unifiedMinPlain) {
308  delete m_unifiedMinPlain;
309  m_unifiedMinPlain = NULL;
310  }
311  if (m_subMaxPlain) {
312  delete m_subMaxPlain;
313  m_subMaxPlain = NULL;
314  }
315  if (m_unifiedMaxPlain) {
316  delete m_unifiedMaxPlain;
317  m_unifiedMaxPlain = NULL;
318  }
319  if (m_subMeanPlain) {
320  delete m_subMeanPlain;
321  m_subMeanPlain = NULL;
322  }
323  if (m_unifiedMeanPlain) {
324  delete m_unifiedMeanPlain;
325  m_unifiedMeanPlain = NULL;
326  }
327  if (m_subMedianPlain) {
328  delete m_subMedianPlain;
329  m_subMedianPlain = NULL;
330  }
331  if (m_unifiedMedianPlain) {
332  delete m_unifiedMedianPlain;
333  m_unifiedMedianPlain = NULL;
334  }
335  if (m_subSampleVariancePlain) {
336  delete m_subSampleVariancePlain;
337  m_subSampleVariancePlain = NULL;
338  }
339  if (m_unifiedSampleVariancePlain) {
340  delete m_unifiedSampleVariancePlain;
341  m_unifiedSampleVariancePlain = NULL;
342  }
343  if (m_subBoxPlain) {
344  delete m_subBoxPlain;
345  m_subBoxPlain = NULL;
346  }
347  if (m_unifiedBoxPlain) {
348  delete m_unifiedBoxPlain;
349  m_unifiedBoxPlain = NULL;
350  }
351 
352  return;
353 }
354 // --------------------------------------------------
355 template <class V, class M>
356 void
357 BaseVectorSequence<V,M>::append(
358  const BaseVectorSequence<V,M>& src,
359  unsigned int initialPos,
360  unsigned int numPos)
361 {
362  UQ_FATAL_TEST_MACRO((src.subSequenceSize() < (initialPos+1)),
363  m_env.worldRank(),
364  "BaseVectorSequence<V,M>::append()",
365  "initialPos is too big");
366 
367  UQ_FATAL_TEST_MACRO((src.subSequenceSize() < (initialPos+numPos)),
368  m_env.worldRank(),
369  "BaseVectorSequence<V,M>::append()",
370  "numPos is too big");
371 
372  this->deleteStoredVectors();
373  unsigned int currentSize = this->subSequenceSize();
374  this->resizeSequence(currentSize+numPos);
375  V tmpVec(src.vectorSpace().zeroVector());
376  for (unsigned int i = 0; i < numPos; ++i) {
377  src.getPositionValues(initialPos+i,tmpVec);
378  this->setPositionValues(currentSize+i,tmpVec);
379  }
380 
381  return;
382 }
383 // --------------------------------------------------
384 template <class V, class M>
385 double
386 BaseVectorSequence<V,M>::subPositionsOfMaximum(
387  const ScalarSequence<double>& subCorrespondingScalarValues,
388  BaseVectorSequence<V,M>& subPositionsOfMaximum)
389 {
390  if (m_env.subDisplayFile()) {
391  *m_env.subDisplayFile() << "Entering BaseVectorSequence<V,M>::subPositionsOfMaximum()"
392  << ": subCorrespondingScalarValues,subSequenceSize() = " << subCorrespondingScalarValues.subSequenceSize()
393  << ", this->subSequenceSize = " << this->subSequenceSize()
394  << std::endl;
395  }
396 
397  UQ_FATAL_TEST_MACRO(subCorrespondingScalarValues.subSequenceSize() != this->subSequenceSize(),
398  m_env.worldRank(),
399  "BaseVectorSequence<V,M>::subPositionsOfMaximum()",
400  "invalid input");
401 
402  double subMaxValue = subCorrespondingScalarValues.subMaxPlain();
403  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
404 
405  unsigned int subNumPos = 0;
406  for (unsigned int i = 0; i < iMax; ++i) {
407  if (subCorrespondingScalarValues[i] == subMaxValue) {
408  subNumPos++;
409  }
410  }
411 
412  V tmpVec(this->vectorSpace().zeroVector());
413  subPositionsOfMaximum.resizeSequence(subNumPos);
414  unsigned int j = 0;
415  for (unsigned int i = 0; i < iMax; ++i) {
416  if (subCorrespondingScalarValues[i] == subMaxValue) {
417  this->getPositionValues (i,tmpVec);
418  subPositionsOfMaximum.setPositionValues(j,tmpVec);
419  j++;
420  }
421  }
422 
423  if (m_env.subDisplayFile()) {
424  *m_env.subDisplayFile() << "Leaving BaseVectorSequence<V,M>::subPositionsOfMaximum()"
425  << std::endl;
426  }
427 
428  return subMaxValue;
429 }
430 // --------------------------------------------------
431 template <class V, class M>
432 double
433 BaseVectorSequence<V,M>::unifiedPositionsOfMaximum( // rr0
434  const ScalarSequence<double>& subCorrespondingScalarValues,
435  BaseVectorSequence<V,M>& unifiedPositionsOfMaximum)
436 {
437  if (m_env.subDisplayFile()) {
438  *m_env.subDisplayFile() << "Entering BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()"
439  << ": subCorrespondingScalarValues,subSequenceSize() = " << subCorrespondingScalarValues.subSequenceSize()
440  << ", this->subSequenceSize = " << this->subSequenceSize()
441  << std::endl;
442  }
443 
444  UQ_FATAL_TEST_MACRO(subCorrespondingScalarValues.subSequenceSize() != this->subSequenceSize(),
445  m_env.worldRank(),
446  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
447  "invalid input");
448 
449  // Compute the max on each process
450  double subMaxValue = subCorrespondingScalarValues.subMaxPlain();
451 
452  //******************************************************************
453  // Get overall max
454  //******************************************************************
455  double unifiedMaxValue;
456  std::vector<double> sendbufPos(1,0.);
457  for (unsigned int i = 0; i < sendbufPos.size(); ++i) {
458  sendbufPos[i] = subMaxValue;
459  }
460  m_env.inter0Comm().Allreduce((void *) &sendbufPos[0], (void *) &unifiedMaxValue, (int) sendbufPos.size(), RawValue_MPI_DOUBLE, RawValue_MPI_MAX,
461  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
462  "failed MPI.Allreduce() for max");
463 
464  // Find number of elements (subNumPos) on each process that attain the
465  // global max. This could be zero.
466  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
467  int subNumPos = 0; // Yes, 'int', due to MPI to be used soon
468  for (unsigned int i = 0; i < iMax; ++i) {
469  if (subCorrespondingScalarValues[i] == unifiedMaxValue) {
470  subNumPos++;
471  }
472  }
473 
474  // Fill up unifiedPositionsOfMaximum with the states that attain maxima
475  // (if they exist)
476  V tmpVec(this->vectorSpace().zeroVector());
477  unifiedPositionsOfMaximum.resizeSequence(subNumPos); // subNumPos could be 0
478  unsigned int j = 0;
479  for (unsigned int i = 0; i < iMax; ++i) {
480  // This 'if' statement is false if subNumPos == 0
481  if (subCorrespondingScalarValues[i] == unifiedMaxValue) {
482  this->getPositionValues (i,tmpVec);
483  unifiedPositionsOfMaximum.setPositionValues(j,tmpVec);
484  j++;
485  }
486  }
487 
488  // Compute total (over all processes) number of elements that attain maxima
489  std::vector<int> auxBuf(1,0);
490  int unifiedNumPos = 0; // Yes, 'int', due to MPI to be used soon
491  auxBuf[0] = subNumPos;
492  m_env.inter0Comm().Allreduce((void *) &auxBuf[0], (void *) &unifiedNumPos, (int) auxBuf.size(), RawValue_MPI_INT, RawValue_MPI_SUM,
493  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
494  "failed MPI.Allreduce() for sum");
495 
496  // Resize positions to hold all elements that attain maxima
497  unifiedPositionsOfMaximum.resizeSequence(unifiedNumPos);
498 
499  //******************************************************************
500  // Use MPI_Gatherv for number of positions
501  //******************************************************************
502  // Gather up *number* of maxima on each chain and store in recvcntsPos
503  unsigned int Np = (unsigned int) m_env.inter0Comm().NumProc();
504 
505  std::vector<int> recvcntsPos(Np,0); // '0' is NOT the correct value for recvcntsPos[0]
506  m_env.inter0Comm().Gather((void *) &subNumPos, 1, RawValue_MPI_INT, (void *) &recvcntsPos[0], (int) 1, RawValue_MPI_INT, 0,
507  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
508  "failed MPI.Gatherv()");
509  if (m_env.inter0Rank() == 0) {
510  UQ_FATAL_TEST_MACRO(recvcntsPos[0] != (int) subNumPos,
511  m_env.worldRank(),
512  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
513  "failed MPI.Gather() result at proc 0 (recvcntsPos[0])");
514  }
515 
516  // Construct offset indices based on number of maxima
517  std::vector<int> displsPos(Np,0);
518  for (unsigned int nodeId = 1; nodeId < Np; ++nodeId) { // Yes, from '1' on
519  displsPos[nodeId] = displsPos[nodeId-1] + recvcntsPos[nodeId-1];
520  }
521  if (m_env.inter0Rank() == 0) {
522  UQ_FATAL_TEST_MACRO(unifiedNumPos != (displsPos[Np - 1] + recvcntsPos[Np - 1]),
523  m_env.worldRank(),
524  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
525  "failed MPI.Gather() result at proc 0 (unifiedNumPos)");
526  }
527 
528  //******************************************************************
529  // Use MPI_Gatherv for number of doubles
530  //******************************************************************
531  // Gather up number of states that attain maxima multiplied by the size of
532  // the state (on each process). Store this in recvcntsDbs.
533  // So recvcntsDbs[i] is the number of states that attain maxima on chain i
534  // multiplied by the size of the state (a state could be a vector of length
535  // > 1).
536  unsigned int dimSize = m_vectorSpace.dimLocal();
537  int subNumDbs = subNumPos * dimSize; // Yes, 'int', due to MPI to be used soon
538  std::vector<int> recvcntsDbs(Np,0); // '0' is NOT the correct value for recvcntsDbs[0]
539  m_env.inter0Comm().Gather((void *) &subNumDbs, 1, RawValue_MPI_INT, (void *) &recvcntsDbs[0], (int) 1, RawValue_MPI_INT, 0,
540  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
541  "failed MPI.Gatherv()");
542  if (m_env.inter0Rank() == 0) {
543  UQ_FATAL_TEST_MACRO(recvcntsDbs[0] != (int) subNumDbs,
544  m_env.worldRank(),
545  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
546  "failed MPI.Gather() result at proc 0 (recvcntsDbs[0])");
547  }
548 
549  // Now gather up the offsets based on the number of states (mulitplied by the
550  // size of the state)
551  std::vector<int> displsDbs(Np,0);
552  for (unsigned int nodeId = 1; nodeId < Np; ++nodeId) { // Yes, from '1' on
553  displsDbs[nodeId] = displsDbs[nodeId-1] + recvcntsDbs[nodeId-1];
554  }
555  if (m_env.inter0Rank() == 0) {
556  UQ_FATAL_TEST_MACRO(((int) (unifiedNumPos*dimSize)) != (displsDbs[Np - 1] + recvcntsDbs[Np - 1]),
557  m_env.worldRank(),
558  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
559  "failed MPI.Gather() result at proc 0 (unifiedNumPos*dimSize)");
560  }
561 
562  //******************************************************************
563  // Prepare counters and buffers for gatherv of maximum positions
564  //******************************************************************
565  // Take all the states on all chains that attain maxima, and put them in a
566  // send buffer ready for MPI
567  std::vector<double> sendbufDbs(subNumDbs,0.);
568  for (unsigned int i = 0; i < (unsigned int) subNumPos; ++i) {
569  unifiedPositionsOfMaximum.getPositionValues(i,tmpVec);
570  for (unsigned int j = 0; j < dimSize; ++j) {
571  sendbufDbs[i*dimSize + j] = tmpVec[j];
572  }
573  }
574 
575  std::vector<double> recvbufDbs(unifiedNumPos * dimSize);
576 
577  // Gather up all states that attain maxima and store then in recvbufDbs
578  m_env.inter0Comm().Gatherv((void *) &sendbufDbs[0], (int) subNumDbs, RawValue_MPI_DOUBLE, (void *) &recvbufDbs[0], (int *) &recvcntsDbs[0], (int *) &displsDbs[0], RawValue_MPI_DOUBLE, 0,
579  "BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()",
580  "failed MPI.Gatherv()");
581 
582  //******************************************************************
583  // Transfer data from 'recvbuf' to 'unifiedPositionsOfMaximum'
584  //******************************************************************
585  // Copy over all the gathered states to the unifiedPositionsOfMaximum
586  // variable on process 0. Process 0 now contains everything.
587  if (m_env.inter0Rank() == (int) 0) {
588  for (unsigned int i = 0; i < (unsigned int) unifiedNumPos; ++i) {
589  for (unsigned int j = 0; j < dimSize; ++j) {
590  tmpVec[j] = recvbufDbs[i*dimSize + j];
591  }
592  unifiedPositionsOfMaximum.setPositionValues(i,tmpVec);
593  }
594  }
595  else {
596  // Process zero has all the states that attain maxima, so let's nuke the
597  // others rather than letting them contain NULL pointers
598  unifiedPositionsOfMaximum.resizeSequence(0);
599  }
600 
601  if (m_env.subDisplayFile()) {
602  *m_env.subDisplayFile() << "Leaving BaseVectorSequence<V,M>::unifiedPositionsOfMaximum()"
603  << std::endl;
604  }
605 
606  return unifiedMaxValue;
607 }
608 // --------------------------------------------------
609 template <class V, class M>
610 void
611 BaseVectorSequence<V,M>::setGaussian(const V& meanVec, const V& stdDevVec)
612 {
613  V gaussianVector(m_vectorSpace.zeroVector());
614  for (unsigned int j = 0; j < this->subSequenceSize(); ++j) {
615  gaussianVector.cwSetGaussian(meanVec,stdDevVec);
616  this->setPositionValues(j,gaussianVector);
617  }
618 
619  this->deleteStoredVectors();
620 
621  return;
622 }
623 // --------------------------------------------------
624 template <class V, class M>
625 void
626 BaseVectorSequence<V,M>::setUniform(const V& aVec, const V& bVec)
627 {
628  V uniformVector(m_vectorSpace.zeroVector());
629  for (unsigned int j = 0; j < this->subSequenceSize(); ++j) {
630  uniformVector.cwSetUniform(aVec,bVec);
631  this->setPositionValues(j,uniformVector);
632  }
633 
634  this->deleteStoredVectors();
635 
636  return;
637 }
638 // --------------------------------------------------
639 template<class V, class M>
640 void
641 BaseVectorSequence<V,M>::computeFilterParams(
642  std::ofstream* passedOfs,
643  unsigned int& initialPos,
644  unsigned int& spacing)
645 {
646  if (m_env.subDisplayFile()) {
647  *m_env.subDisplayFile() << "\n"
648  << "\n-----------------------------------------------------"
649  << "\n Computing filter parameters for chain '" << m_name << "' ..."
650  << "\n-----------------------------------------------------"
651  << "\n"
652  << std::endl;
653  }
654 
655  bool okSituation = ((passedOfs == NULL ) ||
656  ((passedOfs != NULL) && (m_env.subRank() >= 0)));
657  UQ_FATAL_TEST_MACRO(!okSituation,
658  m_env.worldRank(),
659  "BaseVectorSequence<V,M>::computeFilterParams()",
660  "unexpected combination of file pointer and subRank");
661 
662  //initialPos = 0;
663  spacing = 1;
664 
665  if (m_env.subDisplayFile()) {
666  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
667  << "\n Finished computing filter parameters for chain '" << m_name << "'"
668  << ": initialPos = " << initialPos
669  << ", spacing = " << spacing
670  << "\n-----------------------------------------------------"
671  << "\n"
672  << std::endl;
673  }
674 
675  return;
676 }
677 
678 // Protected methods---------------------------------
679 template <class V, class M>
680 void
681 BaseVectorSequence<V,M>::copy(const BaseVectorSequence<V,M>& src)
682 {
683  // FIX ME: should check environments as well ???
684 
685  UQ_FATAL_TEST_MACRO(m_vectorSpace.dimLocal() != src.m_vectorSpace.dimLocal(),
686  m_env.worldRank(),
687  "SequenceOfVectors<V,M>::copy()",
688  "incompatible vector space dimensions");
689 
690  m_name = src.m_name;
691  this->deleteStoredVectors();
692 
693  return;
694 }
695 
696 // --------------------------------------------------
697 // Methods conditionally available ------------------
698 // --------------------------------------------------
699 // --------------------------------------------------
700 
701 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
702 template<class V, class M>
703 void
704 BaseVectorSequence<V,M>::computeStatistics(
705  const SequenceStatisticalOptions& statisticalOptions,
706  std::ofstream* passedOfs)
707 {
708  if (m_env.subDisplayFile()) {
709  *m_env.subDisplayFile() << "\n"
710  << "\n-----------------------------------------------------"
711  << "\n Computing statistics for chain '" << m_name << "' ..."
712  << "\n-----------------------------------------------------"
713  << "\n"
714  << std::endl;
715  }
716 
717  bool okSituation = ((passedOfs == NULL ) ||
718  ((passedOfs != NULL) && (m_env.subRank() >= 0)));
719  UQ_FATAL_TEST_MACRO(!okSituation,
720  m_env.worldRank(),
721  "BaseVectorSequence<V,M>::computeStatistics()",
722  "unexpected combination of file pointer and subRank");
723 
724  int iRC = UQ_OK_RC;
725  struct timeval timevalTmp;
726  iRC = gettimeofday(&timevalTmp, NULL);
727  double tmpRunTime = 0.;
728 
729  // Set initial positions for the computation of chain statistics
730  std::vector<unsigned int> initialPosForStatistics(statisticalOptions.initialDiscardedPortions().size(),0);
731  for (unsigned int i = 0; i < initialPosForStatistics.size(); ++i) {
732  initialPosForStatistics[i] = (unsigned int) (statisticalOptions.initialDiscardedPortions()[i] * (double) (this->subSequenceSize()-1));
733  if (m_env.subDisplayFile()) {
734  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeStatistics()"
735  << ": statisticalOptions.initialDiscardedPortions()[" << i << "] = " << statisticalOptions.initialDiscardedPortions()[i]
736  << ", initialPosForStatistics[" << i << "] = " << initialPosForStatistics[i]
737  << std::endl;
738  }
739  }
740  if (m_env.subDisplayFile()) {
741  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeStatistics(): initial positions for statistics =";
742  for (unsigned int i = 0; i < initialPosForStatistics.size(); ++i) {
743  *m_env.subDisplayFile() << " " << initialPosForStatistics[i];
744  }
745  *m_env.subDisplayFile() << std::endl;
746  }
747 
748  //****************************************************
749  // Compute mean, median, sample std, population std
750  //****************************************************
751  this->computeMeanVars(statisticalOptions,
752  passedOfs,
753  NULL,
754  NULL,
755  NULL,
756  NULL);
757 
758 #ifdef UQ_CODE_HAS_MONITORS
759  if (statisticalOptions.meanMonitorPeriod() != 0) {
760  this->computeMeanEvolution(statisticalOptions,
761  passedOfs);
762  }
763 #endif
764 
765  //****************************************************
766  // Compute variance of sample mean through the 'batch means method' (BMM)
767  //****************************************************
768 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
769  if ((statisticalOptions.bmmRun() ) &&
770  (initialPosForStatistics.size() > 0) &&
771  (statisticalOptions.bmmLengths().size() > 0)) {
772  this->computeBMM(statisticalOptions,
773  initialPosForStatistics,
774  passedOfs);
775  }
776 #endif
777  //****************************************************
778  // Compute FFT of chain, for one parameter only
779  //****************************************************
780 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
781  if ((statisticalOptions.fftCompute() ) &&
782  (initialPosForStatistics.size() > 0)) {
783  this->computeFFT(statisticalOptions,
784  initialPosForStatistics,
785  passedOfs);
786  }
787 #endif
788  //****************************************************
789  // Compute power spectral density (PSD) of chain, for one parameter only
790  //****************************************************
791 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
792  if ((statisticalOptions.psdCompute() ) &&
793  (initialPosForStatistics.size() > 0)) {
794  this->computePSD(statisticalOptions,
795  initialPosForStatistics,
796  passedOfs);
797  }
798 #endif
799  //****************************************************
800  // Compute power spectral density (PSD) of chain at zero frequency
801  //****************************************************
802 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
803  if ((statisticalOptions.psdAtZeroCompute() ) &&
804  (initialPosForStatistics.size() > 0) &&
805  (statisticalOptions.psdAtZeroNumBlocks().size() > 0)) {
806  this->computePSDAtZero(statisticalOptions,
807  initialPosForStatistics,
808  passedOfs);
809  }
810 #endif
811  //****************************************************
812  // Compute Geweke
813  //****************************************************
814 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
815  if ((statisticalOptions.gewekeCompute()) &&
816  (initialPosForStatistics.size() > 0)) {
817  this->computeGeweke(statisticalOptions,
818  initialPosForStatistics,
819  passedOfs);
820  }
821 #endif
822  //****************************************************
823  // Compute mean statistical accuracy
824  //****************************************************
825 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
826  if ((statisticalOptions.meanStaccCompute()) &&
827  (initialPosForStatistics.size() > 0 )) {
828  this->computeMeanStacc(statisticalOptions,
829  initialPosForStatistics,
830  passedOfs);
831  }
832 #endif
833  // Set lags for the computation of chain autocorrelations
834  std::vector<unsigned int> lagsForCorrs(statisticalOptions.autoCorrNumLags(),1);
835  for (unsigned int i = 1; i < lagsForCorrs.size(); ++i) {
836  lagsForCorrs[i] = statisticalOptions.autoCorrSecondLag() + (i-1)*statisticalOptions.autoCorrLagSpacing();
837  }
838 
839  //****************************************************
840  // Compute autocorrelation coefficients via definition
841  //****************************************************
842  if ((statisticalOptions.autoCorrComputeViaDef()) &&
843  (initialPosForStatistics.size() > 0 ) &&
844  (lagsForCorrs.size() > 0 )) {
845  this->computeAutoCorrViaDef(statisticalOptions,
846  initialPosForStatistics,
847  lagsForCorrs,
848  passedOfs);
849  }
850 
851  //****************************************************
852  // Compute autocorrelation coefficients via FFT
853  //****************************************************
854  if ((statisticalOptions.autoCorrComputeViaFft()) &&
855  (initialPosForStatistics.size() > 0 ) &&
856  (lagsForCorrs.size() > 0 )) {
857  this->computeAutoCorrViaFFT(statisticalOptions,
858  initialPosForStatistics,
859  lagsForCorrs,
860  passedOfs);
861  }
862 
863  //****************************************************
864  // Compute histogram and/or cdf stacc and/or Kde
865  //****************************************************
866 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
867  if ((statisticalOptions.histCompute ()) ||
868  (statisticalOptions.cdfStaccCompute()) ||
869  (statisticalOptions.kdeCompute ())) {
870 #else
871  if (statisticalOptions.kdeCompute()) //kemelli deleted single '{' on 5/28/13 - compare to queso-0.45.2
872 #endif
873  this->computeHistCdfstaccKde(statisticalOptions,
874  passedOfs);
875  }
876 
877  //****************************************************
878  // Compute covariance and correlation matrices
879  //****************************************************
880  if ((statisticalOptions.covMatrixCompute ()) ||
881  (statisticalOptions.corrMatrixCompute())) {
882  this->computeCovCorrMatrices(statisticalOptions,
883  passedOfs);
884  }
885 
886  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
887  if (m_env.subDisplayFile()) {
888  *m_env.subDisplayFile() << "All statistics of chain '" << m_name << "'"
889  << " took " << tmpRunTime
890  << " seconds"
891  << std::endl;
892  }
893 
894  if (m_env.subDisplayFile()) {
895  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
896  << "\n Finished computing statistics for chain '" << m_name << "'"
897  << "\n-----------------------------------------------------"
898  << "\n"
899  << std::endl;
900  }
901 
902  return;
903 }
904 //----------------------------------------------
905 template<class V, class M>
906 void
907 BaseVectorSequence<V,M>::computeMeanVars(
908  const SequenceStatisticalOptions& statisticalOptions,
909  std::ofstream* passedOfs,
910  V* subMeanPtr,
911  V* subMedianPtr,
912  V* subSampleVarPtr,
913  V* subPopulVarPtr)
914 {
915  int iRC = UQ_OK_RC;
916  struct timeval timevalTmp;
917  iRC = gettimeofday(&timevalTmp, NULL);
918  double tmpRunTime = 0.;
919 
920  if (m_env.subDisplayFile()) {
921  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
922  << "\nComputing mean, sample variance and population variance"
923  << std::endl;
924  }
925 
926  V subChainMean(m_vectorSpace.zeroVector());
927  this->subMeanExtra(0,
928  this->subSequenceSize(),
929  subChainMean);
930 
931  V subChainMedian(m_vectorSpace.zeroVector());
932  this->subMedianExtra(0,
933  this->subSequenceSize(),
934  subChainMedian);
935 
936  V subChainSampleVariance(m_vectorSpace.zeroVector());
937  this->subSampleVarianceExtra(0,
938  this->subSequenceSize(),
939  subChainMean,
940  subChainSampleVariance);
941 
942  if ((m_env.displayVerbosity() >= 5) && (m_env.subDisplayFile())) {
943  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeMeanVars()"
944  << ": subChainMean.sizeLocal() = " << subChainMean.sizeLocal()
945  << ", subChainMean = " << subChainMean
946  << ", subChainMedian = " << subChainMedian
947  << ", subChainSampleVariance.sizeLocal() = " << subChainSampleVariance.sizeLocal()
948  << ", subChainSampleVariance = " << subChainSampleVariance
949  << std::endl;
950  }
951 
952  V estimatedVarianceOfSampleMean(subChainSampleVariance);
953  estimatedVarianceOfSampleMean /= (double) this->subSequenceSize();
954  bool savedVectorPrintState = estimatedVarianceOfSampleMean.getPrintHorizontally();
955  estimatedVarianceOfSampleMean.setPrintHorizontally(false);
956  if (m_env.subDisplayFile()) {
957  *m_env.subDisplayFile() << "\nEstimated variance of sample mean for the whole chain '" << m_name << "'"
958  << ", under independence assumption:"
959  << "\n"
960  << estimatedVarianceOfSampleMean
961  << std::endl;
962  }
963  estimatedVarianceOfSampleMean.setPrintHorizontally(savedVectorPrintState);
964 
965  V estimatedStdOfSampleMean(estimatedVarianceOfSampleMean);
966  estimatedStdOfSampleMean.cwSqrt();
967  savedVectorPrintState = estimatedStdOfSampleMean.getPrintHorizontally();
968  estimatedStdOfSampleMean.setPrintHorizontally(false);
969  if (m_env.subDisplayFile()) {
970  *m_env.subDisplayFile() << "\nEstimated standard deviation of sample mean for the whole chain '" << m_name << "'"
971  << ", under independence assumption:"
972  << "\n"
973  << estimatedStdOfSampleMean
974  << std::endl;
975  }
976  estimatedStdOfSampleMean.setPrintHorizontally(savedVectorPrintState);
977 
978  V subChainPopulationVariance(m_vectorSpace.zeroVector());
979  this->subPopulationVariance(0,
980  this->subSequenceSize(),
981  subChainMean,
982  subChainPopulationVariance);
983 
984  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
985  if (m_env.subDisplayFile()) {
986  *m_env.subDisplayFile() << "Sub Mean, median, and variances took " << tmpRunTime
987  << " seconds"
988  << std::endl;
989  }
990 
991  if (m_env.subDisplayFile()) {
992  *m_env.subDisplayFile() << "\nSub mean, median, sample std, population std"
993  << std::endl;
994  char line[512];
995  sprintf(line,"%s%4s%s%6s%s%9s%s%9s%s",
996  "Parameter",
997  " ",
998  "Mean",
999  " ",
1000  "Median",
1001  " ",
1002  "SampleStd",
1003  " ",
1004  "Popul.Std");
1005  *m_env.subDisplayFile() << line;
1006 
1007  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1008  sprintf(line,"\n%8.8s%2s%11.4e%2s%11.4e%2s%11.4e%2s%11.4e",
1009  m_vectorSpace.localComponentName(i).c_str(), /*.*/
1010  " ",
1011  subChainMean[i],
1012  " ",
1013  subChainMedian[i],
1014  " ",
1015  std::sqrt(subChainSampleVariance[i]),
1016  " ",
1017  std::sqrt(subChainPopulationVariance[i]));
1018  *m_env.subDisplayFile() << line;
1019  }
1020  *m_env.subDisplayFile() << std::endl;
1021  }
1022 
1023  if (subMeanPtr ) *subMeanPtr = subChainMean;
1024  if (subMedianPtr ) *subMedianPtr = subChainMedian;
1025  if (subSampleVarPtr) *subSampleVarPtr = subChainSampleVariance;
1026  if (subPopulVarPtr ) *subPopulVarPtr = subChainPopulationVariance;
1027 
1028  if (m_env.numSubEnvironments() > 1) {
1029  // Write unified min-max
1030  if (m_vectorSpace.numOfProcsForStorage() == 1) {
1031  V unifiedChainMean(m_vectorSpace.zeroVector());
1032  this->unifiedMeanExtra(0,
1033  this->subSequenceSize(),
1034  unifiedChainMean);
1035 
1036  V unifiedChainMedian(m_vectorSpace.zeroVector());
1037  this->unifiedMedianExtra(0,
1038  this->subSequenceSize(),
1039  unifiedChainMedian);
1040 
1041  V unifiedChainSampleVariance(m_vectorSpace.zeroVector());
1042  this->unifiedSampleVarianceExtra(0,
1043  this->subSequenceSize(),
1044  unifiedChainMean,
1045  unifiedChainSampleVariance);
1046 
1047  V unifiedChainPopulationVariance(m_vectorSpace.zeroVector());
1048  this->unifiedPopulationVariance(0,
1049  this->subSequenceSize(),
1050  unifiedChainMean,
1051  unifiedChainPopulationVariance);
1052 
1053  if (m_env.inter0Rank() == 0) {
1054  if (m_env.subDisplayFile()) {
1055  *m_env.subDisplayFile() << "\nUnif mean, median, sample std, population std"
1056  << std::endl;
1057  char line[512];
1058  sprintf(line,"%s%4s%s%6s%s%9s%s%9s%s",
1059  "Parameter",
1060  " ",
1061  "Mean",
1062  " ",
1063  "Median",
1064  " ",
1065  "SampleStd",
1066  " ",
1067  "Popul.Std");
1068  *m_env.subDisplayFile() << line;
1069 
1070  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1071  sprintf(line,"\n%8.8s%2s%11.4e%2s%11.4e%2s%11.4e%2s%11.4e",
1072  m_vectorSpace.localComponentName(i).c_str(), /*.*/
1073  " ",
1074  unifiedChainMean[i],
1075  " ",
1076  unifiedChainMedian[i],
1077  " ",
1078  std::sqrt(unifiedChainSampleVariance[i]),
1079  " ",
1080  std::sqrt(unifiedChainPopulationVariance[i]));
1081  *m_env.subDisplayFile() << line;
1082  }
1083  *m_env.subDisplayFile() << std::endl;
1084  } // if subDisplayFile
1085  }
1086  }
1087  else {
1088  UQ_FATAL_TEST_MACRO(true,
1089  m_env.worldRank(),
1090  "BaseVectorSequence<V,M>::computeMeanVars()",
1091  "unified min-max writing, parallel vectors not supported yet");
1092  }
1093  } // if numSubEnvs > 1
1094 
1095  if (m_env.subDisplayFile()) {
1096  *m_env.subDisplayFile() << "\nEnded computing mean, sample variance and population variance"
1097  << "\n-----------------------------------------------------"
1098  << std::endl;
1099  }
1100 
1101  return;
1102 }
1103 //--------------------------------------------------------
1104 template<class V, class M>
1105 void
1106 BaseVectorSequence<V,M>::computeAutoCorrViaDef(
1107  const SequenceStatisticalOptions& statisticalOptions,
1108  const std::vector<unsigned int>& initialPosForStatistics,
1109  const std::vector<unsigned int>& lagsForCorrs,
1110  std::ofstream* passedOfs)
1111 {
1112  int iRC = UQ_OK_RC;
1113  struct timeval timevalTmp;
1114  iRC = gettimeofday(&timevalTmp, NULL);
1115  double tmpRunTime = 0.;
1116 
1117  if (m_env.subDisplayFile()) {
1118  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
1119  << "\nComputing autocorrelation coefficients (via def)"
1120  << std::endl;
1121  }
1122 
1123  if (statisticalOptions.autoCorrDisplay() && (m_env.subDisplayFile())) {
1124  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeAutoCorrViaDef(): lags for autocorrelation (via def) = ";
1125  for (unsigned int i = 0; i < lagsForCorrs.size(); ++i) {
1126  *m_env.subDisplayFile() << " " << lagsForCorrs[i];
1127  }
1128  *m_env.subDisplayFile() << std::endl;
1129  }
1130 
1131  TwoDArray<V> _2dArrayOfAutoCorrs(initialPosForStatistics.size(),lagsForCorrs.size());
1132  for (unsigned int i = 0; i < _2dArrayOfAutoCorrs.numRows(); ++i) {
1133  for (unsigned int j = 0; j < _2dArrayOfAutoCorrs.numCols(); ++j) {
1134  _2dArrayOfAutoCorrs.setLocation(i,j,new V(m_vectorSpace.zeroVector()) /*.*/);
1135  }
1136  }
1137  //V corrVec(m_vectorSpace.zeroVector());
1138  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1139  unsigned int initialPos = initialPosForStatistics[initialPosId];
1140  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1141  unsigned int lag = lagsForCorrs[lagId];
1142  this->autoCorrViaDef(initialPos,
1143  this->subSequenceSize()-initialPos,
1144  lag,
1145  _2dArrayOfAutoCorrs(initialPosId,lagId));
1146  //_2dArrayOfAutoCorrs(initialPosId,lagId) = corrVec;
1147  }
1148  }
1149 
1150  // It is not practical to compute the variance of sample mean by computing the autocorrelations via definition for each lag
1151  // The code computes the variance of sample mean by computing the autocorrelations via fft, below, in another routine
1152 
1153  if ((statisticalOptions.autoCorrDisplay()) && (m_env.subDisplayFile())) {
1154  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1155  *m_env.subDisplayFile() << "\nComputed autocorrelation coefficients (via def), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1156  << " (each column corresponds to a different lag)"
1157  << std::endl;
1158  char line[512];
1159  sprintf(line,"%s",
1160  "Parameter");
1161  *m_env.subDisplayFile() << line;
1162  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1163  sprintf(line,"%10s%3d",
1164  " ",
1165  lagsForCorrs[lagId]);
1166  *m_env.subDisplayFile() << line;
1167  }
1168 
1169  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1170  sprintf(line,"\n%9.9s",
1171  m_vectorSpace.localComponentName(i).c_str() /*.*/);
1172  *m_env.subDisplayFile() << line;
1173  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1174  sprintf(line,"%2s%11.4e",
1175  " ",
1176  _2dArrayOfAutoCorrs(initialPosId,lagId)[i]);
1177  *m_env.subDisplayFile() << line;
1178  }
1179  }
1180  *m_env.subDisplayFile() << std::endl;
1181  }
1182  }
1183 
1184  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
1185  if (m_env.subDisplayFile()) {
1186  *m_env.subDisplayFile() << "Chain autocorrelation (via def) took " << tmpRunTime
1187  << " seconds"
1188  << std::endl;
1189  }
1190 
1191  // Write autocorrelations
1192  if (statisticalOptions.autoCorrWrite() && passedOfs) {
1193  std::ofstream& ofsvar = *passedOfs;
1194  ofsvar << m_name << "_corrViaDefLags_sub" << m_env.subIdString() << " = zeros(" << 1
1195  << "," << lagsForCorrs.size()
1196  << ");"
1197  << std::endl;
1198  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1199  ofsvar << m_name << "_corrViaDefLags_sub" << m_env.subIdString() << "(" << 1
1200  << "," << lagId+1
1201  << ") = " << lagsForCorrs[lagId]
1202  << ";"
1203  << std::endl;
1204  }
1205 
1206  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1207  ofsvar << m_name << "_corrViaDefInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << " = zeros(" << this->vectorSizeLocal() /*.*/
1208  << "," << lagsForCorrs.size()
1209  << ");"
1210  << std::endl;
1211  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1212  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1213  ofsvar << m_name << "_corrViaDefInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << "(" << i+1
1214  << "," << lagId+1
1215  << ") = " << _2dArrayOfAutoCorrs(initialPosId,lagId)[i]
1216  << ";"
1217  << std::endl;
1218  }
1219  }
1220  }
1221  }
1222 
1223  return;
1224 }
1225 // --------------------------------------------------
1226 template<class V, class M>
1227 void
1228 BaseVectorSequence<V,M>::computeAutoCorrViaFFT(
1229  const SequenceStatisticalOptions& statisticalOptions,
1230  const std::vector<unsigned int>& initialPosForStatistics,
1231  const std::vector<unsigned int>& lagsForCorrs,
1232  std::ofstream* passedOfs)
1233 {
1234  int iRC = UQ_OK_RC;
1235  struct timeval timevalTmp;
1236  iRC = gettimeofday(&timevalTmp, NULL);
1237  double tmpRunTime = 0.;
1238 
1239  if (m_env.subDisplayFile()) {
1240  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
1241  << "\nComputing autocorrelation coefficients (via fft)"
1242  << std::endl;
1243  }
1244 
1245  if (statisticalOptions.autoCorrDisplay() && (m_env.subDisplayFile())) {
1246  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeAutoCorrViaFFT(): lags for autocorrelation (via fft) = ";
1247  for (unsigned int i = 0; i < lagsForCorrs.size(); ++i) {
1248  *m_env.subDisplayFile() << " " << lagsForCorrs[i];
1249  }
1250  *m_env.subDisplayFile() << std::endl;
1251  }
1252 
1253  TwoDArray<V> _2dArrayOfAutoCorrs(initialPosForStatistics.size(),lagsForCorrs.size());
1254  for (unsigned int i = 0; i < _2dArrayOfAutoCorrs.numRows(); ++i) {
1255  for (unsigned int j = 0; j < _2dArrayOfAutoCorrs.numCols(); ++j) {
1256  _2dArrayOfAutoCorrs.setLocation(i,j,new V(m_vectorSpace.zeroVector()) /*.*/);
1257  }
1258  }
1259  std::vector<V*> corrVecs(lagsForCorrs.size(),NULL);
1260  std::vector<V*> corrSumVecs(initialPosForStatistics.size(),NULL);
1261  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1262  corrSumVecs[initialPosId] = new V(m_vectorSpace.zeroVector()) /*.*/;
1263  unsigned int initialPos = initialPosForStatistics[initialPosId];
1264  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1265  corrVecs[lagId] = new V(m_vectorSpace.zeroVector()) /*.*/;
1266  }
1267  if (m_env.subDisplayFile()) {
1268  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeAutoCorrViaFFT()"
1269  << ": about to call chain.autoCorrViaFft()"
1270  << " with initialPos = " << initialPos
1271  << ", numPos = " << this->subSequenceSize()-initialPos
1272  << ", lagsForCorrs.size() = " << lagsForCorrs.size()
1273  << ", corrVecs.size() = " << corrVecs.size()
1274  << std::endl;
1275  }
1276  this->autoCorrViaFft(initialPos,
1277  this->subSequenceSize()-initialPos, // Use all possible data positions
1278  lagsForCorrs,
1279  corrVecs);
1280  this->autoCorrViaFft(initialPos,
1281  this->subSequenceSize()-initialPos, // Use all possible data positions
1282  (unsigned int) (1.0 * (double) (this->subSequenceSize()-initialPos)), // CHECK
1283  *corrSumVecs[initialPosId]); // Sum of all possibly computable autocorrelations, not only the asked ones in lagsForCorrs
1284  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1285  _2dArrayOfAutoCorrs(initialPosId,lagId) = *(corrVecs[lagId]);
1286  }
1287  }
1288  for (unsigned int j = 0; j < corrVecs.size(); ++j) {
1289  if (corrVecs[j] != NULL) delete corrVecs[j];
1290  }
1291 
1292  if (statisticalOptions.autoCorrDisplay()) {
1293  V subChainMean (m_vectorSpace.zeroVector());
1294  V subChainSampleVariance (m_vectorSpace.zeroVector());
1295  V estimatedVarianceOfSampleMean(m_vectorSpace.zeroVector());
1296  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1297  unsigned int initialPos = initialPosForStatistics[initialPosId];
1298 
1299  this->subMeanExtra(initialPos,
1300  this->subSequenceSize()-initialPos,
1301  subChainMean);
1302 
1303  this->subSampleVarianceExtra(initialPos,
1304  this->subSequenceSize()-initialPos,
1305  subChainMean,
1306  subChainSampleVariance);
1307 
1308  if (m_env.subDisplayFile()) {
1309  *m_env.subDisplayFile() << "\nEstimated variance of sample mean, through autocorrelation (via fft), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1310  << std::endl;
1311  }
1312  estimatedVarianceOfSampleMean.cwSet(-1.); // Yes, '-1' because the autocorrelation at lag 0, which values '+1', is already counted in the sum
1313  estimatedVarianceOfSampleMean += 2.* (*corrSumVecs[initialPosId]);
1314  estimatedVarianceOfSampleMean *= subChainSampleVariance;
1315  estimatedVarianceOfSampleMean /= (double) (this->subSequenceSize() - initialPos);
1316  bool savedVectorPrintState = estimatedVarianceOfSampleMean.getPrintHorizontally();
1317  estimatedVarianceOfSampleMean.setPrintHorizontally(false);
1318  if (m_env.subDisplayFile()) {
1319  *m_env.subDisplayFile() << estimatedVarianceOfSampleMean
1320  << std::endl;
1321  }
1322  estimatedVarianceOfSampleMean.setPrintHorizontally(savedVectorPrintState);
1323 
1324  if (m_env.subDisplayFile()) {
1325  *m_env.subDisplayFile() << "\nComputed autocorrelation coefficients (via fft), for subchain beginning at position " << initialPosForStatistics[initialPosId]
1326  << " (each column corresponds to a different lag)"
1327  << std::endl;
1328 
1329  char line[512];
1330  sprintf(line,"%s",
1331  "Parameter");
1332  *m_env.subDisplayFile() << line;
1333  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1334  sprintf(line,"%10s%3d",
1335  " ",
1336  lagsForCorrs[lagId]);
1337  *m_env.subDisplayFile() << line;
1338  }
1339 
1340  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1341  sprintf(line,"\n%9.9s",
1342  m_vectorSpace.localComponentName(i).c_str() /*.*/);
1343  *m_env.subDisplayFile() << line;
1344  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1345  sprintf(line,"%2s%11.4e",
1346  " ",
1347  _2dArrayOfAutoCorrs(initialPosId,lagId)[i]);
1348  *m_env.subDisplayFile() << line;
1349  }
1350  }
1351  *m_env.subDisplayFile() << std::endl;
1352  }
1353  }
1354  }
1355 
1356  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
1357  if (m_env.subDisplayFile()) {
1358  *m_env.subDisplayFile() << "Chain autocorrelation (via fft) took " << tmpRunTime
1359  << " seconds"
1360  << std::endl;
1361  }
1362 
1363  // Write autocorrelations
1364  if (statisticalOptions.autoCorrWrite() && passedOfs) {
1365  std::ofstream& ofsvar = *passedOfs;
1366  ofsvar << m_name << "_corrViaFftLags_sub" << m_env.subIdString() << " = zeros(" << 1
1367  << "," << lagsForCorrs.size()
1368  << ");"
1369  << std::endl;
1370  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1371  ofsvar << m_name << "_corrViaFftLags_sub" << m_env.subIdString() << "(" << 1
1372  << "," << lagId+1
1373  << ") = " << lagsForCorrs[lagId]
1374  << ";"
1375  << std::endl;
1376  }
1377 
1378  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
1379  ofsvar << m_name << "_corrViaFftInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << " = zeros(" << this->vectorSizeLocal() /*.*/
1380  << "," << lagsForCorrs.size()
1381  << ");"
1382  << std::endl;
1383  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1384  for (unsigned int lagId = 0; lagId < lagsForCorrs.size(); lagId++) {
1385  ofsvar << m_name << "_corrViaFftInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << "(" << i+1
1386  << "," << lagId+1
1387  << ") = " << _2dArrayOfAutoCorrs(initialPosId,lagId)[i]
1388  << ";"
1389  << std::endl;
1390  }
1391  }
1392  }
1393  }
1394 
1395  return;
1396 }
1397 // --------------------------------------------------
1398 template<class V, class M>
1399 void
1400 BaseVectorSequence<V,M>::computeHistCdfstaccKde( // Use the whole chain
1401  const SequenceStatisticalOptions& statisticalOptions,
1402  std::ofstream* passedOfs)
1403 {
1404  if (m_env.subDisplayFile()) {
1405  *m_env.subDisplayFile() << "\n"
1406  << "\n-----------------------------------------------------"
1407  << "\n Computing histogram and/or cdf stacc and/or Kde for chain '" << m_name << "' ..."
1408  << "\n-----------------------------------------------------"
1409  << "\n"
1410  << std::endl;
1411  }
1412 
1413  int iRC = UQ_OK_RC;
1414  struct timeval timevalTmp;
1415 
1416  //****************************************************
1417  // Compute MIN and MAX: for histograms and Kde
1418  //****************************************************
1419  double tmpRunTime = 0.;
1420  iRC = gettimeofday(&timevalTmp, NULL);
1421  if (m_env.subDisplayFile()) {
1422  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
1423  << "\nComputing min and max for histograms and Kde"
1424  << std::endl;
1425  }
1426 
1427  V statsMinPositions(m_vectorSpace.zeroVector());
1428  V statsMaxPositions(m_vectorSpace.zeroVector());
1429  this->subMinMaxExtra(0, // Use the whole chain
1430  this->subSequenceSize(),
1431  statsMinPositions,
1432  statsMaxPositions);
1433 
1434  if (m_env.subDisplayFile()) {
1435  *m_env.subDisplayFile() << "\nComputed min values and max values for chain '" << m_name << "'"
1436  << std::endl;
1437 
1438  char line[512];
1439  sprintf(line,"%s",
1440  "Parameter");
1441  *m_env.subDisplayFile() << line;
1442 
1443  sprintf(line,"%9s%s%9s%s",
1444  " ",
1445  "min",
1446  " ",
1447  "max");
1448  *m_env.subDisplayFile() << line;
1449 
1450  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1451  sprintf(line,"\n%8.8s",
1452  m_vectorSpace.localComponentName(i).c_str() /*.*/);
1453  *m_env.subDisplayFile() << line;
1454 
1455  sprintf(line,"%2s%11.4e%2s%11.4e",
1456  " ",
1457  statsMinPositions[i],
1458  " ",
1459  statsMaxPositions[i]);
1460  *m_env.subDisplayFile() << line;
1461  }
1462  *m_env.subDisplayFile() << std::endl;
1463  }
1464 
1465  V unifiedStatsMinPositions(statsMinPositions);
1466  V unifiedStatsMaxPositions(statsMaxPositions);
1467  if (m_env.numSubEnvironments() > 1) {
1468  // Compute unified min-max
1469  this->unifiedMinMaxExtra(0, // Use the whole chain
1470  this->subSequenceSize(),
1471  unifiedStatsMinPositions,
1472  unifiedStatsMaxPositions);
1473 
1474  // Write unified min-max
1475  if (m_env.subDisplayFile()) {
1476  if (m_vectorSpace.numOfProcsForStorage() == 1) {
1477  if (m_env.inter0Rank() == 0) {
1478  *m_env.subDisplayFile() << "\nComputed unified min values and max values for chain '" << m_name << "'"
1479  << std::endl;
1480 
1481  char line[512];
1482  sprintf(line,"%s",
1483  "Parameter");
1484  *m_env.subDisplayFile() << line;
1485 
1486  sprintf(line,"%9s%s%9s%s",
1487  " ",
1488  "min",
1489  " ",
1490  "max");
1491  *m_env.subDisplayFile() << line;
1492 
1493  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1494  sprintf(line,"\n%8.8s",
1495  m_vectorSpace.localComponentName(i).c_str() /*.*/);
1496  *m_env.subDisplayFile() << line;
1497 
1498  sprintf(line,"%2s%11.4e%2s%11.4e",
1499  " ",
1500  unifiedStatsMinPositions[i],
1501  " ",
1502  unifiedStatsMaxPositions[i]);
1503  *m_env.subDisplayFile() << line;
1504  }
1505  *m_env.subDisplayFile() << std::endl;
1506  }
1507  }
1508  else {
1509  UQ_FATAL_TEST_MACRO(true,
1510  m_env.worldRank(),
1511  "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1512  "unified min-max writing, parallel vectors not supported yet");
1513  }
1514  } // if subDisplayFile
1515  } // if numSubEnvs > 1
1516 
1517  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1518  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
1519  if (m_env.subDisplayFile()) {
1520  *m_env.subDisplayFile() << "Chain min and max took " << tmpRunTime
1521  << " seconds"
1522  << std::endl;
1523  }
1524 
1525  //****************************************************
1526  // Compute histograms
1527  //****************************************************
1528 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
1529  if ((statisticalOptions.histCompute() ) &&
1530  (statisticalOptions.histNumInternalBins() > 0)) {
1531  tmpRunTime = 0.;
1532  iRC = gettimeofday(&timevalTmp, NULL);
1533  if (m_env.subDisplayFile()) {
1534  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
1535  << "\nComputing histograms"
1536  << std::endl;
1537  }
1538 
1539  std::string subCoreName_HistCenters((std::string)( "_HistCenters_sub")+m_env.subIdString());
1540  std::string uniCoreName_HistCenters((std::string)("_unifHistCenters_sub")+m_env.subIdString());
1541  if (m_env.numSubEnvironments() == 1) subCoreName_HistCenters = uniCoreName_HistCenters;
1542 
1543  std::string subCoreName_HistQuantts((std::string)( "_HistQuantts_sub")+m_env.subIdString());
1544  std::string uniCoreName_HistQuantts((std::string)("_unifHistQuantts_sub")+m_env.subIdString());
1545  if (m_env.numSubEnvironments() == 1) subCoreName_HistQuantts = uniCoreName_HistQuantts;
1546 
1547  for (unsigned int i = 0; i < statsMaxPositions.sizeLocal(); ++i) {
1548  statsMaxPositions[i] *= (1. + 1.e-15);
1549  }
1550 
1551  // Compute histograms
1552  std::vector<V*> histCentersForAllBins(statisticalOptions.histNumInternalBins()+2,NULL); // IMPORTANT: +2
1553  std::vector<V*> histQuanttsForAllBins(statisticalOptions.histNumInternalBins()+2,NULL); // IMPORTANT: +2
1554  this->subHistogram(0, // Use the whole chain
1555  statsMinPositions,
1556  statsMaxPositions,
1557  histCentersForAllBins,
1558  histQuanttsForAllBins);
1559 
1560  // Write histograms
1561  if (passedOfs) {
1562  std::ofstream& ofsvar = *passedOfs;
1563  ofsvar << m_name << subCoreName_HistCenters << " = zeros(" << this->vectorSizeLocal() /*.*/
1564  << "," << histCentersForAllBins.size()
1565  << ");"
1566  << std::endl;
1567  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1568  for (unsigned int j = 0; j < histCentersForAllBins.size(); ++j) {
1569  ofsvar << m_name << subCoreName_HistCenters << "(" << i+1
1570  << "," << j+1
1571  << ") = " << (*(histCentersForAllBins[j]))[i]
1572  << ";"
1573  << std::endl;
1574  }
1575  }
1576 
1577  ofsvar << m_name << subCoreName_HistQuantts << " = zeros(" << this->vectorSizeLocal() /*.*/
1578  << "," << histQuanttsForAllBins.size()
1579  << ");"
1580  << std::endl;
1581  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1582  for (unsigned int j = 0; j < histQuanttsForAllBins.size(); ++j) {
1583  ofsvar << m_name << subCoreName_HistQuantts << "(" << i+1
1584  << "," << j+1
1585  << ") = " << (*(histQuanttsForAllBins[j]))[i]
1586  << ";"
1587  << std::endl;
1588  }
1589  }
1590  }
1591 
1592  for (unsigned int i = 0; i < histQuanttsForAllBins.size(); ++i) {
1593  if (histQuanttsForAllBins[i] != NULL) delete histQuanttsForAllBins[i];
1594  }
1595  for (unsigned int i = 0; i < histCentersForAllBins.size(); ++i) {
1596  if (histCentersForAllBins[i] != NULL) delete histCentersForAllBins[i];
1597  }
1598 
1599  std::vector<V*> unifiedHistCentersForAllBins(statisticalOptions.histNumInternalBins()+2,NULL); // IMPORTANT: +2
1600  std::vector<V*> unifiedHistQuanttsForAllBins(statisticalOptions.histNumInternalBins()+2,NULL); // IMPORTANT: +2
1601  if (m_env.numSubEnvironments() > 1) {
1602  // Compute unified histogram
1603  this->unifiedHistogram(0, // Use the whole chain
1604  unifiedStatsMinPositions,
1605  unifiedStatsMaxPositions,
1606  unifiedHistCentersForAllBins,
1607  unifiedHistQuanttsForAllBins);
1608 
1609  // Write unified histogram
1610  if (passedOfs) {
1611  if (m_vectorSpace.numOfProcsForStorage() == 1) {
1612  if (m_env.inter0Rank() == 0) {
1613  std::ofstream& ofsvar = *passedOfs;
1614  ofsvar << m_name << uniCoreName_HistCenters << " = zeros(" << this->vectorSizeLocal() /*.*/
1615  << "," << unifiedHistCentersForAllBins.size()
1616  << ");"
1617  << std::endl;
1618  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1619  for (unsigned int j = 0; j < unifiedHistCentersForAllBins.size(); ++j) {
1620  ofsvar << m_name << uniCoreName_HistCenters << "(" << i+1
1621  << "," << j+1
1622  << ") = " << (*(unifiedHistCentersForAllBins[j]))[i]
1623  << ";"
1624  << std::endl;
1625  }
1626  }
1627 
1628  ofsvar << m_name << uniCoreName_HistQuantts << " = zeros(" << this->vectorSizeLocal() /*.*/
1629  << "," << unifiedHistQuanttsForAllBins.size()
1630  << ");"
1631  << std::endl;
1632  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1633  for (unsigned int j = 0; j < unifiedHistQuanttsForAllBins.size(); ++j) {
1634  ofsvar << m_name << uniCoreName_HistQuantts << "(" << i+1
1635  << "," << j+1
1636  << ") = " << (*(unifiedHistQuanttsForAllBins[j]))[i]
1637  << ";"
1638  << std::endl;
1639  }
1640  }
1641  }
1642  }
1643  else {
1644  UQ_FATAL_TEST_MACRO(true,
1645  m_env.worldRank(),
1646  "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1647  "unified histogram writing, parallel vectors not supported yet");
1648  }
1649  } // if passedOfs
1650 
1651  for (unsigned int i = 0; i < unifiedHistQuanttsForAllBins.size(); ++i) {
1652  if (unifiedHistQuanttsForAllBins[i] != NULL) delete unifiedHistQuanttsForAllBins[i];
1653  }
1654  for (unsigned int i = 0; i < unifiedHistCentersForAllBins.size(); ++i) {
1655  if (unifiedHistCentersForAllBins[i] != NULL) delete unifiedHistCentersForAllBins[i];
1656  }
1657  } // if numSubEnvs > 1
1658 
1659  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1660  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
1661  if (m_env.subDisplayFile()) {
1662  *m_env.subDisplayFile() << "Chain histograms took " << tmpRunTime
1663  << " seconds"
1664  << std::endl;
1665  }
1666  }
1667 #endif
1668 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
1669  //****************************************************
1670  // Compute cdf statistical accuracy
1671  //****************************************************
1672  if ((statisticalOptions.cdfStaccCompute() ) &&
1673  (statisticalOptions.cdfStaccNumEvalPositions() > 0)) {
1674  tmpRunTime = 0.;
1675  iRC = gettimeofday(&timevalTmp, NULL);
1676  if (m_env.subDisplayFile()) {
1677  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
1678  << "\nComputing cdf statistical accuracy"
1679  << std::endl;
1680  }
1681 
1682  std::vector<V*> cdfStaccEvalPositions(statisticalOptions.cdfStaccNumEvalPositions(),NULL);
1683  MiscComputePositionsBetweenMinMax(statsMinPositions,
1684  statsMaxPositions,
1685  cdfStaccEvalPositions);
1686 
1687  std::vector<V*> cdfStaccValues(statisticalOptions.cdfStaccNumEvalPositions(),NULL);
1688  this->subCdfStacc(0, // Use the whole chain
1689  cdfStaccEvalPositions,
1690  cdfStaccValues);
1691 
1692  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1693  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
1694  if (m_env.subDisplayFile()) {
1695  *m_env.subDisplayFile() << "Chain cdf statistical accuracy took " << tmpRunTime
1696  << " seconds"
1697  << std::endl;
1698  }
1699  }
1700 #endif
1701  //****************************************************
1702  // Compute estimations of probability densities
1703  //****************************************************
1704  if ((statisticalOptions.kdeCompute() ) &&
1705  (statisticalOptions.kdeNumEvalPositions() > 0)) {
1706  tmpRunTime = 0.;
1707  iRC = gettimeofday(&timevalTmp, NULL);
1708  if (m_env.subDisplayFile()) {
1709  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
1710  << "\nComputing Kde"
1711  << std::endl;
1712  }
1713 
1714  std::string subCoreName_GaussianKdePositions((std::string)( "_GkdePosits_sub")+m_env.subIdString());
1715  std::string uniCoreName_GaussianKdePositions((std::string)("_unifGkdePosits_sub")+m_env.subIdString());
1716  if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdePositions = uniCoreName_GaussianKdePositions; // avoid temporarily (see '< -1' below)
1717 
1718  std::string subCoreName_GaussianKdeScaleVec ((std::string)( "_GkdeScalev_sub")+m_env.subIdString());
1719  std::string uniCoreName_GaussianKdeScaleVec ((std::string)("_unifGkdeScalev_sub")+m_env.subIdString());
1720  if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdeScaleVec = uniCoreName_GaussianKdeScaleVec; // avoid temporarily (see '< -1' below)
1721 
1722  std::string subCoreName_GaussianKdeValues ((std::string)( "_GkdeValues_sub")+m_env.subIdString());
1723  std::string uniCoreName_GaussianKdeValues ((std::string)("_unifGkdeValues_sub")+m_env.subIdString());
1724  if (m_env.numSubEnvironments() == 1) subCoreName_GaussianKdeValues = uniCoreName_GaussianKdeValues; // avoid temporarily (see '< -1' below)
1725 
1726  V iqrVec(m_vectorSpace.zeroVector());
1727  this->subInterQuantileRange(0, // Use the whole chain
1728  iqrVec);
1729 
1730  V gaussianKdeScaleVec(m_vectorSpace.zeroVector());
1731  this->subScalesForKde(0, // Use the whole chain
1732  iqrVec,
1733  1,
1734  gaussianKdeScaleVec);
1735 
1736  std::vector<V*> kdeEvalPositions(statisticalOptions.kdeNumEvalPositions(),NULL);
1737  MiscComputePositionsBetweenMinMax(statsMinPositions,
1738  statsMaxPositions,
1739  kdeEvalPositions);
1740 
1741  std::vector<V*> gaussianKdeDensities(statisticalOptions.kdeNumEvalPositions(),NULL);
1742  this->subGaussian1dKde(0, // Use the whole chain
1743  gaussianKdeScaleVec,
1744  kdeEvalPositions,
1745  gaussianKdeDensities);
1746 
1747  // Write iqr
1748  if (m_env.subDisplayFile()) {
1749  *m_env.subDisplayFile() << "\nComputed inter quantile ranges for chain '" << m_name << "'"
1750  << std::endl;
1751 
1752  char line[512];
1753  sprintf(line,"%s",
1754  "Parameter");
1755  *m_env.subDisplayFile() << line;
1756 
1757  sprintf(line,"%9s%s",
1758  " ",
1759  "iqr");
1760  *m_env.subDisplayFile() << line;
1761 
1762  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1763  sprintf(line,"\n%8.8s",
1764  m_vectorSpace.localComponentName(i).c_str() /*.*/);
1765  *m_env.subDisplayFile() << line;
1766 
1767  sprintf(line,"%2s%11.4e",
1768  " ",
1769  iqrVec[i]);
1770  *m_env.subDisplayFile() << line;
1771  }
1772  *m_env.subDisplayFile() << std::endl;
1773  }
1774 
1775  // Write Kde
1776  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) { // output debug
1777  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1778  << ", for chain '" << m_name << "'"
1779  << ", about to write sub kde to ofstream with pointer = " << passedOfs
1780  << std::endl;
1781  }
1782  if (passedOfs) {
1783  std::ofstream& ofsvar = *passedOfs;
1784  ofsvar << m_name << subCoreName_GaussianKdePositions << " = zeros(" << this->vectorSizeLocal() /*.*/
1785  << "," << kdeEvalPositions.size()
1786  << ");"
1787  << std::endl;
1788  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1789  for (unsigned int j = 0; j < kdeEvalPositions.size(); ++j) {
1790  ofsvar << m_name << subCoreName_GaussianKdePositions << "(" << i+1
1791  << "," << j+1
1792  << ") = " << (*(kdeEvalPositions[j]))[i]
1793  << ";"
1794  << std::endl;
1795  }
1796  }
1797 
1798  ofsvar << m_name << subCoreName_GaussianKdeScaleVec << " = zeros(" << this->vectorSizeLocal() /*.*/
1799  << ");"
1800  << std::endl;
1801  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1802  ofsvar << m_name << subCoreName_GaussianKdeScaleVec << "(" << i+1
1803  << ") = " << gaussianKdeScaleVec[i]
1804  << ";"
1805  << std::endl;
1806  }
1807 
1808  ofsvar << m_name << subCoreName_GaussianKdeValues << " = zeros(" << this->vectorSizeLocal() /*.*/
1809  << "," << gaussianKdeDensities.size()
1810  << ");"
1811  << std::endl;
1812  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1813  for (unsigned int j = 0; j < gaussianKdeDensities.size(); ++j) {
1814  ofsvar << m_name << subCoreName_GaussianKdeValues << "(" << i+1
1815  << "," << j+1
1816  << ") = " << (*(gaussianKdeDensities[j]))[i]
1817  << ";"
1818  << std::endl;
1819  }
1820  }
1821  }
1822 
1823  for (unsigned int i = 0; i < gaussianKdeDensities.size(); ++i) {
1824  if (gaussianKdeDensities[i] != NULL) delete gaussianKdeDensities[i];
1825  }
1826  for (unsigned int i = 0; i < kdeEvalPositions.size(); ++i) {
1827  if (kdeEvalPositions[i] != NULL) delete kdeEvalPositions[i];
1828  }
1829 
1830  if ((int) m_env.numSubEnvironments() > 1) { // avoid code temporarily
1831  // Compute unified Kde
1832  V unifiedIqrVec(m_vectorSpace.zeroVector());
1833  this->unifiedInterQuantileRange(0, // Use the whole chain
1834  unifiedIqrVec);
1835  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1836 
1837  V unifiedGaussianKdeScaleVec(m_vectorSpace.zeroVector());
1838  this->unifiedScalesForKde(0, // Use the whole chain
1839  unifiedIqrVec,
1840  1,
1841  unifiedGaussianKdeScaleVec);
1842  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1843 
1844  std::vector<V*> unifiedKdeEvalPositions(statisticalOptions.kdeNumEvalPositions(),NULL);
1845  MiscComputePositionsBetweenMinMax(unifiedStatsMinPositions,
1846  unifiedStatsMaxPositions,
1847  unifiedKdeEvalPositions);
1848 
1849  std::vector<V*> unifiedGaussianKdeDensities(statisticalOptions.kdeNumEvalPositions(),NULL);
1850  this->unifiedGaussian1dKde(0, // Use the whole chain
1851  unifiedGaussianKdeScaleVec,
1852  unifiedKdeEvalPositions,
1853  unifiedGaussianKdeDensities);
1854  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1855 
1856  // Write unified iqr
1857  if (m_env.subDisplayFile()) {
1858  if (m_vectorSpace.numOfProcsForStorage() == 1) {
1859  if (m_env.inter0Rank() == 0) {
1860  *m_env.subDisplayFile() << "\nComputed unified inter quantile ranges for chain '" << m_name << "'"
1861  << std::endl;
1862 
1863  char line[512];
1864  sprintf(line,"%s",
1865  "Parameter");
1866  *m_env.subDisplayFile() << line;
1867 
1868  sprintf(line,"%9s%s",
1869  " ",
1870  "iqr");
1871  *m_env.subDisplayFile() << line;
1872 
1873  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1874  sprintf(line,"\n%8.8s",
1875  m_vectorSpace.localComponentName(i).c_str() /*.*/);
1876  *m_env.subDisplayFile() << line;
1877 
1878  sprintf(line,"%2s%11.4e",
1879  " ",
1880  unifiedIqrVec[i]);
1881  *m_env.subDisplayFile() << line;
1882  }
1883  *m_env.subDisplayFile() << std::endl;
1884  }
1885  }
1886  else {
1887  UQ_FATAL_TEST_MACRO(true,
1888  m_env.worldRank(),
1889  "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1890  "unified iqr writing, parallel vectors not supported yet");
1891  }
1892  }
1893 
1894  // Write unified Kde
1895  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) { // output debug
1896  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1897  << ", for chain '" << m_name << "'"
1898  << ", about to write unified kde to ofstream with pointer = " << passedOfs
1899  << std::endl;
1900  }
1901  if (passedOfs) {
1902  if (m_vectorSpace.numOfProcsForStorage() == 1) {
1903  if (m_env.inter0Rank() == 0) {
1904  std::ofstream& ofsvar = *passedOfs;
1905  ofsvar << m_name << uniCoreName_GaussianKdePositions << " = zeros(" << this->vectorSizeLocal() /*.*/
1906  << "," << unifiedKdeEvalPositions.size()
1907  << ");"
1908  << std::endl;
1909  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() > 10)) { // output debug
1910  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeHistCdfstaccKde()"
1911  << ", for chain '" << m_name << "'"
1912  << ": just wrote '... = zeros(.,.);' line to output file, which has pointer = " << passedOfs
1913  << std::endl;
1914  }
1915  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1916  for (unsigned int j = 0; j < unifiedKdeEvalPositions.size(); ++j) {
1917  ofsvar << m_name << uniCoreName_GaussianKdePositions << "(" << i+1
1918  << "," << j+1
1919  << ") = " << (*(unifiedKdeEvalPositions[j]))[i]
1920  << ";"
1921  << std::endl;
1922  }
1923  }
1924 
1925  ofsvar << m_name << uniCoreName_GaussianKdeScaleVec << " = zeros(" << this->vectorSizeLocal() /*.*/
1926  << ");"
1927  << std::endl;
1928  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1929  ofsvar << m_name << uniCoreName_GaussianKdeScaleVec << "(" << i+1
1930  << ") = " << unifiedGaussianKdeScaleVec[i]
1931  << ";"
1932  << std::endl;
1933  }
1934 
1935  ofsvar << m_name << uniCoreName_GaussianKdeValues << " = zeros(" << this->vectorSizeLocal() /*.*/
1936  << "," << unifiedGaussianKdeDensities.size()
1937  << ");"
1938  << std::endl;
1939  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
1940  for (unsigned int j = 0; j < unifiedGaussianKdeDensities.size(); ++j) {
1941  ofsvar << m_name << uniCoreName_GaussianKdeValues << "(" << i+1
1942  << "," << j+1
1943  << ") = " << (*(unifiedGaussianKdeDensities[j]))[i]
1944  << ";"
1945  << std::endl;
1946  }
1947  }
1948  }
1949  }
1950  else {
1951  UQ_FATAL_TEST_MACRO(true,
1952  m_env.worldRank(),
1953  "BaseVectorSequence<V,M>::computeHistCdfstaccKde()",
1954  "unified Kde writing, parallel vectors not supported yet");
1955  }
1956  }
1957 
1958  for (unsigned int i = 0; i < unifiedGaussianKdeDensities.size(); ++i) {
1959  if (unifiedGaussianKdeDensities[i] != NULL) delete unifiedGaussianKdeDensities[i];
1960  }
1961  for (unsigned int i = 0; i < unifiedKdeEvalPositions.size(); ++i) {
1962  if (unifiedKdeEvalPositions[i] != NULL) delete unifiedKdeEvalPositions[i];
1963  }
1964  }
1965 
1966  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1967  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
1968  if (m_env.subDisplayFile()) {
1969  *m_env.subDisplayFile() << "Chain Kde took " << tmpRunTime
1970  << " seconds"
1971  << std::endl;
1972  }
1973  }
1974 
1975  if (m_env.subDisplayFile()) {
1976  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
1977  << "\n Finished computing histogram and/or cdf stacc and/or Kde for chain '" << m_name << "'"
1978  << "\n-----------------------------------------------------"
1979  << "\n"
1980  << std::endl;
1981  }
1982 
1983  return;
1984 }
1985 // --------------------------------------------------
1986 template<class V, class M>
1987 void
1988 BaseVectorSequence<V,M>::computeCovCorrMatrices( // Use the whole chain
1989  const SequenceStatisticalOptions& statisticalOptions,
1990  std::ofstream* passedOfs)
1991 {
1992  if (m_env.subDisplayFile()) {
1993  *m_env.subDisplayFile() << "\n"
1994  << "\n-----------------------------------------------------"
1995  << "\n Computing covariance and correlation matrices for chain '" << m_name << "' ..."
1996  << "\n-----------------------------------------------------"
1997  << "\n"
1998  << std::endl;
1999  }
2000 
2001  //int iRC = UQ_OK_RC;
2002  //struct timeval timevalTmp;
2003  M* covarianceMatrix = new M(m_env,
2004  m_vectorSpace.map(), // number of rows
2005  m_vectorSpace.dimGlobal()); // number of cols
2006  M* correlationMatrix = new M(m_env,
2007  m_vectorSpace.map(), // number of rows
2008  m_vectorSpace.dimGlobal()); // number of cols
2009 
2010  ComputeCovCorrMatricesBetweenVectorSequences(*this,
2011  *this,
2012  this->subSequenceSize(),
2013  *covarianceMatrix,
2014  *correlationMatrix);
2015 
2016  if (m_env.subDisplayFile()) {
2017  if (m_vectorSpace.numOfProcsForStorage() == 1) {
2018  // Only unified covariance matrix is written. And only one processor writes it.
2019  if (m_env.inter0Rank() == 0) {
2020  *m_env.subDisplayFile() << "\nBaseVectorSequence<V,M>::computeCovCorrMatrices"
2021  << ", chain " << m_name
2022  << ": contents of covariance matrix are\n" << *covarianceMatrix
2023  << std::endl;
2024 
2025  *m_env.subDisplayFile() << "\nBaseVectorSequence<V,M>::computeCovCorrMatrices"
2026  << ", chain " << m_name
2027  << ": contents of correlation matrix are\n" << *correlationMatrix
2028  << std::endl;
2029  }
2030  }
2031  else {
2032  UQ_FATAL_TEST_MACRO(true,
2033  m_env.worldRank(),
2034  "BaseVectorSequence<V,M>::computeCovCorrMatrices()",
2035  "parallel vectors not supported yet");
2036  }
2037  }
2038 
2039  delete correlationMatrix;
2040  delete covarianceMatrix;
2041 
2042  if (m_env.subDisplayFile()) {
2043  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2044  << "\n Finished computing covariance and correlation matrices for chain '" << m_name << "'"
2045  << "\n-----------------------------------------------------"
2046  << "\n"
2047  << std::endl;
2048  }
2049 
2050  return;
2051 }
2052 #endif // #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2053 
2054 // --------------------------------------------------
2055 // --------------------------------------------------
2056 // --------------------------------------------------
2057 
2058 #ifdef UQ_CODE_HAS_MONITORS
2059 template<class V, class M>
2060 void
2061 BaseVectorSequence<V,M>::computeMeanEvolution(
2062  const SequenceStatisticalOptions& statisticalOptions,
2063  std::ofstream* passedOfs)
2064 {
2065  int iRC = UQ_OK_RC;
2066  struct timeval timevalTmp;
2067  iRC = gettimeofday(&timevalTmp, NULL);
2068  double tmpRunTime = 0.;
2069 
2070  if (m_env.subDisplayFile()) {
2071  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2072  << "\nComputing mean evolution"
2073  << std::endl;
2074  }
2075 
2076  unsigned int monitorPeriod = statisticalOptions.meanMonitorPeriod();
2077  unsigned int iMin = 0;
2078  unsigned int iMax = (unsigned int) ( ((double) this->subSequenceSize())/((double) monitorPeriod) );
2079 
2080  if (monitorPeriod == 1) {
2081  iMin = 1;
2082  }
2083 
2084  if (m_env.subDisplayFile()) {
2085  *m_env.subDisplayFile() << " Sub sequence size = " << this->subSequenceSize()
2086  << "\n Monitor period = " << monitorPeriod
2087  << "\n Number of monitoring positions = " << iMax
2088  << std::endl;
2089  }
2090 
2091  this->subMeanMonitorAlloc(iMax);
2092  if (m_env.numSubEnvironments() > 1) {
2093  this->subMeanInter0MonitorAlloc(iMax);
2094  this->unifiedMeanMonitorAlloc(iMax);
2095  }
2096 
2097  for (unsigned int i = iMin; i < iMax; ++i) {
2098  unsigned int currentMonitoredFinalPosition = (i+1)*monitorPeriod - 1; // Yes, '-1'
2099  V subMeanVec (m_vectorSpace.zeroVector());
2100  V subMeanCltStd(m_vectorSpace.zeroVector());
2101  this->subMeanMonitorRun(currentMonitoredFinalPosition,
2102  subMeanVec,
2103  subMeanCltStd);
2104  this->subMeanMonitorStore(i,
2105  currentMonitoredFinalPosition,
2106  subMeanVec,
2107  subMeanCltStd);
2108 
2109  if (m_env.numSubEnvironments() > 1) {
2110  V subMeanInter0Mean (m_vectorSpace.zeroVector());
2111  V subMeanInter0Clt95 (m_vectorSpace.zeroVector());
2112  V subMeanInter0Empirical90(m_vectorSpace.zeroVector());
2113  V subMeanInter0Min (m_vectorSpace.zeroVector());
2114  V subMeanInter0Max (m_vectorSpace.zeroVector());
2115  this->subMeanInter0MonitorRun(currentMonitoredFinalPosition,
2116  subMeanInter0Mean,
2117  subMeanInter0Clt95,
2118  subMeanInter0Empirical90,
2119  subMeanInter0Min,
2120  subMeanInter0Max);
2121  this->subMeanInter0MonitorStore(i,
2122  currentMonitoredFinalPosition,
2123  subMeanInter0Mean,
2124  subMeanInter0Clt95,
2125  subMeanInter0Empirical90,
2126  subMeanInter0Min,
2127  subMeanInter0Max);
2128 
2129  V unifiedMeanVec (m_vectorSpace.zeroVector());
2130  V unifiedMeanCltStd(m_vectorSpace.zeroVector());
2131  this->unifiedMeanMonitorRun(currentMonitoredFinalPosition,
2132  unifiedMeanVec,
2133  unifiedMeanCltStd);
2134  this->unifiedMeanMonitorStore(i,
2135  currentMonitoredFinalPosition,
2136  unifiedMeanVec,
2137  unifiedMeanCltStd);
2138  }
2139  }
2140 
2141  if (passedOfs) {
2142  this->subMeanMonitorWrite(*passedOfs);
2143  if (m_env.numSubEnvironments() > 1) {
2144  this->subMeanInter0MonitorWrite(*passedOfs);
2145  this->unifiedMeanMonitorWrite(*passedOfs); // Yes, call 'unified' even though not all nodes might have 'passedOfs != NULL' ????????????? Ernesto
2146  }
2147  }
2148 
2149  this->subMeanMonitorFree();
2150  if (m_env.numSubEnvironments() > 1) {
2151  this->subMeanInter0MonitorFree();
2152  this->unifiedMeanMonitorFree();
2153  }
2154 
2155  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
2156  if (m_env.subDisplayFile()) {
2157  *m_env.subDisplayFile() << "Mean evolution took " << tmpRunTime
2158  << " seconds"
2159  << std::endl;
2160  }
2161 
2162  return;
2163 }
2164 #endif //#ifdef UQ_CODE_HAS_MONITORS
2165 
2166 // --------------------------------------------------
2167 // --------------------------------------------------
2168 // --------------------------------------------------
2169 
2170 
2171 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2172 template<class V, class M>
2173 void
2174 BaseVectorSequence<V,M>::computeBMM(
2175  const SequenceStatisticalOptions& statisticalOptions,
2176  const std::vector<unsigned int>& initialPosForStatistics,
2177  std::ofstream* passedOfs)
2178 {
2179  int iRC = UQ_OK_RC;
2180  struct timeval timevalTmp;
2181  iRC = gettimeofday(&timevalTmp, NULL);
2182  double tmpRunTime = 0.;
2183 
2184  if (m_env.subDisplayFile()) {
2185  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2186  << "\nComputing variance of sample mean through BMM"
2187  << std::endl;
2188  }
2189 
2190  if (m_env.subDisplayFile()) {
2191  *m_env.subDisplayFile() << "In BaseVectorSequence<V,M>::computeBMM(): lengths for batchs in BMM =";
2192  for (unsigned int i = 0; i < statisticalOptions.bmmLengths().size(); ++i) {
2193  *m_env.subDisplayFile() << " " << statisticalOptions.bmmLengths()[i];
2194  }
2195  *m_env.subDisplayFile() << std::endl;
2196  }
2197 
2198  TwoDArray<V> _2dArrayOfBMM(initialPosForStatistics.size(),statisticalOptions.bmmLengths().size());
2199  for (unsigned int i = 0; i < _2dArrayOfBMM.numRows(); ++i) {
2200  for (unsigned int j = 0; j < _2dArrayOfBMM.numCols(); ++j) {
2201  _2dArrayOfBMM.setLocation(i,j,new V(m_vectorSpace.zeroVector()) /*.*/);
2202  }
2203  }
2204  V bmmVec(m_vectorSpace.zeroVector());
2205  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2206  unsigned int initialPos = initialPosForStatistics[initialPosId];
2207  for (unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2208  unsigned int batchLength = statisticalOptions.bmmLengths()[batchLengthId];
2209  this->bmm(initialPos,
2210  batchLength,
2211  bmmVec);
2212  _2dArrayOfBMM(initialPosId,batchLengthId) = bmmVec;
2213  }
2214  }
2215 
2216  if (m_env.subDisplayFile()) {
2217  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2218  *m_env.subDisplayFile() << "\nEstimated variance of sample mean, through batch means method, for subchain beginning at position " << initialPosForStatistics[initialPosId]
2219  << " (each column corresponds to a batch length)"
2220  << std::endl;
2221 
2222  char line[512];
2223  sprintf(line,"%s",
2224  "Parameter");
2225  *m_env.subDisplayFile() << line;
2226  for (unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2227  sprintf(line,"%10s%3d",
2228  " ",
2229  statisticalOptions.bmmLengths()[batchLengthId]);
2230  *m_env.subDisplayFile() << line;
2231  }
2232 
2233  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
2234  sprintf(line,"\n%9.9s",
2235  m_vectorSpace.localComponentName(i).c_str() /*.*/);
2236  *m_env.subDisplayFile() << line;
2237  for (unsigned int batchLengthId = 0; batchLengthId < statisticalOptions.bmmLengths().size(); batchLengthId++) {
2238  sprintf(line,"%2s%11.4e",
2239  " ",
2240  _2dArrayOfBMM(initialPosId,batchLengthId)[i]);
2241  *m_env.subDisplayFile() << line;
2242  }
2243  }
2244  *m_env.subDisplayFile() << std::endl;
2245  }
2246  }
2247 
2248  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
2249  if (m_env.subDisplayFile()) {
2250  *m_env.subDisplayFile() << "Chain BMM took " << tmpRunTime
2251  << " seconds"
2252  << std::endl;
2253  }
2254 
2255  return;
2256 }
2257 // --------------------------------------------------
2258 template<class V, class M>
2259 void
2260 BaseVectorSequence<V,M>::computeFFT(
2261  const SequenceStatisticalOptions& statisticalOptions,
2262  const std::vector<unsigned int>& initialPosForStatistics,
2263  std::ofstream* passedOfs)
2264 {
2265  int iRC = UQ_OK_RC;
2266  struct timeval timevalTmp;
2267  iRC = gettimeofday(&timevalTmp, NULL);
2268  double tmpRunTime = 0.;
2269 
2270  if (m_env.subDisplayFile()) {
2271  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2272  << "\nComputing FFT of chain on parameter of id = " << statisticalOptions.fftParamId()
2273  << std::endl;
2274  }
2275 
2276  std::vector<std::complex<double> > forwardResult(0,std::complex<double>(0.,0.));
2277  std::vector<std::complex<double> > inverseResult(0,std::complex<double>(0.,0.));
2278  Fft<std::complex<double> > fftObj(m_env);
2279  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2280  unsigned int initialPosition = initialPosForStatistics[initialPosId];
2281  this->fftForward(initialPosition,
2282  statisticalOptions.fftSize(),
2283  statisticalOptions.fftParamId(),
2284  forwardResult);
2285 
2286  if (statisticalOptions.fftWrite() && passedOfs) {
2287  std::ofstream& ofsvar = *passedOfs;
2288  ofsvar << m_name << "_fftInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << " = zeros(" << 1
2289  << "," << forwardResult.size()
2290  << ");"
2291  << std::endl;
2292  for (unsigned int j = 0; j < forwardResult.size(); ++j) {
2293  ofsvar << m_name << "_fftInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << "(" << 1
2294  << "," << j+1
2295  << ") = " << forwardResult[j].real()
2296  << " + i*" << forwardResult[j].imag()
2297  << ";"
2298  << std::endl;
2299  }
2300  } // if write
2301 
2302  if (statisticalOptions.fftTestInversion()) {
2303  fftObj.inverse(forwardResult,
2304  statisticalOptions.fftSize(),
2305  inverseResult);
2306  if (statisticalOptions.fftWrite() && passedOfs) {
2307  std::ofstream& ofsvar = *passedOfs;
2308  ofsvar << m_name << "_iftInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << " = zeros(" << 1
2309  << "," << inverseResult.size()
2310  << ");"
2311  << std::endl;
2312  for (unsigned int j = 0; j < inverseResult.size(); ++j) {
2313  ofsvar << m_name << "_iftInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << "(" << 1
2314  << "," << j+1
2315  << ") = " << inverseResult[j].real()
2316  << " + i*" << inverseResult[j].imag()
2317  << ";"
2318  << std::endl;
2319  }
2320  } // if write
2321  }
2322  } // for initialPosId
2323 
2324  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
2325  if (m_env.subDisplayFile()) {
2326  *m_env.subDisplayFile() << "Chain FFT took " << tmpRunTime
2327  << " seconds"
2328  << std::endl;
2329  }
2330 
2331  return;
2332 }
2333 // --------------------------------------------------
2334 template<class V, class M>
2335 void
2336 BaseVectorSequence<V,M>::computePSD(
2337  const SequenceStatisticalOptions& statisticalOptions,
2338  const std::vector<unsigned int>& initialPosForStatistics,
2339  std::ofstream* passedOfs)
2340 {
2341  int iRC = UQ_OK_RC;
2342  struct timeval timevalTmp;
2343  iRC = gettimeofday(&timevalTmp, NULL);
2344  double tmpRunTime = 0.;
2345 
2346  if (m_env.subDisplayFile()) {
2347  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2348  << "\nComputing PSD of chain on parameter of id = " << statisticalOptions.psdParamId()
2349  << std::endl;
2350  }
2351 
2352  std::vector<double> psdResult(0,0.);
2353  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2354  unsigned int initialPosition = initialPosForStatistics[initialPosId];
2355  this->psd(initialPosition,
2356  statisticalOptions.psdNumBlocks(),
2357  statisticalOptions.psdHopSizeRatio(),
2358  statisticalOptions.psdParamId(),
2359  psdResult);
2360 
2361  if (statisticalOptions.psdWrite() && passedOfs) {
2362  std::ofstream& ofsvar = *passedOfs;
2363  ofsvar << m_name << "_psdInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << " = zeros(" << 1
2364  << "," << psdResult.size()
2365  << ");"
2366  << std::endl;
2367  for (unsigned int j = 0; j < psdResult.size(); ++j) {
2368  ofsvar << m_name << "_psdInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << "(" << 1
2369  << "," << j+1
2370  << ") = " << psdResult[j]
2371  << ";"
2372  << std::endl;
2373  }
2374  } // if write
2375  }
2376 
2377  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
2378  if (m_env.subDisplayFile()) {
2379  *m_env.subDisplayFile() << "Chain PSD took " << tmpRunTime
2380  << " seconds"
2381  << std::endl;
2382  }
2383 
2384  return;
2385 }
2386 // --------------------------------------------------
2387 template<class V, class M>
2388 void
2389 BaseVectorSequence<V,M>::computePSDAtZero(
2390  const SequenceStatisticalOptions& statisticalOptions,
2391  const std::vector<unsigned int>& initialPosForStatistics,
2392  std::ofstream* passedOfs)
2393 {
2394  int iRC = UQ_OK_RC;
2395  struct timeval timevalTmp;
2396  iRC = gettimeofday(&timevalTmp, NULL);
2397  double tmpRunTime = 0.;
2398 
2399  if (m_env.subDisplayFile()) {
2400  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2401  << "\nComputing PSD at frequency zero for all parameters"
2402  << std::endl;
2403  }
2404 
2405  TwoDArray<V> _2dArrayOfPSDAtZero(initialPosForStatistics.size(),statisticalOptions.psdAtZeroNumBlocks().size());
2406  for (unsigned int i = 0; i < _2dArrayOfPSDAtZero.numRows(); ++i) {
2407  for (unsigned int j = 0; j < _2dArrayOfPSDAtZero.numCols(); ++j) {
2408  _2dArrayOfPSDAtZero.setLocation(i,j,new V(m_vectorSpace.zeroVector()) /*.*/);
2409  }
2410  }
2411  V psdVec(m_vectorSpace.zeroVector());
2412  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2413  unsigned int initialPosition = initialPosForStatistics[initialPosId];
2414  for (unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2415  unsigned int numBlocks = statisticalOptions.psdAtZeroNumBlocks()[numBlocksId];
2416  this->psdAtZero(initialPosition,
2417  numBlocks,
2418  statisticalOptions.psdAtZeroHopSizeRatio(),
2419  psdVec);
2420  _2dArrayOfPSDAtZero(initialPosId,numBlocksId) = psdVec;
2421  }
2422  }
2423 
2424  // Display PSD at frequency zero
2425  if ((statisticalOptions.psdAtZeroDisplay()) && (m_env.subDisplayFile())) {
2426  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2427  unsigned int initialPos = initialPosForStatistics[initialPosId];
2428  *m_env.subDisplayFile() << "\nComputed PSD at frequency zero for subchain beginning at position " << initialPos
2429  << ", so effective data size = " << this->subSequenceSize() - initialPos
2430  << " (each column corresponds to a number of blocks)"
2431  << std::endl;
2432 
2433  char line[512];
2434  sprintf(line,"%s",
2435  "Parameter");
2436  *m_env.subDisplayFile() << line;
2437  for (unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2438  sprintf(line,"%10s%3d",
2439  " ",
2440  statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]);
2441  *m_env.subDisplayFile() << line;
2442  }
2443 
2444  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
2445  sprintf(line,"\n%9.9s",
2446  m_vectorSpace.localComponentName(i).c_str() /*.*/);
2447  *m_env.subDisplayFile() << line;
2448  for (unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2449  sprintf(line,"%2s%11.4e",
2450  " ",
2451  _2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]);
2452  *m_env.subDisplayFile() << line;
2453  }
2454  }
2455  *m_env.subDisplayFile() << std::endl;
2456  }
2457  }
2458 
2459  // Display estimated variance of sample mean through PSD
2460  if (/*(statisticalOptions.psdAtZeroDisplay()) &&*/ (m_env.subDisplayFile())) {
2461  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2462  unsigned int initialPos = initialPosForStatistics[initialPosId];
2463  *m_env.subDisplayFile() << "\nEstimated variance of sample mean, through psd, for subchain beginning at position " << initialPos
2464  << ", so effective data size = " << this->subSequenceSize() - initialPos
2465  << " (each column corresponds to a number of blocks)"
2466  << std::endl;
2467 
2468  char line[512];
2469  sprintf(line,"%s",
2470  "Parameter");
2471  *m_env.subDisplayFile() << line;
2472  for (unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2473  sprintf(line,"%10s%3d",
2474  " ",
2475  statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]);
2476  *m_env.subDisplayFile() << line;
2477  }
2478 
2479  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
2480  sprintf(line,"\n%9.9s",
2481  m_vectorSpace.localComponentName(i).c_str() /*.*/);
2482  *m_env.subDisplayFile() << line;
2483  for (unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2484  sprintf(line,"%2s%11.4e",
2485  " ",
2486  2.*M_PI*_2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]/(double) (this->subSequenceSize() - initialPos));
2487  *m_env.subDisplayFile() << line;
2488  }
2489  }
2490  *m_env.subDisplayFile() << std::endl;
2491  }
2492  }
2493 
2494  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
2495  if (m_env.subDisplayFile()) {
2496  *m_env.subDisplayFile() << "Chain PSD at frequency zero took " << tmpRunTime
2497  << " seconds"
2498  << std::endl;
2499  }
2500 
2501  // Write PSD at frequency zero
2502  if (statisticalOptions.psdAtZeroWrite() && passedOfs) {
2503  std::ofstream& ofsvar = *passedOfs;
2504  ofsvar << m_name << "_psdAtZeroNumBlocks_sub" << m_env.subIdString() << " = zeros(" << 1
2505  << "," << statisticalOptions.psdAtZeroNumBlocks().size()
2506  << ");"
2507  << std::endl;
2508  for (unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2509  ofsvar << m_name << "_psdAtZeroNumBlocks_sub" << m_env.subIdString() << "(" << 1
2510  << "," << numBlocksId+1
2511  << ") = " << statisticalOptions.psdAtZeroNumBlocks()[numBlocksId]
2512  << ";"
2513  << std::endl;
2514  }
2515 
2516  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2517  ofsvar << m_name << "_psdAtZeroInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << " = zeros(" << this->vectorSizeLocal() /*.*/
2518  << "," << statisticalOptions.psdAtZeroNumBlocks().size()
2519  << ");"
2520  << std::endl;
2521  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
2522  for (unsigned int numBlocksId = 0; numBlocksId < statisticalOptions.psdAtZeroNumBlocks().size(); numBlocksId++) {
2523  ofsvar << m_name << "_psdAtZeroInitPos" << initialPosForStatistics[initialPosId] << "_sub" << m_env.subIdString() << "(" << i+1
2524  << "," << numBlocksId+1
2525  << ") = " << _2dArrayOfPSDAtZero(initialPosId,numBlocksId)[i]
2526  << ";"
2527  << std::endl;
2528  }
2529  }
2530  }
2531  }
2532 
2533  return;
2534 }
2535 // --------------------------------------------------
2536 template<class V, class M>
2537 void
2538 BaseVectorSequence<V,M>::computeGeweke(
2539  const SequenceStatisticalOptions& statisticalOptions,
2540  const std::vector<unsigned int>& initialPosForStatistics,
2541  std::ofstream* passedOfs)
2542 {
2543  int iRC = UQ_OK_RC;
2544  struct timeval timevalTmp;
2545  iRC = gettimeofday(&timevalTmp, NULL);
2546  double tmpRunTime = 0.;
2547 
2548  if (m_env.subDisplayFile()) {
2549  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2550  << "\nComputing Geweke coefficients"
2551  << std::endl;
2552  }
2553 
2554  std::vector<V*> vectorOfGeweke(initialPosForStatistics.size(),NULL);
2555  V gewVec(m_vectorSpace.zeroVector());
2556  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2557  unsigned int initialPosition = initialPosForStatistics[initialPosId];
2558  this->geweke(initialPosition,
2559  statisticalOptions.gewekeNaRatio(),
2560  statisticalOptions.gewekeNbRatio(),
2561  gewVec);
2562  vectorOfGeweke[initialPosId] = new V(gewVec);
2563  }
2564 
2565  if (m_env.subDisplayFile()) {
2566  *m_env.subDisplayFile() << "\nComputed Geweke coefficients with 10% and 50% percentages"
2567  << " (each column corresponds to a different initial position on the full chain)"
2568  << std::endl;
2569 
2570  char line[512];
2571  sprintf(line,"%s",
2572  "Parameter");
2573  *m_env.subDisplayFile() << line;
2574  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2575  sprintf(line,"%10s%3d",
2576  " ",
2577  initialPosForStatistics[initialPosId]);
2578  *m_env.subDisplayFile() << line;
2579  }
2580 
2581  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
2582  sprintf(line,"\n%9.9s",
2583  m_vectorSpace.localComponentName(i).c_str() /*.*/);
2584  *m_env.subDisplayFile() << line;
2585  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2586  sprintf(line,"%2s%11.4e",
2587  " ",
2588  (*(vectorOfGeweke[initialPosId]))[i]);
2589  *m_env.subDisplayFile() << line;
2590  }
2591  }
2592  *m_env.subDisplayFile() << std::endl;
2593  }
2594 
2595  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
2596  if (m_env.subDisplayFile()) {
2597  *m_env.subDisplayFile() << "Chain Geweke took " << tmpRunTime
2598  << " seconds"
2599  << std::endl;
2600  }
2601 
2602  return;
2603 }
2604 // --------------------------------------------------
2605 template<class V, class M>
2606 void
2607 BaseVectorSequence<V,M>::computeMeanStacc(
2608  const SequenceStatisticalOptions& statisticalOptions,
2609  const std::vector<unsigned int>& initialPosForStatistics,
2610  std::ofstream* passedOfs)
2611 {
2612  int iRC = UQ_OK_RC;
2613  struct timeval timevalTmp;
2614  iRC = gettimeofday(&timevalTmp, NULL);
2615  double tmpRunTime = 0.;
2616 
2617  if (m_env.subDisplayFile()) {
2618  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
2619  << "\nComputing mean statistical accuracy"
2620  << std::endl;
2621  }
2622 
2623  std::vector<V*> vectorOfMeanStacc(initialPosForStatistics.size(),NULL);
2624  V meanStaccVec(m_vectorSpace.zeroVector());
2625  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2626  unsigned int initialPosition = initialPosForStatistics[initialPosId];
2627  this->meanStacc(initialPosition,
2628  meanStaccVec);
2629  vectorOfMeanStacc[initialPosId] = new V(meanStaccVec);
2630  }
2631 
2632  if (m_env.subDisplayFile()) {
2633  *m_env.subDisplayFile() << "\nComputed mean statistical accuracy"
2634  << " (each column corresponds to a different initial position on the full chain)"
2635  << std::endl;
2636 
2637  char line[512];
2638  sprintf(line,"%s",
2639  "Parameter");
2640  *m_env.subDisplayFile() << line;
2641  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2642  sprintf(line,"%10s%3d",
2643  " ",
2644  initialPosForStatistics[initialPosId]);
2645  *m_env.subDisplayFile() << line;
2646  }
2647 
2648  for (unsigned int i = 0; i < this->vectorSizeLocal() /*.*/; ++i) {
2649  sprintf(line,"\n%9.9s",
2650  m_vectorSpace.localComponentName(i).c_str() /*.*/);
2651  *m_env.subDisplayFile() << line;
2652  for (unsigned int initialPosId = 0; initialPosId < initialPosForStatistics.size(); initialPosId++) {
2653  sprintf(line,"%2s%11.4e",
2654  " ",
2655  (*(vectorOfMeanStacc[initialPosId]))[i]);
2656  *m_env.subDisplayFile() << line;
2657  }
2658  }
2659  *m_env.subDisplayFile() << std::endl;
2660  }
2661 
2662  tmpRunTime += MiscGetEllapsedSeconds(&timevalTmp);
2663  if (m_env.subDisplayFile()) {
2664  *m_env.subDisplayFile() << "Chain mean statistical accuracy took " << tmpRunTime
2665  << " seconds"
2666  << std::endl;
2667  }
2668 
2669  return;
2670 }
2671 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2672 
2673 // --------------------------------------------------
2674 // Additional methods --------------------------------
2675 // (outside class declaration) ----------------------
2676 //---------------------------------------------------
2677 template <class P_V, class P_M, class Q_V, class Q_M>
2678 void
2679 ComputeCovCorrMatricesBetweenVectorSequences(
2680  const BaseVectorSequence<P_V,P_M>& subPSeq,
2681  const BaseVectorSequence<Q_V,Q_M>& subQSeq,
2682  unsigned int subNumSamples,
2683  P_M& pqCovMatrix,
2684  P_M& pqCorrMatrix)
2685 {
2686  // Check input data consistency
2687  const BaseEnvironment& env = subPSeq.vectorSpace().zeroVector().env();
2688 
2689  bool useOnlyInter0Comm = (subPSeq.vectorSpace().numOfProcsForStorage() == 1) &&
2690  (subQSeq.vectorSpace().numOfProcsForStorage() == 1);
2691 
2692  UQ_FATAL_TEST_MACRO((useOnlyInter0Comm == false),
2693  env.worldRank(),
2694  "ComputeCovCorrMatricesBetweenVectorSequences()",
2695  "parallel vectors not supported yet");
2696 
2697  unsigned int numRowsLocal = subPSeq.vectorSpace().dimLocal();
2698  unsigned int numCols = subQSeq.vectorSpace().dimGlobal();
2699 
2700  UQ_FATAL_TEST_MACRO((numRowsLocal != pqCovMatrix.numRowsLocal()) || (numCols != pqCovMatrix.numCols()),
2701  env.worldRank(),
2702  "ComputeCovCorrMatricesBetweenVectorSequences()",
2703  "inconsistent dimensions for covariance matrix");
2704 
2705  UQ_FATAL_TEST_MACRO((numRowsLocal != pqCorrMatrix.numRowsLocal()) || (numCols != pqCorrMatrix.numCols()),
2706  env.worldRank(),
2707  "ComputeCorrelationBetweenVectorSequences()",
2708  "inconsistent dimensions for correlation matrix");
2709 
2710  UQ_FATAL_TEST_MACRO((subNumSamples > subPSeq.subSequenceSize()) || (subNumSamples > subQSeq.subSequenceSize()),
2711  env.worldRank(),
2712  "ComputeCovCorrMatricesBetweenVectorSequences()",
2713  "subNumSamples is too large");
2714 
2715  // For both P and Q vector sequences: fill them
2716  P_V tmpP(subPSeq.vectorSpace().zeroVector());
2717  Q_V tmpQ(subQSeq.vectorSpace().zeroVector());
2718 
2719  // For both P and Q vector sequences: compute the unified mean
2720  P_V unifiedMeanP(subPSeq.vectorSpace().zeroVector());
2721  subPSeq.unifiedMeanExtra(0,subNumSamples,unifiedMeanP);
2722 
2723  Q_V unifiedMeanQ(subQSeq.vectorSpace().zeroVector());
2724  subQSeq.unifiedMeanExtra(0,subNumSamples,unifiedMeanQ);
2725 
2726  // Compute "sub" covariance matrix
2727  for (unsigned i = 0; i < numRowsLocal; ++i) {
2728  for (unsigned j = 0; j < numCols; ++j) {
2729  pqCovMatrix(i,j) = 0.;
2730  }
2731  }
2732  for (unsigned k = 0; k < subNumSamples; ++k) {
2733  // For both P and Q vector sequences: get the difference (wrt the unified mean) in them
2734  subPSeq.getPositionValues(k,tmpP);
2735  tmpP -= unifiedMeanP;
2736 
2737  subQSeq.getPositionValues(k,tmpQ);
2738  tmpQ -= unifiedMeanQ;
2739 
2740  for (unsigned i = 0; i < numRowsLocal; ++i) {
2741  for (unsigned j = 0; j < numCols; ++j) {
2742  pqCovMatrix(i,j) += tmpP[i]*tmpQ[j];
2743  }
2744  }
2745  }
2746 
2747  // For both P and Q vector sequences: compute the unified variance
2748  P_V unifiedSampleVarianceP(subPSeq.vectorSpace().zeroVector());
2749  subPSeq.unifiedSampleVarianceExtra(0,
2750  subNumSamples,
2751  unifiedMeanP,
2752  unifiedSampleVarianceP);
2753 
2754  Q_V unifiedSampleVarianceQ(subQSeq.vectorSpace().zeroVector());
2755  subQSeq.unifiedSampleVarianceExtra(0,
2756  subNumSamples,
2757  unifiedMeanQ,
2758  unifiedSampleVarianceQ);
2759 
2760  // Compute unified covariance matrix
2761  if (useOnlyInter0Comm) {
2762  if (env.inter0Rank() >= 0) {
2763  unsigned int unifiedNumSamples = 0;
2764  env.inter0Comm().Allreduce((void *) &subNumSamples, (void *) &unifiedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
2765  "ComputeCovCorrMatricesBetweenVectorSequences()",
2766  "failed MPI.Allreduce() for subNumSamples");
2767 
2768  for (unsigned i = 0; i < numRowsLocal; ++i) {
2769  for (unsigned j = 0; j < numCols; ++j) {
2770  double aux = 0.;
2771  env.inter0Comm().Allreduce((void *) &pqCovMatrix(i,j), (void *) &aux, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
2772  "ComputeCovCorrMatricesBetweenVectorSequences()",
2773  "failed MPI.Allreduce() for a matrix position");
2774  pqCovMatrix(i,j) = aux/((double) (unifiedNumSamples-1)); // Yes, '-1' in order to compensate for the 'N-1' denominator factor in the calculations of sample variances above (whose square roots will be used below)
2775  }
2776  }
2777 
2778  for (unsigned i = 0; i < numRowsLocal; ++i) {
2779  for (unsigned j = 0; j < numCols; ++j) {
2780  pqCorrMatrix(i,j) = pqCovMatrix(i,j)/std::sqrt(unifiedSampleVarianceP[i])/std::sqrt(unifiedSampleVarianceQ[j]);
2781  if (((pqCorrMatrix(i,j) + 1.) < -1.e-8) ||
2782  ((pqCorrMatrix(i,j) - 1.) > 1.e-8)) {
2783  if (env.inter0Rank() == 0) {
2784  std::cerr << "In ComputeCovCorrMatricesBetweenVectorSequences()"
2785  << ": worldRank = " << env.worldRank()
2786  << ", i = " << i
2787  << ", j = " << j
2788  << ", pqCorrMatrix(i,j)+1 = " << pqCorrMatrix(i,j)+1.
2789  << ", pqCorrMatrix(i,j)-1 = " << pqCorrMatrix(i,j)-1.
2790  << std::endl;
2791  }
2792  env.inter0Comm().Barrier();
2793  }
2794  UQ_FATAL_TEST_MACRO(((pqCorrMatrix(i,j) + 1.) < -1.e-8) ||
2795  ((pqCorrMatrix(i,j) - 1.) > 1.e-8),
2796  env.worldRank(),
2797  "ComputeCovCorrMatricesBetweenVectorSequences()",
2798  "computed correlation is out of range");
2799  }
2800  }
2801  }
2802  else {
2803  // Node not in the 'inter0' communicator: do nothing extra
2804  }
2805  }
2806  else {
2807  UQ_FATAL_TEST_MACRO((useOnlyInter0Comm == false),
2808  env.worldRank(),
2809  "ComputeCovCorrMatricesBetweenVectorSequences()",
2810  "parallel vectors not supported yet (2)");
2811  }
2812 
2813  return;
2814 }
2815 
2816 } // End namespace QUESO
2817 
2818 template class QUESO::BaseVectorSequence<QUESO::GslVector, QUESO::GslMatrix>;
2819 template void QUESO::ComputeCovCorrMatricesBetweenVectorSequences<QUESO::GslVector, QUESO::GslMatrix, QUESO::GslVector, QUESO::GslMatrix>(QUESO::BaseVectorSequence<QUESO::GslVector, QUESO::GslMatrix> const&, QUESO::BaseVectorSequence<QUESO::GslVector, QUESO::GslMatrix> const&, unsigned int, QUESO::GslMatrix&, QUESO::GslMatrix&);
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:199
QUESO::BaseVectorSequence::clear
void clear()
Reset the values and the size of the sequence of vectors.
Definition: VectorSequence.C:136
QUESO::BaseVectorSequence::append
void append(const BaseVectorSequence< V, M > &src, unsigned int initialPos, unsigned int numPos)
Appends the vector src to this vector.
Definition: VectorSequence.C:357
QUESO::BaseVectorSequence::deleteStoredVectors
void deleteStoredVectors()
Deletes all the stored vectors.
Definition: VectorSequence.C:301
QUESO::BaseVectorSequence::setGaussian
void setGaussian(const V &meanVec, const V &stdDevVec)
Sets the values of the sequence as a Gaussian distribution of mean given by meanVec and standard devi...
Definition: VectorSequence.C:611
QUESO::BaseVectorSequence::unifiedMeanExtra
virtual void unifiedMeanExtra(unsigned int initialPos, unsigned int numPos, V &unifiedMeanVec) const =0
Finds the mean value of the unified sequence of numPos positions starting at position initialPos...
QUESO::Fft
Class for a Fast Fourier Transform (FFT) algorithm.
Definition: Fft.h:66
QUESO::BaseVectorSequence::unifiedSampleVariancePlain
const V & unifiedSampleVariancePlain() const
Finds the variance of a sample of the unified sequence.
Definition: VectorSequence.C:261
QUESO::BaseEnvironment::inter0Rank
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:289
QUESO::BaseVectorSequence::~BaseVectorSequence
virtual ~BaseVectorSequence()
Destructor.
Definition: VectorSequence.C:60
QUESO::BaseVectorSequence::subMaxPlain
const V & subMaxPlain() const
Finds the maximum value of the sub-sequence.
Definition: VectorSequence.C:175
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
RawValue_MPI_DOUBLE
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:48
QUESO::BaseVectorSequence::unifiedMaxPlain
const V & unifiedMaxPlain() const
Finds the maximum value of the unified sequence.
Definition: VectorSequence.C:188
QUESO::BaseVectorSequence::unifiedPositionsOfMaximum
double unifiedPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
Definition: VectorSequence.C:433
QUESO::MpiComm::Allreduce
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
Definition: MpiComm.C:131
QUESO::BaseVectorSequence::m_vectorSpace
const VectorSpace< V, M > & m_vectorSpace
Definition: VectorSequence.h:424
RawValue_MPI_UNSIGNED
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:49
QUESO::BaseEnvironment
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:187
QUESO::BaseVectorSequence::unifiedMinPlain
const V & unifiedMinPlain() const
Finds the minimum value of the unified sequence.
Definition: VectorSequence.C:162
QUESO::BaseVectorSequence::setUniform
void setUniform(const V &aVec, const V &bVec)
Sets the values of the sequence as a uniform distribution between the values given by vectors aVec an...
Definition: VectorSequence.C:626
QUESO::BaseVectorSequence::computeFilterParams
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
Definition: VectorSequence.C:641
QUESO::BaseVectorSequence::setPositionValues
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.
QUESO::BaseVectorSequence::unifiedSampleVarianceExtra
virtual void unifiedSampleVarianceExtra(unsigned int initialPos, unsigned int numPos, const V &unifiedMeanVec, V &unifiedSamVec) const =0
Finds the sample variance of the unified sequence, considering numPos positions starting at position ...
RawValue_MPI_SUM
#define RawValue_MPI_SUM
Definition: MpiComm.h:52
RawValue_MPI_MAX
#define RawValue_MPI_MAX
Definition: MpiComm.h:51
RawValue_MPI_INT
#define RawValue_MPI_INT
Definition: MpiComm.h:47
QUESO::ComputeCovCorrMatricesBetweenVectorSequences
void ComputeCovCorrMatricesBetweenVectorSequences(const BaseVectorSequence< P_V, P_M > &subPSeq, const BaseVectorSequence< Q_V, Q_M > &subQSeq, unsigned int subNumSamples, P_M &pqCovMatrix, P_M &pqCorrMatrix)
Definition: VectorSequence.C:2679
QUESO::BaseVectorSequence::resizeSequence
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
QUESO::BaseVectorSequence::subMedianPlain
const V & subMedianPlain() const
Finds the median value of the sub-sequence.
Definition: VectorSequence.C:225
QUESO::BaseVectorSequence::name
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
Definition: VectorSequence.C:121
QUESO::MpiComm::Barrier
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
Definition: MpiComm.C:143
QUESO::BaseVectorSequence::unifiedBoxPlain
const BoxSubset< V, M > & unifiedBoxPlain() const
Finds a box subset of the unified-sequence (given by the min and max values of the unified sequence c...
Definition: VectorSequence.C:287
QUESO::ScalarSequence::subMaxPlain
const T & subMaxPlain() const
Finds the maximum value of the sub-sequence of scalars.
Definition: ScalarSequence.C:390
QUESO::BaseVectorSequence::unifiedMedianPlain
const V & unifiedMedianPlain() const
Finds the median value of the unified sequence.
Definition: VectorSequence.C:237
QUESO::VectorSpace::dimGlobal
unsigned int dimGlobal() const
Definition: VectorSpace.C:205
QUESO::BaseVectorSequence
Base class for handling vector and array samples (sequence of vectors or arrays). ...
Definition: VectorSequence.h:54
QUESO::BaseVectorSequence::subBoxPlain
const BoxSubset< V, M > & subBoxPlain() const
Finds a box subset of the sub-sequence (given by its min and max values calculated via subMinPlain an...
Definition: VectorSequence.C:273
QUESO::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
QUESO::BaseVectorSequence::subPositionsOfMaximum
double subPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
Definition: VectorSequence.C:386
QUESO::BaseVectorSequence::vectorSizeLocal
unsigned int vectorSizeLocal() const
Local dimension (size) of the vector space.
Definition: VectorSequence.C:100
QUESO::BaseVectorSequence::subSampleVariancePlain
const V & subSampleVariancePlain() const
Finds the variance of a sample of the sub-sequence.
Definition: VectorSequence.C:249
QUESO::BaseVectorSequence::vectorSizeGlobal
unsigned int vectorSizeGlobal() const
Global dimension (size) of the vector space.
Definition: VectorSequence.C:107
QUESO::BoxSubset
Class representing a subset of a vector space shaped like a hypercube.
Definition: BoxSubset.h:41
QUESO::BaseVectorSequence::vectorSpace
const VectorSpace< V, M > & vectorSpace() const
Vector space; access to protected attribute VectorSpace&lt;V,M&gt;&amp; m_vectorSpace.
Definition: VectorSequence.C:114
QUESO::VectorSpace
A class representing a vector space.
Definition: VectorSet.h:46
QUESO::BaseVectorSequence::copy
void copy(const BaseVectorSequence< V, M > &src)
Copies vector sequence src to this.
Definition: VectorSequence.C:681
QUESO::BaseVectorSequence::BaseVectorSequence
BaseVectorSequence(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
Definition: VectorSequence.C:34
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:218
QUESO::ScalarSequence::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Definition: ScalarSequence.C:141
QUESO::BaseVectorSequence::subMinPlain
const V & subMinPlain() const
Finds the minimum value of the sub-sequence.
Definition: VectorSequence.C:149
QUESO::BaseVectorSequence::subMeanPlain
const V & subMeanPlain() const
Finds the mean value of the sub-sequence.
Definition: VectorSequence.C:201
QUESO::VectorSpace::numOfProcsForStorage
unsigned int numOfProcsForStorage() const
Returns total number of processes.
Definition: VectorSpace.C:193
QUESO::BaseEnvironment::inter0Comm
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:295
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:76
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:223
QUESO::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:287
QUESO::BaseVectorSequence::unifiedMeanPlain
const V & unifiedMeanPlain() const
Finds the mean value of the unified sequence.
Definition: VectorSequence.C:213
QUESO::BaseVectorSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize() const
Calculates the size of the unified sequence of vectors.
Definition: VectorSequence.C:70
QUESO::GslMatrix
Class for matrix operations using GSL library.
Definition: GslMatrix.h:47
QUESO::BaseVectorSequence::setName
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
Definition: VectorSequence.C:128
QUESO::BaseVectorSequence::getPositionValues
virtual void getPositionValues(unsigned int posId, V &vec) const =0
Gets the values of the sequence at position posId and stores them at vec. See template specialization...
QUESO::MiscComputePositionsBetweenMinMax
void MiscComputePositionsBetweenMinMax(V minValues, V maxValues, std::vector< V * > &positions)
Definition: Miscellaneous.C:545
Generated on Thu Apr 23 2015 19:26:15 for queso-0.51.1 by  doxygen 1.8.5