Merge remote-tracking branch 'origin/develop'

benchmarks/basic_tests.py

@@ -871,6 +871,15 @@ def mf6_v5_ies_test():
     assert os.path.exists(oe_file)
     pe_file = oe_file.replace(".obs.",".par.")
     assert os.path.exists(pe_file)
+    pcs_file = oe_file.replace(".obs.",".pcs.")
+    assert os.path.exists(pcs_file)
+    df = pd.read_csv(pcs_file,index_col=0)
+    pst_pargp = set(list(pst.parameter_data.pargp.unique()))
+    df_pargp = set(df.index.to_list())
+    d = pst_pargp.symmetric_difference(df_pargp)
+    print(d)
+    assert len(d) == 0,d
+
 
 def mf6_v5_sen_test():
     model_d = "mf6_freyberg"
@@ -959,8 +968,8 @@ def mf6_v5_glm_test():
     #sen_basic_test()
     #salib_verf()
     #tplins1_test()
-    ext_stdcol_test()
-    #mf6_v5_ies_test()
+    #ext_stdcol_test()
+    mf6_v5_ies_test()
     #mf6_v5_sen_test()
     #mf6_v5_opt_stack_test()
     #mf6_v5_glm_test()

src/libs/common/config_os.h

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

src/libs/pestpp_common/EnsembleMethodUtils.cpp

@@ -858,23 +858,98 @@ ParChangeSummarizer::ParChangeSummarizer(ParameterEnsemble *_base_pe_ptr, FileMa
 }
 
 
