Merge pull request #71 from jwhite-usgs/feat_align_ensembles

Feat align ensembles

scripts/build_pestpp_mac.sh

@@ -1,5 +1,7 @@
 #!/bin/sh
-
+realpath() {
+    [[ $1 = /* ]] && echo "$1" || echo "$PWD/${1#./}"
+}
 set -e
 
 full_path=$(realpath $0)

src/libs/common/config_os.h

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

src/libs/common/utilities.cpp

@@ -1142,7 +1142,7 @@ void ExternalCtlFile::read_file()
 
 	getline(f_in, org_header_line);
 	header_line = upper_cp(org_header_line);
-
+	f_rec << "...header line: " << header_line << endl;
 	tokenize(header_line, col_names, delim,false);
 	int hsize = col_names.size();
 	set<string> tset;
@@ -1194,6 +1194,7 @@ void ExternalCtlFile::read_file()
 		strip_ip(next_line, "both");
 		if ((next_line.size() == 0) || (next_line.substr(0,1) == "#"))
 		{
+			f_rec << "... skipping comment line in external file:" << endl << org_next_line << endl;
 			lcount++;
 			continue;
 		}
@@ -1248,6 +1249,7 @@ void ExternalCtlFile::read_file()
 		row_order.push_back(lcount - 1);
 		lcount++;
 	}
+	f_rec << "...read " << lcount << " lines from external file" << endl;
 
 }
 

src/libs/pestpp_common/Ensemble.cpp

@@ -20,6 +20,27 @@ rand_gen_ptr(_rand_gen_ptr)
 {
 }
 
+bool Ensemble::try_align_other_rows(PerformanceLog* performance_log, Ensemble& other)
+{
+	map<string, int> other_real_map = other.get_real_map(), this_real_map = get_real_map();
+	bool need_to_align = false;
+	for (auto item : this_real_map)
+	{
+		//if even one realization name is not common, we are done here
+		if (other_real_map.find(item.first) == other_real_map.end())
+			return false;
+		//if this realization name is not at the same location then we do need to reorder
+		if (item.second != other_real_map[item.first])
+			need_to_align = true;
+	}
+	//if we made it to here and need_to_align is false, then two ensembles are already row aligned
+	if (!need_to_align)
+		return false;
+	performance_log->log_event("Ensemble::try_align_other_rows(): reorderig other ensemble");
+	other.reorder(real_names, vector<string>(),true);
+	return true;
+}
+
 void Ensemble::check_for_dups()
 {
 	vector<string> dups;
@@ -650,14 +671,18 @@ void Ensemble::set_eigen(Eigen::MatrixXd _reals)
 }
 
 
-void Ensemble::reorder(const vector<string> &_real_names, const vector<string> &_var_names)
+void Ensemble::reorder(const vector<string> &_real_names, const vector<string> &_var_names, bool update_org_real_names)
 {
 	//reorder inplace
 	reals = get_eigen(_real_names, _var_names);
 	if (_var_names.size() != 0)
 		var_names = _var_names;
 	if (_real_names.size() != 0)
+	{
 		real_names = _real_names;
+		if (update_org_real_names)
+			org_real_names = _real_names;
+	}
 }
 
 void Ensemble::drop_rows(const vector<int> &row_idxs)

src/libs/pestpp_common/Ensemble.h

@@ -69,7 +69,7 @@ class Ensemble
 
 	Covariance get_diagonal_cov_matrix();
 	pair<Covariance,Covariance> get_empirical_cov_matrices(FileManager* file_manager_ptr);
-	void reorder(const vector<string> &_real_names, const vector<string> &_var_names);
+	void reorder(const vector<string> &_real_names, const vector<string> &_var_names, bool update_org_real_names=false);
 	void drop_rows(const vector<int> &row_idxs);
 	void drop_rows(const vector<string> &drop_names);
 	void drop_cols(const vector<string>& drop_names);
@@ -93,6 +93,8 @@ class Ensemble
 	map<string, int> get_var_map() { return var_map; }
 	map<string, int> get_real_map();
 
+	bool try_align_other_rows(PerformanceLog* performance_log, Ensemble& other);
+
 protected:
 	std::mt19937* rand_gen_ptr;
 	Pest* pest_scenario_ptr;

src/libs/pestpp_common/EnsembleMethodUtils.cpp

@@ -931,24 +931,6 @@ void ParChangeSummarizer::summarize(ParameterEnsemble &pe, int iiter)
 	}
 	cout << endl;
 	frec << endl;
-	vector<string> rnames = pe.get_real_names();
-
-	if (find(rnames.begin(), rnames.end(), "BASE") != rnames.end())
-	{
-		Eigen::VectorXd v = pe.get_real_vector("BASE");
-
-		Parameters pars = pe.get_pest_scenario_ptr()->get_ctl_parameters();
-		pars.update(pe.get_var_names(), eigenvec_2_stlvec(v));
-		pe.get_pest_scenario_ptr()->get_base_par_tran_seq().numeric2ctl_ip(pars);
-		// save parameters to .par file
-		ss.str("");
-		ss << file_manager_ptr->get_base_filename() << "." << iiter << ".par";
-		string filename = ss.str();
-		output_file_writer_ptr->write_par(file_manager_ptr->open_ofile_absolute("par",filename), pars, *(pe.get_pest_scenario_ptr()->get_base_par_tran_seq().get_offset_ptr()),
-			*(pe.get_pest_scenario_ptr()->get_base_par_tran_seq().get_scale_ptr()));
-		file_manager_ptr->close_file("par");
-	}
-
 }
 
 

src/libs/pestpp_common/EnsembleSmoother.cpp

@@ -978,24 +978,7 @@ void IterEnsembleSmoother::initialize()
 		message(0, "control file parameter phi report:");
 		ph.report(true);
 		ph.write(0, 1);
-		save_base_real_par_rei(pest_scenario, _pe, _oe, output_file_writer, file_manager, -1);
-		//ObjectiveFunc obj_func(&(pest_scenario.get_ctl_observations()), &(pest_scenario.get_ctl_observation_info()), &(pest_scenario.get_prior_info()));
-		//Observations obs;
-		//Eigen::VectorXd v = _oe.get_real_vector("BASE");
-		//vector<double> vv;
-		//vv.resize(v.size());
-		//Eigen::VectorXd::Map(&vv[0], v.size()) = v;
-		//obs.update(_oe.get_var_names(), vv);
-
-		//// save parameters to .par file
-		//output_file_writer.write_par(file_manager.open_ofile_ext("base.par"),pars, *(pts.get_offset_ptr()),
-		//	*(pts.get_scale_ptr()));
-		//file_manager.close_file("par");
-
-		//// save new residuals to .rei file
-		//output_file_writer.write_rei(file_manager.open_ofile_ext("base.rei"), 0,
-		//	pest_scenario.get_ctl_observations(),obs,obj_func,pars);
-
+		save_base_real_par_rei(pest_scenario, _pe, _oe, output_file_writer, file_manager, -1);			
 		return;
 	}
 
@@ -1189,6 +1172,39 @@ void IterEnsembleSmoother::initialize()
 		else
 			add_bases();
 
+
+	//now we check to see if we need to try to align the par and obs en
+	//this would only be needed if either of these were not drawn
+	if (!pe_drawn || !oe_drawn)
+	{
+		bool aligned = pe.try_align_other_rows(performance_log, oe);
+		if (aligned)
+		{
+			message(2, "observation ensemble reordered to align rows with parameter ensemble");
+		}
+	}
+
+	//just check to see if common real names are found but are not in the same location
+	map<string, int> pe_map = pe.get_real_map(), oe_map = oe.get_real_map();
+	vector<string> misaligned;
+	for (auto item : pe_map)
+	{
+		if (oe_map.find(item.first) == oe_map.end())
+			continue;
+		if (item.second != oe_map[item.first])
+			misaligned.push_back(item.first);
+	}
+	if (misaligned.size() > 0)
+	{
+		message(1, "WARNING: common realization names shared between the parameter and observation ensembles but they are not in the same row locations, see .rec file for listing");
+		ofstream& frec = file_manager.rec_ofstream();
+		frec << endl <<  "WARNING: the following " << misaligned.size() << " realization names are shared between the parameter and observation ensembles but they are not in the same row locations:" << endl;
+		for (auto ma : misaligned)
+			frec << ma << endl;
+	}
+
+
+
 	message(2, "checking for denormal values in pe");
 	pe.check_for_normal("initial transformed parameter ensemble");
 	ss.str("");

src/libs/pestpp_common/SVDSolver.cpp

@@ -178,9 +178,9 @@ ModelRun SVDSolver::solve(RunManagerAbstract &run_manager, TerminationController
 			termination_ctl.save_state(fout_restart);
 			//write current parameters so we have a backup for restarting
 			RestartController::write_start_parameters_updated(fout_restart, file_manager.build_filename("parb", false));
-			output_file_writer.write_par(file_manager.open_ofile_ext("parb"), best_upgrade_run.get_ctl_pars(), *(par_transform.get_offset_ptr()),
+			output_file_writer.write_par(file_manager.open_ofile_ext("par"), best_upgrade_run.get_ctl_pars(), *(par_transform.get_offset_ptr()),
 				*(par_transform.get_scale_ptr()));
-			file_manager.close_file("parb");
+			file_manager.close_file("par");
 			RestartController::write_finish_parameters_updated(fout_restart, file_manager.build_filename("parb", false));
 
 			// write header for SVD file
@@ -294,7 +294,7 @@ ModelRun SVDSolver::solve(RunManagerAbstract &run_manager, TerminationController
 		file_manager.close_file("par");
 
 		filename.str(""); // reset the stringstream
-		filename << "par" << global_iter_num;
+		filename << global_iter_num << ".par";
 		output_file_writer.write_par(file_manager.open_ofile_ext(filename.str()), best_upgrade_run.get_ctl_pars(), *(par_transform.get_offset_ptr()),
 			*(par_transform.get_scale_ptr()));
 		file_manager.close_file(filename.str());