1 /* _______________________________________________________________________
2
3 PECOS: Parallel Environment for Creation Of Stochastics
4 Copyright (c) 2008, Sandia National Laboratories.
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Pecos directory.
7 _______________________________________________________________________ */
8
9 /** \file pecos_gsg_driver.cpp
10 \brief A driver program for PECOS */
11
12 #include <iostream>
13 #include <fstream>
14 #include <sstream>
15 #include <vector>
16
17 #include "IncrementalSparseGridDriver.hpp"
18 #include "SharedProjectOrthogPolyApproxData.hpp"
19 #include "ProjectOrthogPolyApproximation.hpp"
20 #include "TensorProductDriver.hpp"
21 #include "CubatureDriver.hpp"
22 #include "pecos_data_types.hpp"
23 #include "Teuchos_SerialDenseHelpers.hpp"
24
25 #include "TestFunctions.hpp"
26
27 using namespace std;
28
29 #define POLYTYPE LEGENDRE_ORTHOG
30 #define QUADRULE GAUSS_PATTERSON
31 #define MAX_CHARS_PER_LINE 1000
32 #define BTYPE GLOBAL_PROJECTION_ORTHOGONAL_POLYNOMIAL
33 #define NUMVARS 2
34 #define STARTLEV 1
35 #define NITER 0
36 #define MAXORD 5
37 #define NQOI 1
38 #define VERBOSE 1
39 #define FCNTYPE "gerstner-iso1"
40 #define VARTHRLD 1.e-2
41
42 #ifdef GSGREST
fileExist(const char * fname)43 inline bool fileExist (const char *fname) {
44 std::ifstream in(fname);
45 return in.good();
46 }
47
48 void saveData(const char *fbase, const RealMatrix &evalGrid, const IntVector &evalFlag, const RealMatrix &evalData) ;
49 void loadData(const char *fbase, RealMatrix &evalGrid, IntVector &evalFlag, RealMatrix &evalData) ;
50 RealMatrix getFeval(RealMatrix &evalGrid, IntVector &evalFlag,
51 RealMatrix &evalData, const RealMatrix &var_sets,
52 bool &foundFlag);
53 #else
54 RealMatrix feval(const RealMatrix &dataMat, const int nQoI, String ftype);
55 #endif
56
usage()57 int usage(){
58 printf("usage: pecos_gsg_driver [options]\n");
59 printf(" -h : print out this help message \n");
60 printf(" -d <nvar> : dimensionality of parameter space (default=%d) \n",NUMVARS);
61 printf(" -e <veps> : tolerance for incremental variance threshold (default=%lg) \n",VARTHRLD);
62 printf(" -i <niter> : maximum no. of iterations (default=%d) \n",NITER);
63 printf(" -l <stlev> : starting quadrature level (default=%d) \n",STARTLEV);
64 printf(" -m <mord> : not used - maximum order for shared-poly-data (default=%d) \n",MAXORD);
65 printf(" -n <nqoi> : no. of outputs (default=%d) \n",NQOI);
66 printf(" -p <ptype> : type of basis (default=LEGENDRE) \n");
67 printf(" -t <ftype> : test function (default=%s) \n",FCNTYPE);
68 printf(" -v <verb> : verbosity level (default=%d) \n",VERBOSE);
69 exit(0);
70 return (0);
71 }
72
73 /// A driver program for PECOS.
main(int argc,char * argv[])74 int main(int argc, char* argv[])
75 {
76
77 using namespace Pecos;
78
79 // Define defaults
80 int polyType = POLYTYPE ;
81 int quadRule = QUADRULE ;
82 int nQoI = NQOI ; /* no. of Quantities of Interest (model outputs) */
83 size_t nvar = NUMVARS ; /* dimensionality of input space */
84 unsigned short mOrd = MAXORD ; /* maximum order */
85 unsigned short nIter = NITER ; /* maximum no. of iterations */
86 unsigned short strtlev = STARTLEV ; /* starting quadrature level */
87 unsigned short verb = VERBOSE ; /* verbosity */
88 double varEps = VARTHRLD ;
89 String ftype = String(FCNTYPE);
90 bool extFunc = false ;
91 short btype = (short) BTYPE;
92
93 String pstring, qstring;
94 // Command-line arguments: read user input
95 int c;
96 while ((c=getopt(argc,(char **)argv,"hfd:e:i:l:m:n:p:t:v:"))!=-1){
97 switch (c) {
98 case 'h':
99 usage();
100 break;
101 case 'f':
102 extFunc = true;
103 break;
104 case 'd':
105 nvar = strtol(optarg, (char **)NULL,0);
106 break;
107 case 'e':
108 varEps = strtod(optarg, (char **)NULL);
109 break;
110 case 'i':
111 nIter = strtol(optarg, (char **)NULL,0);
112 break;
113 case 'l':
114 strtlev = strtol(optarg, (char **)NULL,0);
115 break;
116 case 'm':
117 mOrd = strtol(optarg, (char **)NULL,0);
118 break;
119 case 'n':
120 nQoI = strtol(optarg, (char **)NULL,0);
121 break;
122 case 'p':
123 pstring = String(optarg);
124 break;
125 case 't':
126 ftype = String(optarg);
127 break;
128 case 'v':
129 verb = strtol(optarg, (char **)NULL,0);
130 break;
131 default :
132 break;
133 }
134 }
135
136 // Retrieve command-line setup
137 if (pstring == String("LEGENDRE")) {
138 polyType = LEGENDRE_ORTHOG ;
139 quadRule = GAUSS_PATTERSON ;
140 //quadRule = GAUSS_LEGENDRE ;
141 }
142 else if (pstring == String("HERMITE")) {
143 polyType = HERMITE_ORTHOG ;
144 quadRule = GENZ_KEISTER ;
145 }
146
147 if (verb > 2)
148 PCout << "Instantiating IncrementalSparseGridDriver:\n";
149
150 RealVector dimension_pref; // empty -> isotropic
151 short growth_rate = UNRESTRICTED_GROWTH;
152 short refine_cntl = DIMENSION_ADAPTIVE_CONTROL_GENERALIZED;
153
154 // Start
155 // Can either use IntegrationDriver(driver_type) and then assign data or
156 // use IntegrationDriver() and then assign letter.
157 IntegrationDriver int_driver; // empty envelope
158 // assign letter using assign_rep()
159 std::shared_ptr<IncrementalSparseGridDriver> csg_driver =
160 std::make_shared<IncrementalSparseGridDriver>
161 (strtlev, dimension_pref, growth_rate, refine_cntl);
162 int_driver.assign_rep(csg_driver);
163
164 if (verb > 2)
165 PCout << "Instantiating basis...\n";
166
167 std::vector<BasisPolynomial> poly_basis(nvar); // array of envelopes
168 for (int i=0; i<nvar; ++i) {
169 poly_basis[i] = BasisPolynomial(polyType); // envelope operator=
170 poly_basis[i].collocation_rule(quadRule);
171 }
172 csg_driver->initialize_grid(poly_basis);
173
174 if (verb > 2)
175 PCout << " - done\n";
176
177 // Instantiate Pecos Objects
178 if (verb > 2)
179 PCout << "Instantiating pecos objects...\n";
180 ExpansionConfigOptions expcfgopt(INCREMENTAL_SPARSE_GRID, // expsolnapp
181 DEFAULT_BASIS, // expbassus
182 NO_COMBINE, // exp combine type
183 NO_DISCREP, // discrepancy type
184 SILENT_OUTPUT, // output level
185 true, // vbd flag
186 2, // vbd order
187 refine_cntl, // refinement control
188 COVARIANCE_METRIC, // refine metric
189 ACTIVE_EXPANSION_STATS,// refine stats type
190 100, // max refine iter
191 100, // max solver iter
192 1.e-5, // conv tol
193 2); // soft conv limit
194 BasisConfigOptions bcopt;
195 UShortArray aord(mOrd,nvar);
196 SharedBasisApproxData shared_data; // Envelope
197 std::shared_ptr<SharedProjectOrthogPolyApproxData> shared_poly_data =
198 std::make_shared<SharedProjectOrthogPolyApproxData>
199 (BTYPE,aord,nvar,expcfgopt,bcopt); // Letter
200 shared_data.assign_rep(shared_poly_data); // don't increment ref count
201 shared_poly_data->integration_driver_rep(csg_driver);
202 shared_poly_data->polynomial_basis(poly_basis);
203
204 // Instantiate Project poly approx
205 std::vector<BasisApproximation> poly_approx(nQoI); // array of envelopes
206 for ( int iQoI=0; iQoI<nQoI; iQoI++)
207 poly_approx[iQoI].assign_rep
208 (std::make_shared<ProjectOrthogPolyApproximation>(shared_data));
209
210 if (verb > 2)
211 PCout << " - done\n";
212
213 #ifdef GSGREST
214 /* Define saved data*/
215 RealMatrix evalGrid;
216 IntVector evalFlag;
217 RealMatrix evalData(1,nQoI);
218 loadData((char *)"func",evalGrid,evalFlag,evalData);
219 PCout<<evalGrid;
220 PCout<<evalFlag;
221 PCout<<evalData;
222 #endif
223
224 // initial grid and compute reference approximation
225 RealMatrix var_sets;
226 csg_driver->compute_grid(var_sets);
227 int numPts = var_sets.numCols();
228 assert(nvar==var_sets.numRows());
229 if (verb > 1) {
230 PCout<<var_sets<<endl;
231 if (verb > 2) {
232 PCout << "Evaluate function on reference grid, ";
233 PCout << "instantiate SurrogateData and compute coefficients ...\n";
234 }
235 }
236
237 #ifdef GSGREST
238 bool foundFlag ;
239 RealMatrix fev = getFeval(evalGrid,evalFlag,evalData,var_sets,foundFlag);
240 if ( !foundFlag ) {
241 saveData((char *)"func",evalGrid,evalFlag,evalData);
242 return (0);
243 }
244 #else
245 RealMatrix fev = feval(var_sets,nQoI,ftype);
246 #endif
247
248 // Create SurrogateData instances and assign to
249 // ProjectOrthogPolyApproximation instances
250 bool handle = true;
251 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
252 SurrogateData sdi(handle);
253 for( int jCol = 0; jCol < numPts; jCol++) {
254 SurrogateDataVars sdv(nvar,0,0);
255 SurrogateDataResp sdr(1,nvar); // no gradient or hessian
256 sdv.continuous_variables(Teuchos::getCol<int,double>(Teuchos::Copy,var_sets,jCol));
257 sdr.response_function(fev(jCol,iQoI));
258 sdi.push_back(sdv,sdr);
259 }
260 poly_approx[iQoI].surrogate_data(sdi);
261 }
262
263 shared_poly_data->allocate_data();
264 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
265 poly_approx[iQoI].compute_coefficients();
266 poly_approx[iQoI].print_coefficients(PCout,false);
267 }
268
269 if (verb > 2)
270 PCout << " - done\n";
271
272 // start refinement
273 csg_driver->initialize_sets();
274 UShortArraySet a;
275 bool conv_within_tol = false;
276 for ( size_t iter = 0; iter < nIter; iter++) {
277
278 /* Compute base variance */
279 RealVector respVariance(nQoI,0.0);
280 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
281 auto poly_approx_rep = std::static_pointer_cast<PolynomialApproximation>
282 (poly_approx[iQoI].approx_rep());
283 respVariance[iQoI] = poly_approx_rep->variance() ;
284 }
285 Real deltaVar = 0.0;
286 std::vector<short unsigned int> asave;
287
288 a = csg_driver->active_multi_index();
289 if (verb > 1)
290 PCout<<"Refine, iteration: "<<iter+1
291 <<"\n ... starting variance:\n"<<respVariance<<'\n';
292
293 #ifdef GSGREST
294 /* Iterate through all proposed sets and save/exit if some vals
295 are not available */
296 bool foundAllSets = true ;
297 for (UShortArraySet::iterator it=a.begin(); it!=a.end(); ++it) {
298 csg_driver->push_trial_set(*it);
299 if (!shared_poly_data->push_available()) {
300 csg_driver->compute_trial_grid(var_sets);
301 fev = getFeval(evalGrid,evalFlag,evalData,var_sets,foundFlag);
302 if ( !foundFlag ) foundAllSets = false ;
303 } else {
304 csg_driver->push_set();
305 }
306 csg_driver->pop_set();
307 }
308 if ( !foundAllSets ) {
309 saveData((char *)"func",evalGrid,evalFlag,evalData);
310 return (0);
311 }
312 #endif
313
314 int choose = 0;
315 for (UShortArraySet::iterator it=a.begin(); it!=a.end(); ++it) {
316
317 csg_driver->increment_smolyak_multi_index(*it);
318
319 // Update surrogate data
320 numPts = 0;
321 if (shared_poly_data->push_available()) {
322
323 if (verb > 1)
324 PCout<<"Restoring existing index set:\n"<<*it<<endl;
325
326 // Set available -> restore in csg and the rest
327 csg_driver->push_set();
328
329 size_t idxRestore = shared_poly_data->restore_index();
330 shared_poly_data->pre_push_data();
331 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
332 poly_approx[iQoI].push_coefficients();
333 // Also restore the corresponding surrogate data
334 poly_approx[iQoI].surrogate_data().push(idxRestore,true);
335 }
336 shared_poly_data->post_push_data();
337
338 }
339 else {
340
341 // New set -> compute
342 // Create SurrogateData instances and assign to ProjectOrthogPolyApproximation instances
343 csg_driver->compute_trial_grid(var_sets);
344 numPts = var_sets.numCols();
345 if (verb > 1) {
346 PCout<<"Computing new index set:\n"<<*it<<endl;
347 //PCout<<RealMatrix(var_sets,Teuchos::TRANS)<<endl;
348 }
349 #ifdef GSGREST
350 fev = getFeval(evalGrid,evalFlag,evalData,var_sets,foundFlag);
351 if ( !foundFlag ) {
352 saveData((char *)"func",evalGrid,evalFlag,evalData);
353 return (0);
354 }
355 #else
356 fev = feval(var_sets,nQoI,ftype);
357 #endif
358
359 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
360 SurrogateData& sdi = poly_approx[iQoI].surrogate_data();
361 for( int jCol = 0; jCol < numPts; jCol++) {
362 SurrogateDataVars sdv(nvar,0,0);
363 SurrogateDataResp sdr(1,nvar); // no gradient or hessian
364 sdv.continuous_variables(Teuchos::getCol<int,double>(Teuchos::Copy,var_sets,jCol));
365 sdr.response_function(fev(jCol,iQoI));
366 sdi.push_back(sdv,sdr);
367 } // done loop over number of points
368 sdi.pop_count(numPts);
369 } // done loop over QoIs
370
371 shared_poly_data->increment_data();
372 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
373 poly_approx[iQoI].increment_coefficients();
374 }
375 }
376
377 /* Compute (normalized) change in variance */
378 RealVector respVarianceNew(nQoI,0.0);
379 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
380 auto poly_approx_rep = std::static_pointer_cast<PolynomialApproximation>
381 (poly_approx[iQoI].approx_rep());
382 respVarianceNew[iQoI] = poly_approx_rep->variance() ;
383 }
384 if (verb > 1) {
385 PCout<<" ... new variance:\n";
386 for ( int iQoI=0; iQoI<nQoI; iQoI++)
387 PCout<<respVarianceNew[iQoI]<<" ";
388 PCout<<"\n";
389 }
390 respVarianceNew -= respVariance;
391 if (verb > 1) {
392 PCout<<" ... delta variance:\n";
393 for ( int iQoI=0; iQoI<nQoI; iQoI++)
394 PCout<<respVarianceNew[iQoI]<<" ";
395 PCout<<"\n";
396 }
397 Real normChange = respVarianceNew.normFrobenius()/csg_driver->unique_trial_points();
398 if (normChange > deltaVar) {
399 deltaVar = normChange;
400 asave = *it;
401 }
402
403 shared_poly_data->decrement_data();
404 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
405 poly_approx[iQoI].pop_coefficients(true);
406 // Also restore the corresponding surrogate data
407 poly_approx[iQoI].surrogate_data().pop(true);
408 }
409 csg_driver->pop_set(); // reverse order from increment
410
411 } /* End iteration over proposed sets */
412
413 if (verb > 1)
414 PCout<<"Choosing :\n"<<asave
415 <<"\n ... with relative variance: "<<deltaVar<<std::endl ;
416
417 if ( asave.size() > 0 ) {
418 csg_driver->update_sets(asave);
419
420 //need to restore the data
421 size_t idxRestore = shared_poly_data->restore_index();
422 shared_poly_data->pre_push_data();
423 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
424 poly_approx[iQoI].push_coefficients();
425 poly_approx[iQoI].surrogate_data().push(idxRestore,true);
426 }
427 shared_poly_data->post_push_data();
428 csg_driver->update_reference();
429 }
430
431 if ( deltaVar < varEps )
432 { conv_within_tol = true; break; }
433
434 } /* end iteration loop */
435
436 // output sets; implementation above applies best refinement (no revert)
437 csg_driver->finalize_sets(true, conv_within_tol, false);
438
439 // sequence from ApproximationInterface::finalize_approximation():
440
441 // shared pre-finalize:
442 shared_poly_data->pre_finalize_data();
443
444 // per-approximation finalize:
445 for ( int iQoI=0; iQoI<nQoI; iQoI++) {
446 // from Approximation::finalize() called from PecosApproximation::finalize()
447 SurrogateData& sdi = poly_approx[iQoI].surrogate_data();
448 size_t i, num_restore = sdi.popped_sets(); // # of popped trial sets
449 for (i=0; i<num_restore; ++i)
450 sdi.push(shared_poly_data->finalize_index(i),false);
451 sdi.clear_popped();
452 // from PecosApproximation::finalize()
453 poly_approx[iQoI].finalize_coefficients();
454 }
455
456 // shared post-finalize:
457 shared_poly_data->post_finalize_data();
458
459 for ( int iQoI=0; iQoI<nQoI; iQoI++)
460 poly_approx[iQoI].print_coefficients(PCout,false);
461
462 return (0);
463
464 }
465
466 #ifdef GSGREST
467 /* Retrieve feval from saved data */
getFeval(RealMatrix & evalGrid,IntVector & evalFlag,RealMatrix & evalData,const RealMatrix & var_sets,bool & foundFlag)468 RealMatrix getFeval(RealMatrix &evalGrid, IntVector &evalFlag,
469 RealMatrix &evalData, const RealMatrix &var_sets,
470 bool &foundFlag) {
471
472 int numDim = var_sets.numRows(); // Dimensionality
473 int numPts = var_sets.numCols(); // Number of pts
474 int nouts = evalData.numCols(); // Dimensionality of dependent variables
475
476 RealMatrix fev(numPts,nouts);
477 if (evalGrid.numCols() == 0) {
478 /* Empty data, probably first call */
479 evalGrid = var_sets ;
480 evalFlag.size(numPts);
481 evalData.shape(numPts,nouts);
482 foundFlag = false;
483 return fev ;
484 }
485
486 /* check consistency */
487 assert( numDim == evalGrid.numRows());
488 assert( evalGrid.numCols() == evalFlag.length() );
489 assert( evalGrid.numCols() == evalData.numRows() );
490
491 foundFlag = true;
492 /* Loop over all grid points requested */
493 for (size_t i=0; i<numPts; i++) {
494 bool foundPt = false;
495
496 /* Loop over all saved grid points */
497 for (size_t j=0; j<evalGrid.numCols(); j++) {
498
499 /* compute L1 norm */
500 double dsum = 0.0;
501 for (size_t k = 0; k<numDim ; k++ )
502 dsum += fabs(var_sets(k,i)-evalGrid(k,j));
503
504 if ( dsum < 1.e-10 ) {
505
506 if (evalFlag(j) == 1 ) {
507 /* found saved point */
508 foundPt = true ;
509 for (size_t k = 0; k<nouts ; k++ )
510 fev(i,k) = evalData(j,k);
511 }
512
513 }
514
515 if (foundPt) break ;
516
517 } /* end loop over stored points */
518
519 if ( !foundPt ) {
520
521 /* missing point, add it to the list */
522 foundFlag = false;
523 /* Increment Grid */
524 RealMatrix evalGridSave = evalGrid;
525 evalGrid.shapeUninitialized(evalGridSave.numRows(),evalGridSave.numCols()+1);
526 for ( int j1 = 0; j1<evalGridSave.numCols(); j1++ )
527 for ( int i1 = 0; i1<evalGridSave.numRows(); i1++ )
528 evalGrid(i1,j1) = evalGridSave(i1,j1);
529 for ( int i1 = 0; i1<evalGridSave.numRows(); i1++ )
530 evalGrid(i1,evalGridSave.numCols()) = var_sets(i1,i);
531 /* Increment eval flag */
532 IntVector evalFlagSave = evalFlag;
533 evalFlag.size(evalFlagSave.length()+1);
534 for ( int i1 = 0; i1<evalFlagSave.length(); i1++ )
535 evalFlag(i1) = evalFlagSave(i1);
536 evalFlag(evalFlagSave.length()) = 0;
537 /* Increment eval func dimensions */
538 RealMatrix evalDataSave = evalData;
539 evalData.shape(evalDataSave.numRows()+1,evalDataSave.numCols());
540 for ( int j1 = 0; j1<evalDataSave.numCols(); j1++ )
541 for ( int i1 = 0; i1<evalDataSave.numRows(); i1++ )
542 evalData(i1,j1) = evalDataSave(i1,j1);
543 }
544
545 } /* end loop over stored points */
546
547
548 return fev;
549
550 }
551
getMatrixSize(const char * fname,int & nRows,int & nCols)552 void getMatrixSize(const char *fname, int &nRows, int &nCols) {
553
554
555 std::ifstream in(fname);
556
557 if(!in){
558 printf("getMatrixSize() : the requested file %s does not exist !\n",fname) ;
559 exit(1) ;
560 }
561
562 String theLine="";
563
564 // figure out number of lines and columns
565 int ix = 0 ;
566 while(in.good()) {
567
568 getline(in,theLine);
569
570 if ( theLine == "" ) break;
571 if ( theLine.compare(0,1,"#") == 0 ) continue ;
572
573 istringstream s(theLine);
574 int iy = 0 ;
575 double tmp ;
576 while( s >> tmp ) iy++ ;
577
578 if ( ( ix > 0 ) && ( iy != nCols ) )
579 {
580 printf("getMatrixSize() : Error at line %d !!!\n",ix+1) ;
581 printf(" no. of columns should be %d instead of %d\n",nCols,iy) ;
582 exit(1) ;
583 }
584
585 nCols = iy ;
586
587 ix++ ;
588
589 }
590
591 nRows = ix ;
592 in.close();
593
594 return ;
595
596 }
597
getMatrixTrans(const char * fname,RealMatrix & dataIn)598 void getMatrixTrans(const char *fname, RealMatrix &dataIn) {
599
600 int nRows = dataIn.numRows();
601 int nCols = dataIn.numCols();
602
603 if (nRows==0 || nCols==0){
604 printf("getMatrixTrans() : the requested data matrix is empty\n") ;
605 exit(1) ;
606 }
607
608 std::ifstream in(fname);
609
610 if(!in){
611 printf("getMatrixTrans() : the requested file %s does not exist\n",fname) ;
612 exit(1) ;
613 }
614
615 String theLine="";
616 int ix = 0;
617
618 while(in.good()){
619
620 getline(in,theLine);
621
622 if (theLine=="") break;
623 if ( theLine.compare(0, 1, "#") == 0 ) continue ;
624
625 istringstream s(theLine);
626 int iy = 0;
627 Real tmp;
628 while(s >> tmp){
629 dataIn(iy,ix)=tmp;
630 iy++;
631 }
632 if ( iy != nRows ) {
633 printf("getMatrixTrans(): Error at line %d while reading %s; number of columns should be %d\n",
634 ix+1, fname,nRows);
635 exit(1);
636 }
637 ix++;
638 }
639 if ( ix != nCols ) {
640 printf("getMatrixTrans(): Error while reading %s; number of rows should be %d\n",fname,nCols);
641 exit(1);
642 }
643 in.close();
644
645 return;
646
647 }
648
getMatrix(const char * fname,RealMatrix & dataIn)649 void getMatrix(const char *fname, RealMatrix &dataIn) {
650
651 int nRows = dataIn.numRows();
652 int nCols = dataIn.numCols();
653
654 if (nRows==0 || nCols==0){
655 printf("getMatrixTrans() : the requested data matrix is empty\n") ;
656 exit(1) ;
657 }
658
659 std::ifstream in(fname);
660
661 if(!in){
662 printf("getMatrixTrans() : the requested file %s does not exist\n",fname) ;
663 exit(1) ;
664 }
665
666 String theLine="";
667 int ix = 0;
668
669 while(in.good()){
670
671 getline(in,theLine);
672
673 if (theLine=="") break;
674 if ( theLine.compare(0, 1, "#") == 0 ) continue ;
675
676 istringstream s(theLine);
677 int iy = 0;
678 Real tmp;
679 while(s >> tmp){
680 dataIn(ix,iy)=tmp;
681 iy++;
682 }
683 if ( iy != nCols ) {
684 printf("getMatrixTrans(): Error at line %d while reading %s; number of columns should be %d\n",
685 ix+1, fname,nRows);
686 exit(1);
687 }
688 ix++;
689 }
690 if ( ix != nRows ) {
691 printf("getMatrixTrans(): Error while reading %s; number of rows should be %d\n",fname,nCols);
692 exit(1);
693 }
694 in.close();
695
696 return;
697
698 }
699
getIntVector(const char * fname,IntVector & dataIn)700 void getIntVector(const char *fname, IntVector &dataIn) {
701
702 int nvals=dataIn.length();
703
704 if (nvals==0){
705 printf("getIntVector() : the requested data array is empty\n") ;
706 exit(1) ;
707 }
708
709 int nCols=1;
710
711 std::ifstream in(fname);
712
713 if(!in){
714 printf("getIntVector() : the requested file %s does not exist\n",fname) ;
715 exit(1) ;
716 }
717
718 string theLine="";
719 int ix=0;
720
721 while( in.good() ){
722
723 getline(in,theLine);
724
725 if (theLine=="") break;
726
727 istringstream s(theLine);
728 int iy=0;
729 int tmp;
730 s >> tmp;
731 dataIn(ix) = tmp;
732 iy++;
733 if (s>>tmp) {
734 printf("getIntVector(): Error at line %d while reading %s; number of columns should be %d\n",
735 ix+1, fname,nCols);
736 exit(1);
737 }
738 ix++;
739 }
740 if ( ix != nvals ) {
741 printf("getIntVector(): Error while reading %s; number of rows should be %d\n",
742 fname,nvals);
743 exit(1);
744 }
745 in.close();
746
747 return ;
748
749 }
750
751 /* Save feval to files */
saveData(const char * fbase,const RealMatrix & evalGrid,const IntVector & evalFlag,const RealMatrix & evalData)752 void saveData(const char *fbase, const RealMatrix &evalGrid,
753 const IntVector &evalFlag, const RealMatrix &evalData) {
754
755 char fname[100];
756
757 std::ofstream fout;
758
759 /* save grid */
760 sprintf(fname,"%s_grid.dat",fbase);
761 fout.open(fname);
762 for (size_t i=0; i < evalGrid.numCols(); i++ ) {
763 for (size_t j=0; j < evalGrid.numRows(); j++ ) {
764 fout.precision(16);
765 fout << std::scientific << std::setw(24) << evalGrid(j,i) << " ";
766 }
767 fout << std::endl ;
768 }
769 fout.close() ;
770
771 /* save evalFlag */
772 sprintf(fname,"%s_flag.dat",fbase);
773 fout.open(fname);
774 for (size_t i=0; i < evalFlag.length(); i++ )
775 fout << evalFlag(i) << std::endl;
776 fout.close() ;
777
778 /* save output matrix */
779 sprintf(fname,"%s_data.dat",fbase);
780 fout.open(fname);
781 fout.precision(16);
782 for (size_t j=0; j < evalData.numRows(); j++ ) {
783 for (size_t i=0; i < evalData.numCols(); i++ ) {
784 fout << std::scientific << std::setw(24) << evalData(j,i) << " ";
785 }
786 fout << std::endl ;
787 }
788 fout.close() ;
789
790 return ;
791
792 }
793
794 /* Load feval from files */
loadData(const char * fbase,RealMatrix & evalGrid,IntVector & evalFlag,RealMatrix & evalData)795 void loadData(const char *fbase, RealMatrix &evalGrid,
796 IntVector &evalFlag, RealMatrix &evalData) {
797
798 char fname[100];
799
800 /* Check if all files exist */
801 sprintf(fname,"%s_grid.dat",fbase);
802 if ( !fileExist(fname) ) return ;
803 sprintf(fname,"%s_flag.dat",fbase);
804 if ( !fileExist(fname) ) return ;
805 sprintf(fname,"%s_data.dat",fbase);
806 if ( !fileExist(fname) ) return ;
807
808
809 /* load grid */
810 sprintf(fname,"%s_grid.dat",fbase);
811 int nDims, nPts ;
812 getMatrixSize(fname, nPts, nDims) ;
813 evalGrid.shapeUninitialized(nDims,nPts);
814 getMatrixTrans(fname, evalGrid);
815
816 /* load eval flag */
817 evalFlag.sizeUninitialized(nPts);
818 sprintf(fname,"%s_flag.dat",fbase);
819 getIntVector(fname, evalFlag);
820
821 /* load function outputs */
822 int nOuts ;
823 sprintf(fname,"%s_data.dat",fbase);
824 getMatrixSize(fname, nPts, nOuts) ;
825 evalData.shapeUninitialized(nPts,nOuts);
826 assert(nPts == evalGrid.numCols());
827 getMatrix(fname, evalData);
828
829 bool allEvals = true ;
830 for (size_t i=0;i<nPts;i++)
831 if (evalFlag(i) != 1) {
832 PCout<<" Point "<<i<<" is not evaluated!"<<std::endl;
833 allEvals = false;
834 }
835
836 if (!allEvals) {
837 PCout<<" -> Exit !"<<std::endl;
838 std::terminate() ;
839 }
840
841 return ;
842
843 }
844
845 #else
846
feval(const RealMatrix & dataMat,const int nQoI,String ftype)847 RealMatrix feval(const RealMatrix &dataMat, const int nQoI, String ftype)
848 {
849
850 assert(nQoI==1);
851
852 int i, j ;
853
854 int numDim = dataMat.numRows(); // Dimensionality
855 int numPts = dataMat.numCols(); // Number of pts
856
857 /* Count the number of function evaluations; */
858 RealMatrix fev(numPts,nQoI);
859 for (i=0; i<numPts; ++i) {
860 RealVector xIn(numDim);
861 for (j=0; j<numDim; ++j) xIn[j] = dataMat(j,i);
862 //fev(ieval,0)=genz(String("cp1"), xIn);
863 if ( ftype == String("pol2") )
864 fev(i,0) = custPol(ftype, xIn);
865 else if ( ftype == String("gerstner-iso1") )
866 fev(i,0) = gerstner(String("iso1"), xIn);
867 else if ( ftype == String("gerstner-iso2") )
868 fev(i,0) = gerstner(String("iso2"), xIn);
869 else if ( ftype == String("gerstner-iso3") )
870 fev(i,0) = gerstner(String("iso3"), xIn);
871 else if ( ftype == String("gerstner-aniso1") )
872 fev(i,0) = gerstner(String("aniso1"), xIn);
873 else if ( ftype == String("gerstner-aniso2") )
874 fev(i,0) = gerstner(String("aniso2"), xIn);
875 else if ( ftype == String("gerstner-aniso3") )
876 fev(i,0) = gerstner(String("aniso3"), xIn);
877 else
878 throw(std::runtime_error("Pecos::feval() unknown ftype\n"));
879 }
880
881 return fev;
882
883 }
884
885 #endif
886