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/constraints.cpp b/src/libs/pestpp_common/constraints.cpp index f90bc1f0..c860df55 100644 --- a/src/libs/pestpp_common/constraints.cpp +++ b/src/libs/pestpp_common/constraints.cpp @@ -241,7 +241,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 +786,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 +1069,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 +1192,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 +1310,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 +1336,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 +1363,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 +1418,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 +1509,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 +1521,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 +1531,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 +1554,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 +1626,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 +1634,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 +2200,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/sequential_lp.cpp b/src/libs/pestpp_common/sequential_lp.cpp index 309e3a24..fa9a3db8 100644 --- a/src/libs/pestpp_common/sequential_lp.cpp +++ b/src/libs/pestpp_common/sequential_lp.cpp @@ -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/programs/pestpp-opt/pestpp-opt.cpp b/src/programs/pestpp-opt/pestpp-opt.cpp index ddcdd036..8d7ed52a 100644 --- a/src/programs/pestpp-opt/pestpp-opt.cpp +++ b/src/programs/pestpp-opt/pestpp-opt.cpp @@ -365,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();