Skip to content

Commit

Permalink
working on mean shift iters
Browse files Browse the repository at this point in the history
  • Loading branch information
@jtwhite79
jtwhite79 committed Sep 30, 2023
1 parent 48e99f3 commit 6d5f334
Show file tree
Hide file tree
Showing 4 changed files with 124 additions and 49 deletions.
151 changes: 107 additions & 44 deletions src/libs/pestpp_common/EnsembleMethodUtils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3871,26 +3871,27 @@ vector<int> L2PhiHandler::get_idxs_greater_than(double bad_phi, double bad_phi_s
double std = calc_std(&_meas);
vector<int> idxs;
vector<string> names = oe.get_real_names();
double bad_thres = 1.0e+300;
if (bad_phi_sigma > 0) {
bad_thres = mean + (std * bad_phi_sigma);
ofstream &frec = file_manager->rec_ofstream();
frec << "...bad_phi_sigma " << bad_phi_sigma << " and mean " << mean << " yields bad phi threshold of " << bad_thres << endl;
}
else
{
vector<double> meas_vec;
for (auto & m : _meas)
meas_vec.push_back(m.second);
sort(meas_vec.begin(),meas_vec.end());
double qval = -1. * bad_phi_sigma / 100.;
int qidx = (int)(qval * (double)meas_vec.size());
qidx = min(0,qidx);
qidx = max(((int)meas_vec.size()) - 2,qidx);
bad_thres = meas_vec[qidx];
ofstream& frec = file_manager->rec_ofstream();
frec << "...bad_phi_sigma quantile value " << qval << " yields bad phi threshold of " << bad_thres << endl;
double bad_thres = std::numeric_limits<double>::max();
if (bad_phi_sigma < std::numeric_limits<double>::max()) {
if (bad_phi_sigma > 0) {
bad_thres = min(std::numeric_limits<double>::max(), mean + (std * bad_phi_sigma));
ofstream &frec = file_manager->rec_ofstream();
frec << "...bad_phi_sigma " << bad_phi_sigma << " and mean " << mean << " yields bad phi threshold of "
<< bad_thres << endl;
} else {
vector<double> meas_vec;
for (auto &m : _meas)
meas_vec.push_back(m.second);
sort(meas_vec.begin(), meas_vec.end());
double qval = -1. * bad_phi_sigma / 100.;
int qidx = (int) (qval * (double) meas_vec.size());
qidx = min(0, qidx);
qidx = max(((int) meas_vec.size()) - 2, qidx);
bad_thres = min(std::numeric_limits<double>::max(),meas_vec[qidx]);
ofstream &frec = file_manager->rec_ofstream();
frec << "...bad_phi_sigma quantile value " << qval << " yields bad phi threshold of " << bad_thres << endl;

}
}

for (int i = 0; i < names.size(); i++) {
Expand Down Expand Up @@ -5623,7 +5624,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing)
reg_factor = pest_scenario.get_pestpp_options().get_ies_reg_factor();
message(1, "using reg_factor: ", reg_factor);
double bad_phi = pest_scenario.get_pestpp_options().get_ies_bad_phi();
if (bad_phi < 1.0e+30)
if (bad_phi < std::numeric_limits<double>::max())
message(1, "using bad_phi: ", bad_phi);

int num_reals = pest_scenario.get_pestpp_options().get_ies_num_reals();
Expand Down Expand Up @@ -7314,7 +7315,6 @@ bool EnsembleMethod::solve(bool use_mda, vector<double> inflation_factors, vecto
local_subset_size = 4;
message(2,ss.str());
}

}

if ((use_subset) && (local_subset_size > pe.shape().first))
Expand Down Expand Up @@ -7500,7 +7500,7 @@ bool EnsembleMethod::solve(bool use_mda, vector<double> inflation_factors, vecto

vector<map<int, int>> real_run_ids_lams;
int best_idx = -1;
double best_mean = 1.0e+30, best_std = 1.0e+30; // todo (Ayman): read those from input
double best_mean = 1.0e+300, best_std = 1.0e+300; // todo (Ayman): read those from input
double mean, std;

message(0, "running upgrade ensembles");
Expand Down Expand Up @@ -7598,7 +7598,67 @@ bool EnsembleMethod::solve(bool use_mda, vector<double> inflation_factors, vecto


//subset stuff here
if ((best_idx != -1) && (use_subset) && (local_subset_size < pe.shape().first))

if (iter <= pest_scenario.get_pestpp_options().get_ies_n_iter_mean()) {
performance_log->log_event("this is a mean-only iteration");
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());

performance_log->log_event("getting current parameter ensemble");
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);
}

vector<double> mean_vec = pe_best.get_mean_stl_var_vector();
//TODO: save par file for mean vector
pe_best = pe.zeros_like(0);
if (verbose_level > 2)
{
performance_log->log_event("saving 'best' parameter ensemble that was used to form mean vector);
//TODO:...
//pe_best.to_csv()
}
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);
}
performance_log->log_event("running mean-shifted prior realizations");
run_ensemble_util(performance_log,frec,new_pe,oe_lam_best,run_mgr_ptr);
pe_lams[best_idx] = new_pe;
//make sure we dont try to process the subset stuff below
local_subset_size = pe.shape().first;
}
else if ((best_idx != -1) && (use_subset) && (local_subset_size < pe.shape().first))
{
double acc_phi = last_best_mean * acc_fac;
Expand Down Expand Up @@ -7672,7 +7732,6 @@ bool EnsembleMethod::solve(bool use_mda, vector<double> inflation_factors, vecto
//need to work out which par and obs en real names to run - some may have failed during subset testing...
ObservationEnsemble remaining_oe_lam = oe;//copy

ParameterEnsemble remaining_pe_lam = pe_lams[best_idx];
remaining_pe_lam.set_fixed_info(pe.get_fixed_info());
Expand Down Expand Up @@ -7700,23 +7759,7 @@ bool EnsembleMethod::solve(bool use_mda, vector<double> inflation_factors, vecto
//remove any failed runs from subset testing
if (missing.size() > 0)
remaining_pe_lam.drop_rows(missing);
for (auto& pe_filename : pe_filenames)
{
performance_log->log_event("removing upgrade ensemble '" + pe_filename + "'");
try
{
remove(pe_filename.c_str());
}
catch (exception& e)
{
message(2, "error removing upgrade ensemble: '" + pe_filename + "' :" + e.what());
}
catch (...)
{
message(2, "error removing upgrade ensemble: '" + pe_filename);
}

}
remove_external_pe_filenames(pe_filenames);
}
Expand Down Expand Up @@ -7907,13 +7950,33 @@ bool EnsembleMethod::solve(bool use_mda, vector<double> inflation_factors, vecto
last_best_lam = new_lam;
}
save_ensembles("rejected",cycle,pe_lams[best_idx],oe_lam_best);


}
//report_and_save();
report_and_save();
return true;
}
void EnsembleMethod::remove_external_pe_filenames(vector<string>& pe_filenames)
{
for (auto& pe_filename : pe_filenames)
{
performance_log->log_event("removing upgrade ensemble '" + pe_filename + "'");
try
{
remove(pe_filename.c_str());
}
catch (exception& e)
{
message(2, "error removing upgrade ensemble: '" + pe_filename + "' :" + e.what());
}
catch (...)
{
message(2, "error removing upgrade ensemble: '" + pe_filename);
}
}
}
//template<typename T, typename A>
//void EnsembleMethod::message(int level, const string& _message, vector<T, A> _extras, bool echo)
//{
Expand Down
2 changes: 2 additions & 0 deletions src/libs/pestpp_common/EnsembleMethodUtils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -499,5 +499,7 @@ class EnsembleMethod

