Merge remote-tracking branch 'upstream/develop'

README.md

@@ -2,6 +2,8 @@
   <img src="documentation/pestpplogo.png" alt="pestpplogo image">
 </p>
 
+
+
 # PEST++
 
 ## a Software Suite for Parameter Estimation, Uncertainty Analysis, Management Optimization and Sensitivity Analysis

documentation/pestpp_users_manual.md

@@ -1,13 +1,13 @@
 
  <img src="0d3cb7750c90b712af04ea3a51c8ecb968d784cc.png" style="width:6.26806in;height:1.68194in" alt="A close up of a purple sign Description automatically generated" />
 
-# <a id='s1' />Version 5.1.18
+# <a id='s1' />Version 5.1.20
 
 <img src="0e14ec9848f78a9809081572ca785af9990c2d38.png" style="width:6.26806in;height:3.05972in" />
 
 PEST++ Development Team
 
-July 2022
+August 2022
 
 # <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.1.18](#s1)
+- [Version 5.1.20](#s1)
 - [Acknowledgements](#s2)
 - [Preface](#s3)
 - [License](#s4)

src/libs/common/config_os.h

@@ -2,7 +2,7 @@
 #define CONFIG_OS_H_
 
 
-#define PESTPP_VERSION "5.1.18";
+#define PESTPP_VERSION "5.1.20";
 
 #if defined(_WIN32) || defined(_WIN64)
 #define OS_WIN

src/libs/pestpp_common/Ensemble.cpp

@@ -2173,6 +2173,7 @@ map<int,int> ParameterEnsemble::add_runs(RunManagerAbstract *run_mgr_ptr,const v
 	{
 		par_transform.active_ctl2model_ip(pars);
 	}
+    //map<string, pair<string, double>> tied_items = par_transform.get_tied_ptr()->get_items();
 	Parameters pars_real = pars;
 	Eigen::VectorXd evec;
 	vector<double> svec;
@@ -2937,6 +2938,10 @@ void ParameterEnsemble::to_binary_ordered(string file_name)
 		throw runtime_error("error opening file for binary parameter ensemble:" + file_name);
 	}
 
+	bool has_tied = false;
+	if (par_transform.get_tied_ptr()->get_items().size() > 0)
+	    has_tied = true;
+
 	//vector<string> vnames = var_names;
 	vector<string> vnames = pest_scenario_ptr->get_ctl_ordered_par_names();
 	//vnames.insert(vnames.end(), fixed_names.begin(), fixed_names.end());
@@ -3010,6 +3015,11 @@ void ParameterEnsemble::to_binary_ordered(string file_name)
 			par_transform.model2ctl_ip(pars);
 		else if (tstat == transStatus::NUM)
 			par_transform.numeric2ctl_ip(pars);
+        else if ((has_tied) && (tstat == transStatus::CTL))
+        {
+            par_transform.active_ctl2model_ip(pars);
+            par_transform.model2ctl_ip(pars);
+        }
 		replace_fixed(real_names[irow], pars);
 
 		for (int jcol = 0; jcol<n_var; ++jcol)
@@ -3128,6 +3138,9 @@ void ParameterEnsemble::to_dense_ordered(string file_name)
 
 	//write matrix
 	n = 0;
+	bool has_tied = false;
+	if (par_transform.get_tied_ptr()->get_items().size() > 0)
+	    has_tied = true;
 	//map<string, double>::const_iterator found_pi_par;
 	//map<string, double>::const_iterator not_found_pi_par;
 	//icount = row_idxs + 1 + col_idxs * self.shape[0]
@@ -3145,6 +3158,11 @@ void ParameterEnsemble::to_dense_ordered(string file_name)
 			par_transform.model2ctl_ip(pars);
 		else if (tstat == transStatus::NUM)
 			par_transform.numeric2ctl_ip(pars);
+        else if ((has_tied) && (tstat == transStatus::CTL))
+        {
+            par_transform.active_ctl2model_ip(pars);
+            par_transform.model2ctl_ip(pars);
+        }
 		replace_fixed(real_names[irow], pars);
 		string name = real_names[irow];
 		tmp = name.size();
@@ -3331,7 +3349,9 @@ void ParameterEnsemble::to_csv_by_vars(ofstream &csv, bool write_header)
 	}
 	if (write_header)
 		csv << endl;
-
+    bool has_tied = false;
+    if (par_transform.get_tied_ptr()->get_items().size() > 0)
+        has_tied = true;
 	//get the pars and transform to be in sync with ensemble trans status
 	/*Parameters pars = pest_scenario_ptr->get_ctl_parameters();
 	if (tstat == transStatus::NUM)
@@ -3358,6 +3378,11 @@ void ParameterEnsemble::to_csv_by_vars(ofstream &csv, bool write_header)
 		{
 			par_transform.active_ctl2model_ip(pars);
 		}
+        else if ((has_tied) && (tstat == transStatus::CTL))
+        {
+            par_transform.active_ctl2model_ip(pars);
+            par_transform.model2ctl_ip(pars);
+        }
 		ireal = real_map[rname];
 		pars.update_without_clear(var_names, reals.row(ireal));
 		if (tstat == transStatus::MODEL)
@@ -3404,6 +3429,10 @@ void ParameterEnsemble::to_csv_by_reals(ofstream &csv, bool write_header)
 	}*/
 
 	//for (int ireal = 0; ireal < reals.rows(); ireal++)
+	bool has_tied = false;
+	if (par_transform.get_tied_ptr()->get_items().size() > 0)
+	    has_tied = true;
+
 	int ireal = 0;
 	map<string, int> real_map;
 	for (int i = 0; i < real_names.size(); i++)
@@ -3435,7 +3464,14 @@ void ParameterEnsemble::to_csv_by_reals(ofstream &csv, bool write_header)
 			par_transform.model2ctl_ip(pars);
 		else if (tstat == transStatus::NUM)
 			par_transform.numeric2ctl_ip(pars);
+		// here we need to check for tied parameters...
+		else if ((has_tied) && (tstat == transStatus::CTL))
+        {
+            par_transform.active_ctl2model_ip(pars);
+            par_transform.model2ctl_ip(pars);
+        }
 		replace_fixed(real_names[ireal],pars);
+
 		for (auto &name : names)
 			csv << ',' << pars[name];
 		csv << endl;

src/libs/pestpp_common/EnsembleMethodUtils.cpp

@@ -3883,7 +3883,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing)
 	if ((obs_restart_csv.size() == 0) && (!pe_drawn) && (oe_drawn))
 	{
 		vector<string> rnames = pe.get_real_names();
-		oe_base.set_real_names(rnames);
+		oe_base.set_real_names(rnames,true);
 		message(2, "resetting obs + noise ensemble real names to parameter ensemble real names");
 	}
 

src/libs/pestpp_common/MOEA.cpp

@@ -1358,7 +1358,8 @@ void MOEA::update_archive_nsga(ObservationEnsemble& _op, ParameterEnsemble& _dp)
 		op_archive.keep_rows(keep);
 		dp_archive.keep_rows(keep);
 	}
-
+	dp_archive.reset_org_real_names();
+	op_archive.reset_org_real_names();
 	save_populations(dp_archive, op_archive, "archive");
 }
 
