diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index c1379895..2654f3c2 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -2492,16 +2492,16 @@ double L2PhiHandler::calc_std(map *phi_map) double L2PhiHandler::get_representative_phi(phiType pt) { - //if (pest_scenario->get_pestpp_options().get_ies_n_iter_mean() < 0) + //if (pest_scenario->get_pestpp_options().get_ies_n_iter_reinflate() < 0) bool use_min = false; - for (auto& fac : pest_scenario->get_pestpp_options().get_ies_n_iter_mean()) - { - if (fac < 0) - { - use_min = true; - break; - } - } +// for (auto& fac : pest_scenario->get_pestpp_options().get_ies_n_iter_reinflate()) +// { +// if (fac < 0) +// { +// use_min = true; +// break; +// } +// } if (use_min) { return get_min(pt); @@ -4317,7 +4317,7 @@ bool EnsembleMethod::should_terminate(int current_n_iter_mean) message(1, "phiredstp: ", phiredstp); message(1, "nphistp: ", nphistp); message(1, "nphinored (also used for consecutive bad lambda cycles): ", nphinored); - //int n_mean_iter = pest_scenario.get_pestpp_options().get_ies_n_iter_mean(); + //int n_mean_iter = pest_scenario.get_pestpp_options().get_ies_n_iter_reinflate(); vector::iterator begin_idx = best_mean_phis.begin(); //if ((current_n_iter_mean > 0) && (best_mean_phis.size() > current_n_iter_mean)) // begin_idx = best_mean_phis.end() - (current_n_iter_mean+1); //bc of prior phi and then adding the mean shift to the list @@ -4910,7 +4910,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) consec_bad_lambda_cycles = 0; reinflate_to_minphi_real = false; bool use_min = false; - for (auto& fac : pest_scenario.get_pestpp_options().get_ies_n_iter_mean()) + for (auto& fac : pest_scenario.get_pestpp_options().get_ies_n_iter_reinflate()) { if (fac < 0) { @@ -4923,7 +4923,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing) { message(2,"n_iter_reinflate < 0, using min-phi real for re-inflation, resetting n_iter_reinflate to positive"); reinflate_to_minphi_real = true; - //pest_scenario.get_pestpp_options_ptr()->set_ies_n_iter_mean(-1 * pest_scenario.get_pestpp_options().get_ies_n_iter_mean()); + //pest_scenario.get_pestpp_options_ptr()->set_ies_n_iter_reinflate(-1 * pest_scenario.get_pestpp_options().get_ies_n_iter_reinflate()); } lam_mults = pest_scenario.get_pestpp_options().get_ies_lam_mults(); if (lam_mults.size() == 0) @@ -7438,35 +7438,35 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto return true; } -void EnsembleMethod::reset_par_ensemble_to_prior_mean(){ +void EnsembleMethod::reset_par_ensemble_to_prior_mean(double reinflate_factor){ + string min_phi_name = ""; - if (true) - { - map pmap = ph.get_phi_map(L2PhiHandler::phiType::ACTUAL); - double min_phi = numeric_limits::max(); + //find the min phi real... + map pmap = ph.get_phi_map(L2PhiHandler::phiType::ACTUAL); + double min_phi = numeric_limits::max(); - for (auto& item : pmap) - { - if (item.second < min_phi) - { - min_phi = item.second; - min_phi_name = item.first; - } - } - vector real_names = oe.get_real_names(); - int idx = distance(real_names.begin(),find(real_names.begin(),real_names.end(),min_phi_name)); - if (idx == real_names.size()) + for (auto& item : pmap) + { + if (item.second < min_phi) { - throw_em_error("min phi realization name not found in obs ensemble"); + min_phi = item.second; + min_phi_name = item.first; } - real_names = pe.get_real_names(); - min_phi_name = real_names[idx]; - } + vector real_names = oe.get_real_names(); + int idx = distance(real_names.begin(),find(real_names.begin(),real_names.end(),min_phi_name)); + if (idx == real_names.size()) + { + throw_em_error("min phi realization name not found in obs ensemble"); + } + real_names = pe.get_real_names(); + min_phi_name = real_names[idx]; + message(0,"resetting current parameter ensemble to prior ensemble with current ensemble mean"); + message(1,"reinflation factor:",reinflate_factor); performance_log->log_event("getting prior parameter ensemble mean-centered anomalies"); Eigen::MatrixXd anoms = pe_base.get_eigen_anomalies(pe_base.get_real_names(), pe.get_var_names(), pest_scenario.get_pestpp_options().get_ies_center_on()); - + anoms = anoms * reinflate_factor; performance_log->log_event("getting current parameter ensemble mean vector"); vector mean_vec = pe.get_mean_stl_var_vector(); Eigen::VectorXd offset(mean_vec.size()); @@ -7478,7 +7478,6 @@ void EnsembleMethod::reset_par_ensemble_to_prior_mean(){ message(2,"using min-phi realization for offset"); } - performance_log->log_event("adding offset to anomalies"); for (int i = 0; i < mean_vec.size(); i++) { diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.h b/src/libs/pestpp_common/EnsembleMethodUtils.h index fcb1575f..0b27bb13 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.h +++ b/src/libs/pestpp_common/EnsembleMethodUtils.h @@ -499,7 +499,7 @@ class EnsembleMethod double get_lambda(); - void reset_par_ensemble_to_prior_mean(); + void reset_par_ensemble_to_prior_mean(double reinflate_factor); }; #endif diff --git a/src/libs/pestpp_common/EnsembleSmoother.cpp b/src/libs/pestpp_common/EnsembleSmoother.cpp index 96038c14..5fc06523 100644 --- a/src/libs/pestpp_common/EnsembleSmoother.cpp +++ b/src/libs/pestpp_common/EnsembleSmoother.cpp @@ -36,11 +36,13 @@ void IterEnsembleSmoother::iterate_2_solution() ofstream &frec = file_manager.rec_ofstream(); bool accept; - vector n_iter_mean = pest_scenario.get_pestpp_options().get_ies_n_iter_mean(); + vector n_iter_reinflate = pest_scenario.get_pestpp_options().get_ies_n_iter_reinflate(); + vector reinflate_factor = pest_scenario.get_pestpp_options().get_ies_reinflate_factor(); int iters_since_reinflate = 0; - int n_iter_mean_idx = 0; - int current_n_iter_mean = n_iter_mean[n_iter_mean_idx]; + int n_iter_reinflate_idx = 0; + int current_n_iter_mean = n_iter_reinflate[n_iter_reinflate_idx]; + double current_reinflate_factor = reinflate_factor[n_iter_reinflate_idx]; int solution_iter = 0; for (int i = 0; i < pest_scenario.get_control_info().noptmax; i++) { @@ -76,23 +78,29 @@ void IterEnsembleSmoother::iterate_2_solution() else consec_bad_lambda_cycles++; - //if ((n_iter_mean > 0) && (solution_iter % n_iter_mean == 0)) + //if ((n_iter_reinflate > 0) && (solution_iter % n_iter_reinflate == 0)) iters_since_reinflate++; if ((current_n_iter_mean != 0) && (iters_since_reinflate >= current_n_iter_mean)) { + message(2,"incrementing iteration count for reinflation cycle"); iter++; - reset_par_ensemble_to_prior_mean(); + + reset_par_ensemble_to_prior_mean(current_reinflate_factor); iters_since_reinflate = 0; - n_iter_mean_idx++; - if (n_iter_mean.size() > n_iter_mean_idx) + n_iter_reinflate_idx++; + if (reinflate_factor.size() > n_iter_reinflate_idx) + { + current_reinflate_factor = reinflate_factor[n_iter_reinflate_idx]; + } + if (n_iter_reinflate.size() > n_iter_reinflate_idx) { - current_n_iter_mean = n_iter_mean[n_iter_mean_idx]; + current_n_iter_mean = n_iter_reinflate[n_iter_reinflate_idx]; } } if (should_terminate(current_n_iter_mean)) { - //if (iter > pest_scenario.get_pestpp_options().get_ies_n_iter_mean()) { + //if (iter > pest_scenario.get_pestpp_options().get_ies_n_iter_reinflate()) { if (current_n_iter_mean == 0) { break; } diff --git a/src/libs/pestpp_common/pest_data_structs.cpp b/src/libs/pestpp_common/pest_data_structs.cpp index 56878d5e..4c6d4283 100644 --- a/src/libs/pestpp_common/pest_data_structs.cpp +++ b/src/libs/pestpp_common/pest_data_structs.cpp @@ -1185,12 +1185,23 @@ bool PestppOptions::assign_ies_value_by_key(const string& key, const string& val { passed_args.insert("IES_N_ITER_MEAN"); passed_args.insert("IES_N_ITER_REINFLATE"); - ies_n_iter_mean.clear(); + ies_n_iter_reinflate.clear(); vector tok; tokenize(value, tok, ","); for (const auto& fac : tok) { - ies_n_iter_mean.push_back(convert_cp(fac)); + ies_n_iter_reinflate.push_back(convert_cp(fac)); + } + return true; + } + else if (key == "IES_REINFLATE_FACTOR") + { + vector tok; + tokenize(value, tok, ","); + ies_reinflate_factor.clear(); + for (const auto& fac : tok) + { + ies_reinflate_factor.push_back(convert_cp(fac)); } return true; } @@ -1844,7 +1855,11 @@ void PestppOptions::summary(ostream& os) const 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_reinflate: " << endl; - for (auto v : ies_n_iter_mean) + for (auto v : ies_n_iter_reinflate) + os << v << ","; + os << endl; + os << "ies_reinflate_factor: " << endl; + for (auto v : ies_reinflate_factor) os << v << ","; os << endl; os << "ies_updatebyreals: " << ies_updatebyreals << endl; @@ -2033,11 +2048,12 @@ 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(vector{0}); + set_ies_n_iter_reinflate(vector < int > {0}); + set_ies_reinflate_factor(vector < double > {1.0}); + set_ies_updatebyreals(false); set_save_dense(false); - // 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 858086a9..5756970c 100644 --- a/src/libs/pestpp_common/pest_data_structs.h +++ b/src/libs/pestpp_common/pest_data_structs.h @@ -550,8 +550,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(vector _n_iter_mean) {ies_n_iter_mean = _n_iter_mean;} - vector get_ies_n_iter_mean() const {return ies_n_iter_mean;} + void set_ies_n_iter_reinflate(vector _n_iter_reinflate) { ies_n_iter_reinflate = _n_iter_reinflate;} + vector get_ies_n_iter_reinflate() const {return ies_n_iter_reinflate;} + void set_ies_reinflate_factor(vector reinflate_factor) { ies_reinflate_factor = reinflate_factor;} + vector get_ies_reinflate_factor() const {return ies_reinflate_factor;} + string get_gsa_method() const { return gsa_method; } @@ -828,7 +831,8 @@ class PestppOptions { bool ies_localizer_forgive_missing; string ies_phi_fractions_file; bool ies_phi_factors_by_real; - vector ies_n_iter_mean; + vector ies_n_iter_reinflate; + vector ies_reinflate_factor; bool ies_updatebyreals;