diff --git a/README.md b/README.md index ecc17f71..d824192c 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,8 @@ pestpplogo image

+ + # PEST++ ## a Software Suite for Parameter Estimation, Uncertainty Analysis, Management Optimization and Sensitivity Analysis diff --git a/documentation/pestpp_users_guide_v5.1.18.docx b/documentation/pestpp_users_guide_v5.1.20.docx similarity index 72% rename from documentation/pestpp_users_guide_v5.1.18.docx rename to documentation/pestpp_users_guide_v5.1.20.docx index 9409c38f..91ca95ca 100644 Binary files a/documentation/pestpp_users_guide_v5.1.18.docx and b/documentation/pestpp_users_guide_v5.1.20.docx differ diff --git a/documentation/pestpp_users_manual.md b/documentation/pestpp_users_manual.md index 92e092b9..2821063f 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.1.18 +# Version 5.1.20 PEST++ Development Team -July 2022 +August 2022 # Acknowledgements @@ -70,7 +70,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI # Table of Contents -- [Version 5.1.18](#s1) +- [Version 5.1.20](#s1) - [Acknowledgements](#s2) - [Preface](#s3) - [License](#s4) diff --git a/src/libs/common/config_os.h b/src/libs/common/config_os.h index 9e318fb2..a5f742d7 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.1.18"; +#define PESTPP_VERSION "5.1.20"; #if defined(_WIN32) || defined(_WIN64) #define OS_WIN diff --git a/src/libs/pestpp_common/Ensemble.cpp b/src/libs/pestpp_common/Ensemble.cpp index 64bc68ef..7e971b7c 100644 --- a/src/libs/pestpp_common/Ensemble.cpp +++ b/src/libs/pestpp_common/Ensemble.cpp @@ -2173,6 +2173,7 @@ map ParameterEnsemble::add_runs(RunManagerAbstract *run_mgr_ptr,const v { par_transform.active_ctl2model_ip(pars); } + //map> tied_items = par_transform.get_tied_ptr()->get_items(); Parameters pars_real = pars; Eigen::VectorXd evec; vector svec; @@ -2937,6 +2938,10 @@ void ParameterEnsemble::to_binary_ordered(string file_name) throw runtime_error("error opening file for binary parameter ensemble:" + file_name); } + bool has_tied = false; + if (par_transform.get_tied_ptr()->get_items().size() > 0) + has_tied = true; + //vector vnames = var_names; vector vnames = pest_scenario_ptr->get_ctl_ordered_par_names(); //vnames.insert(vnames.end(), fixed_names.begin(), fixed_names.end()); @@ -3010,6 +3015,11 @@ void ParameterEnsemble::to_binary_ordered(string file_name) par_transform.model2ctl_ip(pars); else if (tstat == transStatus::NUM) par_transform.numeric2ctl_ip(pars); + else if ((has_tied) && (tstat == transStatus::CTL)) + { + par_transform.active_ctl2model_ip(pars); + par_transform.model2ctl_ip(pars); + } replace_fixed(real_names[irow], pars); for (int jcol = 0; jcolget_items().size() > 0) + has_tied = true; //map::const_iterator found_pi_par; //map::const_iterator not_found_pi_par; //icount = row_idxs + 1 + col_idxs * self.shape[0] @@ -3145,6 +3158,11 @@ void ParameterEnsemble::to_dense_ordered(string file_name) par_transform.model2ctl_ip(pars); else if (tstat == transStatus::NUM) par_transform.numeric2ctl_ip(pars); + else if ((has_tied) && (tstat == transStatus::CTL)) + { + par_transform.active_ctl2model_ip(pars); + par_transform.model2ctl_ip(pars); + } replace_fixed(real_names[irow], pars); string name = real_names[irow]; tmp = name.size(); @@ -3331,7 +3349,9 @@ void ParameterEnsemble::to_csv_by_vars(ofstream &csv, bool write_header) } if (write_header) csv << endl; - + bool has_tied = false; + if (par_transform.get_tied_ptr()->get_items().size() > 0) + has_tied = true; //get the pars and transform to be in sync with ensemble trans status /*Parameters pars = pest_scenario_ptr->get_ctl_parameters(); if (tstat == transStatus::NUM) @@ -3358,6 +3378,11 @@ void ParameterEnsemble::to_csv_by_vars(ofstream &csv, bool write_header) { par_transform.active_ctl2model_ip(pars); } + else if ((has_tied) && (tstat == transStatus::CTL)) + { + par_transform.active_ctl2model_ip(pars); + par_transform.model2ctl_ip(pars); + } ireal = real_map[rname]; pars.update_without_clear(var_names, reals.row(ireal)); if (tstat == transStatus::MODEL) @@ -3404,6 +3429,10 @@ void ParameterEnsemble::to_csv_by_reals(ofstream &csv, bool write_header) }*/ //for (int ireal = 0; ireal < reals.rows(); ireal++) + bool has_tied = false; + if (par_transform.get_tied_ptr()->get_items().size() > 0) + has_tied = true; + int ireal = 0; map real_map; for (int i = 0; i < real_names.size(); i++) @@ -3435,7 +3464,14 @@ void ParameterEnsemble::to_csv_by_reals(ofstream &csv, bool write_header) par_transform.model2ctl_ip(pars); else if (tstat == transStatus::NUM) par_transform.numeric2ctl_ip(pars); + // here we need to check for tied parameters... + else if ((has_tied) && (tstat == transStatus::CTL)) + { + par_transform.active_ctl2model_ip(pars); + par_transform.model2ctl_ip(pars); + } replace_fixed(real_names[ireal],pars); + for (auto &name : names) csv << ',' << pars[name]; csv << endl; diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index 6019afeb..3c2291d0 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -3883,7 +3883,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) if ((obs_restart_csv.size() == 0) && (!pe_drawn) && (oe_drawn)) { vector rnames = pe.get_real_names(); - oe_base.set_real_names(rnames); + oe_base.set_real_names(rnames,true); message(2, "resetting obs + noise ensemble real names to parameter ensemble real names"); } diff --git a/src/libs/pestpp_common/MOEA.cpp b/src/libs/pestpp_common/MOEA.cpp index 8f383c4b..b8a1ec27 100644 --- a/src/libs/pestpp_common/MOEA.cpp +++ b/src/libs/pestpp_common/MOEA.cpp @@ -1358,7 +1358,8 @@ void MOEA::update_archive_nsga(ObservationEnsemble& _op, ParameterEnsemble& _dp) op_archive.keep_rows(keep); dp_archive.keep_rows(keep); } - + dp_archive.reset_org_real_names(); + op_archive.reset_org_real_names(); save_populations(dp_archive, op_archive, "archive"); } @@ -1417,7 +1418,8 @@ void MOEA::update_archive_spea(ObservationEnsemble& _op, ParameterEnsemble& _dp) op_archive.keep_rows(keep); dp_archive.keep_rows(keep); } - + dp_archive.reset_org_real_names(); + op_archive.reset_org_real_names(); save_populations(dp_archive, op_archive, "archive"); } } @@ -3776,7 +3778,7 @@ ParameterEnsemble MOEA::generate_sbx_population(int num_members, ParameterEnsemb return tmp_dp; } -void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, string tag) +void MOEA::save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag) { stringstream ss; @@ -3791,16 +3793,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri if (pest_scenario.get_pestpp_options().get_save_binary()) { ss << ".jcb"; - dp.to_binary(ss.str()); + _dp.to_binary(ss.str()); } else { ss << ".csv"; - dp.to_csv(ss.str()); + _dp.to_csv(ss.str()); } string name = ss.str(); ss.str(""); - ss << " saved decision variable population of size " << dp.shape().first << " X " << dp.shape().second << " to '" << name << "'"; + ss << " saved decision variable population of size " << _dp.shape().first << " X " << _dp.shape().second << " to '" << name << "'"; message(1, ss.str()); ss.str(""); if ((save_every > 0) && (iter % save_every == 0)) @@ -3814,16 +3816,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri if (pest_scenario.get_pestpp_options().get_save_binary()) { ss << ".jcb"; - dp.to_binary(ss.str()); + _dp.to_binary(ss.str()); } else { ss << ".csv"; - dp.to_csv(ss.str()); + _dp.to_csv(ss.str()); } string name = ss.str(); ss.str(""); - ss << " saved generation-specific decision variable population of size " << dp.shape().first << " X " << dp.shape().second << " to '" << name << "'"; + ss << " saved generation-specific decision variable population of size " << _dp.shape().first << " X " << _dp.shape().second << " to '" << name << "'"; message(1, ss.str()); } @@ -3838,16 +3840,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri if (pest_scenario.get_pestpp_options().get_save_binary()) { ss << ".jcb"; - op.to_binary(ss.str()); + _op.to_binary(ss.str()); } else { ss << ".csv"; - op.to_csv(ss.str()); + _op.to_csv(ss.str()); } name = ss.str(); ss.str(""); - ss << " saved observation population of size " << op.shape().first << " X " << op.shape().second << " to '" << name << "'"; + ss << " saved observation population of size " << _op.shape().first << " X " << _op.shape().second << " to '" << name << "'"; message(1, ss.str()); if ((save_every > 0) && (iter % save_every == 0)) @@ -3862,16 +3864,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri if (pest_scenario.get_pestpp_options().get_save_binary()) { ss << ".jcb"; - op.to_binary(ss.str()); + _op.to_binary(ss.str()); } else { ss << ".csv"; - op.to_csv(ss.str()); + _op.to_csv(ss.str()); } name = ss.str(); ss.str(""); - ss << " saved generation-specific observation population of size " << op.shape().first << " X " << op.shape().second << " to '" << name << "'"; + ss << " saved generation-specific observation population of size " << _op.shape().first << " X " << _op.shape().second << " to '" << name << "'"; message(1, ss.str()); } diff --git a/src/libs/pestpp_common/MOEA.h b/src/libs/pestpp_common/MOEA.h index 73da48ac..4355d0bd 100644 --- a/src/libs/pestpp_common/MOEA.h +++ b/src/libs/pestpp_common/MOEA.h @@ -84,8 +84,6 @@ class ParetoObjectives void fill_domination_containers(map>& _member_struct, map>&solutions_dominated_map, map& num_dominating_map, bool dup_as_dom=false); - - bool compare_two_nsga(string& first, string& second); bool compare_two_spea(string& first, string& second); @@ -222,7 +220,7 @@ class MOEA string get_new_member_name(string tag = string()); - void save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, string tag = string()); + void save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag = string()); void gauss_mutation_ip(ParameterEnsemble& _dp); pair sbx(double probability, double eta_m, int idx1, int idx2); pair sbx_new(double crossover_probability, double di, Eigen::VectorXd& parent1, diff --git a/src/libs/pestpp_common/constraints.cpp b/src/libs/pestpp_common/constraints.cpp index 8312167a..9651dea3 100644 --- a/src/libs/pestpp_common/constraints.cpp +++ b/src/libs/pestpp_common/constraints.cpp @@ -1649,7 +1649,8 @@ double Constraints::get_probit() double Constraints::get_probit(double _risk) { - /* the probit function estimate - needed for the fosm-basd chance constraints*/ + /* the probit function estimate + * - needed for the fosm-basd chance constraints*/ double output = sqrt(2.0) * ErfInv2((2.0 * _risk) - 1.0); return output; } @@ -2221,9 +2222,13 @@ void Constraints::postsolve_obs_constraints_report(Observations& old_obs, Observ This really only applies to the simplex solution because we can, with the linear assumption, project what the resulting constraint values should be*/ + int width = 20; + for(auto& n : ctl_ord_obs_constraint_names) + width = max(width,(int)n.size()); + width = width + 2; ofstream& f_rec = file_mgr_ptr->rec_ofstream(); f_rec << endl << endl << " " << tag << " constraint/objective information at end of iteration " << iter << endl << endl; - f_rec << setw(20) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25); + f_rec << setw(width) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25); if (status_map.size() > 0) f_rec << "simplex status"; if (price_map.size() > 0) @@ -2247,7 +2252,7 @@ void Constraints::postsolve_obs_constraints_report(Observations& old_obs, Observ { string name = ctl_ord_obs_constraint_names[i]; sim_val = new_obs[name]; - f_rec << setw(20) << left << name; + f_rec << setw(width) << left << name; f_rec << setw(15) << right << constraint_sense_name[name]; f_rec << setw(15) << constraints_obs.get_rec(name); if (status_map.size() > 0) @@ -2287,8 +2292,12 @@ void Constraints::postsolve_pi_constraints_report(Parameters& old_pars, Paramete ofstream &f_rec = file_mgr_ptr->rec_ofstream(); if (num_pi_constraints() > 0) { + int width = 20; + for (auto& n : ctl_ord_pi_constraint_names) + width = max(width,(int)n.size()); + width = width + 2; f_rec << endl << endl << " --- prior information constraint summary at end of iteration " << iter << " --- " << endl << endl; - f_rec << setw(20) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25); + f_rec << setw(width) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25); if (status_map.size() > 0) f_rec << "simplex status"; if (price_map.size() > 0) @@ -2301,7 +2310,7 @@ void Constraints::postsolve_pi_constraints_report(Parameters& old_pars, Paramete PriorInformationRec pi_rec = constraints_pi.get_pi_rec(name); pair cur_sim_resid = pi_rec.calc_sim_and_resid(old_pars); pair new_sim_resid = pi_rec.calc_sim_and_resid(new_pars); - f_rec << setw(20) << left << name; + f_rec << setw(width) << left << name; f_rec << setw(15) << right << constraint_sense_name[name]; f_rec << setw(15) << pi_rec.get_obs_value(); if (status_map.size() > 0) diff --git a/src/libs/pestpp_common/sequential_lp.cpp b/src/libs/pestpp_common/sequential_lp.cpp index 98bae0d1..a2f2e788 100644 --- a/src/libs/pestpp_common/sequential_lp.cpp +++ b/src/libs/pestpp_common/sequential_lp.cpp @@ -185,10 +185,14 @@ map sequentialLP::get_out_of_bounds_dec_vars(Parameters &upgrade_ pair sequentialLP::postsolve_decision_var_report(Parameters &upgrade_pars) { + int width = 12; + for(auto& n : dv_names) + width = max(width,(int)n.size()); + width = width + 2; ofstream &f_rec = file_mgr_ptr->rec_ofstream(); const double *reduced_cost = model.getReducedCost(); f_rec << endl << endl << " decision variable information at end of SLP iteration " << slp_iter << endl << endl; - f_rec << setw(20) << left << "name" << right << setw(15) << "current" << setw(15) << "new"; + f_rec << setw(width) << left << "name" << right << setw(15) << "current" << setw(15) << "new"; f_rec << setw(15) << "objfunc coef" << setw(15) << "cur contrib" << setw(15) << "new contrib" << setw(15) << "reduced cost"; f_rec << setw(25) << "simplex status" << endl; string name; @@ -206,7 +210,7 @@ pair sequentialLP::postsolve_decision_var_report(Parameters &upgr new_val = upgrade_pars[name]; actual_pars.update_rec(name, new_val); - f_rec << setw(20) << left << name; + f_rec << setw(width) << left << name; f_rec << setw(15) << right << cur_val; f_rec << setw(15) << new_val; f_rec << setw(15) << obj_coef; @@ -667,8 +671,8 @@ void sequentialLP::solve() cout << " ---------------------------------" << endl << endl << endl; iter_presolve(); - iter_solve(); - iter_postsolve(); + iter_solve(); + iter_postsolve(); if (terminate) break; slp_iter++; if (slp_iter > pest_scenario.get_control_info().noptmax) @@ -841,6 +845,7 @@ void sequentialLP::iter_postsolve() { iter_obj_values.push_back(cur_new_obj.first); obj_best = cur_new_obj.second; + best_pars.update_without_clear(dv_names, upgrade_pars.get_data_vec(dv_names)); } iter_obj_values.push_back(cur_new_obj.second); double obj_func_change = abs(cur_new_obj.first - cur_new_obj.second) / abs(max(max(cur_new_obj.first,cur_new_obj.second),1.0)); @@ -1107,10 +1112,11 @@ void sequentialLP::iter_presolve() bool init_obs = false; if (slp_iter == 1) init_obs = true; - + bool success = jco.build_runs(current_pars, current_constraints_sim, names_to_run, par_trans, pest_scenario.get_base_group_info(), pest_scenario.get_ctl_parameter_info(), *run_mgr_ptr, out_of_bounds,false,init_obs); + if (!success) { const set failed = jco.get_failed_parameter_names(); @@ -1128,6 +1134,7 @@ void sequentialLP::iter_presolve() set failed = run_mgr_ptr->get_failed_run_ids(); //process the remaining responses + success = jco.process_runs(par_trans, pest_scenario.get_base_group_info(), *run_mgr_ptr, *null_prior, false,false); if (!success) @@ -1149,8 +1156,13 @@ void sequentialLP::iter_presolve() if (init_obs) { - //Observations temp_obs; + //Observations temp_obs + //Parameters temp = current_pars; run_mgr_ptr->get_run(0, current_pars, current_constraints_sim, false); +// for (auto& n : dv_names) { +// current_pars.update_rec(n, temp.get_rec(n)); +// } + par_trans.model2ctl_ip(current_pars); } }