@@ -1417,7 +1418,8 @@ void MOEA::update_archive_spea(ObservationEnsemble& _op, ParameterEnsemble& _dp)
 			op_archive.keep_rows(keep);
 			dp_archive.keep_rows(keep);
 		}
-
+        dp_archive.reset_org_real_names();
+        op_archive.reset_org_real_names();
 		save_populations(dp_archive, op_archive, "archive");
 	}
 }
@@ -3776,7 +3778,7 @@ ParameterEnsemble MOEA::generate_sbx_population(int num_members, ParameterEnsemb
 	return tmp_dp;
 }
 
-void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, string tag)
+void MOEA::save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag)
 {
 
 	stringstream ss;
@@ -3791,16 +3793,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri
 	if (pest_scenario.get_pestpp_options().get_save_binary())
 	{
 		ss << ".jcb";
-		dp.to_binary(ss.str());
+		_dp.to_binary(ss.str());
 	}
 	else
 	{
 		ss << ".csv";
-		dp.to_csv(ss.str());
+		_dp.to_csv(ss.str());
 	}
 	string name = ss.str();
 	ss.str("");
-	ss << " saved decision variable population of size " << dp.shape().first << " X " << dp.shape().second << " to '" << name << "'";
+	ss << " saved decision variable population of size " << _dp.shape().first << " X " << _dp.shape().second << " to '" << name << "'";
 	message(1, ss.str());
 	ss.str("");
 	if ((save_every > 0) && (iter % save_every == 0))
@@ -3814,16 +3816,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri
 		if (pest_scenario.get_pestpp_options().get_save_binary())
 		{
 			ss << ".jcb";
-			dp.to_binary(ss.str());
+			_dp.to_binary(ss.str());
 		}
 		else
 		{
 			ss << ".csv";
-			dp.to_csv(ss.str());
+			_dp.to_csv(ss.str());
 		}
 		string name = ss.str();
 		ss.str("");
-		ss << " saved generation-specific decision variable population of size " << dp.shape().first << " X " << dp.shape().second << " to '" << name << "'";
+		ss << " saved generation-specific decision variable population of size " << _dp.shape().first << " X " << _dp.shape().second << " to '" << name << "'";
 		message(1, ss.str());
 	}
 
@@ -3838,16 +3840,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri
 	if (pest_scenario.get_pestpp_options().get_save_binary())
 	{
 		ss << ".jcb";
-		op.to_binary(ss.str());
+		_op.to_binary(ss.str());
 	}
 	else
 	{
 		ss << ".csv";
-		op.to_csv(ss.str());
+		_op.to_csv(ss.str());
 	}
 	name = ss.str();
 	ss.str("");
-	ss << " saved observation population of size " << op.shape().first << " X " << op.shape().second << " to '" << name << "'";
+	ss << " saved observation population of size " << _op.shape().first << " X " << _op.shape().second << " to '" << name << "'";
 	message(1, ss.str());
 
 	if ((save_every > 0) && (iter % save_every == 0))
