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

Generated on Tue Jun 5 2018 19:48:54 for queso-0.57.1 by  doxygen 1.8.5