-void ParChangeSummarizer::summarize(ParameterEnsemble &pe, int iiter)
+void ParChangeSummarizer::summarize(ParameterEnsemble &pe, int iiter, string filename)
 {
-
-	pair<map<string, double>, map<string, double>> moments = pe.get_moment_maps();
-	init_moments = base_pe_ptr->get_moment_maps(pe.get_real_names());
+
 	stringstream ss;
 	ofstream &frec = file_manager_ptr->rec_ofstream();
-	ss << endl << "   parameter group percent change summmary" << endl;
+	ss << endl << "   ---  Parameter Group Change Summmary  ---    " << endl;
 	ss << "   (compared to the initial ensemble using active realizations)" << endl;
 	cout << ss.str();
 	frec << ss.str();
 	ss.str("");
 	ss << setw(15) << "group" << setw(12) << "mean change" << setw(12) << "std change" << setw(18) << "num at/near bnds" << setw(16) << "% at/near bnds" << endl;
 	cout << ss.str();
+	frec << ss.str();	
+	update(pe);
+	vector<string> grp_names;// = base_pe_ptr->get_pest_scenario().get_ctl_ordered_par_group_names();
+	vector<pair<double, string>> mean_pairs;
+	for (auto m : mean_change)
+	{
+		mean_pairs.push_back(pair<double, string>(abs(m.second), m.first));
+	}
+	sort(mean_pairs.begin(), mean_pairs.end());
+	//for (auto m : mean_pairs)
+	for (int i = mean_pairs.size() - 1; i >= 0; i--)
+	{
+		grp_names.push_back(mean_pairs[i].second);
+	}
+
+	int i = 0;
+	for (auto &grp_name : grp_names)
+	{
+		double mean_diff = mean_change[grp_name];
+		double std_diff = std_change[grp_name];
+		int num_out = num_at_bounds[grp_name];
+		int percent_out = percent_at_bounds[grp_name];
+		ss.str("");
+		ss << setw(15) << pest_utils::lower_cp(grp_name) << setw(12) << mean_diff * 100.0 << setw(12) << std_diff * 100.0 << setw(18);
+		ss << num_out << setw(16) << setprecision(2) << percent_out << endl;
+		if (i < 15)
+			cout << ss.str();
+		frec << ss.str();
+		i++;
+	}
+
+	ss.str("");
+	ss << "    Note: parameter change summary sorted according to abs 'mean change'." << endl;
+	cout << ss.str();
 	frec << ss.str();
+	if (grp_names.size() > 15)
+	{
+		ss.str("");
+		ss << "    Note: Only the first 15 parameter groups shown, see rec file for full listing" << endl;
+		cout << ss.str();
+	}
+
+	cout << endl;
+	frec << endl;
+
+	if (filename.size() > 0)
+		write_to_csv(filename);
+
+}
+
+void ParChangeSummarizer::write_to_csv(string& filename)
+{
+	ofstream f(filename);
+	if (f.bad())
+		throw runtime_error("ParChangeSummarizer::write_to_csv() error opening file " + filename);
+
+	f << "group,mean_change,std_change,num_at_near_bounds,percent_at_near_bounds" << endl;
+	for (auto grp_name : base_pe_ptr->get_pest_scenario_ptr()->get_ctl_ordered_par_group_names())
+	{
+		f << pest_utils::lower_cp(grp_name) << "," << mean_change[grp_name] << "," << std_change[grp_name] << ",";
+		f << num_at_bounds[grp_name] << "," << percent_at_bounds[grp_name] << endl;
+	}
+	f.close();
+	file_manager_ptr->rec_ofstream() << "...saved parameter change summary to " << filename << endl;
+	cout << "...saved parameter change summary to " << filename << endl;
+
+}
+
+
+void ParChangeSummarizer::update(ParameterEnsemble& pe)
+{
+	mean_change.clear();
+	std_change.clear();
+	num_at_bounds.clear();
+	percent_at_bounds.clear();
+	pair<map<string, double>, map<string, double>> moments = pe.get_moment_maps();
+	init_moments = base_pe_ptr->get_moment_maps(pe.get_real_names());
+
 	double mean_diff = 0.0, std_diff = 0.0;
-	double dsize, value1, value2,v;
+	double dsize, value1, value2, v;
 	vector<string> pnames = pe.get_var_names();
 	Parameters lb = pe.get_pest_scenario_ptr()->get_ctl_parameter_info().get_low_bnd(pnames);
 	pe.get_pest_scenario_ptr()->get_base_par_tran_seq().active_ctl2numeric_ip(lb);
@@ -884,25 +959,23 @@ void ParChangeSummarizer::summarize(ParameterEnsemble &pe, int iiter)
 	map<string, int> idx_map;
 	for (int i = 0; i < pnames.size(); i++)
 		idx_map[pnames[i]] = i;
-	int num_out,num_pars;
+	int num_out, num_pars;
 	int num_reals = pe.get_real_names().size();
 	Eigen::ArrayXd arr;
-	for (auto &grp_name : grp_names)
+	for (auto& grp_name : grp_names)
 	{
 		mean_diff = 0.0, std_diff = 0.0;
 		num_pars = pargp2par_map[grp_name].size();
 		num_out = 0;
-		for (auto & par_name : pargp2par_map[grp_name])
+		for (auto& par_name : pargp2par_map[grp_name])
 		{
-
 			arr = pe.get_eigen_ptr()->col(idx_map[par_name]).array();
 			for (int i = 0; i < num_reals; i++)
 			{
 				v = arr[i];
 				if ((v > (ub[par_name] * 1.01)) || (v < (lb[par_name] * 0.99)))
 					num_out++;
 			}
-
 			value1 = init_moments.first[par_name];
 			value2 = value1 - moments.first[par_name];
 			if ((value1 != 0.0) && (value2 != 0.0))
@@ -917,20 +990,17 @@ void ParChangeSummarizer::summarize(ParameterEnsemble &pe, int iiter)
 			mean_diff = mean_diff / dsize;
 		if (std_diff != 0.0)
 			std_diff = std_diff / dsize;
-
-		ss.str("");
+
 		double percent_out = 0;
 		if (num_pars > 0)
 			percent_out = double(num_out) / double(num_pars * num_reals) * 100;
-		//ss << setw(15) << "group" << setw(12) << "mean change" << setw(12) << "std change" << setw(18) << "num at/near bnds" << setw(16) << "% at/near bnds" << endl;
-		ss << setw(15) << pest_utils::lower_cp(grp_name) << setw(12) << mean_diff * 100.0 << setw(12) << std_diff * 100.0 << setw(18);
-		ss << num_out << setw(16) << setprecision(2) << percent_out << endl;
-		cout << ss.str();
-		frec << ss.str();
+
+		mean_change[grp_name] = mean_diff;
+		std_change[grp_name] = std_diff;
+		num_at_bounds[grp_name] = num_out;
+		percent_at_bounds[grp_name] = percent_out;
 
 	}
-	cout << endl;
-	frec << endl;
 }
 
 

src/libs/pestpp_common/EnsembleMethodUtils.h

@@ -105,15 +105,22 @@ class ParChangeSummarizer
 public:
 	ParChangeSummarizer() { ; }
 	ParChangeSummarizer(ParameterEnsemble *_base_pe_ptr, FileManager *_file_manager_ptr, OutputFileWriter* _output_file_writer_ptr);
-	void summarize(ParameterEnsemble &pe, int iiter);
+	void summarize(ParameterEnsemble &pe, int iiter, string filename = string());
+
 
 private:
 	ParameterEnsemble * base_pe_ptr;
 	FileManager *file_manager_ptr;
 	OutputFileWriter* output_file_writer_ptr;
 	map<string, set<string>> pargp2par_map;
 	pair<map<string,double>, map<string, double>> init_moments;
+	map<string, double> mean_change;
+	map<string, double> std_change;
+	map<string, int> num_at_bounds;
+	map<string, int> percent_at_bounds;
 
+	void update(ParameterEnsemble& pe);
+	void write_to_csv(string& filename);
 
 };
 