@@ -3862,16 +3864,16 @@ void MOEA::save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, stri
 		if (pest_scenario.get_pestpp_options().get_save_binary())
 		{
 			ss << ".jcb";
-			op.to_binary(ss.str());
+			_op.to_binary(ss.str());
 		}
 		else
 		{
 			ss << ".csv";
-			op.to_csv(ss.str());
+			_op.to_csv(ss.str());
 		}
 		name = ss.str();
 		ss.str("");
-		ss << " saved generation-specific observation population of size " << op.shape().first << " X " << op.shape().second << " to '" << name << "'";
+		ss << " saved generation-specific observation population of size " << _op.shape().first << " X " << _op.shape().second << " to '" << name << "'";
 		message(1, ss.str());
 	}
 

src/libs/pestpp_common/MOEA.h

@@ -84,8 +84,6 @@ class ParetoObjectives
 	void fill_domination_containers(map<string, map<string, double>>& _member_struct, map<string,
 		vector<string>>&solutions_dominated_map, map<string, int>& num_dominating_map, bool dup_as_dom=false);
 
-
-
 	bool compare_two_nsga(string& first, string& second);
 	bool compare_two_spea(string& first, string& second);
 
@@ -222,7 +220,7 @@ class MOEA
 
 	string get_new_member_name(string tag = string());
 
-	void save_populations(ParameterEnsemble& dp, ObservationEnsemble& op, string tag = string());
+	void save_populations(ParameterEnsemble& _dp, ObservationEnsemble& _op, string tag = string());
 	void gauss_mutation_ip(ParameterEnsemble& _dp);
 	pair<Eigen::VectorXd, Eigen::VectorXd> sbx(double probability, double eta_m, int idx1, int idx2);
 	pair<Eigen::VectorXd, Eigen::VectorXd> sbx_new(double crossover_probability, double di, Eigen::VectorXd& parent1,

src/libs/pestpp_common/constraints.cpp

@@ -1649,7 +1649,8 @@ double Constraints::get_probit()
 
 double Constraints::get_probit(double _risk)
 {
-	/* the probit function estimate - needed for the fosm-basd chance constraints*/
+	/* the probit function estimate
+	 * - needed for the fosm-basd chance constraints*/
 	double output = sqrt(2.0) * ErfInv2((2.0 * _risk) - 1.0);
 	return output;
 }
@@ -2221,9 +2222,13 @@ void Constraints::postsolve_obs_constraints_report(Observations& old_obs, Observ
 	This really only applies to the simplex solution because we can, with the linear assumption, 
 	project what the resulting constraint values should be*/
 
+	int width = 20;
+	for(auto& n : ctl_ord_obs_constraint_names)
+	    width = max(width,(int)n.size());
+	width = width + 2;
 	ofstream& f_rec = file_mgr_ptr->rec_ofstream();
 	f_rec << endl << endl << "     " << tag << " constraint/objective information at end of iteration " << iter << endl << endl;
-	f_rec << setw(20) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25);
+	f_rec << setw(width) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25);
 	if (status_map.size() > 0)
 		f_rec << "simplex status";
 	if (price_map.size() > 0)
@@ -2247,7 +2252,7 @@ void Constraints::postsolve_obs_constraints_report(Observations& old_obs, Observ
 	{
 		string name = ctl_ord_obs_constraint_names[i];
 		sim_val = new_obs[name];
-		f_rec << setw(20) << left << name;
+		f_rec << setw(width) << left << name;
 		f_rec << setw(15) << right << constraint_sense_name[name];
 		f_rec << setw(15) << constraints_obs.get_rec(name);
 		if (status_map.size() > 0)
@@ -2287,8 +2292,12 @@ void Constraints::postsolve_pi_constraints_report(Parameters& old_pars, Paramete
 	ofstream &f_rec = file_mgr_ptr->rec_ofstream();
 	if (num_pi_constraints() > 0)
 	{
+	    int width = 20;
+	    for (auto& n : ctl_ord_pi_constraint_names)
+	        width = max(width,(int)n.size());
+	    width = width + 2;
 		f_rec << endl << endl << " --- prior information constraint summary at end of iteration " << iter << " --- " << endl << endl;
-		f_rec << setw(20) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25);
+		f_rec << setw(width) << left << "name" << right << setw(15) << "sense" << setw(15) << "required" << setw(25);
 		if (status_map.size() > 0)
 			f_rec << "simplex status";
 		if (price_map.size() > 0)
@@ -2301,7 +2310,7 @@ void Constraints::postsolve_pi_constraints_report(Parameters& old_pars, Paramete
 			PriorInformationRec pi_rec = constraints_pi.get_pi_rec(name);
 			pair<double,double> cur_sim_resid = pi_rec.calc_sim_and_resid(old_pars);
 			pair<double,double> new_sim_resid = pi_rec.calc_sim_and_resid(new_pars);
-			f_rec << setw(20) << left << name;
+			f_rec << setw(width) << left << name;
 			f_rec << setw(15) << right << constraint_sense_name[name];
 			f_rec << setw(15) << pi_rec.get_obs_value();
 			if (status_map.size() > 0)