queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
Miscellaneous.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-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <cstring>
26 #include <queso/Defines.h>
27 #include <queso/Miscellaneous.h>
28 #include <queso/GslVector.h>
29 #include <queso/GslMatrix.h>
30 #include <sys/time.h>
31 #include <iostream>
32 #include <fstream>
33 #include <libgen.h>
34 #include <sys/stat.h>
35 #include <cmath>
36 
37 namespace QUESO {
38 
39 void
41  const std::string& inputString,
42  std::vector<double>& outputDoubles)
43 {
44  //std::cout << "In MiscReadDoublesFromString()"
45  // << ": inputString = " << inputString
46  // << std::endl;
47  outputDoubles.clear();
48  bool aDoubleIsBeingRead = false;
49  std::string::size_type positionOfFirstChar = 0;
50  std::string::size_type numberOfChars = 0;
51  for (std::string::size_type i = 0; i < inputString.size(); ++i) {
52  queso_require_not_equal_to_msg(inputString[i], '\0', "character '\0' should not be found!");
53  if (inputString[i] == ' ') {
54  if (aDoubleIsBeingRead == true) {
55  // We just finished reading the current string/double. Convert string to double now.
56  char tmpVar[numberOfChars+1];
57  for (std::string::size_type j = 0; j < numberOfChars; ++j) {
58  tmpVar[j] = inputString[positionOfFirstChar+j];
59  }
60  tmpVar[numberOfChars] = '\0';
61  outputDoubles.push_back(strtod(tmpVar,NULL));
62 
63  // Continue loop
64  aDoubleIsBeingRead = false;
65  positionOfFirstChar = 0;
66  numberOfChars = 0;
67  }
68  }
69  else {
70  if (aDoubleIsBeingRead == false) {
71  aDoubleIsBeingRead = true;
72  positionOfFirstChar = i;
73  }
74  numberOfChars++;
75  }
76  } // for
77  if (aDoubleIsBeingRead == true) {
78  // We just finished reading the current string/double. Convert string to double now.
79  char tmpVar[numberOfChars+1];
80  for (std::string::size_type j = 0; j < numberOfChars; ++j) {
81  tmpVar[j] = inputString[positionOfFirstChar+j];
82  }
83  tmpVar[numberOfChars] = '\0';
84  outputDoubles.push_back(strtod(tmpVar,NULL));
85  }
86  std::vector<double>(outputDoubles).swap(outputDoubles);
87 
88  return;
89 }
90 
91 void
93  const std::string& inputString,
94  std::vector<std::string>& outputWords)
95 {
96  //std::cout << "In MiscReadWordsFromString()"
97  // << ": inputString = " << inputString
98  // << std::endl;
99  outputWords.clear();
100  bool aWordIsBeingRead = false;
101  std::string::size_type positionOfFirstChar = 0;
102  std::string::size_type numberOfChars = 0;
103  for (std::string::size_type i = 0; i < inputString.size(); ++i) {
104  queso_require_not_equal_to_msg(inputString[i], '\0', "character '\0' should not be found!");
105  if (inputString[i] == ' ') {
106  if (aWordIsBeingRead == true) {
107  // We just finished reading the current string/word.
108  char tmpVar[numberOfChars+1];
109  for (std::string::size_type j = 0; j < numberOfChars; ++j) {
110  tmpVar[j] = inputString[positionOfFirstChar+j];
111  }
112  tmpVar[numberOfChars] = '\0';
113  outputWords.push_back(tmpVar);
114 
115  // Continue loop
116  aWordIsBeingRead = false;
117  positionOfFirstChar = 0;
118  numberOfChars = 0;
119  }
120  }
121  else {
122  if (aWordIsBeingRead == false) {
123  aWordIsBeingRead = true;
124  positionOfFirstChar = i;
125  }
126  numberOfChars++;
127  }
128  } // for
129  if (aWordIsBeingRead == true) {
130  // We just finished reading the current string/word.
131  char tmpVar[numberOfChars+1];
132  for (std::string::size_type j = 0; j < numberOfChars; ++j) {
133  tmpVar[j] = inputString[positionOfFirstChar+j];
134  }
135  tmpVar[numberOfChars] = '\0';
136  outputWords.push_back(tmpVar);
137  }
138  std::vector<std::string>(outputWords).swap(outputWords);
139 
140  return;
141 }
142 
143 //void
144 //MiscExtractDoubleFromString(
145 // std::string& inputString,
146 // double& outputDouble)
147 //{
148 // return;
149 //}
150 
151 //void
152 //MiscExtractWordFromString(
153 // std::string& inputString,
154 // std::string& outputWord)
155 //{
156 // return;
157 //}
158 
159 int
161  std::ifstream& ifs,
162  std::string& termString,
163  double* termValue)
164 {
165  int iRC = UQ_OK_RC;
166 
167  ifs >> termString;
168  if ((ifs.rdstate() & std::ifstream::failbit)) {
170  }
171  else if (termValue) {
172  if (termString == std::string("inf")) {
173  *termValue = INFINITY;
174  }
175  else if (termString == std::string("-inf")) {
176  *termValue = -INFINITY;
177  }
178  else if (termString == std::string("nan")) {
179  *termValue = nan("");
180  }
181  else {
182  *termValue = strtod(termString.c_str(),NULL);
183  }
184  }
185  //if (!iRC) std::cout << "Read termString = " << termString << std::endl;
186 
187  return iRC;
188 }
189 
190 int
192  std::ifstream& ifs,
193  std::string& termString,
194  double* termValue,
195  bool& endOfLineAchieved)
196 {
197  int iRC = UQ_OK_RC;
198  endOfLineAchieved = false;
199 
200  char c = ' ';
201  while (c == ' ') {
202  ifs.get(c);
203  if ((ifs.rdstate() & std::ifstream::failbit)) {
205  break;
206  }
207  };
208 
209  char term[512];
210  unsigned int pos = 0;
211 
212  if (!iRC) {
213  while ((pos < 512) && (c != '\n') && (c != '\0') && (c != ' ')) {
214  term[pos++] = c;
215  if ((ifs.rdstate() & std::ifstream::failbit)) {
217  break;
218  }
219  ifs.get(c);
220  };
221  }
222 
223  if (!iRC) {
224  if (c == '\n') endOfLineAchieved = true;
225  term[pos] = '\0';
226  termString = term;
227  //std::cout << "Read chars = " << termString << std::endl;
228  if (termValue) {
229  if (termString == std::string("inf")) {
230  *termValue = INFINITY;
231  }
232  else if (termString == std::string("-inf")) {
233  *termValue = -INFINITY;
234  }
235  else if (termString == std::string("nan")) {
236  *termValue = nan("");
237  }
238  else {
239  *termValue = strtod(termString.c_str(),NULL);
240  }
241  }
242  }
243 
244  return iRC;
245 }
246 
247 double
249  double a,
250  double b,
251  const RngBase* rngObject)
252 {
253  double result = 0.;
254  if (a < 1.) {
255  result = MiscGammar(1.+a,b,rngObject)*std::pow( rngObject->uniformSample(),1./a );
256  }
257  else {
258  double d = a-1./3.;
259  double c = 1./std::sqrt(9.*d);
260  double x = 0.;
261  double w = 0.;
262  while (1) {
263  while (1) {
264  x = rngObject->gaussianSample(1.);
265  w = 1.+c*x;
266  if (w > 0.) break;
267  }
268  w = std::pow(w,3.);
269  double u = rngObject->uniformSample();
270  double compValue = 1.-0.0331*std::pow(x,4.);
271  if (u < compValue) break;
272  compValue = 0.5*std::pow(x,2.)+d*(1.-w+log(w));
273  if (log(u) < compValue) break;
274  }
275  result = b*d*w;
276  }
277 
278  return result;
279 }
280 
281 double
282 MiscGetEllapsedSeconds(struct timeval *timeval0)
283 {
284  double result = 0.;
285 
286  struct timeval timevalNow;
287  /*int iRC;*/
288  /*iRC = */gettimeofday(&timevalNow, NULL);
289 
290  result = (double) (timevalNow.tv_sec - timeval0->tv_sec );
291  result *= 1.e+6;
292  result += (double) (timevalNow.tv_usec - timeval0->tv_usec);
293  result *= 1.e-6;
294 
295  return result;
296 }
297 
298 double MiscHammingWindow(unsigned int N, unsigned int j)
299 {
300  double angle = 2.*M_PI*((double) j)/((double) N);
301  double result = 0.53836 - 0.46164*cos(angle);
302 
303  return result;
304 }
305 
306 double MiscGaussianDensity(double x, double mu, double sigma)
307 {
308  double sigma2 = sigma*sigma;
309  double diff = x-mu;
310 
311  return (1./std::sqrt(2*M_PI*sigma2))*std::exp(-.5*diff*diff/sigma2);
312 }
313 
314 unsigned int MiscUintDebugMessage(
315  unsigned int value,
316  const char* message)
317 {
318  if (message) {
319  std::cout << "Passing in MiscUintDebugMessage(), value = " << value << ", message = " << message << std::endl;
320  }
321  return value;
322 }
323 
325  int value,
326  const char* message)
327 {
328  if (message) {
329  std::cout << "Passing in MiscIntDebugMessage(), value = " << value << ", message = " << message << std::endl;
330  }
331  return value;
332 }
333 
335  double value,
336  const char* message)
337 {
338  if (message) {
339  std::cout << "Passing in MiscDoubleDebugMessage(), value = " << value << ", message = " << message << std::endl;
340  }
341  return value;
342 }
343 
357 
358 // ------------------------------------------
359 // Following routines borrowed from libGRVY
360 // ------------------------------------------
361 
362 int CheckFilePath(const char *pathname)
363  {
364 
365  // verify parent directories in path exist (and create if not).
366 
367 #ifdef HAVE_GRVY
368  return(grvy_check_file_path(pathname));
369 #else
370 
371  const int MAX_DEPTH = 50;
372 
373  char *pathlocal;
374  char *parents;
375  char *dirstring;
376  char *token;
377  int depth = 0;
378 
379  // Save a copy of pathname and look for the parent directories.
380 
381  pathlocal = strdup(pathname);
382  dirstring = strdup(pathname);
383  parents = dirname(pathlocal);
384 
385  if(strcmp(parents,".") == 0)
386  {
387  free(pathlocal);
388  free(dirstring);
389  return 0;
390  }
391 
392  // Deal with the possibility of an absolute path being provided
393 
394  bool abs_path = false;
395 
396  std::string leading_char("");
397  std::string path_to_check;
398 
399  if(strncmp(parents,"/",1) == 0)
400  {
401  leading_char = "/";
402  abs_path = true;
403  }
404 
405  // Verify existence of top-level directory
406 
407  if( (token = strtok(parents,"/")) != NULL )
408  {
409  path_to_check += leading_char + token;
410 
411  if ( GRVY_CheckDir(path_to_check.c_str()) )
412  {
413  free(pathlocal);
414  free(dirstring);
415  return -1;
416  }
417 
418  // Now, search for any remaining parent directories.
419 
420  if(abs_path)
421  sprintf(dirstring,"/%s",token);
422  else
423  sprintf(dirstring,"%s",token);
424 
425  while ( (token = strtok(0,"/")) && (depth < MAX_DEPTH) )
426  {
427  dirstring = strcat(dirstring,"/");
428 
429  if(GRVY_CheckDir(strcat(dirstring,token)))
430  {
431  free(pathlocal);
432  free(dirstring);
433  return -1;
434  }
435  depth++;
436  };
437 
438  if(depth >= MAX_DEPTH )
439  {
440  std::cerr << __func__ << ": error - Max directory depth exceeded, limit = " << MAX_DEPTH << std::endl;
441  free(pathlocal);
442  free(dirstring);
443  return -1;
444  }
445  }
446 
447  // Clean Up
448  free(pathlocal);
449  free(dirstring);
450 
451  return 0;
452 #endif
453  }
454 
455 
456 int GRVY_CheckDir(const char *dirname)
457 {
458  struct stat st;
459 
460  if(stat(dirname,&st) != 0)
461  {
462  if( mkdir(dirname,0700) != 0 )
463  {
464  std::cerr << __func__ << ": error - unable to create directory " << dirname << std::endl;
465  return -1;
466  }
467  }
468  else if (!S_ISDIR(st.st_mode))
469  {
470  std::cerr << __func__ << ": error - entry exists, but is not a directory " << dirname << std::endl;
471  return -1;
472  }
473 
474  return 0;
475 }
476 
477 template <class T>
478 bool
479 MiscCheckForSameValueInAllNodes(T& inputValue, // Yes, 'not' const
480  double acceptableTreshold,
481  const MpiComm& comm,
482  const char* whereString)
483 {
484  // Filter out those nodes that should not participate
485  if (comm.MyPID() < 0) return true;
486 
487  double localValue = (double) inputValue;
488  double sumValue = 0.;
489  comm.Allreduce<double>(&localValue, &sumValue, (int) 1, RawValue_MPI_SUM,
490  whereString,
491  "failed MPI on 'sumValue' inside MiscCheckForSameValueInAllNodes()");
492 
493  double totalNumNodes = (double) comm.NumProc();
494  double testValue = fabs(1. - localValue/(sumValue/totalNumNodes));
495  unsigned int boolSum = 0;
496 #if 1
497  unsigned int boolResult = 0;
498  if (testValue > acceptableTreshold) boolResult = 1;
499  comm.Allreduce<unsigned int>(&boolResult, &boolSum, (int) 1, RawValue_MPI_SUM,
500  whereString,
501  "failed MPI on 'boolSum' inside MiscCheckForSameValueInAllNodes()");
502 
503  if (boolSum > 0) {
504  comm.Barrier();
505  for (int i = 0; i < comm.NumProc(); ++i) {
506  if (i == comm.MyPID()) {
507  std::cerr << "WARNING, "
508  << whereString
509  << ", inside MiscCheckForSameValueInAllNodes()"
510  << ", rank (in this communicator) = " << i
511  << ": boolSum = " << boolSum
512  << ", localValue = " << localValue
513  << ", sumValue = " << sumValue
514  << ", totalNumNodes = " << totalNumNodes
515  << ", avgValue = " << (sumValue/totalNumNodes)
516  << ", relativeTest = " << testValue
517  << std::endl;
518  }
519  comm.Barrier();
520  }
521  comm.Barrier();
522 
523  comm.Bcast((void *) &localValue, (int) 1, RawValue_MPI_DOUBLE, 0,
524  whereString,
525  "failed MPI on 'boolSum' inside MiscCheckForSameValueInAllNodes()");
526  inputValue = localValue; // IMPORTANT
527  }
528 #else
529  queso_require_less_equal_msg(testValue, acceptableTreshold, "not all nodes have the same value inside MiscCheckForSameValueInAllNodes()");
530 #endif
531 
532  return (boolSum == 0);
533 }
534 
535 template <class V>
536 void
538  V maxValues,
539  std::vector<V*>& positions)
540 {
541  double factor = 0.5;
542  switch (positions.size()) {
543  case 0:
544  // Do nothing
545  break;
546 
547  case 1:
548  positions[0] = new V((1. - factor) * minValues + factor * maxValues);
549  break;
550 
551  default:
552  for (unsigned int i = 0; i < positions.size(); ++i) {
553  factor = ((double) i)/(((double) positions.size()) - 1.);
554  positions[i] = new V((1. - factor) * minValues + factor * maxValues);
555  }
556  break;
557  }
558 
559  return;
560 }
561 
562 template <class V1,class V2>
563 void
564 MiscCheckTheParallelEnvironment(const V1& vec1, const V2& vec2)
565 {
566  const BaseEnvironment& env = vec1.env();
567 
568  if (env.numSubEnvironments() == (unsigned int) env.fullComm().NumProc()) {
569  queso_require_equal_to_msg(env.subRank(), 0, "there should exist only one processor per sub environment");
570  queso_require_equal_to_msg(vec1.numOfProcsForStorage(), 1,
571  "only 1 processor (per sub environment) should be necessary for the storage of a parameter vector");
572  queso_require_equal_to_msg(vec2.numOfProcsForStorage(), 1,
573  "only 1 processor (per sub environment) should be necessary for the storage of a parameter vector");
574  }
575  else if (env.numSubEnvironments() < (unsigned int) env.fullComm().NumProc()) {
576  queso_require_equal_to_msg(env.fullComm().NumProc()%env.numSubEnvironments(), 0, "total number of processors should be a multiple of the number of sub environments");
577  unsigned int numProcsPerSubEnvironment = env.fullComm().NumProc()/env.numSubEnvironments();
578  queso_require_equal_to_msg(env.subComm().NumProc(), (int) numProcsPerSubEnvironment, "inconsistent number of processors per sub environment");
579  if ((vec1.numOfProcsForStorage() == 1) &&
580  (vec2.numOfProcsForStorage() == 1)) {
581  // Ok
582  }
583  else if ((vec1.numOfProcsForStorage() == numProcsPerSubEnvironment) &&
584  (vec2.numOfProcsForStorage() == numProcsPerSubEnvironment)) {
585  queso_error_msg("parallel vectors are not supported yet");
586  }
587  else {
588  queso_error_msg("number of processors required for a vector storage should be equal to either 1 or to the number of processors in the sub environment");
589  }
590  }
591  else {
592  queso_error_msg("number of processors per sub environment is less than 1!");
593  }
594 
595  return;
596 }
597 
598 } // End namespace QUESO
599 
600 template void QUESO::MiscCheckTheParallelEnvironment<QUESO::GslVector, QUESO::GslVector>(QUESO::GslVector const&, QUESO::GslVector const&);
601 template bool QUESO::MiscCheckForSameValueInAllNodes<bool>(bool&, double, QUESO::MpiComm const&, char const*);
602 template bool QUESO::MiscCheckForSameValueInAllNodes<double>(double&, double, QUESO::MpiComm const&, char const*);
virtual double gaussianSample(double stdDev) const =0
Samples a value from a Gaussian distribution with standard deviation given by stdDev.
int MiscReadStringAndDoubleFromFile(std::ifstream &ifs, std::string &termString, double *termValue)
int GRVY_CheckDir(const char *dirname)
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
Definition: MpiComm.C:174
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:198
const int UQ_OK_RC
Definition: Defines.h:92
double MiscHammingWindow(unsigned int N, unsigned int j)
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix queso_require_not_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the absence of an options input file"))
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
double MiscGammar(double a, double b, const RngBase *rngObject)
int MyPID() const
Return my process ID.
Definition: MpiComm.C:124
Class for random number generation (base class for either GSL or Boost RNG).
Definition: RngBase.h:45
double MiscDoubleDebugMessage(double value, const char *message)
int CheckFilePath(const char *path)
double MiscGaussianDensity(double x, double mu, double sigma)
int MiscIntDebugMessage(int value, const char *message)
const int UQ_FAILED_READING_FILE_RC
Definition: Defines.h:101
const MpiComm & subComm() const
Access function for each sub-environment&#39;s communicator.
Definition: Environment.C:293
unsigned int MiscUintDebugMessage(unsigned int value, const char *message)
void MiscReadDoublesFromString(const std::string &inputString, std::vector< double > &outputDoubles)
Definition: Miscellaneous.C:40
Class for vector operations using GSL library.
Definition: GslVector.h:49
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
Definition: MpiComm.C:143
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"))
void MiscComputePositionsBetweenMinMax(V minValues, V maxValues, std::vector< V * > &positions)
The QUESO MPI Communicator Class.
Definition: MpiComm.h:203
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:133
bool MiscCheckForSameValueInAllNodes(T &inputValue, double acceptableTreshold, const MpiComm &comm, const char *whereString)
const MpiComm & fullComm() const
Access function for the communicator that was passed to QUESO&#39;s environment.
Definition: Environment.C:274
int MiscReadCharsAndDoubleFromFile(std::ifstream &ifs, std::string &termString, double *termValue, bool &endOfLineAchieved)
void MiscReadWordsFromString(const std::string &inputString, std::vector< std::string > &outputWords)
Definition: Miscellaneous.C:92
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:335
double MiscGetEllapsedSeconds(struct timeval *timeval0)
void MiscCheckTheParallelEnvironment(const V1 &vec1, const V2 &vec2)
void Bcast(void *buffer, int count, RawType_MPI_Datatype datatype, int root, const char *whereMsg, const char *whatMsg) const
Broadcast values from the root process to the slave processes.
Definition: MpiComm.C:191

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