diff --git a/benchmarks/basic_tests.py b/benchmarks/basic_tests.py index 331db4445..a9eec28e5 100644 --- a/benchmarks/basic_tests.py +++ b/benchmarks/basic_tests.py @@ -187,6 +187,8 @@ def sweep_forgive_test(): assert diff.max().max() == 0.0 + + def inv_regul_test(): model_d = "ies_10par_xsec" @@ -1437,11 +1439,56 @@ def run(): pyemu.os_utils.start_workers(t_d, exe_path, pst_name, num_workers=15, worker_root=model_d, port=4004) +def sweep_bin_test(): + + model_d = "ies_10par_xsec" + t_d = os.path.join(model_d,"template") + m_d = os.path.join(model_d,"master_sweep_bin") + if os.path.exists(m_d): + shutil.rmtree(m_d) + pst = pyemu.Pst(os.path.join(t_d,"pest.pst")) + pe = pyemu.ParameterEnsemble.from_uniform_draw(pst,num_reals=50)#.loc[:,pst.par_names[:2]] + + pe.to_csv(os.path.join(t_d,"sweep_in.csv")) + pe._df.index = pe.index.map(str) + print(pe.index) + pe.to_dense(os.path.join(t_d,"sweep_in.bin")) + pst.pestpp_options["ies_par_en"] = "sweep_in.csv" + pst.pestpp_options["sweep_forgive"] = True + pst.pestpp_options["sweep_parameter_file"] = "sweep_in.bin" + pst.control_data.noptmax = -1 + pst.pestpp_options.pop("ies_num_reals",None) + pst.write(os.path.join(t_d,"pest_forgive.pst")) + pst.pestpp_options["sweep_output_file"] = "sweep_out.bin" + pst.pestpp_options["sweep_chunk"] = 9 + pst.pestpp_options["ies_include_base"] = False + pst.write(os.path.join(t_d,"pest_forgive.pst")) + m_d = os.path.join(model_d,"master_sweep_bin_base") + pyemu.os_utils.start_workers(t_d, exe_path, "pest_forgive.pst", 10, master_dir=m_d, + worker_root=model_d,port=port) + df1 = pd.read_csv(os.path.join(m_d, "pest_forgive.0.obs.csv"),index_col=0) + assert df1.shape[0] == pe.shape[0] + m_d = os.path.join(model_d, "master_sweep_bin") + pyemu.os_utils.start_workers(t_d, exe_path.replace("-ies", "-swp"), "pest_forgive.pst", 10, master_dir=m_d, + worker_root=model_d, port=port) + df2 = pyemu.Matrix.from_binary(os.path.join(m_d,"sweep_out.bin")).to_dataframe() + print(df2) + print(df1) + assert df2.shape == df1.shape + diff = (df1.values - df2.values) + print(diff) + print(diff.max()) + print(np.abs(diff).max()) + assert np.abs(diff).max() < 1e-7 + + + if __name__ == "__main__": #run() #mf6_v5_ies_test() #prep_ends() - mf6_v5_sen_test() + sweep_bin_test() + #mf6_v5_sen_test() #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-glm.exe"),os.path.join("..","bin","win","pestpp-glm.exe")) #shutil.copy2(os.path.join("..", "exe", "windows", "x64", "Debug", "pestpp-ies.exe"), # os.path.join("..", "bin", "win", "pestpp-ies.exe")) diff --git a/benchmarks/ies_10par_xsec/template/10par_xsec.hds b/benchmarks/ies_10par_xsec/template/10par_xsec.hds index fa6fa39b2..e69de29bb 100644 --- a/benchmarks/ies_10par_xsec/template/10par_xsec.hds +++ b/benchmarks/ies_10par_xsec/template/10par_xsec.hds @@ -1,2 +0,0 @@ - 1.000 1.200 1.400 1.600 1.800 2.000 2.200 2.400 2.600 2.800 - 1.000 1.400 1.800 2.200 2.600 3.000 3.400 3.800 4.200 4.600 diff --git a/benchmarks/ies_10par_xsec/template/hk_Layer_1.ref b/benchmarks/ies_10par_xsec/template/hk_Layer_1.ref index 078687bf9..e69de29bb 100644 --- a/benchmarks/ies_10par_xsec/template/hk_Layer_1.ref +++ b/benchmarks/ies_10par_xsec/template/hk_Layer_1.ref @@ -1 +0,0 @@ - 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 diff --git a/benchmarks/ies_10par_xsec/template/nested/really/deep/hk_Layer_1.ref b/benchmarks/ies_10par_xsec/template/nested/really/deep/hk_Layer_1.ref index 078687bf9..e69de29bb 100644 --- a/benchmarks/ies_10par_xsec/template/nested/really/deep/hk_Layer_1.ref +++ b/benchmarks/ies_10par_xsec/template/nested/really/deep/hk_Layer_1.ref @@ -1 +0,0 @@ - 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 2.500000000000 diff --git a/benchmarks/ies_10par_xsec/template/strt_Layer_1.ref b/benchmarks/ies_10par_xsec/template/strt_Layer_1.ref index 5fd4b9532..e69de29bb 100644 --- a/benchmarks/ies_10par_xsec/template/strt_Layer_1.ref +++ b/benchmarks/ies_10par_xsec/template/strt_Layer_1.ref @@ -1 +0,0 @@ - 1.00000000000000 1.000000E+00 1.000000E+00 1.000000E+00 1.000000E+00 1.000000E+00 1.000000E+00 1.000000E+00 1.000000E+00 1.000000E+00 diff --git a/benchmarks/mf6_freyberg/template/freyberg6_run_sen.pst b/benchmarks/mf6_freyberg/template/freyberg6_run_sen.pst index a9dfbedba..3d589ae67 100644 --- a/benchmarks/mf6_freyberg/template/freyberg6_run_sen.pst +++ b/benchmarks/mf6_freyberg/template/freyberg6_run_sen.pst @@ -1,13 +1,13 @@ pcf version=2 * control data keyword -pestmode estimation -noptmax 0 -svdmode 1 -maxsing 10000000 -eigthresh 1e-06 -eigwrite 1 -ies_par_en prior.jcb -tie_by_group True +pestmode estimation +noptmax 0 +svdmode 1 +maxsing 10000000 +eigthresh 1e-06 +eigwrite 1 +ies_par_en prior.jcb +tie_by_group true * parameter groups external freyberg6_run_sen.pargrp_data.csv * parameter data external diff --git a/documentation/pestpp_users_guide_v5.2.9.docx b/documentation/pestpp_users_guide_v5.2.10.docx similarity index 72% rename from documentation/pestpp_users_guide_v5.2.9.docx rename to documentation/pestpp_users_guide_v5.2.10.docx index 8abd4153d..07c660ec1 100644 Binary files a/documentation/pestpp_users_guide_v5.2.9.docx and b/documentation/pestpp_users_guide_v5.2.10.docx differ diff --git a/documentation/pestpp_users_manual.md b/documentation/pestpp_users_manual.md index 5f08b6602..747ea558b 100644 --- a/documentation/pestpp_users_manual.md +++ b/documentation/pestpp_users_manual.md @@ -1,13 +1,13 @@ A close up of a purple sign Description automatically generated -# Version 5.2.9 +# Version 5.2.10 PEST++ Development Team -Jan 2024 +April 2024 # Acknowledgements @@ -70,7 +70,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI # Table of Contents -- [Version 5.2.9](#s1) +- [Version 5.2.10](#s1) - [Acknowledgements](#s2) - [Preface](#s3) - [License](#s4) @@ -3475,6 +3475,8 @@ Optionally, users can forego the use of measurement noise realizations through t Still another option to cope with the need to define weights and noise that are non-commensurate is to employ the *external* control file section format for observation data (e.g., *\* observation data external*) and in the external observation data file(s), supply the column *standard_deviation* as values of expected noise and supply the *weight* column in the external files as the values necessary to form a balanced composite objective function. In this way, PESTPP-IES will generate realizations of observation noise using the covariance matrix implied by the *standard_deviation* entries and will use the weights during the upgrade calculation process and during objective function calculations. +Given all the issues that have been encountered by users with the way noise and weights interact, starting inversion 5.2.10, if a user does not indicate to use noise relizations either explicitly (through the *ies_no_noise* option) or implicitly (through the ies_observation_ensemble option) or by providing an observation noise representation (through either an *obscov* or by supplying a *standard_deviation* for at least one non-zero weighted observations), PESTPP-IES will automatically default to a no-noise setting; this default choice will be echoed to the rec and stdout at runtime. + ### 9.1.7 Regularization The term “regularization” generally describes numerical devices that are employed to achieve uniqueness of ill-posed inverse problems. Regularization is thus fundamental to the notion of calibration. If regularization is properly implemented, the calibrated parameter field is one of minimized error variance. In most calibration contexts, regularization eliminates from the calibrated parameter field any heterogeneity that is not supported by direct measurements of hydraulic properties or by the calibration dataset. @@ -3489,7 +3491,7 @@ As with PEST(\_HP) and PESTPP-GLM, the SVD truncation controls (i.e., MAXSING an ### 9.1.8 Base Realization -Optionally (and by default), the parameter ensemble used by PESTPP-IES can include a “base realization”. Parameter values that comprise this realization are those which are listed in the “parameter data” section of the PEST control file which PESTPP-IES reads on commencement of execution. Ideally, this set of parameter values are those of minimum pre-calibration error variance; that is, they comprise the expected values of parameters from an expert knowledge point of view. PESTPP-IES pairs this realization with an observation dataset that has no “manufactured” measurement noise associated with it; this dataset is comprised of measurements that appear in the “observation data” section of the PEST control file. +Optionally (and by default), the parameter ensemble used by PESTPP-IES can include a “base realization”. Parameter values that comprise this realization are those which are listed in the “parameter data” section of the PEST control file which PESTPP-IES reads on commencement of execution. Ideally, this set of parameter values are those of minimum pre-calibration error variance; that is, they comprise the expected values of parameters from an expert knowledge point of view. PESTPP-IES pairs this realization with an observation dataset that has no “manufactured” measurement noise associated with it; this dataset is comprised of measurements that appear in the “observation data” section of the PEST control file. After completion of PESTPP-IES, this “base” realization is the approximate minimum error variance parameter set and is as close a result to what standard regularized least-squares parameter estimation will yield. If a user is interested in using only one parameter set from a PESTPP-IES analysis, they should always and only ever select the base realization for additional analysis. By monitoring the fate of the parameter set comprising this base realization, a user can witness the effect that the ensemble smoother process has on a non-random parameter field. Any bias that is introduced to this parameter field, or any incredulous heterogeneity that is introduced to this parameter field, is also presumably introduced to the other random parameter fields which comprise the ensemble. Inspection of this field can aid a modeller in assessing whether, in his/her opinion, the parameter field ensemble that emerges from the ensemble smoother process comprises a legitimate sample of the posterior parameter probability distribution. @@ -3719,7 +3721,7 @@ Where a single parameter field is being estimated in a standard model calibratio The situation is different, however, when many parameter fields comprising an ensemble are adjusted, such as is done by PESTPP-IES. For an ensemble smoother, the testing of parameter upgrades and the filling of a Jacobian matrix are the same operation. The cost of this operation would become prohibitive if a trial-and-error lambda search procedure accompanied the adjustment of each parameter realization comprising an ensemble. Hence a different strategy must be adopted. This strategy is to undertake lambda testing for only a limited number of parameter fields, and to then use the best lambda that emerges from that testing process when upgrading the rest of them. This limited number of realizations used to evaluate objective function values during upgrade testing is referred to as the “subset”. Options for how to pick the subset are available through the *ies_subset_how()* argument; valid choices for *ies_subset_how()* are “first” (the first *ies_subset_size()* realizations), “last” (the last *ies_subset_size()* realizations), “random” (randomly select *ies_subset_size()* realizations for each iteration), or “phi_based” (select *ies_subset_size()* realizations across the previous composite objective function distribution). Note that if the “base” parameter realization is present, it is always included in the selected subset as the objective function calculated for this realization has special significance. -The number of realizations that comprise the ensemble subset used for lambda testing is set by the value of the *ies_subset_size()* control variable. During each iteration of the ensemble smoother process, values of the Marquardt lambda used for testing realization upgrades are determined by applying a set of multipliers to the best lambda found during the previous iteration. These multipliers are provided through the *ies_lambda_mults()* control variable. A comma separated list of multipliers should be supplied by the user as arguments to this keyword; at least one of these multipliers should be less than 1.0 while, or course, one of them should be greater than 1.0. Line search factors (otherwise known as scale factors) that are applied to each of these lambdas can also be supplied. If so, this is done through the *lambda_scale_fac()* control variable, the same variable that is used by PESTPP-GLM. As for *ies_lambda_mults()*, scale factors should be supplied as a comma-separated list of numbers spanning a range from below 1.0 to greater than 1.0. The total number of model runs required to test parameter upgrades during a given iteration is thus *ies_subset_size()* times the number of multipliers supplied with the *ies_lambda_mults()* control variable times the number of factors supplied with the *lambda_scale_fac()* control variable. +The number of realizations that comprise the ensemble subset used for lambda testing is set by the value of the *ies_subset_size()* control variable, which is given a default value of -10, which indicates to use 10% of the current number of realizations for subset testing; experience has shown this is a pretty good default value in most cases. During each iteration of the ensemble smoother process, values of the Marquardt lambda used for testing realization upgrades are determined by applying a set of multipliers to the best lambda found during the previous iteration. These multipliers are provided through the *ies_lambda_mults()* control variable. A comma separated list of multipliers should be supplied by the user as arguments to this keyword; at least one of these multipliers should be less than 1.0 while, or course, one of them should be greater than 1.0. Line search factors (otherwise known as scale factors) that are applied to each of these lambdas can also be supplied. If so, this is done through the *lambda_scale_fac()* control variable, the same variable that is used by PESTPP-GLM. As for *ies_lambda_mults()*, scale factors should be supplied as a comma-separated list of numbers spanning a range from below 1.0 to greater than 1.0. The total number of model runs required to test parameter upgrades during a given iteration is thus *ies_subset_size()* times the number of multipliers supplied with the *ies_lambda_mults()* control variable times the number of factors supplied with the *lambda_scale_fac()* control variable. The value of the Marquardt lambda to use during the first iteration of the ensemble smoother process can be supplied through the *ies_initial_lambda()* control variable. Lambda multipliers supplied through *ies_lambda_mults()* are applied to this value during the first iteration of this process. The PESTPP-IES default value for *ies_initial_lambda()* is $10^{floor\left( \log_{10}\frac{\mu_{Փ}}{2n} \right)}$ where *μ*Փ is the mean of objective functions achieved using realizations comprising the initial ensemble, and *n* is the number of non-zero-weighted observations featured in the “observation data” section of the PEST control file. @@ -4165,13 +4167,13 @@ On most occasions of PESTPP-SWP usage, model runs are conducted in parallel. Use As usual, variables which control how PESTPP-SWP operates must be placed in a PEST control file whose name is supplied on its command line; these variables should appear on lines that begin with the “++” character string. -PESTPP-SWP is directed to a CSV or JCO/JCB input file through the value supplied for its *sweep_parameter_csv_file()* control variable; the type of file is recognized by its extension. A CSV file must have as many columns as there are parameters featured in the PEST control file, plus one extra column on the left. The first column is reserved for the user-supplied realization name. Parameter and realizations names are provided in a JCB or JCO file according to the respective protocols of these files. +PESTPP-SWP is directed to a CSV, JCO/JCB, or BIN (dense binary) input file through the value supplied for its *sweep_parameter file()* control variable; the type of file is recognized by its extension. A CSV file must have as many columns as there are parameters featured in the PEST control file, plus one extra column on the left. The first column is reserved for the user-supplied realization name. Parameter and realizations names are provided in a JCB or JCO file according to the respective protocols of these files. A dense binary file can be written by pyEMU and is an efficient binary format that allows for chunks of runs to be read (compared to the jacobian formats that require reading the entire file at once). The number of realizations contained in a user-prepared PESTPP-SWP input file depends on the number of parameter sets for which model runs are required. These realizations can be named according to the user’s taste. PESTPP-SWP carries out one model run for each realization. -Parameter names provided in a CSV of JCO/JCB file must correspond to those that are featured in a PEST control file. +Parameter names provided in the parameter file must correspond to those that are featured in a PEST control file. -If the *sweep_parameter_csv_file()* control variable does not appear in the PEST control file that is cited on the PESTPP-SWP command line, PESTPP-SWP assumes an input filename of *sweep_in.csv*. +If the *sweep_parameter_file()* control variable does not appear in the PEST control file that is cited on the PESTPP-SWP command line, PESTPP-SWP assumes an input filename of *sweep_in.csv*. Note the following. @@ -4183,7 +4185,7 @@ Note the following. PESTPP-SWP can fill in values for fixed and tied parameters if these are missing from its input file. Actually, it can provide values for other missing parameters as well if the *sweep_forgive()* control variable is set to *true*. These missing values are taken from the PEST control file which is read by PESTPP-SWP. -PESTPP-SWP writes a model output file whose name is provided through the *sweep_output_csv_file()* control variable. If this variable is not provided, PESTPP-SWP employs the name *sweep_out.csv* for its output file. This fiile contains input and output run identifiers, objective function information, and the simulated value for each control file “observation” for each run. Note that, by default, PESTPP-SWP will not calculate the contribution of prior information equations to the total objective function reported in the sweep output file. This can be overridden with the *sweep_include_regul_phi* option. +PESTPP-SWP writes a model output file whose name is provided through the *sweep_output_file()* control variable. If this variable is not provided, PESTPP-SWP employs the name *sweep_out.csv* for its output file, again the file extension matters: CSV, JCO/JCB, or BIN. This file contains input and output run identifiers, objective function information, and the simulated value for each control file “observation” for each run. Note that, by default, PESTPP-SWP will not calculate the contribution of prior information equations to the total objective function reported in the sweep output file. This can be overridden with the *sweep_include_regul_phi* option. The control variable *sweep_chunk()* pertains to parallelization of model runs. Runs are done in bundles of size *N*, where *N* is the value supplied for this variable. (A chunk of 500 is the default). This number should be chosen wisely. It should be a multiple of the number of agents that PESTPP-SWP can use for carrying out model runs. @@ -5042,6 +5044,7 @@ Variables discussed in section 5.3.6 that control parallel run management are no | *Mou_simplex_factors(0.5,.0.7,0.8)* | double | Backtracking points to test along each reflected simplex individual. | | *Mou_simplex_mutation(false)* | boolean | Flag to add guassian mutation to the reflected simplex individuals. Default is false | | *Mou_multigen_population* | boolean | Flag to retain and reuse all members across all generations when evaluating dominance and feasibility. This can result in re-enforcing the preference for feasible solutions, which can help with highly nonlinearly constrained problems. However, with a large population and many generations, this option can slow down the execution time of PESTPP-MOU since it must dominance sort a much large number of members. Default is false. Note the option activate automatically when using chances with “all’ chance points when chances are reused across generations. | +| *Mou_chance_schedule()* | text | A two column ascii file that defines when to re-evaluate chances. Generations not listed are set to false. This can be useful to have more granular control regarding chance evaluation. | Table 13.2. PESTPP-MOU specific control arguments. PESTPP-MOU shares many other control arguments with PESTPP-OPT diff --git a/scripts/build_pestpp_win_no_ifort.bat b/scripts/build_pestpp_win_no_ifort.bat index 5fb346164..9948a134b 100644 --- a/scripts/build_pestpp_win_no_ifort.bat +++ b/scripts/build_pestpp_win_no_ifort.bat @@ -18,6 +18,7 @@ copy /y *.zip ..\ cd .. rem call "C:\Program Files (x86)\Microsoft Visual Studio\2019\Community\VC\Auxiliary\Build\vcvarsall.bat" x64 +call "C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvarsall.bat" x64 rmdir /Q /S build mkdir build cd build diff --git a/src/libs/common/common.vcxproj b/src/libs/common/common.vcxproj index 78f2911a3..8b03ce990 100644 --- a/src/libs/common/common.vcxproj +++ b/src/libs/common/common.vcxproj @@ -36,38 +36,38 @@ StaticLibrary true - v142 + v143 NotSet StaticLibrary true - v142 + v143 NotSet StaticLibrary true - v142 + v143 NotSet DynamicLibrary true - v142 + v143 NotSet StaticLibrary false - v142 + v143 true NotSet StaticLibrary false - v142 + v143 true NotSet diff --git a/src/libs/common/config_os.h b/src/libs/common/config_os.h index f4595371d..15d9c7d41 100644 --- a/src/libs/common/config_os.h +++ b/src/libs/common/config_os.h @@ -2,7 +2,7 @@ #define CONFIG_OS_H_ -#define PESTPP_VERSION "5.2.9"; +#define PESTPP_VERSION "5.2.10"; #if defined(_WIN32) || defined(_WIN64) #define OS_WIN diff --git a/src/libs/common/utilities.cpp b/src/libs/common/utilities.cpp index 405a1ebfe..222878f73 100644 --- a/src/libs/common/utilities.cpp +++ b/src/libs/common/utilities.cpp @@ -851,7 +851,250 @@ bool parse_string_arg_to_bool(string arg) } } +vector read_dense_binary_col_names(ifstream& in,int n_col) +{ + //first read the names of the columns + stringstream ss; + vector col_name_sizes; + int name_size = 0; + for (int i = 0; i < n_col; i++) + { + in.read((char*)&(name_size), sizeof(name_size)); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format error reading size column name size for column number " << i; + throw runtime_error(ss.str()); + } + + col_name_sizes.push_back(name_size); + } + int i = 0; + string name; + vector col_names; + for (auto col_name_size : col_name_sizes) + { + char* col_name = new char[col_name_size]; + in.read(col_name, col_name_size); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format error reading column name for column number " << i << ", size " << col_name_size; + throw runtime_error(ss.str()); + } + name = string(col_name, col_name_size); + pest_utils::strip_ip(name); + pest_utils::upper_ip(name); + col_names.push_back(name); + i++; + delete[] col_name; + } + return col_names; + +} + +bool read_dense_binary_records(ifstream& in,int n_records,int n_col,vector& row_names, vector>& rec_vecs) +{ + stringstream ss; + int i, name_size; + streampos current_pos = in.tellg(); + in.seekg(0,std::ios::end); + streampos end = in.tellg(); + in.seekg(current_pos,std::ios::beg); + bool success = true; + string name; + double data; + vector rec; + rec_vecs.clear(); + row_names.clear(); + i = 0; + while (true) + { + //finished + //if ((in.bad()) || (in.eof())) + if ((i > 0) && (!in.good())) + { + break; + } + if (in.tellg() == end) { + break; + } + in.read((char*)&(name_size), sizeof(name_size)); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error reading row name size for row number " << i << "...continuing"; + cout << ss.str(); + success = false; + break; + } + char* row_name = new char[name_size]; + in.read(row_name, name_size); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error reading row name for row number " << i << "...continuing"; + cout << ss.str(); + success = false; + break; + } + name = string(row_name, name_size); + delete[] row_name; + pest_utils::strip_ip(name); + pest_utils::upper_ip(name); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error skipping values for row " << i << "...continuing "; + cout << ss.str(); + success = false; + break; + } + + rec.clear(); + rec.resize(n_col,0); + for (int j = 0; j < n_col; j++) + { + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error reading row,col value " << i << "," << j << "...continuing "; + cout << ss.str(); + success = false; + break; + } + in.read((char*)&(data), sizeof(data)); + //rec.push_back(data); + rec[j] = data; + + } + if (in.eof()) + break; + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error skipping values for row " << i << "...continuing "; + cout << ss.str(); + success = false; + break; + } + row_names.push_back(name); + rec_vecs.push_back(rec); + + i++; + if (i >= n_records) + { + break; + } + + } + return success; +} + +vector read_dense_binary_remaining_row_names(ifstream& in,const vector& col_names) +{ + stringstream ss; + int name_size = 0; + string name; + int i = 0; + double data = -1.; + vector row_names; + // record current position in file + streampos current_pos = in.tellg(); + in.seekg(0,std::ios::end); + streampos end = in.tellg(); + in.seekg(current_pos,std::ios::beg); + + while (true) + { + //finished + //if ((in.bad()) || (in.eof())) + if ((i > 0) && (!in.good())) + { + break; + } + if (in.tellg() == end) { + break; + } + in.read((char*)&(name_size), sizeof(name_size)); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error reading row name size for row number " << i << "...continuing"; + cout << ss.str(); + break; + } + char* row_name = new char[name_size]; + in.read(row_name, name_size); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error reading row name for row number " << i << "...continuing"; + cout << ss.str(); + break; + } + name = string(row_name, name_size); + delete[] row_name; + pest_utils::strip_ip(name); + pest_utils::upper_ip(name); + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error skipping values for row " << i << "...continuing "; + cout << ss.str(); + break; + } + + //skip the values + //in.seekg(col_names.size() * sizeof(double), ios_base::cur); + char* rest_of_line = new char[sizeof(double) * col_names.size()]; + in.read(rest_of_line, sizeof(double)* col_names.size()); + delete[] rest_of_line; + if (in.eof()) + break; + if (!in.good()) + { + ss.str(""); + ss << "read_dense_binary(), dense format incomplete record: error skipping values for row " << i << "...continuing "; + cout << ss.str(); + break; + } + row_names.push_back(name); + i++; + } + in.seekg(current_pos,std::ios::beg); + return row_names; +} + +void read_binary_matrix_header(ifstream& in, int& tmp1, int& tmp2, int& tmp3) +{ + stringstream ss; + + if (!in.good()) + { + ss.str(""); + ss << "read_binary_matrix_header - stream is not good"; + throw runtime_error(ss.str()); + } + in.read((char*)&tmp1, sizeof(tmp1)); + in.read((char*)&tmp2, sizeof(tmp2)); + in.read((char*)&tmp3, sizeof(tmp3)); + if (!in.good()) + { + ss.str(""); + ss << "read_binary_matrix_header - stream is not good"; + throw runtime_error(ss.str()); + } + + + +} + +bool is_dense_binary_matrix(int tmp1, int tmp2, int tmp3) +{ + return ((tmp1 == 0) && (tmp2 < 0) && (tmp3 < 0) && (tmp2 == tmp3)); +} void read_dense_binary(const string& filename, vector& row_names, vector& col_names, Eigen::MatrixXd& matrix) { @@ -874,125 +1117,24 @@ void read_dense_binary(const string& filename, vector& row_names, vector int n_obs_and_pi; int i, j; double data; - //char* col_name; - //char* row_name; - // read header - in.read((char*)&n_par, sizeof(n_par)); - in.read((char*)&n_obs_and_pi, sizeof(n_obs_and_pi)); - in.read((char*)&n_nonzero, sizeof(n_nonzero)); + read_binary_matrix_header(in,n_par,n_obs_and_pi,n_nonzero); - if ((n_par == 0) && (n_obs_and_pi < 0) && (n_nonzero < 0) && (n_obs_and_pi == n_nonzero)) + //if ((n_par == 0) && (n_obs_and_pi < 0) && (n_nonzero < 0) && (n_obs_and_pi == n_nonzero)) + if (is_dense_binary_matrix(n_par,n_obs_and_pi,n_nonzero)) { n_obs_and_pi *= -1; cout << "reading 'dense' format matrix with " << n_obs_and_pi << " columns" << endl; - //first read the names of the columns - vector col_name_sizes; - int name_size = 0; - for (int i = 0; i < n_obs_and_pi; i++) - { - in.read((char*)&(name_size), sizeof(name_size)); - if (!in.good()) - { - ss.str(""); - ss << "read_dense_binary(), dense format error reading size column name size for column number " << i; - throw runtime_error(ss.str()); - } + //first read the names of the columns and the rows + col_names = read_dense_binary_col_names(in,n_obs_and_pi); + row_names = read_dense_binary_remaining_row_names(in,col_names); - col_name_sizes.push_back(name_size); - } - int i = 0; - string name; - for (auto col_name_size : col_name_sizes) - { - char* col_name = new char[col_name_size]; - in.read(col_name, col_name_size); - if (!in.good()) - { - ss.str(""); - ss << "read_dense_binary(), dense format error reading column name for column number " << i << ", size " << col_name_size; - throw runtime_error(ss.str()); - } - name = string(col_name, col_name_size); - pest_utils::strip_ip(name); - pest_utils::upper_ip(name); - col_names.push_back(name); - i++; - delete[] col_name; - } - i = 0; - double data = -1.; - // record current position in file - streampos begin_rows = in.tellg(); - in.seekg(0,std::ios::end); - streampos end = in.tellg(); - in.seekg(begin_rows,std::ios::beg); - //read the row names so we can dimension the matrix - while (true) - { - //finished - //if ((in.bad()) || (in.eof())) - if ((i > 0) && (!in.good())) - { - break; - } - if (in.tellg() == end) { - break; - } - in.read((char*)&(name_size), sizeof(name_size)); - if (!in.good()) - { - ss.str(""); - ss << "read_dense_binary(), dense format incomplete record: error reading row name size for row number " << i << "...continuing"; - cout << ss.str(); - break; - } - char* row_name = new char[name_size]; - in.read(row_name, name_size); - if (!in.good()) - { - ss.str(""); - ss << "read_dense_binary(), dense format incomplete record: error reading row name for row number " << i << "...continuing"; - cout << ss.str(); - break; - } - name = string(row_name, name_size); - delete[] row_name; - pest_utils::strip_ip(name); - pest_utils::upper_ip(name); - if (!in.good()) - { - ss.str(""); - ss << "read_dense_binary(), dense format incomplete record: error skipping values for row " << i << "...continuing "; - cout << ss.str(); - break; - } - - //skip the values - //in.seekg(col_names.size() * sizeof(double), ios_base::cur); - char* rest_of_line = new char[sizeof(double) * col_names.size()]; - in.read(rest_of_line, sizeof(double)* col_names.size()); - delete[] rest_of_line; - if (in.eof()) - break; - if (!in.good()) - { - ss.str(""); - ss << "read_dense_binary(), dense format incomplete record: error skipping values for row " << i << "...continuing "; - cout << ss.str(); - break; - } - row_names.push_back(name); - i++; - } in.close(); in.open(filename.c_str(), ifstream::binary); //resize the matrix now that we know big it should be matrix.resize(row_names.size(), col_names.size()); - //seek back to the first row - in.seekg(begin_rows, ios_base::beg); for (int i=0;i& row_names, vector } } } + else + { + throw runtime_error("binary matrix header values do not indicate a dense matrix format"); + } } void read_binary_matrix_header(const string& filename, int& tmp1, int& tmp2, int& tmp3) @@ -1332,6 +1478,109 @@ void save_binary_extfmt(const string &filename, const vector &row_names, jout.close(); } +void save_dense_binary(ofstream& out,const string& row_name,Eigen::VectorXd& data) +{ + if (!out.good()) + { + throw runtime_error("save_dense_binary(): stream not good"); + } + int tmp; + double d; + tmp = row_name.size(); + char *real_name = new char[tmp]; + out.write((char *) &tmp, sizeof(tmp)); + pest_utils::string_to_fortran_char(row_name, real_name, tmp); + out.write(real_name, tmp); + delete[] real_name; + for (int jcol = 0; jcol < data.size(); ++jcol) + { + d = data(jcol); + out.write((char*)&(d), sizeof(d)); + } + + if (!out.good()) + { + throw runtime_error("save_dense_binary(): stream not good"); + } +} + + +void save_dense_binary(ofstream& out,const vector& row_names,Eigen::MatrixXd& data) +{ + if (!out.good()) + { + throw runtime_error("save_dense_binary(): stream not good"); + } + streampos current_pos = out.tellp(); + if (current_pos == 0) + { + throw runtime_error("save_dense_binary(): stream is uninitialized"); + } + int tmp; + double d; + Eigen::VectorXd row; + string name; + for (int irow=0;irow& col_names) +{ + if (!out.good()) + { + throw runtime_error("prep_save_dense_binary(): stream not good"); + } + int tmp = 0; + out.write((char*)&tmp, sizeof(tmp)); + int n_var = col_names.size(); + int n = -1 * n_var; + out.write((char*)&n, sizeof(n)); + out.write((char*)&n, sizeof(n)); + for (vector::const_iterator b = col_names.begin(), e = col_names.end(); + b != e; ++b) + { + string name = pest_utils::lower_cp(*b); + tmp = name.size(); + out.write((char*)&tmp, sizeof(tmp)); + //mx = max(tmp, mx); + } + if (!out.good()) + { + throw runtime_error("prep_save_dense_binary(): stream not good"); + } + for (vector::const_iterator b = col_names.begin(), e = col_names.end(); + b != e; ++b) + { + string name = pest_utils::lower_cp(*b); + char* par_name = new char[name.size()]; + pest_utils::string_to_fortran_char(name, par_name, name.size()); + out.write(par_name, name.size()); + delete[] par_name; + } + if (!out.good()) + { + throw runtime_error("prep_save_dense_binary(): stream not good"); + } +} void save_binary_orgfmt(const string &filename, const vector &row_names, const vector &col_names, const Eigen::SparseMatrix &matrix) { diff --git a/src/libs/common/utilities.h b/src/libs/common/utilities.h index 3e7f14996..f97121725 100644 --- a/src/libs/common/utilities.h +++ b/src/libs/common/utilities.h @@ -299,12 +299,19 @@ class thread_flag // void read_binary_matrix_header(const string& filename, int& tmp1, int& tmp2, int& tmp3); +void read_binary_matrix_header(ifstream& in, int& tmp1, int& tmp2, int& tmp3); +vector read_dense_binary_remaining_row_names(ifstream& in,const vector& col_names); +vector read_dense_binary_col_names(ifstream& in,int n_col); +bool is_dense_binary_matrix(int tmp1, int tmp2, int tmp3); +bool read_dense_binary_records(ifstream& in,int n_records, int n_col,vector& row_names, vector>& rec_vecs); void read_dense_binary(const string& filename, vector& row_names, vector& col_names, Eigen::MatrixXd& matrix); bool read_binary(const string &filename, vector &row_names, vector &col_names, Eigen::SparseMatrix &matrix); bool read_binary(const string &filename, vector &row_names, vector &col_names, Eigen::MatrixXd &matrix); - +void prep_save_dense_binary(ofstream& out,const vector& col_names); +void save_dense_binary(ofstream& out,const vector& row_names,Eigen::MatrixXd& data); +void save_dense_binary(ofstream& out,const string& row_name,Eigen::VectorXd& data); void save_binary(const string &filename, const vector &row_names, const vector &col_names, const Eigen::SparseMatrix &matrix); void save_binary_extfmt(const string &filename,const vector &row_names, const vector &col_names, const Eigen::SparseMatrix &matrix); void save_binary_orgfmt(const string &filename, const vector &row_names, const vector &col_names, const Eigen::SparseMatrix &matrix); diff --git a/src/libs/opt/libopt/libopt.vcxproj b/src/libs/opt/libopt/libopt.vcxproj index f3edffbc8..8118afa5a 100644 --- a/src/libs/opt/libopt/libopt.vcxproj +++ b/src/libs/opt/libopt/libopt.vcxproj @@ -259,26 +259,26 @@ StaticLibrary true - v142 + v143 Unicode StaticLibrary false - v142 + v143 true Unicode StaticLibrary true - v142 + v143 Unicode StaticLibrary false - v142 + v143 true Unicode diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index ea9a87d14..fe1181c24 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -5674,11 +5674,35 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) message(1, "not using prior parameter covariance matrix scaling"); } - if ((oe_base.shape().first > 0) && (cycle != NetPackage::NULL_DA_CYCLE)) - { - message(2, "using previously initialized observation (simulated output) ensemble"); - oe_drawn = false; - } + //check to see if any explicit obs noise options are set + bool reset_to_nonoise = true; + if (ppo->get_passed_args().find("IES_NO_NOISE") != ppo->get_passed_args().end()) + { + reset_to_nonoise = false; + } + else if (!ppo->get_ies_obs_csv().empty()) + reset_to_nonoise = false; + else if (!ppo->get_obscov_filename().empty()) + reset_to_nonoise = false; + else + { + map obs_std = pest_scenario.get_ext_file_double_map("observation data external", "standard_deviation"); + if (obs_std.size() > 0) + reset_to_nonoise = false; + } + if (reset_to_nonoise) + { + ss.str(""); + ss << "NOTE: no obs-noise-specific options have been passed, resetting to `ies_no_noise` to true"; + message(0,ss.str()); + ppo->set_ies_no_noise(true); + } + + + if ((oe_base.shape().first > 0) && (cycle != NetPackage::NULL_DA_CYCLE)) { + message(2, "using previously initialized observation (simulated output) ensemble"); + oe_drawn = false; + } else oe_drawn = initialize_oe(obscov); diff --git a/src/libs/pestpp_common/Jacobian_1to1.cpp b/src/libs/pestpp_common/Jacobian_1to1.cpp index a488d57a6..0b7a2d885 100644 --- a/src/libs/pestpp_common/Jacobian_1to1.cpp +++ b/src/libs/pestpp_common/Jacobian_1to1.cpp @@ -394,7 +394,8 @@ bool Jacobian_1to1::forward_diff(const string &par_name, double base_derivative_ // perturb derivative parameters double incr = derivative_inc(par_name, group_info, base_derivative_val, false); - if (incr == 0.0) return false; + if (incr == 0.0) + return false; new_par_val = new_par[par_name] = base_derivative_val + incr; // try forward derivative out_of_bound_forward = out_of_bounds(new_par, par_info_ptr); diff --git a/src/libs/pestpp_common/MOEA.cpp b/src/libs/pestpp_common/MOEA.cpp index edd96cdf1..9f825f162 100644 --- a/src/libs/pestpp_common/MOEA.cpp +++ b/src/libs/pestpp_common/MOEA.cpp @@ -1337,7 +1337,7 @@ void MOEA::sanity_checks() if ((ppo->get_mou_population_size() < error_min_members) && (population_dv_file.size() == 0)) { ss.str(""); - ss << "population_size < " << error_min_members << ", this is redic, increaing to " << warn_min_members; + ss << "population_size < " << error_min_members << ", this is redic, increasing to " << warn_min_members; warnings.push_back(ss.str()); ppo->set_ies_num_reals(warn_min_members); } @@ -1543,7 +1543,9 @@ void MOEA::queue_chance_runs(ParameterEnsemble& _dp) pest_scenario.get_base_par_tran_seq().ctl2numeric_ip(pars); Observations obs = pest_scenario.get_ctl_observations(); //if this is the first iter and no restart - + ss.str(""); + ss << "queuing chance runs for generation " << iter; + message(1,ss.str()); if (chancepoints == chancePoints::SINGLE) { //dont use the _dp, use the class attr dp and op here @@ -2070,6 +2072,7 @@ void MOEA::initialize() n_adaptive_dvs++; } constraints.initialize(dv_names, numeric_limits::max()); + constraints.initial_report(); if (pest_scenario.get_control_info().noptmax == 0) @@ -2291,7 +2294,7 @@ void MOEA::initialize() { //this can be done, but we need to make sure the appropriate chance restart //args were supplied: base_jacobian or obs_stack - throw_moea_error("chance constraints not yet supported with restart"); + throw_moea_error("chance constraints/objectives not yet supported with restart"); } //since mou reqs strict linking of realization names, let's see if we can find an intersection set @@ -3075,6 +3078,7 @@ void MOEA::initialize_population_schedule() } in.close(); } + ofstream& frec = file_manager.rec_ofstream(); frec << "...population schedule: generation,population size:" << endl; for (int i=0;i& ctl_ord_dec_var_names, double _dbl_ } probit_val = get_probit(); + initialize_chance_schedule(f_rec); //if the std weight options was selected, use it - it overrides all other options if (std_weights) { @@ -1202,10 +1203,10 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl } Eigen::MatrixXd missing_dv_mat = pe.get_eigen(missing, dvnames); Eigen::VectorXd missing_dv_vec, stack_dv_vec; - double dist, min_dist, max_dist; + double dist, min_dist, max_dist,cutoff; string min_real_name,missing_real_name; vector distances; - vector factors,temp; + vector factors,temp, temp2; vector dreal_names; Eigen::MatrixXd shifts; @@ -1242,6 +1243,13 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl for (auto& d : distances) temp.push_back(d / max_dist); max_dist = *max_element(temp.begin(),temp.end()); + temp2 = temp; + sort(temp2.begin(),temp2.end()); + cutoff = temp2[temp2.size()-1]; + if (temp2.size() > 5) + { + cutoff = temp2[4]; + } factors.clear(); factor_sum = 0.0; for (auto& t : temp) @@ -1249,6 +1257,10 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl if (t == 0) factor = 10000.0; + else if (t > cutoff) + { + factor = 0.0; + } else factor = 1.0/t; factor_sum += factor; @@ -2272,11 +2284,92 @@ bool Constraints::should_update_chance(int iter) else return true; } - else if ((iter) % pest_scenario.get_pestpp_options().get_opt_recalc_fosm_every() == 0) + //else if ((iter) % pest_scenario.get_pestpp_options().get_opt_recalc_fosm_every() == 0) + else if (chance_schedule[iter]) return true; return false; } +void Constraints::initialize_chance_schedule(ofstream& frec) +{ + stringstream ss; + chance_schedule.clear(); + + string fname = pest_scenario.get_pestpp_options().get_opt_chance_schedule(); + string line; + vector tokens; + int lcount = 0, gen; + bool should_eval = false; + int recalc_every = pest_scenario.get_pestpp_options().get_opt_recalc_fosm_every(); + for (int i=0;i 0) + { + ss.str(""); + ss << "...reading population schedule from file '" << fname << "' " << endl; + frec << ss.str(); + cout << ss.str(); + ifstream in(fname); + if (in.bad()) + { + throw runtime_error("error opening opt_chance_schedule file '"+fname+"'"); + } + while (getline(in,line)) + { + lcount++; + tokens.clear(); + pest_utils::tokenize(line,tokens,"\t ,"); + if (tokens.size() < 2) + { + ss.str(""); + ss << "opt_chance_schedule file '" << fname << "' line " << lcount << " needs at least two entries"; + throw runtime_error(ss.str()); + } + try + { + gen = stoi(tokens[0]); + } + catch (...) + { + ss.str(""); + ss << "error casting '" << tokens[0] << "' to generation integer on line " << lcount << " in opt_chance_schedule file"; + throw runtime_error(ss.str()); + } + + try + { + should_eval = pest_utils::parse_string_arg_to_bool(tokens[1]); + } + catch (...) + { + ss.str(""); + ss << "error casting '" << tokens[1] << "' to bool on line " << lcount << " in opt_chance_schedule file"; + throw runtime_error(ss.str()); + } + chance_schedule[gen] = should_eval; + } + in.close(); + } + + if (!chance_schedule[0]) + { + ss.str(""); + ss << "...chances must be evaluated during the first generation, resetting chance_schedule" << endl; + frec << ss.str(); + cout << ss.str(); + chance_schedule[0] = true; + } + + frec << "...chance schedule: generation,should_eval_chances:" << endl; + for (int i=0;i status_map, map price_map) diff --git a/src/libs/pestpp_common/constraints.h b/src/libs/pestpp_common/constraints.h index 71c220516..d5c8df702 100644 --- a/src/libs/pestpp_common/constraints.h +++ b/src/libs/pestpp_common/constraints.h @@ -181,6 +181,7 @@ class Constraints double risk; double probit_val; double dbl_max; + map chance_schedule; Covariance obscov; @@ -255,5 +256,9 @@ class Constraints pair,vector> get_working_set(Parameters& par_and_dec_vars, Observations& constraints_sim, bool do_shift, double working_set_tol=0.1); void augment_constraint_mat_with_pi(Mat& mat, vector& pi_names); + + void initialize_chance_schedule(ofstream& frec); + + }; #endif diff --git a/src/libs/pestpp_common/covariance.cpp b/src/libs/pestpp_common/covariance.cpp index 491f31981..3d76a3137 100644 --- a/src/libs/pestpp_common/covariance.cpp +++ b/src/libs/pestpp_common/covariance.cpp @@ -1540,7 +1540,8 @@ void Covariance::from_parameter_bounds(ofstream& frec, const vector &par void Covariance::from_parameter_bounds(Pest &pest_scenario, ofstream& frec) { map par_std = pest_scenario.get_ext_file_double_map("parameter data external", "standard_deviation"); - if (par_std.size() > 0) + const ParameterInfo pi = pest_scenario.get_ctl_parameter_info(); + if (par_std.size() > 0) { frec << "Note: the following parameters have 'standard_deviation' defined - this will be used" << endl; frec << " instead of bounds for the prior parameter covariance matrix : " << endl; @@ -1551,8 +1552,12 @@ void Covariance::from_parameter_bounds(Pest &pest_scenario, ofstream& frec) { if (par_std[pname] <= 0.0) { - frec << "Warning: parameter " << pname << " 'standard_deviation' less than or equal to zero, using bounds instead" << endl; - remove.push_back(pname); + ParameterRec::TRAN_TYPE t = pi.get_parameter_rec_ptr(pname)->tranform_type; + if ((t != ParameterRec::TRAN_TYPE::FIXED) && (t != ParameterRec::TRAN_TYPE::TIED)) { + frec << "Warning: parameter " << pname + << " 'standard_deviation' less than or equal to zero, using bounds instead" << endl; + remove.push_back(pname); + } } else frec << pname << ' ' << par_std[pname] << endl; @@ -1658,6 +1663,7 @@ void Covariance::from_observation_weights(Pest &pest_scenario, ofstream& frec) { map obs_std = pest_scenario.get_ext_file_double_map("observation data external", "standard_deviation"); vector remove; + const ObservationInfo oi = pest_scenario.get_ctl_observation_info(); if (obs_std.size() > 0) { frec << "Note: the following observations have 'standard_deviation' defined - this will be used" << endl; @@ -1666,8 +1672,7 @@ void Covariance::from_observation_weights(Pest &pest_scenario, ofstream& frec) { if (obs_std.find(oname) != obs_std.end()) { - if (obs_std[oname] <= 0.0) - { + if ((obs_std[oname] <= 0.0) && (oi.get_weight(oname) > 0.0)) { frec << "Warning: observation " << oname << " 'standard_deviation' less than or equal to zero, using weight" << endl; remove.push_back(oname); } diff --git a/src/libs/pestpp_common/pest_data_structs.cpp b/src/libs/pestpp_common/pest_data_structs.cpp index 7b3e0ba92..ef28e3568 100644 --- a/src/libs/pestpp_common/pest_data_structs.cpp +++ b/src/libs/pestpp_common/pest_data_structs.cpp @@ -554,18 +554,21 @@ PestppOptions::ARG_STATUS PestppOptions::assign_value_by_key(string key, const s //convert_ip(value, condor_submit_file); condor_submit_file = org_value; } - else if ((key == "SWEEP_PARAMETER_CSV_FILE") || (key == "SWEEP_PAR_CSV")) + else if ((key == "SWEEP_PARAMETER_CSV_FILE") || (key == "SWEEP_PAR_CSV") || (key == "SWEEP_PARAMETER_FILE")) { passed_args.insert("SWEEP_PARAMETER_CSV_FILE"); passed_args.insert("SWEEP_PAR_CSV"); - - //convert_ip(org_value, sweep_parameter_csv_file); + passed_args.insert("SWEEP_PARAMETER_FILE"); + + + //convert_ip(org_value, sweep_parameter_csv_file); sweep_parameter_csv_file = org_value; } - else if ((key == "SWEEP_OUTPUT_CSV_FILE") || (key == "SWEEP_OBS_CSV")) + else if ((key == "SWEEP_OUTPUT_CSV_FILE") || (key == "SWEEP_OBS_CSV") || (key == "SWEEP_OUTPUT_FILE")) { passed_args.insert("SWEEP_OUTPUT_CSV_FILE"); passed_args.insert("SWEEP_OBS_CSV"); + passed_args.insert("SWEEP_OUTPUT_FILE"); //convert_ip(org_value, sweep_output_csv_file); sweep_output_csv_file = org_value; @@ -1482,6 +1485,12 @@ bool PestppOptions::assign_mou_value_by_key(const string& key, const string& val mou_population_schedule = org_value; return true; } + else if (key == "OPT_CHANCE_SCHEDULE") + { + opt_chance_schedule = org_value; + return true; + } + else if (key == "MOU_SIMPLEX_REFLECTIONS") { convert_ip(value, mou_simplex_reflections); @@ -1707,6 +1716,8 @@ void PestppOptions::summary(ostream& os) const os << "opt_iter_tol: " << opt_iter_tol << endl; os << "opt_recalc_fosm_every: " << opt_recalc_fosm_every << endl; os << "opt_chance_points: " << opt_chance_points << endl; + os << "opt_chance_schedule: " << opt_chance_schedule << endl; + @@ -1905,6 +1916,7 @@ void PestppOptions::set_defaults() set_opt_par_stack(""); set_opt_obs_stack(""); set_opt_chance_points("SINGLE"); + set_opt_chance_schedule(""); set_sqp_dv_en(""); @@ -1945,7 +1957,7 @@ void PestppOptions::set_defaults() set_ies_lam_mults(vector()); set_ies_init_lam(0.0); set_ies_use_approx(true); - set_ies_subset_size(4); + set_ies_subset_size(-10); set_ies_reg_factor(0.0); set_ies_verbose_level(1); set_ies_use_prior_scaling(false); diff --git a/src/libs/pestpp_common/pest_data_structs.h b/src/libs/pestpp_common/pest_data_structs.h index 815f6a70e..a569a3232 100644 --- a/src/libs/pestpp_common/pest_data_structs.h +++ b/src/libs/pestpp_common/pest_data_structs.h @@ -363,6 +363,8 @@ class PestppOptions { void set_opt_obs_stack(string _stack) { opt_obs_stack = _stack; } string get_opt_chance_points() const { return opt_chance_points; } void set_opt_chance_points(string chance_points) { opt_chance_points = chance_points; } + string get_opt_chance_schedule() const {return opt_chance_schedule;} + void set_opt_chance_schedule(string fname) {opt_chance_schedule = fname;} string get_sqp_dv_en()const { return sqp_dv_en; } @@ -726,6 +728,7 @@ class PestppOptions { string opt_par_stack; string opt_obs_stack; string opt_chance_points; + string opt_chance_schedule; string sqp_dv_en; string sqp_obs_restart_en; diff --git a/src/libs/pestpp_common/pestpp_common.vcxproj b/src/libs/pestpp_common/pestpp_common.vcxproj index dc8112a3c..515fb8e0a 100644 --- a/src/libs/pestpp_common/pestpp_common.vcxproj +++ b/src/libs/pestpp_common/pestpp_common.vcxproj @@ -110,38 +110,38 @@ StaticLibrary true - v142 + v143 MultiByte StaticLibrary true - v142 + v143 MultiByte StaticLibrary true - v142 + v143 MultiByte DynamicLibrary true - v142 + v143 MultiByte StaticLibrary false - v142 + v143 true MultiByte StaticLibrary false - v142 + v143 true MultiByte No diff --git a/src/libs/run_managers/abstract_base/RunStorage.cpp b/src/libs/run_managers/abstract_base/RunStorage.cpp index 3d38d3da8..3c217bad5 100644 --- a/src/libs/run_managers/abstract_base/RunStorage.cpp +++ b/src/libs/run_managers/abstract_base/RunStorage.cpp @@ -56,7 +56,7 @@ void RunStorage::reset(const vector &_par_names, const vector &_ buf_stream.close(); } - if ((pest_utils::check_exist_in(filename)) || (pest_utils::check_exist_out(filename))) + /*if ((pest_utils::check_exist_in(filename)) || (pest_utils::check_exist_out(filename))) { int flag = remove(filename.c_str()); if (flag != 0) @@ -64,7 +64,7 @@ void RunStorage::reset(const vector &_par_names, const vector &_ throw runtime_error("RunStorage::reset(): error removing existing file '"+filename+"'"); } - } + }*/ buf_stream.open(filename.c_str(), ios_base::out | ios_base::binary); if (!buf_stream) diff --git a/src/libs/run_managers/abstract_base/abstract_base.vcxproj b/src/libs/run_managers/abstract_base/abstract_base.vcxproj index 5913f1ce7..48c57b353 100644 --- a/src/libs/run_managers/abstract_base/abstract_base.vcxproj +++ b/src/libs/run_managers/abstract_base/abstract_base.vcxproj @@ -27,26 +27,26 @@ StaticLibrary true - v142 + v143 MultiByte StaticLibrary true - v142 + v143 MultiByte StaticLibrary false - v142 + v143 true MultiByte StaticLibrary false - v142 + v143 true MultiByte diff --git a/src/libs/run_managers/external/external.vcxproj b/src/libs/run_managers/external/external.vcxproj index bf6b7746e..47b754bc4 100644 --- a/src/libs/run_managers/external/external.vcxproj +++ b/src/libs/run_managers/external/external.vcxproj @@ -27,26 +27,26 @@ StaticLibrary true - v142 + v143 MultiByte StaticLibrary true - v142 + v143 MultiByte StaticLibrary false - v142 + v143 true MultiByte StaticLibrary false - v142 + v143 true MultiByte diff --git a/src/libs/run_managers/serial/serial.vcxproj b/src/libs/run_managers/serial/serial.vcxproj index 1e95a1b42..71bc55e77 100644 --- a/src/libs/run_managers/serial/serial.vcxproj +++ b/src/libs/run_managers/serial/serial.vcxproj @@ -27,26 +27,26 @@ StaticLibrary true - v142 + v143 MultiByte StaticLibrary true - v142 + v143 MultiByte StaticLibrary false - v142 + v143 true MultiByte StaticLibrary false - v142 + v143 true MultiByte diff --git a/src/libs/run_managers/wrappers/wrappers.vcxproj b/src/libs/run_managers/wrappers/wrappers.vcxproj index 627ca80f7..e2918a6b0 100644 --- a/src/libs/run_managers/wrappers/wrappers.vcxproj +++ b/src/libs/run_managers/wrappers/wrappers.vcxproj @@ -27,26 +27,26 @@ StaticLibrary true - v142 + v143 MultiByte StaticLibrary true - v142 + v143 MultiByte StaticLibrary false - v142 + v143 true MultiByte StaticLibrary false - v142 + v143 true MultiByte diff --git a/src/libs/run_managers/yamr/RunManagerPanther.cpp b/src/libs/run_managers/yamr/RunManagerPanther.cpp index 1151269ba..f30305976 100644 --- a/src/libs/run_managers/yamr/RunManagerPanther.cpp +++ b/src/libs/run_managers/yamr/RunManagerPanther.cpp @@ -603,8 +603,8 @@ RunManagerAbstract::RUN_UNTIL_COND RunManagerPanther::run_until(RUN_UNTIL_COND c } if (open_file_trans_streams.size() == 0) { - listen(); - break; + if (!listen()) + break; } int q = pest_utils::quit_file_found(); if ((q == 1) || (q == 2) || (q == 4)) diff --git a/src/libs/run_managers/yamr/yamr.vcxproj b/src/libs/run_managers/yamr/yamr.vcxproj index 9acf27907..237ee409e 100644 --- a/src/libs/run_managers/yamr/yamr.vcxproj +++ b/src/libs/run_managers/yamr/yamr.vcxproj @@ -45,39 +45,39 @@ StaticLibrary true MultiByte - v142 + v143 StaticLibrary true MultiByte - v142 + v143 StaticLibrary true MultiByte - v142 + v143 DynamicLibrary true MultiByte - v142 + v143 StaticLibrary false true MultiByte - v142 + v143 StaticLibrary false true MultiByte - v142 + v143 No diff --git a/src/programs/gsa/gsa.vcxproj b/src/programs/gsa/gsa.vcxproj index c055d8eea..8e0511c4e 100644 --- a/src/programs/gsa/gsa.vcxproj +++ b/src/programs/gsa/gsa.vcxproj @@ -60,75 +60,75 @@ Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application true - v142 + v143 MultiByte Application false - v142 + v143 true MultiByte Application false - v142 + v143 true MultiByte Application false - v142 + v143 true MultiByte diff --git a/src/programs/pest++/pest++.vcxproj b/src/programs/pest++/pest++.vcxproj index 15cb84cdf..0c4bcbc96 100644 --- a/src/programs/pest++/pest++.vcxproj +++ b/src/programs/pest++/pest++.vcxproj @@ -38,28 +38,28 @@ true NotSet Sequential - v142 + v143 Application true NotSet Sequential - v142 + v143 Application true NotSet No - v142 + v143 Application true NotSet Sequential - v142 + v143 Application @@ -67,7 +67,7 @@ true NotSet Sequential - v142 + v143 Application @@ -75,7 +75,7 @@ true NotSet No - v142 + v143 diff --git a/src/programs/pestpp-da/pestpp-da.cpp b/src/programs/pestpp-da/pestpp-da.cpp index 99435de0c..294239480 100644 --- a/src/programs/pestpp-da/pestpp-da.cpp +++ b/src/programs/pestpp-da/pestpp-da.cpp @@ -370,7 +370,7 @@ int main(int argc, char* argv[]) } - if (obs_cycle_info.find(*icycle) != obs_cycle_info.end()) + if ((obs_cycle_info.find(*icycle) != obs_cycle_info.end()) || (weight_cycle_info.find(*icycle) != weight_cycle_info.end())) { performance_log.log_event("updating obs using da obs cycle table info"); map cycle_map = obs_cycle_info[*icycle]; @@ -388,15 +388,22 @@ int main(int argc, char* argv[]) } else { - childPest.get_ctl_observations_4_mod().update_rec(tbl_obs_name, cycle_map[tbl_obs_name]); - //check if this obs is in this cycle's weight info - if (weight_cycle_map.find(tbl_obs_name) != weight_cycle_map.end()) - { - oi.set_weight(tbl_obs_name, weight_cycle_map[tbl_obs_name]); - //pest_scenario.get_observation_info_ptr()->set_weight(tbl_obs_name, weight_cycle_map[tbl_obs_name]); - } + childPest.get_ctl_observations_4_mod().update_rec(tbl_obs_name, cycle_map[tbl_obs_name]); } } + for (auto tbl_obs_name : weights_in_tbl) + { + if (weight_cycle_map.find(tbl_obs_name) == weight_cycle_map.end()) + { + oi.set_weight(tbl_obs_name, 0.0); + } + else + { + oi.set_weight(tbl_obs_name, weight_cycle_map[tbl_obs_name]); + } + } + + childPest.set_observation_info(oi); } @@ -504,6 +511,7 @@ int main(int argc, char* argv[]) pest_scenario.get_pestpp_options().get_additional_ins_delimiters(), pest_scenario.get_pestpp_options().get_num_tpl_ins_threads(), pest_scenario.get_pestpp_options().get_tpl_force_decimal()); + run_manager_ptr->initialize(pest_scenario.get_ctl_parameters(), pest_scenario.get_ctl_observations()); } //generate a parent ensemble which includes all parameters across all cycles @@ -619,7 +627,8 @@ int main(int argc, char* argv[]) } //check for entries in the obs cycle table //cout << childPest.get_ctl_observations().get_rec("GAGE_1") << ", " << pest_scenario.get_observation_info_ptr()->get_weight("GAGE_1") << endl; - if (obs_cycle_info.find(*icycle) != obs_cycle_info.end()) + if ((obs_cycle_info.find(*icycle) != obs_cycle_info.end()) || (false)) + { performance_log.log_event("updating obs using da obs cycle table info"); map cycle_map = obs_cycle_info[*icycle]; diff --git a/src/programs/pestpp-da/pestpp-da.vcxproj b/src/programs/pestpp-da/pestpp-da.vcxproj index 32a526177..4517c233c 100644 --- a/src/programs/pestpp-da/pestpp-da.vcxproj +++ b/src/programs/pestpp-da/pestpp-da.vcxproj @@ -29,26 +29,26 @@ Application true - v142 + v143 Unicode Application false - v142 + v143 true Unicode Application true - v142 + v143 NotSet Application false - v142 + v143 true Unicode diff --git a/src/programs/pestpp-ies/pestpp-ies.vcxproj b/src/programs/pestpp-ies/pestpp-ies.vcxproj index f453d8e77..76b7cf1bf 100644 --- a/src/programs/pestpp-ies/pestpp-ies.vcxproj +++ b/src/programs/pestpp-ies/pestpp-ies.vcxproj @@ -28,26 +28,26 @@ Application true - v142 + v143 Unicode Application false - v142 + v143 true Unicode Application true - v142 + v143 NotSet Application false - v142 + v143 true Unicode diff --git a/src/programs/pestpp-mou/pestpp-mou.vcxproj b/src/programs/pestpp-mou/pestpp-mou.vcxproj index a6611387e..aaaf8fb17 100644 --- a/src/programs/pestpp-mou/pestpp-mou.vcxproj +++ b/src/programs/pestpp-mou/pestpp-mou.vcxproj @@ -28,26 +28,26 @@ Application true - v142 + v143 Unicode Application false - v142 + v143 true Unicode Application true - v142 + v143 NotSet Application false - v142 + v143 true Unicode diff --git a/src/programs/pestpp-opt/pestpp-opt.vcxproj b/src/programs/pestpp-opt/pestpp-opt.vcxproj index b76ca7cc8..a4f7190eb 100644 --- a/src/programs/pestpp-opt/pestpp-opt.vcxproj +++ b/src/programs/pestpp-opt/pestpp-opt.vcxproj @@ -28,26 +28,26 @@ Application true - v142 + v143 Unicode Application false - v142 + v143 true Unicode Application true - v142 + v143 NotSet Application false - v142 + v143 true NotSet diff --git a/src/programs/pestpp-sqp/pestpp-sqp.vcxproj b/src/programs/pestpp-sqp/pestpp-sqp.vcxproj index 16c9d5635..16c5abbc1 100644 --- a/src/programs/pestpp-sqp/pestpp-sqp.vcxproj +++ b/src/programs/pestpp-sqp/pestpp-sqp.vcxproj @@ -28,26 +28,26 @@ Application true - v142 + v143 Unicode Application false - v142 + v143 true Unicode Application true - v142 + v143 NotSet Application false - v142 + v143 true Unicode diff --git a/src/programs/sweep/pestpp-swp.vcxproj b/src/programs/sweep/pestpp-swp.vcxproj index fe7cde36d..b28b6cb8b 100644 --- a/src/programs/sweep/pestpp-swp.vcxproj +++ b/src/programs/sweep/pestpp-swp.vcxproj @@ -28,26 +28,26 @@ Application true - v142 + v143 Unicode Application false - v142 + v143 true Unicode Application true - v142 + v143 NotSet Application false - v142 + v143 true Unicode diff --git a/src/programs/sweep/sweep.cpp b/src/programs/sweep/sweep.cpp index 1b47f6734..4740ed386 100644 --- a/src/programs/sweep/sweep.cpp +++ b/src/programs/sweep/sweep.cpp @@ -116,6 +116,67 @@ map prepare_parameter_csv(Parameters pars, ifstream &csv, bool forgi return header_info; } +map prepare_parameter_dense_binary(Parameters pars, ifstream &in, bool forgive, vector& header_tokens) +{ + stringstream ss; + if (!in.good()) + { + throw runtime_error("ifstream not good in prepare_parameter_dense_binary()"); + } + + //process the header + //any missing header labels will be marked to ignore those columns later + int tmp1,n_col,tmp3; + read_binary_matrix_header(in,tmp1,n_col,tmp3); + if (!is_dense_binary_matrix(tmp1,n_col,tmp3)) + { + throw runtime_error("prepare_parameter_dense_binary() file does not contain a dense binary matrix"); + } + n_col *= -1; + header_tokens.clear(); + header_tokens = read_dense_binary_col_names(in,n_col); + ss.str(""); + ss << "..." << header_tokens.size() << " valid parameter entries found in dense binary parameter file (" + << n_col << " total columns)" << endl; + cout << ss.str(); + // check for parameter names that in the pest control file but that are missing from the csv file + vector missing_names; + string name; + set stokens(header_tokens.begin(),header_tokens.end()); + for (auto &p : pars) + if (stokens.find(p.first) == stokens.end()) + missing_names.push_back(p.first); + + if (missing_names.size() > 0) + { + stringstream ss; + ss << " the following pest control file parameter names were not found in the parameter csv file:" << endl; + for (auto &n : missing_names) ss << n << endl; + if (!forgive) + throw runtime_error(ss.str()); + else + cout << ss.str() << endl << "continuing anyway..." << endl; + } + + //build up a list of idxs to use + vector ctl_pnames = pars.get_keys(); + unordered_set s_pnames(ctl_pnames.begin(), ctl_pnames.end()); + unordered_set::iterator end = s_pnames.end(); + ctl_pnames.resize(0); + vector header_idxs; + map header_info; + for (int i = 0; i < header_tokens.size(); i++) + { + //if (find(ctl_pnames.begin(), ctl_pnames.end(), header_tokens[i]) != ctl_pnames.end()) + if (s_pnames.find(header_tokens[i]) != end) + { + //header_idxs.push_back(i); + header_info[header_tokens[i]] = i; + } + } + return header_info; +} + //pair,vector> load_parameters_from_csv(map &header_info, ifstream &csv, int chunk, const Parameters &ctl_pars, vector &run_ids, vector &sweep_pars) void load_parameters_from_csv(map& header_info, ifstream& csv, int chunk, const Parameters& ctl_pars, vector& run_ids, vector& sweep_pars) @@ -204,86 +265,171 @@ void load_parameters_from_csv(map& header_info, ifstream& csv, int //return pair,vector> (run_ids,sweep_pars); } +void load_parameters_from_dense_binary(map& header_info, ifstream& in, int chunk, const Parameters& ctl_pars, vector& run_ids, vector& sweep_pars) -void prep_sweep_output_file(Pest &pest_scenario, ofstream &csv) { - //ofstream csv(pest_scenario.get_pestpp_options().get_sweep_output_csv_file()); - csv.open(pest_scenario.get_pestpp_options().get_sweep_output_csv_file()); - if (!csv.good()) - { - throw runtime_error("could not open sweep_output_csv_file for writing: " + - pest_scenario.get_pestpp_options().get_sweep_output_csv_file()); - } - csv << setprecision(numeric_limits::digits10); - csv << "run_id,input_run_id,failed_flag"; - csv << ",phi,meas_phi,regul_phi"; - for (auto &ogrp : pest_scenario.get_ctl_ordered_obs_group_names()) - { - csv << ',' << pest_utils::lower_cp(ogrp); - } - for (auto &oname : pest_scenario.get_ctl_ordered_obs_names()) - csv << ',' << pest_utils::lower_cp(oname); - csv << endl; - csv.flush(); - //return csv; - + cout << endl; + //process each parameter value line in the csv file + streampos current_pos = in.tellg(); + in.seekg(0,std::ios::beg); + int tmp1,n_col,tmp3; + read_binary_matrix_header(in,tmp1,n_col,tmp3); + if (!is_dense_binary_matrix(tmp1,n_col,tmp3)) + { + throw runtime_error("prepare_parameter_dense_binary() file does not contain a dense binary matrix"); + } + n_col *= -1; + in.seekg(current_pos,std::ios::beg); + + + vector> rec_vecs; + + run_ids.clear(); + sweep_pars.clear(); + + bool success = read_dense_binary_records(in,chunk,n_col,run_ids,rec_vecs); + if (!success) + { + throw runtime_error("error processing dense binary file parameter records"); + } + sweep_pars.reserve(rec_vecs.size()); + + double val; + string line; + vector tokens,names; + vector vals; + Parameters pars = ctl_pars; + string run_id; + + char c; + + for (int irec=0;irec < run_ids.size();irec++) + { + //for (int i = 0; i < header_tokens.size(); i++) + for (auto hi : header_info) + { + pars.update_rec(hi.first, rec_vecs[irec][hi.second]); + //} + } + //pars.update_without_clear(names, vals); + sweep_pars.push_back(pars); + //sweep_pars.emplace_back(pars); + //run_ids.push_back(run_id); + + if (pars.size() > 10000) + cout << irec << "\r" << flush; + } + //csv.close(); + //return pair,vector> (run_ids,sweep_pars); } +void prep_sweep_output_file(Pest &pest_scenario, ofstream &out, bool& is_binary) +{ + //ofstream out(pest_scenario.get_pestpp_options().get_sweep_output_csv_file()); + vector obs_names = pest_scenario.get_ctl_ordered_obs_names(); + string filename = pest_scenario.get_pestpp_options().get_sweep_output_csv_file(); + string obs_ext = filename.substr(filename.size() - 3, filename.size()); + pest_utils::lower_ip(obs_ext); + if (obs_ext.compare("bin") == 0) + { + out.open(filename,ios::binary); + if (!out.good()) { + throw runtime_error("could not open sweep_output_file for writing: " + + filename); + } + prep_save_dense_binary(out,obs_names); + is_binary = true; + + } + else { + out.open(filename); + -void process_sweep_runs(ofstream &csv, Pest &pest_scenario, RunManagerAbstract* run_manager_ptr, vector run_ids, vector listed_run_ids, ObjectiveFunc obj_func,int total_runs_done) + if (!out.good()) { + throw runtime_error("could not open sweep_output_file for writing: " + + filename); + } + out << setprecision(numeric_limits::digits10); + out << "run_id,input_run_id,failed_flag"; + out << ",phi,meas_phi,regul_phi"; + for (auto &ogrp: pest_scenario.get_ctl_ordered_obs_group_names()) { + out << ',' << pest_utils::lower_cp(ogrp); + } + for (auto &oname: pest_scenario.get_ctl_ordered_obs_names()) + out << ',' << pest_utils::lower_cp(oname); + out << endl; + out.flush(); + is_binary = false; + //return out; + } +} + +void process_sweep_runs(ofstream &out, Pest &pest_scenario, RunManagerAbstract* run_manager_ptr, vector run_ids, vector listed_run_ids, + ObjectiveFunc obj_func,int total_runs_done,bool is_binary_output) { Parameters pars; Observations obs; - double fail_val = -1.0E+10; + double fail_val = -1.0E+100; int run_id; string listed_run_id; + bool success; //for (auto &run_id : run_ids) + + vector obs_names = pest_scenario.get_ctl_ordered_obs_names(); + vector obs_group_names = pest_scenario.get_ctl_ordered_obs_group_names(); + Eigen::VectorXd vec(obs_names.size()); + vec.setConstant(fail_val); for (int i = 0;i get_run(run_id, pars, obs)) - { - PhiData phi_data = obj_func.phi_report(obs, pars, *(pest_scenario.get_regul_scheme_ptr())); - - csv << ",0"; - - csv << ',' << phi_data.total(); - csv << ',' << phi_data.meas; - csv << ',' << phi_data.regul; - for (auto &obs_grp : pest_scenario.get_ctl_ordered_obs_group_names()) - { - csv << ',' << phi_data.group_phi.at(obs_grp); - } - for (auto oname : pest_scenario.get_ctl_ordered_obs_names()) - { - csv << ',' << obs[oname]; - } - csv << endl; - } - //if the run bombed - else - { - csv << ",1"; - csv << ",,,"; - for (auto &ogrp : pest_scenario.get_ctl_ordered_obs_group_names()) - { - csv << ','; - } - for (int i = 0; i < pest_scenario.get_ctl_ordered_obs_names().size(); i++) - { - csv << ',' << fail_val; - } - csv << endl; - } + success = run_manager_ptr->get_run(run_id, pars, obs); + if (is_binary_output) + { + vec.setConstant(fail_val); + if (success) + { + vec = obs.get_data_eigen_vec(obs_names); + } + save_dense_binary(out,listed_run_id,vec); + } + else { + out << run_id + total_runs_done; + out << ',' << listed_run_id; + // if the run was successful + if (success) { + PhiData phi_data = obj_func.phi_report(obs, pars, *(pest_scenario.get_regul_scheme_ptr())); + + out << ",0"; + + out << ',' << phi_data.total(); + out << ',' << phi_data.meas; + out << ',' << phi_data.regul; + for (auto &obs_grp: obs_group_names) { + out << ',' << phi_data.group_phi.at(obs_grp); + } + for (auto &oname: obs_names) { + out << ',' << obs[oname]; + } + out << endl; + } + //if the run bombed + else { + out << ",1"; + out << ",,,"; + for (auto &ogrp: obs_group_names) { + out << ','; + } + for (int i = 0; i < obs_names.size(); i++) { + out << ',' << fail_val; + } + out << endl; + } + } } } - int main(int argc, char* argv[]) { #ifndef _DEBUG @@ -393,7 +539,7 @@ int main(int argc, char* argv[]) if (!restart_flag || save_restart_rec_header) { fout_rec << " pestpp-swp.exe - a parameteric sweep utility" << endl << "for PEST(++) datasets " << endl << endl; - fout_rec << " by the PEST++ developement team" << endl << endl << endl; + fout_rec << " by the PEST++ development team" << endl << endl << endl; fout_rec << endl; fout_rec << endl << endl << "version: " << version << endl; fout_rec << "binary compiled on " << __DATE__ << " at " << __TIME__ << endl << endl; @@ -423,9 +569,9 @@ int main(int argc, char* argv[]) } catch (PestError e) { - cerr << "Error prococessing control file: " << filename << endl << endl; + cerr << "Error processing control file: " << filename << endl << endl; cerr << e.what() << endl << endl; - fout_rec << "Error prococessing control file: " << filename << endl << endl; + fout_rec << "Error processing control file: " << filename << endl << endl; fout_rec << e.what() << endl << endl; fout_rec.close(); throw(e); @@ -469,7 +615,7 @@ int main(int argc, char* argv[]) ifstream par_stream(par_csv_file); if (!par_stream.good()) { - throw runtime_error("could not open parameter csv file " + par_csv_file); + throw runtime_error("could not open parameter sweep file " + par_csv_file); } RunManagerAbstract *run_manager_ptr; @@ -532,10 +678,12 @@ int main(int argc, char* argv[]) if ((par_ext.compare("jcb") == 0) || (par_ext.compare("jco") == 0)) { cout << " --- binary jco-type file detected for par_csv" << endl; + fout_rec << " --- binary jco-type file detected for par_csv" << endl; use_jco = true; Jacobian jco(file_manager); jco.read(par_csv_file); cout << jco.get_matrix_ptr()->rows() << " runs found in binary jco-type file" << endl; + fout_rec << jco.get_matrix_ptr()->rows() << " runs found in binary jco-type file" << endl; //check that the jco is compatible with the control file vector names = jco.get_base_numeric_par_names(); jco_col_names = jco.get_sim_obs_names(); @@ -557,15 +705,35 @@ int main(int argc, char* argv[]) jco_mat = jco.get_matrix(jco.get_sim_obs_names(), pest_scenario.get_ctl_ordered_par_names()).toDense(); } - else + else if (par_ext.compare("csv") == 0) { + cout << " --- csv file detected for par_csv" << endl; + fout_rec << " --- csv file detected for par_csv" << endl; header_info = prepare_parameter_csv(pest_scenario.get_ctl_parameters(), par_stream, pest_scenario.get_pestpp_options().get_sweep_forgive()); } + else if (par_ext.compare("bin")==0) + { + cout << " --- dense binary file detected for par_csv" << endl; + fout_rec << " --- dense binary file detected for par_csv" << endl; + vector col_names; + header_info = prepare_parameter_dense_binary(pest_scenario.get_ctl_parameters(),par_stream, + pest_scenario.get_pestpp_options().get_sweep_forgive(),col_names); + vector all_row_names = read_dense_binary_remaining_row_names(par_stream,col_names); + cout << "..." << all_row_names.size() << " parameter sets found in dense binary file" << endl; + fout_rec << "..." << all_row_names.size() << " parameter sets found in dense binary file" << endl; + + + } + else + { + throw runtime_error("unrecognized parameter sweep input file extension: '"+par_ext+"'"); + } // prepare the output file ofstream obs_stream; - prep_sweep_output_file(pest_scenario,obs_stream); + bool is_binary_output = false; + prep_sweep_output_file(pest_scenario,obs_stream,is_binary_output); int chunk = pest_scenario.get_pestpp_options().get_sweep_chunk(); //pair,vector> sweep_par_info; @@ -609,21 +777,45 @@ int main(int argc, char* argv[]) } else { - try { - performance_log.log_event("starting to read parameter csv file"); - load_parameters_from_csv(header_info, par_stream, chunk, pest_scenario.get_ctl_parameters(), run_ids,sweep_pars); - performance_log.log_event("finished reading parameter csv file"); - } - catch (exception &e) - { - stringstream ss; - ss << "error processing parameter csv file: " << e.what(); - performance_log.log_event(ss.str()); - fout_rec << endl << ss.str() << endl; - fout_rec.close(); + if (par_ext.compare("bin")==0) + { + try { + performance_log.log_event("starting to read parameter dense binary file"); + load_parameters_from_dense_binary(header_info, par_stream, chunk, pest_scenario.get_ctl_parameters(), run_ids,sweep_pars); + performance_log.log_event("finished reading parameter dense binary file"); + } + catch (exception &e) + { + stringstream ss; + ss << "error processing parameter dense binary file: " << e.what(); + performance_log.log_event(ss.str()); + fout_rec << endl << ss.str() << endl; + fout_rec.close(); + + throw runtime_error(ss.str()); + } + + } + else + { + try { + performance_log.log_event("starting to read parameter csv file"); + load_parameters_from_csv(header_info, par_stream, chunk, pest_scenario.get_ctl_parameters(), run_ids,sweep_pars); + performance_log.log_event("finished reading parameter csv file"); + } + catch (exception &e) + { + stringstream ss; + ss << "error processing parameter csv file: " << e.what(); + performance_log.log_event(ss.str()); + fout_rec << endl << ss.str() << endl; + fout_rec.close(); + + throw runtime_error(ss.str()); + } + + } - throw runtime_error(ss.str()); - } } cout << "done" << endl; // if there are no parameters to run, we are done @@ -643,7 +835,6 @@ int main(int argc, char* argv[]) { //Parameters temp = base_trans_seq.active_ctl2model_cp(par); irun_ids.push_back(run_manager_ptr->add_run(base_trans_seq.active_ctl2model_cp(par))); - } //make some runs @@ -653,7 +844,7 @@ int main(int argc, char* argv[]) //process the runs cout << "processing runs..."; //process_sweep_runs(obs_stream, pest_scenario, run_manager_ptr, run_ids, obj_func,total_runs_done); - process_sweep_runs(obs_stream, pest_scenario, run_manager_ptr,irun_ids, run_ids, obj_func, total_runs_done); + process_sweep_runs(obs_stream, pest_scenario, run_manager_ptr,irun_ids, run_ids, obj_func, total_runs_done,is_binary_output); cout << "done" << endl; total_runs_done += run_ids.size(); @@ -679,6 +870,7 @@ int main(int argc, char* argv[]) // clean up obs_stream.close(); + par_stream.close(); delete run_manager_ptr; string case_name = file_manager.get_base_filename();