src/libs/pestpp_common/EnsembleSmoother.cpp

@@ -1774,7 +1774,9 @@ void IterEnsembleSmoother::iterate_2_solution()
 		ph.write(iter, run_mgr_ptr->get_total_runs());
 		if (pest_scenario.get_pestpp_options().get_ies_save_rescov())
 			ph.save_residual_cov(oe,iter);
-		pcs.summarize(pe,iter);
+		ss.str("");
+		ss << file_manager.get_base_filename() << "." << iter << ".pcs.csv";
+		pcs.summarize(pe,iter,ss.str());
 
 
 		if (accept)

src/libs/pestpp_common/OutputFileWriter.cpp

@@ -632,17 +632,35 @@ void OutputFileWriter::obs_report(ostream &os, const Observations &obs, const Ob
 	//ObservationInfo oi = pest_scenario.get_ctl_observation_info();
 	//for(vector<string>::const_iterator b = obs_name_vec.begin(),
 	//	e = obs_name_vec.end(); b!=e; ++b)
-	for (auto &b : obs_name_vec)
+	if (obs_name_vec.size() < 100000)
 	{
-		obs_val = obs.get_rec(b);
-		sim_val = sim.get_rec(b);
-		os << " " << setw(20) << lower_cp(b)
-			<< " " << setw(12) << lower_cp(oi.get_observation_rec_ptr(b)->group)
-			<< " " << showpoint << setw(20) << obs_val
-			<< " " << showpoint << setw(20) << sim_val
-			<< " " << showpoint << setw(20) << obs_val - sim_val
-			<< " " << showpoint << setw(20) << oi.get_observation_rec_ptr(b)->weight << endl;
+		for (auto& b : obs_name_vec)
+		{
+			obs_val = obs.get_rec(b);
+			sim_val = sim.get_rec(b);
+			os << " " << setw(20) << lower_cp(b)
+				<< " " << setw(12) << lower_cp(oi.get_observation_rec_ptr(b)->group)
+				<< " " << showpoint << setw(20) << obs_val
+				<< " " << showpoint << setw(20) << sim_val
+				<< " " << showpoint << setw(20) << obs_val - sim_val
+				<< " " << showpoint << setw(20) << oi.get_observation_rec_ptr(b)->weight << endl;
+		}
 	}
+	else
+	{
+		for (auto& b : obs_name_vec)
+		{
+			obs_val = obs.get_rec(b);
+			sim_val = sim.get_rec(b);
+			os << " " << lower_cp(b) 
+				<< " " << lower_cp(oi.get_observation_rec_ptr(b)->group)
+				<< " " << obs_val
+				<< " " << sim_val
+				<< " " << obs_val - sim_val
+				<< " " << oi.get_observation_rec_ptr(b)->weight << endl;
+		}
+	}
+
 
 }