diff --git a/documentation/media/image7.png b/documentation/media/image7.png index cae621ec1..319e4c644 100644 Binary files a/documentation/media/image7.png and b/documentation/media/image7.png differ diff --git a/documentation/media/image8.png b/documentation/media/image8.png new file mode 100644 index 000000000..cae621ec1 Binary files /dev/null and b/documentation/media/image8.png differ diff --git a/documentation/pestpp_users_guide_v5.2.12.docx b/documentation/pestpp_users_guide_v5.2.13.docx similarity index 70% rename from documentation/pestpp_users_guide_v5.2.12.docx rename to documentation/pestpp_users_guide_v5.2.13.docx index 39ec256d9..e49b24cd4 100644 Binary files a/documentation/pestpp_users_guide_v5.2.12.docx and b/documentation/pestpp_users_guide_v5.2.13.docx differ diff --git a/documentation/pestpp_users_manual.md b/documentation/pestpp_users_manual.md index 8b6850cda..c09da4f45 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.2.12 +# Version 5.2.13 PEST++ Development Team -June 2024 +October 2024 # Acknowledgements @@ -70,7 +70,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI # Table of Contents -- [Version 5.2.12](#s1) +- [Version 5.2.13](#s1) - [Acknowledgements](#s2) - [Preface](#s3) - [License](#s4) @@ -253,6 +253,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI - [9.2.8 Reporting ](#s13-2-8) - [9.2.9 Termination Criteria, Objective Functions, and Upgrade Acceptance ](#s13-2-9) - [9.2.10 Internal Weight Adjustment ](#s13-2-10) + - [9.2.11 Selective Updates ](#s13-2-11) - [9.3 PESTPP-IES Output Files](#s13-3) - [9.3.1 CSV Output Files](#s13-3-1) - [9.3.2 Non-CSV Output Files](#s13-3-2) @@ -3821,6 +3822,14 @@ When using a weight ensemble and the multi-modal solution process, users may als If users want to have more fine-grained control of the weight adjustment, option are available in both pyEMU and the PEST utilities. +### 9.2.11 Selective Updates + +In highly nonlinear settings, some realizations may show an increase in phi across iterations, while the majority of realizations shows decreases – see Figure 9.4 for an example. This indicates that, because of nonlinearity, the optimal parameter ensemble update that is effective at reducing phi for the ensemble as a whole is not ideal for all realizations. To combat this problem, PESTPP-IES will, by default (as of version 5.2.13), only update realizations that meet the phi reduction criteria. This is identical to the partial upgrade process that is triggered automatically of the mean phi of the entire ensemble does not meet the phi reduction criteria. This behaviour is controlled with the *ies_update_by_reals* option, which is true by default. + +A graph of a graph Description automatically generated with medium confidence + +Figure 9.4 – The ensemble upgrade from iteration 2 to iteration 3 shows that some realizations have an increase in phi values, while most realizations show a decrease in phi. + ## 9.3 PESTPP-IES Output Files ### 9.3.1 CSV Output Files @@ -4152,7 +4161,12 @@ Note also that the number of control variables may change with time. Refer to th ies_n_iter_mean int -The number of mean-shift iterations to use. Default is 0- +The number of mean-shift iterations to use. Default is 0, which is indicates not to do the prior-mean-shifting process + + +ies_update_by_reals +bool +Flag to indicate whether or not to update each realization according to its phi reduction. @@ -4722,7 +4736,7 @@ In this way, the string-based cycle values allow users to apply sophisticated ru Although PESTPP-DA is a tool designed for flexible sequential and batch data assimilation, the generalized nature of the cycle concept, in concert with the observation and weight cycle tables, also provides a range of other functionality. In this way, the cycle concept can be thought of as an outer iteration process. For example, users can undertake the advanced “direct predictive hypothesis testing” analysis (e.g., Moore et al., 2010) with PESTPP-DA by constructing a generic weight cycle table where each cycle includes increasing weight on a control file observation quantity that represents a simulated outcome of interest. For example, assume a model has been constructed to simulate surface-water/groundwater exchange (SGE) along an important river reach. Further assume that the simulated SGE along this reach is included in the control file as an observation. To test the hypothesis that the SGE for this reach could be zero, users should set the observation value quantity in the control file to 0.0 and set the weight to 1.0 (this weight will not be used but simply activates this quantity in the PESTPP-DA cycle process). Now users can construct a weight cycle table. Let’s use 10 cycles. For the historic observations that are being assimilated, the entries for all cycles in the weight cycle table for these observations should be identical to the weights in the control file. The entries for the SGE “observation” in the weight cycle table should slow increase from 0.0 in the first cycle to a value large enough to dominate the objective function in the last cycle. Conceptually, during each PESTPP-DA “cycle”, a (iterative) ensemble smoother formulation will be used to minimize the objective function, but as cycles progress, the desire to force the SGE towards zero increasingly features in the objective function. In this way, the compatibility between the fitting the historic observations and the ability to make SGE be zero is directly tested. If the ability to fit the past observations is maintained while also making the simulated SGE zero, then one cannot reject the hypothesis that the SGE could be zero on the basis of compatibility with historic observations. This technique is very similar to “pareto mode” in PEST(\_HP), except here, we can take advantage of the computational efficiency of the iterative ensemble solver in PESTPP-DA. Figure 12.XXX depicts the results of such an analysis -Chart, scatter chart Description automatically generated +Chart, scatter chart Description automatically generated Figure 12.XXX. Results of a direct predictive hypothesis testing analysis where the relation between fitting historic observations and a desire to make surface-water/groundwater exchange (SGE) zero is evaluated. The ensemble-based pareto trade-off between these two quantities shows that simulating an SGE of zero is not compatible with the historic observations. diff --git a/src/libs/common/config_os.h b/src/libs/common/config_os.h index ace10f9a8..353586c33 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.2.12"; +#define PESTPP_VERSION "5.2.13"; #if defined(_WIN32) || defined(_WIN64) #define OS_WIN diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp index 58d5b1c7b..727af4e0d 100644 --- a/src/libs/pestpp_common/EnsembleMethodUtils.cpp +++ b/src/libs/pestpp_common/EnsembleMethodUtils.cpp @@ -7089,12 +7089,17 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto continue; } + if (pest_scenario.get_pestpp_options().get_ies_updatebyreals()) + { + message(1, "updating realizations with reduced phi"); + update_reals_by_phi(pe_lams[i], oe_lams[i],subset_idxs); + } + ph.update(oe_lams[i], pe_lams[i], weights); message(0, "phi summary for lambda, scale fac:", vals, echo); ph.report(echo); - mean = ph.get_mean(L2PhiHandler::phiType::COMPOSITE); std = ph.get_std(L2PhiHandler::phiType::COMPOSITE); ph.write_lambda(iter,oe_lams[i].shape().first,last_best_lam,last_best_mean, @@ -7110,8 +7115,8 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto } if (best_idx == -1) { - message(0, "WARNING: unsuccessful upgrade testing, resetting lambda to 10000.0 and returning to upgrade calculations"); - last_best_lam = 10000.0; + message(0, "WARNING: unsuccessful upgrade testing, multiplying lambda by 10000.0 and returning to upgrade calculations"); + last_best_lam *= 10000.0; return false; } @@ -7121,7 +7126,6 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto if ((best_idx != -1) && (use_subset) && (local_subset_size < pe.shape().first)) { - double acc_phi = last_best_mean * acc_fac; if (pest_scenario.get_pestpp_options().get_ies_debug_high_subset_phi()) @@ -7348,9 +7352,15 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto message(0, "updating parameter ensemble"); performance_log->log_event("updating parameter ensemble"); last_best_mean = best_mean; - - pe = pe_lams[best_idx]; - oe = oe_lam_best; + if (pest_scenario.get_pestpp_options().get_ies_updatebyreals()) + { + message(0, "only updating realizations with reduced phi"); + update_reals_by_phi(pe_lams[best_idx], oe_lam_best); + } + else { + pe = pe_lams[best_idx]; + oe = oe_lam_best; + } if (best_std < last_best_std * acc_fac) { double new_lam = lam_vals[best_idx] * lam_dec; @@ -7368,11 +7378,10 @@ bool EnsembleMethod::solve(bool use_mda, vector inflation_factors, vecto else { //message(0, "not updating parameter ensemble"); - if (!use_mda) + if ((!use_mda) && (!pest_scenario.get_pestpp_options().get_ies_updatebyreals())) { message(0, "only updating realizations with reduced phi"); update_reals_by_phi(pe_lams[best_idx], oe_lam_best); - } ph.update(oe, pe, weights); //re-check phi diff --git a/src/libs/pestpp_common/pest_data_structs.cpp b/src/libs/pestpp_common/pest_data_structs.cpp index 3bea7f892..cfad3779b 100644 --- a/src/libs/pestpp_common/pest_data_structs.cpp +++ b/src/libs/pestpp_common/pest_data_structs.cpp @@ -1186,6 +1186,11 @@ bool PestppOptions::assign_ies_value_by_key(const string& key, const string& val convert_ip(value,ies_n_iter_mean); return true; } + else if (key == "IES_UPDATE_BY_REALS") + { + ies_updatebyreals = pest_utils::parse_string_arg_to_bool(value); + return true; + } return false; @@ -1824,9 +1829,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_mean: " << ies_n_iter_mean << endl; + os << "ies_updatebyreals: " << ies_updatebyreals << endl; + - os << endl << "pestpp-sen options: " << endl; + os << endl << "pestpp-sen options: " << endl; os << "gsa_method: " << gsa_method << endl; os << "gsa_morris_pooled_obs: " << gsa_morris_pooled_obs << endl; os << "gsa_morris_obs_sen: " << gsa_morris_obs_sen << endl; @@ -2010,6 +2017,7 @@ void PestppOptions::set_defaults() set_ies_phi_fractions_files(""); set_ies_phi_factors_by_real(false); set_ies_n_iter_mean(0); + set_ies_updatebyreals(true); // DA parameters diff --git a/src/libs/pestpp_common/pest_data_structs.h b/src/libs/pestpp_common/pest_data_structs.h index 0201bf379..b51331d09 100644 --- a/src/libs/pestpp_common/pest_data_structs.h +++ b/src/libs/pestpp_common/pest_data_structs.h @@ -492,6 +492,8 @@ class PestppOptions { void set_ies_phi_fractions_files(string _file) {ies_phi_fractions_file = _file;} bool get_ies_phi_factors_by_real() const {return ies_phi_factors_by_real;} void set_ies_phi_factors_by_real(bool _flag) {ies_phi_factors_by_real = _flag;} + bool get_ies_updatebyreals() const {return ies_updatebyreals;} + void set_ies_updatebyreals(bool _flag) {ies_updatebyreals = _flag;} @@ -823,6 +825,7 @@ class PestppOptions { string ies_phi_fractions_file; bool ies_phi_factors_by_real; int ies_n_iter_mean; + bool ies_updatebyreals; // Data Assimilation parameters diff --git a/src/libs/run_managers/abstract_base/model_interface.cpp b/src/libs/run_managers/abstract_base/model_interface.cpp index dfa314b4c..2702ea73b 100644 --- a/src/libs/run_managers/abstract_base/model_interface.cpp +++ b/src/libs/run_managers/abstract_base/model_interface.cpp @@ -718,6 +718,10 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ cout << "exit_code: " << exitcode << endl; throw std::runtime_error("GetExitCodeProcess() returned error status for command: " + cmd_string); } + if (should_echo) + { + cout << "...exit_code: " << exitcode << endl; + } break; } //else cout << exitcode << "...still waiting for command " << cmd_string << endl; @@ -765,16 +769,22 @@ void ModelInterface::run(pest_utils::thread_flag* terminate, pest_utils::thread_ //check if process is still active int status = 0; pid_t exit_code = waitpid(command_pid, &status, WNOHANG); - //if the process ended, break + + //if the process ended, break if ((exit_code == -1) || (status != 0)) { + cout << "exit_code: " << exit_code << endl; + cout << "status: " << status << endl; finished->set(true); - cout << "exit_code: " << exit_code << endl; - cout << "status: " << status << endl; throw std::runtime_error("waitpid() returned error status for command: " + cmd_string); } else if (exit_code != 0) { + if (should_echo) + { + cout << "...exit_code: " << exit_code << endl; + cout << "...status: " << status << endl; + } break; } //check for termination flag