Skip to content

Commit

Permalink
weaving in option from pestpp-opt to not use stack in anomaly form if…
Browse files Browse the repository at this point in the history 
… possible - depends on if stack was evaluated for the current iteration
  • Loading branch information
@jtwhite79
jtwhite79 committed Jul 26, 2024
1 parent b07a526 commit 830a8ac
Show file tree
Hide file tree
Showing 5 changed files with 86 additions and 39 deletions.
6 changes: 3 additions & 3 deletions src/libs/common/utilities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
}
}
Expand Down
59 changes: 39 additions & 20 deletions src/libs/pestpp_common/constraints.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,7 +241,7 @@ void Constraints::initialize(vector<string>& 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);
Expand Down Expand Up @@ -786,7 +786,10 @@ void Constraints::initialize(vector<string>& 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"))
{
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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();
Expand All @@ -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) {
Expand All @@ -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();
Expand Down Expand Up @@ -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<string> names = shifted_obs.get_keys();
Expand Down Expand Up @@ -1506,7 +1509,8 @@ ObservationEnsemble Constraints::get_stack_mean(map<string, ObservationEnsemble>
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;

Expand All @@ -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);
Expand All @@ -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<string, int> 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<string, double>, map<string, double>> mm = _stack_oe.get_moment_maps();
map<string, double> 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);
Expand All @@ -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];
Expand Down Expand Up @@ -1607,15 +1626,15 @@ vector<double> Constraints::get_constraint_residual_vec(Observations& sim)
return residuals_vec;
}

pair<vector<double>,vector<double>> Constraints::get_constraint_bound_vectors(Parameters& current_pars, Observations& current_obs)
pair<vector<double>,vector<double>> 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.
These are needed for the simplex solve*/
vector<double> 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
Expand Down Expand Up @@ -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)
Expand Down
11 changes: 5 additions & 6 deletions src/libs/pestpp_common/constraints.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>,vector<double>> get_constraint_bound_vectors(Parameters& current_pars, Observations &current_obs);
pair<vector<double>,vector<double>> get_constraint_bound_vectors(Parameters& current_pars, Observations &current_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<string> get_fosm_par_names();

Expand Down Expand Up @@ -151,10 +150,10 @@ class Constraints
static pair<ConstraintSense, string> 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());

Expand Down Expand Up @@ -194,7 +193,7 @@ class Constraints
map<string, ObservationEnsemble> stack_oe_map;
map<string, Parameters> 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;
Expand Down
25 changes: 21 additions & 4 deletions src/libs/pestpp_common/sequential_lp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>, vector<double>> bounds = constraints.get_constraint_bound_vectors(current_pars, current_constraints_sim);
}
pair<vector<double>, vector<double>> 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);
Expand All @@ -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);
Expand Down Expand Up @@ -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
Expand Down
24 changes: 18 additions & 6 deletions src/programs/pestpp-opt/pestpp-opt.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down

0 comments on commit 830a8ac

Please sign in to comment.