Merge pull request #315 from jtwhite79/feat_upgradebyreals

Feat upgradebyreals

documentation/pestpp_users_manual.md

@@ -1,13 +1,13 @@
 
  <img src="./media/image1.png" style="width:6.26806in;height:1.68194in" alt="A close up of a purple sign Description automatically generated" />
 
-# <a id='s1' />Version 5.2.12
+# <a id='s1' />Version 5.2.13
 
 <img src="./media/image2.png" style="width:6.26806in;height:3.05972in" />
 
 PEST++ Development Team
 
-June 2024
+October 2024
 
 # <a id='s2' />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.
 
+### <a id='s13-2-11' />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.
+
+<img src="./media/image7.png" style="width:6.26806in;height:3.20347in" alt="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.
+
 ## <a id='s13-3' />9.3 PESTPP-IES Output Files
 
 ### <a id='s13-3-1' />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
 <tr class="odd">
 <td><em>ies_n_iter_mean</em></td>
 <td>int</td>
-<td>The number of mean-shift iterations to use. Default is 0-</td>
+<td>The number of mean-shift iterations to use. Default is 0, which is indicates not to do the prior-mean-shifting process</td>
+</tr>
+<tr class="even">
+<td><em>ies_update_by_reals</em></td>
+<td>bool</td>
+<td>Flag to indicate whether or not to update each realization according to its phi reduction.</td>
 </tr>
 </tbody>
 </table>
@@ -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
 
-<img src="./media/image7.png" style="width:6.26806in;height:6.29514in" alt="Chart, scatter chart Description automatically generated" />
+<img src="./media/image8.png" style="width:6.26806in;height:6.29514in" alt="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.
 

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

src/libs/pestpp_common/EnsembleMethodUtils.cpp

@@ -7089,12 +7089,17 @@ bool EnsembleMethod::solve(bool use_mda, vector<double> 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<double> 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<double> 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<double> 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<double> 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

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

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

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