void prep_drop_violations();

void remove_external_pe_filenames(vector<string>& pe_filenames);

};
#endif
14 changes: 10 additions & 4 deletions src/libs/pestpp_common/pest_data_structs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1177,8 +1177,12 @@ bool PestppOptions::assign_ies_value_by_key(const string& key, const string& val
ies_phi_factors_by_real = pest_utils::parse_string_arg_to_bool(value);
return true;
}

return false;
else if (key == "IES_N_ITER_MEAN")
{
convert_ip(value,ies_n_iter_mean);
return true;
}
return false;
}

bool PestppOptions::assign_da_value_by_key(const string& key, const string& value, const string& org_value)
Expand Down Expand Up @@ -1805,6 +1809,7 @@ void PestppOptions::summary(ostream& os) const
os << "ies_localizer_forgive_extra: " << ies_localizer_forgive_missing << endl;
os << "ies_phi_factors_file: " << ies_phi_fractions_file << endl;
os << "ies_phi_factors_by_real: " << ies_phi_factors_by_real << endl;
os << "ies_n_iter_mean: " << ies_n_iter_mean << endl;


os << endl << "pestpp-sen options: " << endl;
Expand Down Expand Up @@ -1945,8 +1950,8 @@ void PestppOptions::set_defaults()
set_ies_verbose_level(1);
set_ies_use_prior_scaling(false);
set_ies_num_reals(50);
set_ies_bad_phi(1.0e+300);
set_ies_bad_phi_sigma(1.0e+300);
set_ies_bad_phi(std::numeric_limits<double>::max());
set_ies_bad_phi_sigma(std::numeric_limits<double>::max());
set_ies_include_base(true);
set_ies_use_empirical_prior(false);
set_ies_group_draws(true);
Expand Down Expand Up @@ -1988,6 +1993,7 @@ void PestppOptions::set_defaults()
set_ies_localizer_forgive_missing(false);
set_ies_phi_fractions_files("");
set_ies_phi_factors_by_real(false);
set_ies_n_iter_mean(0);

// DA parameters
//set_da_use_ies(false);
Expand Down
6 changes: 5 additions & 1 deletion src/libs/pestpp_common/pest_data_structs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -540,8 +540,11 @@ class PestppOptions {
void set_ies_multimodal_alpha(double _flag) { ies_multimodal_alpha = _flag; }
void set_ensemble_output_precision(int prec) { ensemble_output_precision = prec;}
int get_ensemble_output_precision() const {return ensemble_output_precision;}
void set_ies_n_iter_mean(int _n_iter_mean) {ies_n_iter_mean = _n_iter_mean;}
int get_ies_n_iter_mean() const {return ies_n_iter_mean;}

string get_gsa_method() const { return gsa_method; }

string get_gsa_method() const { return gsa_method; }
void set_gsa_method(string _m) { gsa_method = _m; }
bool get_gsa_morris_pooled_obs() const { return gsa_morris_pooled_obs; }
void set_gsa_morris_pooled_obs(bool _flag) {gsa_morris_pooled_obs = _flag; }
Expand Down Expand Up @@ -812,6 +815,7 @@ class PestppOptions {
bool ies_localizer_forgive_missing;
string ies_phi_fractions_file;
bool ies_phi_factors_by_real;
int ies_n_iter_mean;


// Data Assimilation parameters
Expand Down

0 comments on commit 6d5f334

Please sign in to comment.