diff --git a/benchmarks/basic_tests.py b/benchmarks/basic_tests.py index 0ae729ea..fb10f86b 100644 --- a/benchmarks/basic_tests.py +++ b/benchmarks/basic_tests.py @@ -1485,7 +1485,14 @@ def sweep_bin_test(): #def fail_test(): # raise Exception("fail please") +def invest(): + pst = pyemu.Pst(os.path.join("PESTPPTest","PestPilotPointTest.pst")) + + + + if __name__ == "__main__": + invest() #run() #mf6_v5_ies_test() #prep_ends() @@ -1529,7 +1536,7 @@ def sweep_bin_test(): #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-ies.exe"),os.path.join("..","bin","win","pestpp-ies.exe")) #tplins1_test() - fr_timeout_test() + #fr_timeout_test() #mf6_v5_ies_test() #mf6_v5_sen_test() diff --git a/documentation/pestpp_users_guide_v5.2.10.docx b/documentation/pestpp_users_guide_v5.2.12.docx similarity index 71% rename from documentation/pestpp_users_guide_v5.2.10.docx rename to documentation/pestpp_users_guide_v5.2.12.docx index 07c660ec..39ec256d 100644 Binary files a/documentation/pestpp_users_guide_v5.2.10.docx and b/documentation/pestpp_users_guide_v5.2.12.docx differ diff --git a/documentation/pestpp_users_manual.md b/documentation/pestpp_users_manual.md index 747ea558..8b6850cd 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.10 +# Version 5.2.12 PEST++ Development Team -April 2024 +June 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.10](#s1) +- [Version 5.2.12](#s1) - [Acknowledgements](#s2) - [Preface](#s3) - [License](#s4) @@ -3505,12 +3505,16 @@ As is usual practice in the PEST++ suite, if a parameter is designated as log-tr ### 9.1.10 Inequality Observations -PESTPP-IES introduces a special observation type that is not available in PESTPP-GLM or PEST, but resembles constraints supported by PESTPP-OPT. This is the “one way observation” type, which is synonymous with the inequality constraints used in PESTPP-OPT. For observations of this type, a residual is zero unless model outputs are either greater than or less than its “measured” value; the user specifies which of these apply. This reflects the nature of some types of measurements. However, their use is broader than this. “Greater than” and “less than” observations can comprise a powerful mechanism for inserting “soft knowledge” into the history-matching process. +PESTPP-IES introduces a special observation type that is not available in PESTPP-GLM or PEST, but resembles constraints supported by PESTPP-OPT. This is the “one way observation” type, which is synonymous with the inequality constraints used in PESTPP-OPT. For observations of this type, a residual is zero unless model outputs are either greater than or less than its “measured” value; the user specifies which of these apply. This reflects the nature of some types of measurements. However, their use is broader than this. “Greater than” and “less than” observations can comprise a powerful mechanism for inserting “soft knowledge” into the history-matching process, as well as accommodating increasingly imprecise data within the data assimilation process. If an observation belongs to an observation group whose name begins with the string “g\_” or “greater\_”, this observation is a “greater than” observation. No objective function penalty is incurred if the modelled value of the pertinent quantity is greater than its observed value as listed in the “observation data” section of the PEST control file. However, if the model-calculated value is less than the observed value, the objective function penalty is calculated in the usual manner, that is as the squared residual times the squared weight. That is, these observations are treated as quadratic penalty inequality constraints. Similarly, if an observation belongs to an observation group whose name begins with the string “l\_” or “less\_”, then this observation is a “less than” observation. No objective function penalty is incurred if the modelled value of the pertinent quantity is less than the measured value listed in the “observation data” section of the PEST control file. However, if the model-calculated value is greater than the measured value, the objective function penalty is calculated in the usual manner, that is as the squared residual times the squared weight. +As of version 5.2.12, users can describe inequality observations with the version 2 control file by adding an optional “less_than” and/or “greater_than” column to the \* observation data section of the control file. In this usecase, only observations that have a non-null/non-nan value for “lower_bound” or “upper_bound” and that have a non-zero weight are used as inequality observations. Note, and this getting fancy, users can supply both less_than and greater_than value for the same non-zero weighted observation, which then triggers that observation to be treated as a so-called “range observation”, employing a double-sided equality to constrain the acceptable range of that observation. This usecase is synonymous with a uniform error model, in contrast to the standard gaussian error model that is typically used. Note that if less_than and/or greater_than are supplied for a non-zero weighted observation, if the observation group meets the standard inequality definitions (stated above), then the observation value is used preferentially. + +Note that observation noise is never added to inequality type observations and even if noise realizations are supplied through an external file, PESTPP-IES will ignore the noise values. + Users that wish to (very) strongly enforce inequality conditions are advised to avoid using extremely high weights for inequalities, and instead are encouraged to use the “drop_violations” functionality described later. ### 9.1.11 Localization diff --git a/src/libs/common/config_os.h b/src/libs/common/config_os.h index eb876fd0..ace10f9a 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.11"; +#define PESTPP_VERSION "5.2.12"; #if defined(_WIN32) || defined(_WIN64) #define OS_WIN diff --git a/src/libs/common/utilities.cpp b/src/libs/common/utilities.cpp index 3f716824..c49398b9 100644 --- a/src/libs/common/utilities.cpp +++ b/src/libs/common/utilities.cpp @@ -636,14 +636,14 @@ void read_res(string& res_filename, Observations& obs) if (extra.size() > 0) { stringstream ss; - ss << "extra obs found res file...ignoring: "; - int i = 0; + ss << extra.size() << " extra obs found res file...ignoring: "; + /*int i = 0; for (auto& n : extra) { ss << n << ' '; i++; if (i % 5 == 0) ss << endl; - } + }*/ cout << ss.str(); } } diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index aa260072..58d5b1c7 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -2029,7 +2029,10 @@ L2PhiHandler::L2PhiHandler(Pest *_pest_scenario, FileManager *_file_manager, string og; double weight; const ObservationInfo* oi = pest_scenario->get_ctl_observation_info_ptr(); - for (auto &oname : pest_scenario->get_ctl_ordered_obs_names()) + map extfile_gt = pest_scenario->get_ext_file_double_map("observation data external","greater_than"); + map extfile_lt = pest_scenario->get_ext_file_double_map("observation data external","less_than"); + + for (auto &oname : pest_scenario->get_ctl_ordered_obs_names()) { og = oi->get_group(oname); weight = oi->get_weight(oname); @@ -2047,7 +2050,23 @@ L2PhiHandler::L2PhiHandler(Pest *_pest_scenario, FileManager *_file_manager, { gt_obs_names.push_back(oname); } - } + else + { + if ((extfile_gt.find(oname) != extfile_gt.end()) && + (extfile_lt.find(oname) != extfile_lt.end())) + { + double_obs_bounds[oname] = pair(extfile_gt.at(oname),extfile_lt.at(oname)); + } + else if (extfile_gt.find(oname) != extfile_gt.end()) + { + gt_obs_bounds[oname] = extfile_gt.at(oname); + } + else if (extfile_lt.find(oname) != extfile_lt.end()) + { + lt_obs_bounds[oname] = extfile_lt.at(oname); + } + } + } //save the org reg factor and org q vector org_reg_factor = pest_scenario->get_pestpp_options().get_ies_reg_factor(); @@ -2064,9 +2083,12 @@ L2PhiHandler::L2PhiHandler(Pest *_pest_scenario, FileManager *_file_manager, { prepare_csv(file_manager->open_ofile_ext(tag+"phi.actual.csv"), oreal_names); prepare_csv(file_manager->open_ofile_ext(tag+"phi.meas.csv"), oreal_names); + //prepare_csv(file_manager->open_ofile_ext(tag+"phi.noise.csv"), oreal_names); prepare_csv(file_manager->open_ofile_ext(tag+"phi.composite.csv"), oreal_names); prepare_csv(file_manager->open_ofile_ext(tag+"phi.regul.csv"), preal_names); prepare_group_csv(file_manager->open_ofile_ext(tag+"phi.group.csv")); + prepare_lambda_csv(file_manager->open_ofile_ext(tag+"phi.lambda.csv")); + } } @@ -2074,11 +2096,12 @@ L2PhiHandler::L2PhiHandler(Pest *_pest_scenario, FileManager *_file_manager, Eigen::MatrixXd L2PhiHandler::get_obs_resid(ObservationEnsemble &oe, bool apply_ineq) { vector names = oe_base->get_var_names(); - Eigen::MatrixXd resid = oe.get_eigen(vector(),names) - + Eigen::MatrixXd oe_vals = oe.get_eigen(vector(),names); + Eigen::MatrixXd resid = oe_vals - oe_base->get_eigen(oe.get_real_names(), vector()); if (apply_ineq) - apply_ineq_constraints(resid,names); + apply_ineq_constraints(resid,oe_vals, names); return resid; } @@ -2171,9 +2194,10 @@ Eigen::MatrixXd L2PhiHandler::get_obs_resid_subset(ObservationEnsemble &oe, bool real_names = oe.get_real_names(); } vector names = oe.get_var_names(); - Eigen::MatrixXd resid = oe.get_eigen(real_names,vector()) - oe_base->get_eigen(real_names, names); + Eigen::MatrixXd oe_vals = oe.get_eigen(real_names,vector()); + Eigen::MatrixXd resid = oe_vals - oe_base->get_eigen(real_names, names); if (apply_ineq) - apply_ineq_constraints(resid, names); + apply_ineq_constraints(resid, oe_vals,names); return resid; } @@ -2206,10 +2230,39 @@ Eigen::MatrixXd L2PhiHandler::get_actual_obs_resid(ObservationEnsemble &oe) ovals.transposeInPlace(); for (int i = 0; i < resid.rows(); i++) resid.row(i) = oe_vals.row(i) - ovals; - apply_ineq_constraints(resid, act_obs_names); + apply_ineq_constraints(resid, oe_vals, act_obs_names); return resid; } +//Eigen::MatrixXd L2PhiHandler::get_noise_resid(ObservationEnsemble &oe, bool apply_ineq) +//{ +// //vector act_obs_names = pest_scenario->get_ctl_ordered_nz_obs_names(); +// vector act_obs_names = oe_base->get_var_names(); +// +// Eigen::MatrixXd oe_vals = oe.get_eigen(vector(),act_obs_names); +// Eigen::MatrixXd noise_vals = oe_base->get_eigen(oe.get_real_names(), vector()); +// +// Observations obs = pest_scenario->get_ctl_observations(); +// Eigen::MatrixXd ovals = obs.get_data_eigen_vec(act_obs_names); +// ovals.transposeInPlace(); +// for (int i = 0; i < oe_vals.rows(); i++) +// { +// cout << "ovals " << ovals << endl; +// cout << "o1 " << oe_vals.row(i) << endl; +// cout << "n1 " << noise_vals.row(i) << endl; +// +// oe_vals.row(i) = oe_vals.row(i) - ovals; +// noise_vals.row(i) = noise_vals.row(i) - ovals; +// cout << "o2 " << oe_vals.row(i) << endl; +// cout << "n2 " << noise_vals.row(i) << endl << endl; +// } +// +// Eigen::MatrixXd resid = oe_vals - noise_vals; +// if (apply_ineq) +// apply_ineq_constraints(resid,oe_vals, act_obs_names); +// return resid; +//} + Eigen::VectorXd L2PhiHandler::get_q_vector() { @@ -2283,6 +2336,13 @@ void L2PhiHandler::update(ObservationEnsemble & oe, ParameterEnsemble & pe) meas[pv.first] = pv.second.sum(); } + +// noise.clear(); +// map noise_map = calc_noise(oe, q); +// for (auto &pv : noise_map) { +// noise[pv.first] = pv.second.sum(); +// } + if (org_reg_factor != 0.0) { par_group_idx_map.clear(); @@ -2348,6 +2408,14 @@ void L2PhiHandler::update(ObservationEnsemble & oe, ParameterEnsemble & pe, Obse meas[pv.first] = pv.second.sum(); } +// noise.clear(); +// map noise_map = calc_noise(oe, weights); +// for (auto &pv : noise_map) +// { +// noise[pv.first] = pv.second.sum(); +// +// } + if (org_reg_factor != 0.0) { par_group_idx_map.clear(); @@ -2438,6 +2506,8 @@ map* L2PhiHandler::get_phi_map_ptr(L2PhiHandler::phiType pt) return &meas; case L2PhiHandler::phiType::REGUL: return ®ul; + case L2PhiHandler::phiType::NOISE: + return &noise; } throw runtime_error("PhiHandler::get_phi_map() didn't find a phi map..."); } @@ -2454,6 +2524,8 @@ map L2PhiHandler::get_phi_map(L2PhiHandler::phiType pt) return meas; case L2PhiHandler::phiType::REGUL: return regul; + case L2PhiHandler::phiType::NOISE: + return noise; } throw runtime_error("PhiHandler::get_phi_map() didn't find a phi map..."); } @@ -2553,6 +2625,9 @@ string L2PhiHandler::get_summary_string(L2PhiHandler::phiType pt) case L2PhiHandler::phiType::COMPOSITE: typ = "composite"; break; + case L2PhiHandler::phiType::NOISE: + typ = "noise"; + break; } ss << setw(15) << typ << setw(15) << stats["mean"] << setw(15) << stats["std"] << setw(15) << stats["min"] << setw(15) << stats["max"] << endl; return ss.str(); @@ -2663,7 +2738,7 @@ void L2PhiHandler::report_group(bool echo) { ss << " --- observation group phi summary --- " << endl; ss << " (computed using 'actual' phi)" << endl; ss << " (sorted by mean phi)" << endl; - ss << left << setw(len) << "group" << right << setw(10) << "mean" << setw(10) << "std"; + ss << left << setw(len) << "group" << right << setw(6) << "count" << setw(10) << "mean" << setw(10) << "std"; ss << setw(10) << "min" << setw(10) << "max"; ss << setw(10) << "percent" << setw(10) << "std" << endl; //<< setw(10) << "min " << setw(10) << "max " << endl; f << ss.str(); @@ -2680,12 +2755,22 @@ void L2PhiHandler::report_group(bool echo) { sort(pairs.begin(),pairs.end(),cmp_pair); c = 0; + int nzc = 0; string g; for (auto& pair : pairs) { g = pair.first; + nzc = 0; + for (auto& n : oi_ptr->observations) + { + if ((n.second.weight > 0) && (n.second.group == g)) + { + nzc++; + } + } ss.str(""); ss << left << setw(len) << pest_utils::lower_cp(g) << " "; + ss << right << setw(5) << nzc << " "; ss << right << setw(9) << setprecision(3) << mn_map[g] << " "; ss << setw(9) << setprecision(3) << std_map[g] << " "; ss << setw(9) << setprecision(3) << mmn_map[g] << " "; @@ -2756,7 +2841,11 @@ void L2PhiHandler::report(bool echo, bool group_report) s = get_summary_string(L2PhiHandler::phiType::ACTUAL); f << s; if (echo) - cout << s; + cout << s; +// s = get_summary_string(L2PhiHandler::phiType::NOISE); +// f << s; +// if (echo) +// cout << s; } else { @@ -2772,6 +2861,11 @@ void L2PhiHandler::report(bool echo, bool group_report) f << s; if (echo) cout << s; +// s = get_summary_string(L2PhiHandler::phiType::NOISE); +// f << s; +// if (echo) +// cout << s; + } } @@ -2821,12 +2915,15 @@ void L2PhiHandler::write(int iter_num, int total_runs, bool write_group) { write_csv(iter_num, total_runs, file_manager->get_ofstream(tag+"phi.actual.csv"), phiType::ACTUAL,oreal_names); write_csv(iter_num, total_runs, file_manager->get_ofstream(tag+"phi.meas.csv"), phiType::MEAS, oreal_names); - if (pest_scenario->get_pestpp_options().get_ies_reg_factor() != 0.0) + //write_csv(iter_num, total_runs, file_manager->get_ofstream(tag+"phi.noise.csv"), phiType::NOISE, oreal_names); + + if (pest_scenario->get_pestpp_options().get_ies_reg_factor() != 0.0) { write_csv(iter_num, total_runs, file_manager->get_ofstream(tag+"phi.regul.csv"), phiType::REGUL, preal_names); } write_csv(iter_num, total_runs, file_manager->get_ofstream(tag+"phi.composite.csv"), phiType::COMPOSITE, oreal_names); - if (write_group) + + if (write_group) write_group_csv(iter_num, total_runs, file_manager->get_ofstream(tag+"phi.group.csv")); } @@ -2894,6 +2991,16 @@ void L2PhiHandler::write_group_csv(int iter_num, int total_runs, ofstream &csv, } } +void L2PhiHandler::write_lambda(int iteration,int num_reals,double current_lambda,double current_comp_mean_phi, + double current_comp_std_phi,double lambda_mult, + double lambda, double comp_mean_phi,double comp_std_phi) +{ + ofstream& csv = file_manager->get_ofstream(tag+"phi.lambda.csv"); + csv << iteration << ',' << num_reals << ',' << current_lambda << ',' << current_comp_mean_phi << ',', + csv << current_comp_std_phi << ',' << lambda_mult << ','; + csv << lambda << ',' << comp_mean_phi << ',' << comp_std_phi << endl; +} + void L2PhiHandler::prepare_group_csv(ofstream &csv, vector extra) { csv << "iteration,total_runs,obs_realization,par_realization"; @@ -2909,6 +3016,13 @@ void L2PhiHandler::prepare_group_csv(ofstream &csv, vector extra) csv << endl; } +void L2PhiHandler::prepare_lambda_csv(ofstream &csv) +{ + csv << "iteration,num_reals,current_lambda,current_comp_mean_phi,current_comp_std_phi,lambda_scale_fac,lambda,comp_mean_phi,comp_std_phi" << endl; +} + + + vector L2PhiHandler::get_idxs_greater_than(double bad_phi, double bad_phi_sigma, ObservationEnsemble &oe, ObservationEnsemble& weights) { //todo: handle weights ensemble here... @@ -2956,17 +3070,17 @@ vector L2PhiHandler::get_idxs_greater_than(double bad_phi, double bad_phi_s return idxs; } -map L2PhiHandler::get_meas_phi(ObservationEnsemble& oe, Eigen::VectorXd& q_vec) -{ - map meas_map = calc_meas(oe,q_vec); - map meas_phi; - for (auto& m : meas_map) - { - meas_phi[m.first] = m.second.sum(); - } - return meas_phi; - -} +//map L2PhiHandler::get_meas_phi(ObservationEnsemble& oe, Eigen::VectorXd& q_vec) +//{ +// map meas_map = calc_meas(oe,q_vec); +// map meas_phi; +// for (auto& m : meas_map) +// { +// meas_phi[m.first] = m.second.sum(); +// } +// return meas_phi; +// +//} //void PhiThread::work(int thread_id, Eigen::MatrixXd wmat, Eigen::MatrixXd resid, vector oe_real_names, map>& phi_map) void upgrade_thread_function_phi(int id, Eigen::MatrixXd& wmat, Eigen::MatrixXd& resid, vector oe_real_names, map>& phi_map, PhiThread& worker, exception_ptr& eptr) { try { @@ -3171,6 +3285,36 @@ map L2PhiHandler::calc_meas(ObservationEnsemble & oe, E return phi_map; } +//map L2PhiHandler::calc_noise(ObservationEnsemble & oe, Eigen::VectorXd& q_vec) +//{ +// map phi_map; +// Eigen::VectorXd oe_base_vec, oe_vec, diff, w_vec; +// vector act_obs_names = pest_scenario->get_ctl_ordered_nz_obs_names(); +// +// if (act_obs_names.size() == 0) +// { +// for (auto name : oe.get_real_names()) +// phi_map[name] = Eigen::VectorXd(); +// return phi_map; +// } +// vector base_real_names = oe_base->get_real_names(), oe_real_names = oe.get_real_names(); +// set bset(base_real_names.begin(),base_real_names.end()); +// set::iterator end = bset.end(); +// string rname; +// Eigen::MatrixXd resid = get_noise_resid(oe); +// assert(oe_real_names.size() == resid.rows()); +// resid = resid.array().rowwise() * q_vec.transpose().array(); +// resid = resid.array().cwiseProduct(resid.array()); +// for (int i = 0; i L2PhiHandler::calc_meas(ObservationEnsemble & oe, ObservationEnsemble& weights) { map phi_map; @@ -3203,6 +3347,39 @@ map L2PhiHandler::calc_meas(ObservationEnsemble & oe, O } +//map L2PhiHandler::calc_noise(ObservationEnsemble & oe, ObservationEnsemble& weights) +//{ +// map phi_map; +// Eigen::VectorXd oe_base_vec, oe_vec, diff, w_vec; +// vector act_obs_names = pest_scenario->get_ctl_ordered_nz_obs_names(); +// +// if (act_obs_names.size() == 0) +// { +// for (auto name : oe.get_real_names()) +// phi_map[name] = Eigen::VectorXd(); +// return phi_map; +// } +// vector base_real_names = oe_base->get_real_names(), oe_real_names = oe.get_real_names(); +// set bset(base_real_names.begin(),base_real_names.end()); +// set::iterator end = bset.end(); +// string rname; +// Eigen::MatrixXd resid = get_noise_resid(oe); +// Eigen::MatrixXd wmat = weights.get_eigen(oe.get_real_names(),oe_base->get_var_names()); +// assert(oe_real_names.size() == resid.rows()); +// resid = resid.array() * wmat.array(); +// resid = resid.array().cwiseProduct(resid.array()); +// for (int i = 0; i L2PhiHandler::calc_regul(ParameterEnsemble & pe) { map phi_map; @@ -3292,13 +3469,15 @@ vector L2PhiHandler::get_violating_realizations(ObservationEnsemble& oe, } -void L2PhiHandler::apply_ineq_constraints(Eigen::MatrixXd &resid, vector &names) +void L2PhiHandler::apply_ineq_constraints(Eigen::MatrixXd &resid, Eigen::MatrixXd &sim_vals,vector &names) { //vector names = oe_base->get_var_names(); //vector lt_names = get_lt_obs_names(), gt_names = get_gt_obs_names(); //vector act_obs_names = pest_scenario->get_ctl_ordered_nz_obs_names(); - if ((lt_obs_names.size() == 0) && (gt_obs_names.size() == 0)) + if ((lt_obs_names.empty()) && (gt_obs_names.empty()) && + (lt_obs_bounds.empty()) && (gt_obs_bounds.empty()) && + (double_obs_bounds.empty())) return; assert(names.size() == resid.cols()); @@ -3308,7 +3487,9 @@ void L2PhiHandler::apply_ineq_constraints(Eigen::MatrixXd &resid, vector lt_vals[n] = obs.get_rec(n); for (auto &n : gt_obs_names) gt_vals[n] = obs.get_rec(n); - if ((lt_vals.size() == 0) && (gt_vals.size() == 0)) + if ((lt_vals.empty()) && (gt_vals.empty()) && + (lt_obs_bounds.empty()) && (gt_obs_bounds.empty()) && + (double_obs_bounds.empty())) return; map idxs; //for (int i = 0; i < act_obs_names.size(); i++) @@ -3316,10 +3497,27 @@ void L2PhiHandler::apply_ineq_constraints(Eigen::MatrixXd &resid, vector for (int i = 0; i < names.size(); i++) idxs[names[i]] = i; int idx; - double val; - Eigen::VectorXd col; + double val,val2; + Eigen::VectorXd col, scol; + + for (auto const iv : double_obs_bounds) + { + idx = idxs[iv.first]; + col = resid.col(idx); + scol = sim_vals.col(idx); + val = iv.second.first; + val2 = iv.second.second; - for (auto iv : lt_vals) + //cout << resid.col(idx) << endl; + for (int i = 0; i < resid.rows(); i++) + col(i) = ((scol(i) > val) && (scol(i) < val2)) ? 0.0 : col(i); + //cout << resid.col(idx) << endl; + cout << col << endl << endl; + resid.col(idx) = col; + //cout << resid.col(idx) << endl; + } + + for (auto const iv : lt_vals) { idx = idxs[iv.first]; col = resid.col(idx); @@ -3332,7 +3530,25 @@ void L2PhiHandler::apply_ineq_constraints(Eigen::MatrixXd &resid, vector //cout << resid.col(idx) << endl; } - for (auto iv : gt_vals) + + for (auto const iv : lt_obs_bounds) + { + idx = idxs[iv.first]; + col = resid.col(idx); + scol = sim_vals.col(idx); + val = iv.second; + //cout << resid.col(idx) << endl; + for (int i = 0; i < resid.rows(); i++) + col(i) = (scol(i) < val) ? 0.0 : col(i); + //cout << resid.col(idx) << endl; + cout << col << endl << endl; + resid.col(idx) = col; + //cout << resid.col(idx) << endl; + } + + + //Eigen::MatrixXd temp = resid; + for (auto const iv : gt_vals) { idx = idxs[iv.first]; col = resid.col(idx); @@ -3341,9 +3557,21 @@ void L2PhiHandler::apply_ineq_constraints(Eigen::MatrixXd &resid, vector col(i) = (col(i) > 0.0) ? 0.0 : col(i); resid.col(idx) = col; } + + for (auto const iv : gt_obs_bounds) + { + idx = idxs[iv.first]; + col = resid.col(idx); + scol = sim_vals.col(idx); + val = iv.second; + for (int i = 0; i < resid.rows(); i++) + col(i) = (scol(i) > val) ? 0.0 : col(i); + resid.col(idx) = col; + } } + map L2PhiHandler::calc_actual(ObservationEnsemble & oe, Eigen::VectorXd &q_vec) { map phi_map; @@ -3434,16 +3662,18 @@ void ParChangeSummarizer::summarize(ParameterEnsemble &pe, string filename) { update(pe); vector grp_names;// = base_pe_ptr->get_pest_scenario().get_ctl_ordered_par_group_names(); - vector> mean_pairs; - for (auto m : mean_change) - { - mean_pairs.push_back(pair(abs(m.second), m.first)); + vector> pairs; + //for (auto m : mean_change) + for (auto& m : percent_at_lbound) + { + //mean_pairs.push_back(pair(abs(m.second), m.first)); + pairs.push_back(pair(m.second+percent_at_ubound.at(m.first),m.first)); } - sort(mean_pairs.begin(), mean_pairs.end()); + sort(pairs.begin(), pairs.end()); //for (auto m : mean_pairs) - for (int i = mean_pairs.size() - 1; i >= 0; i--) + for (int i = pairs.size() - 1; i >= 0; i--) { - grp_names.push_back(mean_pairs[i].second); + grp_names.push_back(pairs[i].second); } int mxlen = 15; for (auto& g : grp_names) @@ -3489,7 +3719,7 @@ void ParChangeSummarizer::summarize(ParameterEnsemble &pe, string filename) } ss.str(""); - ss << " Note: parameter change summary sorted according to abs 'mean change'." << endl; + ss << " Note: parameter change summary sorted according to percent at bounds." << endl; //ss << " 'n CV decr' is the number of parameters with current CV less " << cv_dec_threshold*100.0 << "% of the initial CV" << endl; ss << " Note: the parameter change statistics implicitly include the effect of " << endl; ss << " realizations that have failed or have been dropped." << endl; @@ -4548,6 +4778,46 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) { message(1, "greater_than inequality defined for observations: ", ph.get_gt_obs_names().size()); } + map t; + t = ph.get_lt_obs_bounds(); + if (!t.empty()) + { + ss.str(""); + ss << "less_than inequality defined through 'less_than' data for observations:" << endl; + for (const auto it : t) + { + ss << it.first << "," << it.second << endl; + } + ss << endl; + message(1,ss.str()); + } + t = ph.get_gt_obs_bounds(); + if (!t.empty()) + { + ss.str(""); + ss << "greater_than inequality defined through 'greater_than' data for observations:" << endl; + for (const auto it : t) + { + ss << it.first << "," << it.second << endl; + } + ss << endl; + message(1,ss.str()); + } + t.clear(); + map> tt = ph.get_double_obs_bounds(); + if (!tt.empty()) + { + ss.str(""); + ss << "double inequality defined through 'greater_than' and 'less_than' data for observations:" << endl; + for (const auto it : tt) + { + ss << it.first << "," << it.second.first << " to " << it.second.second << endl; + } + ss << endl; + message(1,ss.str()); + } + + message(1, "running control file parameter values"); vector failed_idxs = run_ensemble(_pe, _oe,vector(),cycle); @@ -5033,7 +5303,47 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) { message(1, "greater_than inequality defined for observations: ", ph.get_gt_obs_names().size()); } - message(1, "running mean parameter values"); + map t; + t = ph.get_lt_obs_bounds(); + if (!t.empty()) + { + ss.str(""); + ss << "less_than inequality defined through 'less_than' data for observations:" << endl; + for (const auto it : t) + { + ss << it.first << "," << it.second << endl; + } + ss << endl; + message(1,ss.str()); + } + t = ph.get_gt_obs_bounds(); + if (!t.empty()) + { + ss.str(""); + ss << "greater_than inequality defined through 'greater_than' data for observations:" << endl; + for (const auto it : t) + { + ss << it.first << "," << it.second << endl; + } + ss << endl; + message(1,ss.str()); + } + t.clear(); + map> tt = ph.get_double_obs_bounds(); + if (!tt.empty()) + { + ss.str(""); + ss << "double inequality defined through 'greater_than' and 'less_than' data for observations:" << endl; + for (const auto it : tt) + { + ss << it.first << "," << it.second.first << " to " << it.second.second << endl; + } + ss << endl; + message(1,ss.str()); + } + + + message(1, "running mean parameter values"); vector failed_idxs = run_ensemble(_pe, _oe,vector(),cycle); if (failed_idxs.size() != 0) @@ -5198,8 +5508,64 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) { message(1, "greater_than inequality defined for observations: ", ph.get_gt_obs_names().size()); } + map t; + t = ph.get_lt_obs_bounds(); + if (!t.empty()) + { + ss.str(""); + if (t.size() < 50) { + ss << "less_than inequality defined through 'less_than' data for observations:" << endl; + for (const auto it: t) { + ss << it.first << "," << it.second << endl; + } + ss << endl; + } + else + { + ss << "less_than inequality defined through 'less_than' data for " << t.size() << " observations" << endl; + } + message(1,ss.str()); + } + t = ph.get_gt_obs_bounds(); + if (!t.empty()) + { + ss.str(""); + if (t.size() < 50) { + ss << "greater_than inequality defined through 'greater_than' data for observations:" << endl; + for (const auto it: t) { + ss << it.first << "," << it.second << endl; + } + ss << endl; + } + message(1,ss.str()); + } + else + { + ss << "greater_than inequality defined through 'greater_than' data for " << t.size() << " observations" << endl; - ph.update(oe, pe, weights); + } + t.clear(); + map> tt = ph.get_double_obs_bounds(); + if (!tt.empty()) + { + ss.str(""); + if (tt.size() < 50) { + ss << "double inequality defined through 'greater_than' and 'less_than' data for observations:" << endl; + for (const auto it: tt) { + ss << it.first << "," << it.second.first << " to " << it.second.second << endl; + } + ss << endl; + } + else + { + ss << "double inequality defined through 'greater_than' and 'less_than' data for " << tt.size() << " observations" << endl; + + } + message(1,ss.str()); + } + + + ph.update(oe, pe, weights); message(0, "pre-drop initial phi summary"); ph.report(true); @@ -5286,18 +5652,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) } } - if ((ppo->get_ies_phi_fractions_file().size() > 0) || - (ppo->get_obscov_filename().size() > 0) || - (in_conflict.size() > 0)) - { - string filename = file_manager.get_base_filename() + ".adjusted.obs_data.csv"; - ofstream f_obs(filename); - if (f_obs.bad()) - throw_em_error("error opening: " + filename); - output_file_writer.scenario_obs_csv(f_obs); - f_obs.close(); - } drop_bad_reals(pe, oe); @@ -5340,6 +5695,18 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) } message(1, "saved adjusted weight ensemble to ", ss.str()); } + if ((ppo->get_ies_phi_fractions_file().size() > 0) || + (ppo->get_obscov_filename().size() > 0) || + (in_conflict.size() > 0)) + { + string filename = file_manager.get_base_filename() + ".adjusted.obs_data.csv"; + ofstream f_obs(filename); + if (f_obs.bad()) + throw_em_error("error opening: " + filename); + output_file_writer.scenario_obs_csv(f_obs); + f_obs.close(); + + } performance_log->log_event("calc initial phi"); @@ -5885,7 +6252,7 @@ void EnsembleMethod::adjust_weights_single(map>& group_to_ message(1,"original mean phi: ",cur_mean_phi); //if the current is really low, just return and the traps in initialize() will catch it. - if (cur_mean_phi < 1.0e-10) + if (cur_mean_phi < 1.0e-30) { performance_log->log_event("mean phi too low - returning"); return; @@ -6684,6 +7051,9 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto } } + double acc_fac = pest_scenario.get_pestpp_options().get_ies_accept_phi_fac(); + double lam_inc = pest_scenario.get_pestpp_options().get_ies_lambda_inc_fac(); + double lam_dec = pest_scenario.get_pestpp_options().get_ies_lambda_dec_fac(); ObservationEnsemble oe_lam_best(&pest_scenario); bool echo = false; if (verbose_level > 1) @@ -6724,8 +7094,12 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto message(0, "phi summary for lambda, scale fac:", vals, echo); ph.report(echo); + mean = ph.get_mean(L2PhiHandler::phiType::COMPOSITE); std = ph.get_std(L2PhiHandler::phiType::COMPOSITE); + ph.write_lambda(iter,oe_lams[i].shape().first,last_best_lam,last_best_mean, + last_best_std, + scale_vals[i],lam_vals[i],mean,std); if (mean < best_mean) { oe_lam_best = oe_lams[i]; @@ -6741,100 +7115,9 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto return false; } - double acc_fac = pest_scenario.get_pestpp_options().get_ies_accept_phi_fac(); - double lam_inc = pest_scenario.get_pestpp_options().get_ies_lambda_inc_fac(); - double lam_dec = pest_scenario.get_pestpp_options().get_ies_lambda_dec_fac(); - - - /*if (iter <= pest_scenario.get_pestpp_options().get_ies_n_iter_mean()) { - message(1,"processing mean-only upgrade"); - - message(0, "phi summary for best lambda, scale fac: ", vector({ lam_vals[best_idx],scale_vals[best_idx] })); - ph.update(oe_lams[best_idx], pe_lams[best_idx],weights); - ph.report(true,false); - - performance_log->log_event("getting prior parameter ensemble mean-centered anomalies"); - Eigen::MatrixXd anoms = pe_base.get_eigen_anomalies(pe.get_real_names(), pe.get_var_names()); - for (int i = 0; i < pe_lams.size(); i++) - { - if (i == best_idx) - continue; - pe_lams[i] = pe.zeros_like(0); - } - - ParameterEnsemble pe_best = pe_lams[best_idx]; - pe_best.set_fixed_info(pe.get_fixed_info()); - - if (pe_filenames.size() > 0) - { - performance_log->log_event("'ies_upgrades_in_memory' is 'false', loading 'best' parameter ensemble from file '" + pe_filenames[best_idx] + "'"); - - pe_best.from_binary(pe_filenames[best_idx]); - pe_best.transform_ip(ParameterEnsemble::transStatus::NUM); - //remove any failed runs from subset testing - remove_external_pe_filenames(pe_filenames); - } - performance_log->log_event("getting best-lambda parameter ensemble mean vector"); - vector mean_vec = pe_best.get_mean_stl_var_vector(); - Parameters pars = pest_scenario.get_ctl_parameters(); - pe.get_par_transform().ctl2numeric_ip(pars); - pars.update_without_clear(pe_best.get_var_names(),mean_vec); - pe.get_par_transform().numeric2ctl_ip(pars); - ss.str(""); - ss << file_manager.get_base_filename() << "." << iter << ".best.par"; - ofstream pfile(ss.str()); - if (!pfile) - { - throw_em_error("error opening best-lambda mean par file '" + ss.str() + "' for writing"); - } - output_file_writer.write_par(pfile,pars,*pe.get_par_transform().get_offset_ptr(),*pe.get_par_transform().get_scale_ptr()); - pfile.close(); - message(2,"saved best-lambda par ensemble mean vector to ",ss.str()); - pe_best = pe.zeros_like(0); - if (verbose_level > 2) - { - performance_log->log_event("saving best-lambda par ensemble"); - ss.str(""); - ss << file_manager.get_base_filename() << "." << iter << ".best.par"; - if (pest_scenario.get_pestpp_options().get_save_binary()) - { - ss << ".jcb"; - pe_best.to_binary(ss.str()); - } - else - { - ss << ".csv"; - pe_best.to_csv(ss.str()); - } - } - - performance_log->log_event("adding mean to anomalies"); - for (int i = 0; i < mean_vec.size(); i++) - { - anoms.col(i) = anoms.col(i).array() + mean_vec[i]; - } - performance_log->log_event("forming new parameter ensemble of mean-shifted prior realizations"); - ParameterEnsemble new_pe = ParameterEnsemble(&pest_scenario,&rand_gen,anoms,pe.get_real_names(),pe.get_var_names()); - - new_pe.set_trans_status(pe.get_trans_status()); - new_pe.set_fixed_info(pe.get_fixed_info()); - if (pest_scenario.get_pestpp_options().get_ies_enforce_bounds()) { - new_pe.enforce_bounds(performance_log, false); - } - oe_lam_best = oe; //copy - - message(0,"running mean-shifted prior realizations: ",new_pe.shape().first); - ss.str(""); - ss << "iteration:" << iter; - vector temp; - run_ensemble_util(performance_log,frec,new_pe,oe_lam_best,run_mgr_ptr,false,temp,NetPackage::NULL_DA_CYCLE, ss.str()); - pe_lams[best_idx] = new_pe; - //make sure we dont try to process the subset stuff below - local_subset_size = pe.shape().first; - }*/ if ((best_idx != -1) && (use_subset) && (local_subset_size < pe.shape().first)) { @@ -8430,7 +8713,7 @@ vector EnsembleMethod::detect_prior_data_conflict(bool save) use_stat_dist = false; double smn, sstd, omn, ostd, dist; - double sd = pest_scenario.get_pestpp_options().get_ies_pdc_sigma_distance(); + double sd = abs(pest_scenario.get_pestpp_options().get_ies_pdc_sigma_distance()); int oe_nr = oe.shape().first; int oe_base_nr = oe_base.shape().first; Eigen::VectorXd t; diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.h b/src/libs/pestpp_common/EnsembleMethodUtils.h index 3d001c8d..93c06b17 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.h +++ b/src/libs/pestpp_common/EnsembleMethodUtils.h @@ -66,7 +66,7 @@ class L2PhiHandler { public: - enum phiType { MEAS, COMPOSITE, REGUL, ACTUAL }; + enum phiType { MEAS, COMPOSITE, REGUL, ACTUAL, NOISE }; L2PhiHandler() { ; } L2PhiHandler(Pest *_pest_scenario, FileManager *_file_manager, ObservationEnsemble *_oe_base, ParameterEnsemble *_pe_base, @@ -88,6 +88,9 @@ class L2PhiHandler void write(int iter_num, int total_runs, bool write_group = true); void write_group(int iter_num, int total_runs, vector extra); + //csv << "iteration,num_reals,current_lambda,accept_phi,lambda_mult,lambda_scale_fac,lambda,meas_phi" << endl; + void write_lambda(int iteration,int num_reals,double current_lambda,double current_comp_mean_phi,double current_comp_std_phi, + double lambda_mult,double lambda, double comp_mean_phi, double comp_std_phi); vector get_idxs_greater_than(double bad_phi, double bad_phi_sigma, ObservationEnsemble &oe, ObservationEnsemble& weights); Eigen::MatrixXd get_obs_resid(ObservationEnsemble &oe, bool apply_ineq=true); @@ -96,15 +99,19 @@ class L2PhiHandler Eigen::MatrixXd get_par_resid(ParameterEnsemble &pe); Eigen::MatrixXd get_par_resid_subset(ParameterEnsemble &pe,vector real_names=vector()); Eigen::MatrixXd get_actual_obs_resid(ObservationEnsemble &oe); + //Eigen::MatrixXd get_noise_resid(ObservationEnsemble &oe, bool apply_ineq=true); Eigen::VectorXd get_q_vector(); vector get_lt_obs_names() { return lt_obs_names; } vector get_gt_obs_names() { return gt_obs_names; } + map get_lt_obs_bounds() {return lt_obs_bounds;} + map get_gt_obs_bounds() {return gt_obs_bounds;} + map> get_double_obs_bounds() {return double_obs_bounds;} - void apply_ineq_constraints(Eigen::MatrixXd &resid, vector &names); + void apply_ineq_constraints(Eigen::MatrixXd &resid, Eigen::MatrixXd &sim_vals, vector &names); void save_residual_cov(ObservationEnsemble& oe, int iter); - map get_meas_phi(ObservationEnsemble& oe, Eigen::VectorXd& q_vec); + //map get_meas_phi(ObservationEnsemble& oe, Eigen::VectorXd& q_vec); map> get_swr_real_map(ObservationEnsemble& oe, ObservationEnsemble& weights,phiType ptype=phiType::MEAS); @@ -120,10 +127,12 @@ class L2PhiHandler string get_summary_header(); void prepare_csv(ofstream &csv,vector &names); void prepare_group_csv(ofstream &csv, vector extra = vector()); + void prepare_lambda_csv(ofstream &csv); map calc_meas(ObservationEnsemble &oe, Eigen::VectorXd& q_vec); map calc_meas(ObservationEnsemble &oe, ObservationEnsemble& weights); - + //map calc_noise(ObservationEnsemble &oe, ObservationEnsemble& weights); + //map calc_noise(ObservationEnsemble &oe, Eigen::VectorXd& q_vec); map calc_regul(ParameterEnsemble &pe);// , double _reg_fac); map calc_actual(ObservationEnsemble &oe, Eigen::VectorXd& q_vec); map calc_actual(ObservationEnsemble & oe, ObservationEnsemble& weights); @@ -147,9 +156,13 @@ class L2PhiHandler map regul; map composite; map actual; + map noise; vector lt_obs_names; vector gt_obs_names; + map lt_obs_bounds; + map gt_obs_bounds; + map> double_obs_bounds; map> obs_group_idx_map; map> par_group_idx_map; diff --git a/src/libs/pestpp_common/Pest.cpp b/src/libs/pestpp_common/Pest.cpp index 738dc8eb..60b5a5e2 100644 --- a/src/libs/pestpp_common/Pest.cpp +++ b/src/libs/pestpp_common/Pest.cpp @@ -676,7 +676,14 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec) { convert_ip(tokens[0], num_tpl_file); if (tokens.size() >= 5) { - convert_ip(tokens[4], control_info.numcom); + try { + convert_ip(tokens[4], control_info.numcom); + } + catch (...) + { + cout << "WARNING: error parsing '" << tokens[4] <<"' to numcom option...continuing" << endl; + control_info.numcom = 0; + } } else { control_info.numcom = 0; diff --git a/src/libs/pestpp_common/Pest.h b/src/libs/pestpp_common/Pest.h index b82ab652..ef2b089e 100644 --- a/src/libs/pestpp_common/Pest.h +++ b/src/libs/pestpp_common/Pest.h @@ -122,7 +122,7 @@ class Pest { protected: //this is the list of external file cols that have meaning... - set efile_keep_cols{ "standard_deviation", "obsnme","parnme","name", "upper_bound","lower_bound", "cycle", "state_par_link","drop_violations" }; + set efile_keep_cols{ "standard_deviation", "obsnme","parnme","name", "upper_bound","lower_bound", "cycle", "state_par_link","drop_violations","greater_than","less_than" }; int n_adj_par = 0; string prior_info_string; ControlInfo control_info; diff --git a/src/libs/pestpp_common/SVDSolver.cpp b/src/libs/pestpp_common/SVDSolver.cpp index 73794145..0d6c1c09 100644 --- a/src/libs/pestpp_common/SVDSolver.cpp +++ b/src/libs/pestpp_common/SVDSolver.cpp @@ -610,12 +610,17 @@ void SVDSolver::calc_upgrade_vec(double i_lambda, Parameters &prev_frozen_active if (upgrade_active_ctl_pars.get_notnormal_keys().size() > 0) { stringstream ss; - for (auto &n : upgrade_active_ctl_pars.get_notnormal_keys()) + ss << "not normal floating point in upgrade pars :" << endl; + for (auto &n : upgrade_active_ctl_pars.get_notnormal_keys()) ss << n << ','; - string message = "not normal floating point in upgrade pars: " + ss.str(); - file_manager.rec_ofstream() << message << endl; - cout << message; - throw runtime_error(message); + + ss << "(this usually happens when the SVD truncation is too aggressive)" << endl; + ss << "(try increasing eigthresh)" << endl << endl; + + file_manager.rec_ofstream() << ss.str(); + + cout << ss.str(); + throw runtime_error(ss.str()); } performance_log->log_event("commencing check of parameter bounds"); @@ -629,7 +634,7 @@ void SVDSolver::calc_upgrade_vec(double i_lambda, Parameters &prev_frozen_active pair ctl_info = pest_scenario.enforce_par_limits(performance_log, notfrozen_upgrade_active_ctl_pars, base_run_active_ctl_pars, true, true); ss.str(""); - ss << "change limit/bound enforement for lambda " << i_lambda << ": " << ctl_info.first << ", scaling factor: " << ctl_info.second; + ss << "change limit/bound enforcement for lambda " << i_lambda << ": " << ctl_info.first << ", scaling factor: " << ctl_info.second; performance_log->log_event(ss.str()); if (ctl_info.second < 0.01) { diff --git a/src/libs/pestpp_common/constraints.cpp b/src/libs/pestpp_common/constraints.cpp index f90bc1f0..e2958de7 100644 --- a/src/libs/pestpp_common/constraints.cpp +++ b/src/libs/pestpp_common/constraints.cpp @@ -127,6 +127,7 @@ void OptObjFunc::throw_optobjfunc_error(string message) void OptObjFunc::initialize(vector _constraint_names, vector _dv_names) { + stringstream ss; //initialize the objective function obj_func_str = pest_scenario.get_pestpp_options().get_opt_obj_func(); obj_sense = (pest_scenario.get_pestpp_options().get_opt_direction() == 1) ? "minimize" : "maximize"; @@ -143,6 +144,10 @@ void OptObjFunc::initialize(vector _constraint_names, vector _dv use_obj_obs = true; obj_obs = obj_func_str; //check + ss.str(""); + ss << "...objective function defined by observation '" << obj_func_str << "'" << endl; + cout << ss.str(); + f_rec << ss.str(); set names(constraint_names.begin(), constraint_names.end()); if (names.find(obj_obs) != names.end()) { @@ -162,25 +167,39 @@ void OptObjFunc::initialize(vector _constraint_names, vector _dv { if (obj_func_str.size() == 0) { - f_rec << " warning: no ++opt_objective_function-->forming a generic objective function (1.0 coef for each decision var)" << endl; - for (auto& name : dv_names) + f_rec << " note: no ++opt_objective_function-->forming a generic objective function (1.0 coef for each decision var)" << endl; + cout << " note: no ++opt_objective_function-->forming a generic objective function (1.0 coef for each decision var)" << endl; + for (auto& name : dv_names) obj_func_coef_map[name] = 1.0; } //or if it is a prior info equation else if (pest_scenario.get_prior_info().find(obj_func_str) != pest_scenario.get_prior_info().end()) { + ss.str(""); + ss << "...objective function defined by prior information equation '" << obj_func_str << "'" << endl; + cout << ss.str(); + f_rec << ss.str(); obj_func_coef_map = pest_scenario.get_prior_info().get_pi_rec(obj_func_str).get_atom_factors(); //throw_sequentialLP_error("prior-information-based objective function not implemented"); } else { //check if this obj_str is a filename - ifstream if_obj(obj_func_str); + obj_func_str = pest_scenario.get_pestpp_options().get_org_opt_obj_func(); + ss.str(""); + ss << "...objective function defined by 2-column external file '" << obj_func_str << "'" << endl; + cout << ss.str(); + f_rec << ss.str(); + if (!pest_utils::check_exist_in(obj_func_str)) + { + throw_optobjfunc_error("unable to open objective function file '"+obj_func_str+"' for reading"); + } + /*ifstream if_obj(obj_func_str); if (!if_obj.good()) throw_optobjfunc_error("unrecognized ++opt_objective_function arg: " + obj_func_str); - else - obj_func_coef_map = pest_utils::read_twocol_ascii_to_map(obj_func_str); + else*/ + obj_func_coef_map = pest_utils::read_twocol_ascii_to_map(obj_func_str); } @@ -241,7 +260,7 @@ void Constraints::initialize(vector& ctl_ord_dec_var_names, double _dbl_ std_weights = pest_scenario.get_pestpp_options().get_opt_std_weights(); if ((!std_weights) && ((stack_size > 0) || (par_stack_name.size() > 0) || (obs_stack_name.size() > 0))) use_fosm = false; - //initialize the stack constainers (ensemble class instances) + //initialize the stack containers (ensemble class instances) stack_pe.set_pest_scenario(&pest_scenario); stack_pe.set_rand_gen(&rand_gen); stack_oe.set_pest_scenario_ptr(&pest_scenario); @@ -786,7 +805,10 @@ void Constraints::initialize(vector& ctl_ord_dec_var_names, double _dbl_ missing.push_back(name); if (missing.size() > 0) throw_constraints_error("obs stack missing the following constraints: ", missing); - + if (stack_size == 0) + { + stack_size = stack_oe.shape().first; + } //a restart with "all" if ((stack_pe_map.size() > 0) && (pest_scenario.get_pestpp_options().get_opt_chance_points() == "ALL")) { @@ -1066,11 +1088,11 @@ double Constraints::get_max_constraint_change(Observations& current_obs, Observa return max_abs_constraint_change; } -//Observations Constraints::get_chance_shifted_constraints() +//Observations Constraints::get_stack_shifted_chance_constraints() //{ // /* one version of this method that doesnt take any args just use the pointer to the // current sim constraint values*/ -// return get_chance_shifted_constraints(*current_constraints_sim_ptr); +// return get_stack_shifted_chance_constraints(*current_constraints_sim_ptr); //} double Constraints::get_sum_of_violations(Parameters& pars, Observations& obs) @@ -1189,7 +1211,7 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl { real_vec = shifted_oe.get_real_vector(real_name); sim.update_without_clear(onames, real_vec); - sim_shifted = get_chance_shifted_constraints(sim, stack_oe_map[real_name],risk_map[real_name], true); + sim_shifted = get_stack_shifted_chance_constraints(sim, stack_oe_map[real_name], risk_map[real_name],true, false); stack_oe_map[real_name].fill_moment_maps(stack_mean,stack_std); shifted_oe.replace(real_map[real_name], sim_shifted); for (auto& constraint : ctl_ord_obs_constraint_names) { @@ -1307,7 +1329,7 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl continue; } //this call uses the class stack_oe attribute; - sim_shifted = get_chance_shifted_constraints(sim, stack_oe_map[dreal_names[ii]],_risk,true); + sim_shifted = get_stack_shifted_chance_constraints(sim, stack_oe_map[dreal_names[ii]], _risk, true, false); shifts.row(ii) = sim_shifted.get_data_eigen_vec(shifted_oe.get_var_names()) * factors[ii]; } real_vec = shifts.colwise().sum(); @@ -1333,7 +1355,7 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl real_vec = shifted_oe.get_real_vector(missing[i]); sim.update_without_clear(onames, real_vec); //this call uses the class stack_oe attribute; - sim_shifted = get_chance_shifted_constraints(sim, stack_oe_map[min_real_name],_risk,true); + sim_shifted = get_stack_shifted_chance_constraints(sim, stack_oe_map[min_real_name], _risk, true, false); //shifted_oe.replace(real_map.at(missing[i]), sim_shifted); stack_oe_map[min_real_name].fill_moment_maps(stack_mean,stack_std); for (auto& constraint : ctl_ord_obs_constraint_names) { @@ -1360,7 +1382,7 @@ Observations Constraints::get_chance_shifted_constraints(Observations& current_o -Observations Constraints::get_chance_shifted_constraints(Observations& current_obs, double _risk) +Observations Constraints::get_chance_shifted_constraints(Observations& current_obs, double _risk, bool use_stack_anomalies) { /* get the simulated constraint values with the chance shift applied*/ ofstream& f_rec = file_mgr_ptr->rec_ofstream(); @@ -1415,11 +1437,11 @@ Observations Constraints::get_chance_shifted_constraints(Observations& current_o if (stack_oe_map.size() > 0) { ObservationEnsemble _mean_stack = get_stack_mean(stack_oe_map); - shifted_obs = get_chance_shifted_constraints(current_obs, _mean_stack, _risk); + shifted_obs = get_stack_shifted_chance_constraints(current_obs, _mean_stack, _risk, false, false); } else - shifted_obs = get_chance_shifted_constraints(current_obs, stack_oe, _risk); + shifted_obs = get_stack_shifted_chance_constraints(current_obs, stack_oe, _risk,false, use_stack_anomalies); } vector names = shifted_obs.get_keys(); @@ -1506,7 +1528,8 @@ ObservationEnsemble Constraints::get_stack_mean(map return _mean_stack; } -Observations Constraints::get_chance_shifted_constraints(Observations& current_obs, ObservationEnsemble& _stack_oe, double _risk, bool full_obs) +Observations Constraints::get_stack_shifted_chance_constraints(Observations& current_obs, ObservationEnsemble& _stack_oe, + double _risk, bool full_obs, bool use_stack_anomalies) { double old_constraint_val, required_val, pt_offset, new_constraint_val; @@ -1517,7 +1540,7 @@ Observations Constraints::get_chance_shifted_constraints(Observations& current_o if (_stack_oe.shape().first < 3) throw_constraints_error("too few (<3) stack members, cannot continue with stack-based chance constraints/objectives"); int cur_num_reals = _stack_oe.shape().first; - //get the inde (which realization number) represents the risk value according to constraint sense + //get the index (which realization number) represents the risk value according to constraint sense // //the "less than" realization index is just the risk value times the number of reals int lt_idx = int(_risk * cur_num_reals); @@ -1527,14 +1550,21 @@ Observations Constraints::get_chance_shifted_constraints(Observations& current_o int eq_idx = int(0.5 * cur_num_reals); //get the mean-centered anomalies - we want to subtract off the mean //in case these stack values are being re-used from a previous iteration - Eigen::MatrixXd anom = _stack_oe.get_eigen_anomalies(); - //get the map of realization name to index location in the stack + //Eigen::MatrixXd anom = _stack_oe.get_eigen_anomalies(); + Eigen::MatrixXd anom; + if (use_stack_anomalies) + anom = _stack_oe.get_eigen_anomalies(); + else + anom = _stack_oe.get_eigen(); + + //get the map of realization name to index location in the stack _stack_oe.update_var_map(); map var_map = _stack_oe.get_var_map(); - //get the mean and stdev summary containters, summarized by observation (e.g. constraint) name + //get the mean and stdev summary containers, summarized by observation (e.g. constraint) name //pair, map> mm = _stack_oe.get_moment_maps(); map mean_map, std_map; _stack_oe.fill_moment_maps(mean_map, std_map); + Eigen::VectorXd cvec; for (auto& name : ctl_ord_obs_constraint_names) { old_constraint_val = current_obs.get_rec(name); @@ -1543,10 +1573,18 @@ Observations Constraints::get_chance_shifted_constraints(Observations& current_o // the realized values of this stack are the anomalies added to the //current constraint value - this assumes the current value //is the mean of the stack distribution - Eigen::VectorXd cvec = anom.col(var_map[name]).array() + old_constraint_val; + if (use_stack_anomalies) + cvec = anom.col(var_map[name]).array() + old_constraint_val; + else + cvec = anom.col(var_map[name]).array();// + old_constraint_val; //now sort the anomolies + current (mean) value vector sort(cvec.data(), cvec.data() + cvec.size()); - //set the stdev container info - this isnt used in + /*cout << cvec << endl << endl; + Eigen::VectorXd temp = _stack_oe.get_var_vector(name); + sort(temp.data(),temp.data() + temp.size()); + cout << temp << endl << endl;*/ + + //set the stdev container info - this isnt used in //calculations but gets reported prior_constraint_stdev[name] = std_map[name]; post_constraint_stdev[name] = std_map[name]; @@ -1607,7 +1645,7 @@ vector Constraints::get_constraint_residual_vec(Observations& sim) return residuals_vec; } -pair,vector> Constraints::get_constraint_bound_vectors(Parameters& current_pars, Observations& current_obs) +pair,vector> Constraints::get_constraint_bound_vectors(Parameters& current_pars, Observations& current_obs, bool use_stack_anomalies) { /* get the upper and lower bound constraint vectors. For less than constraints, the lower bound is set to double max, for greater than constraints, the upper bound is set to double max. @@ -1615,7 +1653,7 @@ pair,vector> Constraints::get_constraint_bound_vectors(Pa vector residuals; if (use_chance) { - Observations current_constraints_chance = get_chance_shifted_constraints(current_obs); + Observations current_constraints_chance = get_chance_shifted_constraints(current_obs, risk, use_stack_anomalies); residuals = get_constraint_residual_vec(current_constraints_chance); } else @@ -2181,7 +2219,7 @@ void Constraints::stack_summary(int iter, Observations& shifted_obs, bool echo, if (!use_chance) return; - //Observations current_constraints_chance = get_chance_shifted_constraints(current_obs); + //Observations current_constraints_chance = get_stack_shifted_chance_constraints(current_obs); int nsize = 20; for (auto o : ctl_ord_obs_constraint_names) diff --git a/src/libs/pestpp_common/constraints.h b/src/libs/pestpp_common/constraints.h index d5c8df70..0b22c535 100644 --- a/src/libs/pestpp_common/constraints.h +++ b/src/libs/pestpp_common/constraints.h @@ -90,12 +90,11 @@ class Constraints //void risk_shift_obs_en_ip(ObservationEnsemble& oe); //get the (risk-shifted) distance to constraints (upper and lower bounds) - pair,vector> get_constraint_bound_vectors(Parameters& current_pars, Observations ¤t_obs); + pair,vector> get_constraint_bound_vectors(Parameters& current_pars, Observations ¤t_obs, bool use_stack_anomalies); //setters void set_jco(Jacobian_1to1& _jco) { jco = _jco; } - void set_initial_constraints_sim(Observations _initial_constraints_sim) { initial_constraints_sim = _initial_constraints_sim; } - + //get fosm-based parameter names vector get_fosm_par_names(); @@ -151,10 +150,10 @@ class Constraints static pair get_sense_from_group_name(const string& name); //get risk-shifted simulated constraint values using current_constraints_sim_ptr - //Observations get_chance_shifted_constraints(); + //Observations get_stack_shifted_chance_constraints(); //get risk-shifted simulated constraint values using _constraints_sim arg Observations get_chance_shifted_constraints(Observations& _constraints_sim); - Observations get_chance_shifted_constraints(Observations& _constraints_sim, double _risk); + Observations get_chance_shifted_constraints(Observations& _constraints_sim, double _risk, bool use_stack_anomalies=true); ObservationEnsemble get_chance_shifted_constraints(ParameterEnsemble& pe, ObservationEnsemble& oe, int gen, string risk_obj = string(), string opt_member=string()); @@ -194,7 +193,7 @@ class Constraints map stack_oe_map; map stack_pe_map; - Observations get_chance_shifted_constraints(Observations& current_obs, ObservationEnsemble& _stack_oe, double _risk, bool full_obs = false); + Observations get_stack_shifted_chance_constraints(Observations& current_obs, ObservationEnsemble& _stack_oe, double _risk, bool full_obs, bool use_stack_anomalies); PriorInformation* null_prior = new PriorInformation(); PriorInformation constraints_pi; diff --git a/src/libs/pestpp_common/pest_data_structs.cpp b/src/libs/pestpp_common/pest_data_structs.cpp index 11de6594..3bea7f89 100644 --- a/src/libs/pestpp_common/pest_data_structs.cpp +++ b/src/libs/pestpp_common/pest_data_structs.cpp @@ -679,6 +679,7 @@ PestppOptions::ARG_STATUS PestppOptions::assign_value_by_key(string key, const s passed_args.insert("OPT_OBJ_FUNC"); passed_args.insert("OPT_OBJECTIVE_FUNCTION"); convert_ip(value,opt_obj_func); + org_opt_obj_func = org_value; } else if (key == "OPT_COIN_LOG") { @@ -1185,6 +1186,8 @@ bool PestppOptions::assign_ies_value_by_key(const string& key, const string& val convert_ip(value,ies_n_iter_mean); return true; } + + return false; } @@ -1899,6 +1902,7 @@ void PestppOptions::set_defaults() set_opt_obj_func(""); + set_org_opt_obj_func(""); set_opt_coin_log(true); set_opt_skip_final(false); set_opt_std_weights(false); @@ -2007,6 +2011,7 @@ void PestppOptions::set_defaults() set_ies_phi_factors_by_real(false); set_ies_n_iter_mean(0); + // DA parameters //set_da_use_ies(false); set_da_par_cycle_table(""); diff --git a/src/libs/pestpp_common/pest_data_structs.h b/src/libs/pestpp_common/pest_data_structs.h index 0236ef2e..0201bf37 100644 --- a/src/libs/pestpp_common/pest_data_structs.h +++ b/src/libs/pestpp_common/pest_data_structs.h @@ -329,7 +329,10 @@ class PestppOptions { string get_opt_obj_func()const { return opt_obj_func; } void set_opt_obj_func(string _opt_obj_func) { opt_obj_func = _opt_obj_func; } - bool get_opt_coin_log()const { return opt_coin_log; } + string get_org_opt_obj_func()const { return org_opt_obj_func; } + void set_org_opt_obj_func(string _org_opt_obj_func) { org_opt_obj_func = _org_opt_obj_func; } + + bool get_opt_coin_log()const { return opt_coin_log; } void set_opt_coin_log(bool _log) { opt_coin_log = _log; } bool get_opt_skip_final()const { return opt_skip_final; } void set_opt_skip_final(bool _skip_final) { opt_skip_final = _skip_final; } @@ -712,6 +715,7 @@ class PestppOptions { bool de_dither_f; string opt_obj_func; + string org_opt_obj_func; bool opt_coin_log; bool opt_skip_final; vector opt_dec_var_groups; diff --git a/src/libs/pestpp_common/sequential_lp.cpp b/src/libs/pestpp_common/sequential_lp.cpp index cf30b41f..fa9a3db8 100644 --- a/src/libs/pestpp_common/sequential_lp.cpp +++ b/src/libs/pestpp_common/sequential_lp.cpp @@ -224,7 +224,7 @@ pair sequentialLP::postsolve_decision_var_report(Parameters &upgr } stringstream ss; ss << slp_iter << ".par"; - + par_trans.active_ctl2ctl_ip(actual_pars); of_wr.write_par(file_mgr_ptr->open_ofile_ext(ss.str()),actual_pars,*par_trans.get_offset_ptr(),*par_trans.get_scale_ptr()); file_mgr_ptr->close_file(ss.str()); of_wr.write_par(file_mgr_ptr->open_ofile_ext("par"), actual_pars, *par_trans.get_offset_ptr(), *par_trans.get_scale_ptr()); @@ -510,10 +510,27 @@ void sequentialLP::iter_solve() //convert Jacobian_1to1 to CoinPackedMatrix cout << " --- forming LP model --- " << endl; CoinPackedMatrix matrix = jacobian_to_coinpackedmatrix(); - + stringstream ss; build_dec_var_bounds(); + bool use_stack_anamolies = true; + if ((constraints.get_use_chance()) && (!constraints.get_use_fosm())) + { + if (constraints.should_update_chance(slp_iter)) + { + ss.str(""); + ss << "...using direct stack simulated results in chance calculations"; + use_stack_anamolies = false; + } + else + { + ss.str(""); + ss << "...using stack anomalies and current simulated constraint values in chance calculations"; + } + f_rec << ss.str() << endl; + cout << ss.str() << endl; - pair, vector> bounds = constraints.get_constraint_bound_vectors(current_pars, current_constraints_sim); + } + pair, vector> bounds = constraints.get_constraint_bound_vectors(current_pars, current_constraints_sim,use_stack_anamolies); constraints.presolve_report(slp_iter,current_pars, current_constraints_sim); //load the linear simplex model //model.loadProblem(matrix, dec_var_lb, dec_var_ub, ctl_ord_obj_func_coefs, constraint_lb, constraint_ub); @@ -531,7 +548,7 @@ void sequentialLP::iter_solve() //if maximum ++opt_coin_loglev, then also write iteration specific mps files if (pest_scenario.get_pestpp_options().get_opt_coin_log()) { - stringstream ss; + ss.str(""); ss << slp_iter << ".mps"; string mps_name = file_mgr_ptr->build_filename(ss.str()); model.writeMps(mps_name.c_str(),0,1); @@ -1172,7 +1189,7 @@ void sequentialLP::iter_presolve() stringstream ss; ss << slp_iter << ".jcb"; string rspmat_file = file_mgr_ptr->build_filename(ss.str()); - f_rec << endl << "saving iteration " << slp_iter << " reponse matrix to file: " << rspmat_file << endl; + f_rec << endl << "saving iteration " << slp_iter << " response matrix to file: " << rspmat_file << endl; jco.save(ss.str()); //check for failed runs diff --git a/src/libs/run_managers/yamr/RunManagerPanther.cpp b/src/libs/run_managers/yamr/RunManagerPanther.cpp index f3030597..5aa36051 100644 --- a/src/libs/run_managers/yamr/RunManagerPanther.cpp +++ b/src/libs/run_managers/yamr/RunManagerPanther.cpp @@ -612,7 +612,7 @@ RunManagerAbstract::RUN_UNTIL_COND RunManagerPanther::run_until(RUN_UNTIL_COND c cout << "'pest.stp' found" << endl; kill_all_active_runs(); ss.str(""); - ss << "'pest.stp' found, closing all reminaing open file transfers"; + ss << "'pest.stp' found, closing all remaining open file transfers"; report(ss.str(),true); for (auto& m : open_file_socket_map) { string fname = m.second; diff --git a/src/programs/gsa/main.cpp b/src/programs/gsa/main.cpp index 84e204e6..b1e250c6 100644 --- a/src/programs/gsa/main.cpp +++ b/src/programs/gsa/main.cpp @@ -94,7 +94,9 @@ int main(int argc, char* argv[]) } catch (PestError e) { - cerr << "Error prococessing control file: " << ctl_file << endl << endl; + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; + cerr << "Error prococessing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw(e); } diff --git a/src/programs/pest++/pest++.cpp b/src/programs/pest++/pest++.cpp index 4ead90de..5920850a 100644 --- a/src/programs/pest++/pest++.cpp +++ b/src/programs/pest++/pest++.cpp @@ -104,8 +104,10 @@ int main(int argc, char* argv[]) ctl_file = file_manager.build_filename("pst"); yam_agent.process_ctl_file(ctl_file); } - catch (PestError e) + catch (exception &e) { + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; cerr << "Error prococessing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw(e); @@ -203,9 +205,11 @@ int main(int argc, char* argv[]) performance_log.log_event("finished processing control file"); #ifndef _DEBUG } - catch (PestError e) + catch (exception &e) { - cerr << "Error prococessing control file: " << filename << endl << endl; + fout_rec << "Error processing control file: " << filename << endl << endl; + fout_rec << e.what() << endl << endl; + cerr << "Error processing control file: " << filename << endl << endl; cerr << e.what() << endl << endl; throw(e); } diff --git a/src/programs/pestpp-da/pestpp-da.cpp b/src/programs/pestpp-da/pestpp-da.cpp index 29423948..13c4bef0 100644 --- a/src/programs/pestpp-da/pestpp-da.cpp +++ b/src/programs/pestpp-da/pestpp-da.cpp @@ -123,8 +123,10 @@ int main(int argc, char* argv[]) yam_agent.process_ctl_file(ctl_file); } - catch (PestError e) + catch (exception &e) { + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; cerr << "Error processing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw(e); @@ -197,13 +199,13 @@ int main(int argc, char* argv[]) pest_scenario.assign_da_cycles(fout_rec); performance_log.log_event("finished processing control file"); } - catch (PestError e) + catch (exception &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); } pest_scenario.check_inputs(fout_rec); diff --git a/src/programs/pestpp-ies/pestpp-ies.cpp b/src/programs/pestpp-ies/pestpp-ies.cpp index cce7cdb8..59031a29 100644 --- a/src/programs/pestpp-ies/pestpp-ies.cpp +++ b/src/programs/pestpp-ies/pestpp-ies.cpp @@ -81,7 +81,9 @@ int main(int argc, char* argv[]) yam_agent.process_ctl_file(ctl_file); } - catch (PestError e) { + catch (exception &e) { + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; cerr << "Error processing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw (e); @@ -154,8 +156,18 @@ int main(int argc, char* argv[]) Pest pest_scenario; //try { performance_log.log_event("starting to process control file"); - pest_scenario.process_ctl_file(file_manager.open_ifile_ext("pst"), file_manager.build_filename("pst"), - fout_rec); + try { + pest_scenario.process_ctl_file(file_manager.open_ifile_ext("pst"), file_manager.build_filename("pst"), + fout_rec); + } + catch (exception &e) + { + fout_rec << "Error processing control file: " << file_manager.build_filename("pst") << endl << endl; + fout_rec << e.what() << endl << endl; + cerr << "Error processing control file: " << file_manager.build_filename("pst") << endl << endl; + cerr << e.what() << endl << endl; + throw(e); + } file_manager.close_file("pst"); performance_log.log_event("finished processing control file"); /*} @@ -308,6 +320,7 @@ int main(int argc, char* argv[]) catch (exception &e) { cout << "Error condition prevents further execution: " << endl << e.what() << endl; + //cout << "press enter to continue" << endl; //char buf[256]; //OperSys::gets_s(buf, sizeof(buf)); diff --git a/src/programs/pestpp-mou/pestpp-mou.cpp b/src/programs/pestpp-mou/pestpp-mou.cpp index 23e2a3e0..c294a684 100644 --- a/src/programs/pestpp-mou/pestpp-mou.cpp +++ b/src/programs/pestpp-mou/pestpp-mou.cpp @@ -100,6 +100,8 @@ int main(int argc, char* argv[]) } catch (PestError e) { + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; cerr << "Error processing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw(e); @@ -180,13 +182,13 @@ int main(int argc, char* argv[]) file_manager.close_file("pst"); performance_log.log_event("finished processing control file"); } - catch (PestError e) + catch (exception &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); } pest_scenario.check_inputs(fout_rec); diff --git a/src/programs/pestpp-opt/pestpp-opt.cpp b/src/programs/pestpp-opt/pestpp-opt.cpp index b548215e..8d7ed52a 100644 --- a/src/programs/pestpp-opt/pestpp-opt.cpp +++ b/src/programs/pestpp-opt/pestpp-opt.cpp @@ -96,6 +96,8 @@ int main(int argc, char* argv[]) } catch (PestError e) { + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; cerr << "Error prococessing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw(e); @@ -168,9 +170,9 @@ int main(int argc, char* argv[]) performance_log.log_event("finished processing control file"); #ifndef _DEBUG } - catch (PestError e) + catch (exception &e) { - cerr << "Error prococessing control file: " << filename << endl << endl; + cerr << "Error processing control file: " << filename << endl << endl; cerr << e.what() << endl << endl; throw(e); } @@ -363,12 +365,24 @@ int main(int argc, char* argv[]) } fout_rec << " ----- Starting Optimization Iterations ---- " << endl << endl; - - Covariance parcov; - parcov.try_from(pest_scenario, file_manager); - sequentialLP slp(pest_scenario, run_manager_ptr, parcov, &file_manager, output_file_writer, performance_log); - slp.solve(); - fout_rec << "Number of forward model runs performed during optimization: " << run_manager_ptr->get_total_runs() << endl; + try { + Covariance parcov; + parcov.try_from(pest_scenario, file_manager); + sequentialLP slp(pest_scenario, run_manager_ptr, parcov, &file_manager, output_file_writer, + performance_log); + slp.solve(); + } + catch (exception &e) + { + fout_rec << "ERROR: " << e.what() << endl; + throw runtime_error(e.what()); + + } + catch (...) + { + throw runtime_error("ERROR in sLP process"); + } + fout_rec << "Number of forward model runs performed during optimization: " << run_manager_ptr->get_total_runs() << endl; } // clean up //fout_rec.close(); diff --git a/src/programs/pestpp-sqp/pestpp-sqp.cpp b/src/programs/pestpp-sqp/pestpp-sqp.cpp index 6deafe34..d2182411 100644 --- a/src/programs/pestpp-sqp/pestpp-sqp.cpp +++ b/src/programs/pestpp-sqp/pestpp-sqp.cpp @@ -84,8 +84,10 @@ int main(int argc, char* argv[]) yam_agent.process_ctl_file(ctl_file); } - catch (PestError e) + catch (exception &e) { + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; cerr << "Error processing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw(e); @@ -170,11 +172,11 @@ int main(int argc, char* argv[]) file_manager.close_file("pst"); performance_log.log_event("finished processing control file"); } - catch (PestError e) + catch (exception &e) { - cerr << "Error prococessing control file: " << cmdline.ctl_file_name << endl << endl; + cerr << "Error processing control file: " << cmdline.ctl_file_name << endl << endl; cerr << e.what() << endl << endl; - fout_rec << "Error prococessing control file: " << cmdline.ctl_file_name << endl << endl; + fout_rec << "Error processing control file: " << cmdline.ctl_file_name << endl << endl; fout_rec << e.what() << endl << endl; fout_rec.close(); throw(e); diff --git a/src/programs/sweep/sweep.cpp b/src/programs/sweep/sweep.cpp index c96b9333..1e10ebdf 100644 --- a/src/programs/sweep/sweep.cpp +++ b/src/programs/sweep/sweep.cpp @@ -493,8 +493,10 @@ int main(int argc, char* argv[]) yam_agent.process_ctl_file(ctl_file); } - catch (PestError e) + catch (exception &e) { + frec << "Error processing control file: " << ctl_file << endl << endl; + frec << e.what() << endl << endl; cerr << "Error prococessing control file: " << ctl_file << endl << endl; cerr << e.what() << endl << endl; throw(e); @@ -567,7 +569,7 @@ int main(int argc, char* argv[]) file_manager.close_file("pst"); performance_log.log_event("finished processing control file"); } - catch (PestError e) + catch (exception &e) { cerr << "Error processing control file: " << filename << endl << endl; cerr << e.what() << endl << endl;