Merge branch 'feat_sqp' into develop

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.17
+# <a id='s1' />Version 5.1.18
 
 <img src="0e14ec9848f78a9809081572ca785af9990c2d38.png" style="width:6.26806in;height:3.05972in" />
 
 PEST++ Development Team
 
-June 2022
+July 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.17](#s1)
+- [Version 5.1.18](#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.17";
+#define PESTPP_VERSION "5.1.18";
 
 #if defined(_WIN32) || defined(_WIN64)
 #define OS_WIN

src/libs/common/utilities.cpp

@@ -377,7 +377,7 @@ map<string, double> read_twocol_ascii_to_map(string filename, int header_lines,
 {
 	map<string, double> result;
 	ifstream fin(filename);
-	if (!fin.good())
+    if (!fin.good())
 		throw runtime_error("could not open file " + filename + " for reading");
 	string line;
 	double value;
@@ -406,7 +406,7 @@ vector<string> read_onecol_ascii_to_vector(std::string filename)
 {
 	vector<string> result;
 	ifstream fin(filename);
-	if (!fin.good())
+    if (!fin.good())
 		throw runtime_error("could not open file " + filename + " for reading");
 	string line;
 	vector<string> tokens;
@@ -430,7 +430,7 @@ void read_res(string& res_filename, Observations& obs)
 {
 	map<string, double> result;
 	ifstream fin(res_filename);
-	if (!fin.good())
+    if (!fin.good())
 		throw runtime_error("could not open residuals file " + res_filename + " for reading");
 	vector<string> tokens;
 	string line, name;
@@ -720,7 +720,7 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 	stringstream ss;
 	ifstream in;
 	in.open(filename.c_str(), ifstream::binary);
-	if (in.bad())
+    if (!in.good())
 	{
 		ss.str("");
 		ss << "read_dense_binary() error opening binary file " << filename << " for reading";
@@ -754,7 +754,7 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 		for (int i = 0; i < n_obs_and_pi; i++)
 		{
 			in.read((char*)&(name_size), sizeof(name_size));
-			if (in.bad())
+            if (!in.good())
 			{
 				ss.str("");
 				ss << "read_dense_binary(), dense format error reading size column name size for column number " << i;
@@ -769,7 +769,7 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 		{
 			char* col_name = new char[col_name_size];
 			in.read(col_name, col_name_size);
-			if (in.bad())
+            if (!in.good())
 			{
 				ss.str("");
 				ss << "read_dense_binary(), dense format error reading column name for column number " << i << ", size " << col_name_size;
@@ -790,13 +790,14 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 		while (true)
 		{
 			//finished
-			if ((in.bad()) || (in.eof()))
+			//if ((in.bad()) || (in.eof()))
+			if ((i > 0) && (!in.good()))
 			{
 				break;
 			}
 
 			in.read((char*)&(name_size), sizeof(name_size));
-			if (in.bad())
+            if (!in.good())
 			{
 				ss.str("");
 				ss << "read_dense_binary(), dense format incomplete record: error reading row name size for row number " << i << "...continuing";
@@ -805,7 +806,7 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 			}
 			char* row_name = new char[name_size];
 			in.read(row_name, name_size);
-			if (in.bad())
+            if (!in.good())
 			{
 				ss.str("");
 				ss << "read_dense_binary(), dense format incomplete record: error reading row name for row number " << i << "...continuing";
@@ -815,7 +816,7 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 			name = string(row_name, name_size);
 			pest_utils::strip_ip(name);
 			pest_utils::upper_ip(name);
-			if (in.bad())
+            if (!in.good())
 			{
 				ss.str("");
 				ss << "read_dense_binary(), dense format incomplete record: error skipping values for row  " << i << "...continuing ";
@@ -828,14 +829,15 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 			in.read(row_name, sizeof(double)* col_names.size());
 			if (in.eof())
 				break;
-			if (in.bad())
+            if (!in.good())
 			{
 				ss.str("");
 				ss << "read_dense_binary(), dense format incomplete record: error skipping values for row  " << i << "...continuing ";
 				cout << ss.str();
 				break;
 			}
 			row_names.push_back(name);
+            i++;
 		}
 
 		in.close();
@@ -852,7 +854,7 @@ void read_dense_binary(const string& filename, vector<string>& row_names, vector
 			in.seekg(sizeof(int) + row_names[i].size(), ios_base::cur);
 			for (int j = 0; j < col_names.size(); j++)
 			{
-				if (in.bad())
+                if (!in.good())
 				{
 					ss.str("");
 					ss << "read_dense_binary(), dense format incomplete record: error reading row,col value  " << i << "," << j << "...continuing ";
@@ -871,7 +873,7 @@ void read_binary_matrix_header(const string& filename, int& tmp1, int& tmp2, int
 	stringstream ss;
 	ifstream in;
 	in.open(filename.c_str(), ifstream::binary);
-	if (in.bad())
+    if (!in.good())
 	{
 		ss.str("");
 		ss << "pest_utils::read_binary_matrix_header() error opening binary file " << filename << " for reading";
@@ -880,7 +882,7 @@ void read_binary_matrix_header(const string& filename, int& tmp1, int& tmp2, int
 	in.read((char*)&tmp1, sizeof(tmp1));
 	in.read((char*)&tmp2, sizeof(tmp2));
 	in.read((char*)&tmp3, sizeof(tmp3));
-	if (in.bad())
+    if (!in.good())
 	{
 		ss.str("");
 		ss << "pest_utils::read_binary_matrix_header() error header from binary file " << filename << " for reading";
@@ -929,7 +931,7 @@ bool read_binary(const string &filename, vector<string> &row_names, vector<strin
 	stringstream ss;
 	ifstream in;
 	in.open(filename.c_str(), ifstream::binary);
-	if (in.bad())
+	if (!in.good())
 	{
 		ss.str("");
 		ss << "pest_utils::read_binary() error opening binary file " << filename << " for reading";
@@ -1376,7 +1378,7 @@ void ExternalCtlFile::read_file(ofstream& f_rec)
 		throw_externalctrlfile_error("filename '" + filename + "' not found");
 	}
 	ifstream f_in(filename);
-	if (!f_in.good())
+    if (!f_in.good())
 	{
 		throw_externalctrlfile_error("error opening filename '" + filename + "' for reading");
 	}

src/libs/pestpp_common/Ensemble.cpp

@@ -422,7 +422,7 @@ pair<Covariance,Covariance> Ensemble::get_empirical_cov_matrices(FileManager* fi
 	}
 	double scale = (num_reals / ((num_reals - 1.) * (num_reals - 1.) * (num_reals - 1.))) * wij_sum;
 	scale = scale / demon;
-	cout << "optimal residual covariance matrix shrinkage factor: " << scale << endl;
+	//cout << "optimal residual covariance matrix shrinkage factor: " << scale << endl;
 	file_manager_ptr->rec_ofstream() << "optimal residual covariance matrix shrinkage factor : " << scale << endl;
 
 	Covariance rcov_diag;
@@ -1977,6 +1977,7 @@ ParameterEnsemble::ParameterEnsemble(Pest *_pest_scenario_ptr, std::mt19937* _ra
 	reals = _reals;
 	var_names = _var_names;
 	real_names = _real_names;
+	org_real_names = _real_names;
 	tstat = transStatus::CTL;
 	set_fixed_names();
 }
@@ -2258,7 +2259,7 @@ void ParameterEnsemble::from_binary(string file_name, bool forgive)
 	map<string,int> header_info = Ensemble::from_binary(file_name, names, false);
 	unordered_set<string>svar_names(var_names.begin(), var_names.end());
 	vector<string> missing;
-	for (auto& name : pest_scenario_ptr->get_ctl_ordered_adj_par_names())
+	for (auto& name : names)
 	{
 		if (svar_names.find(name) == svar_names.end())
 		{
@@ -2282,6 +2283,25 @@ void ParameterEnsemble::from_binary(string file_name, bool forgive)
 			throw_ensemble_error("from_binary() error: the following adjustable parameter names in the control file are not in the binary parameter ensemble file:", missing);
 		}
 	}
+	missing.clear();
+	svar_names.clear();
+    names = pest_scenario_ptr->get_ctl_ordered_par_names();
+	svar_names.insert(names.begin(),names.end());
+	unordered_set<string>::iterator send = svar_names.end();
+	for (auto& name: var_names)
+    {
+	    if (svar_names.find(name) == send)
+        {
+	        missing.push_back(name);
+        }
+    }
+	if (missing.size() > 0)
+    {
+	    drop_cols(missing);
+    }
+
+
+
 
 	prep_par_ensemble_after_read(header_info);
 }
@@ -3735,6 +3755,23 @@ void ObservationEnsemble::from_binary(string file_name)
 	if (missing.size() > 0)
 		throw_ensemble_error("from_binary() error: the following non-zero-weighted obs names in the control file are not in the binary obs ensemble file:", missing);
 	names = pest_scenario_ptr->get_ctl_ordered_obs_names();
+	missing.clear();
+	svar_names.clear();
+	svar_names.insert(names.begin(),names.end());
+	unordered_set<string>::iterator send = svar_names.end();
+
+    for (auto& name : var_names)
+    {
+        if (svar_names.find(name) == send)
+        {
+            missing.push_back(name);
+        }
+    }
+    if (missing.size() > 0)
+    {
+        //drop these extra vars
+        drop_cols(missing);
+    }
 	if (var_names.size() < names.size())
 	{
 		update_var_map();

src/libs/pestpp_common/EnsembleMethodUtils.cpp

@@ -4133,6 +4133,10 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing)
         vector<int> failed = run_ensemble(pe, oe, vector<int>(), cycle);
         if (pe.shape().first == 0)
             throw_em_error("all realizations failed during initial evaluation");
+        if (pest_scenario.get_pestpp_options().get_ies_debug_fail_remainder())
+        {
+            failed.push_back(0);
+        }
         if (failed.size() > 0)
         {
             ss.str("");
@@ -4200,7 +4204,7 @@ void EnsembleMethod::initialize(int cycle, bool run, bool use_existing)
 		ss << " with the prior simulated ensemble." << endl;
 		message(0, ss.str());
 
-		cout << "...see rec file or " << file_manager.get_base_filename() << ".pdc.csv" << "for listing of conflicted observations" << endl << endl;
+		cout << "...see rec file or " << file_manager.get_base_filename() << ".pdc.csv" << " for listing of conflicted observations" << endl << endl;
 		ofstream& frec = file_manager.rec_ofstream();
 		frec << endl << "...conflicted observations: " << endl;
 		for (auto oname : in_conflict)

src/libs/pestpp_common/Jacobian_1to1.cpp

@@ -296,11 +296,12 @@ bool Jacobian_1to1::process_runs(ParamTransformSeq &par_transform,
 			run_manager.get_model_parameters(par_run.second[0], run_list.back().ctl_pars);
 			bool success = run_manager.get_observations_vec(par_run.second[0], run_list.back().obs_vec);
 			run_list.back().numeric_derivative_par = cur_numeric_par_value;
-			/*if ((debug_fail) && (i_run == 1))
+			if (debug_fail)
 			{
 				file_manager.rec_ofstream() << "NOTE: 'GLM_DEBUG_DER_FAIL' is true, failing jco run for parameter '" << cur_par_name << "'" << endl;
 				success = false;
-			}*/
+				debug_fail = false;
+			}
 			if (success)
 			{
 				par_transform.model2ctl_ip(run_list.back().ctl_pars);

src/libs/pestpp_common/Regularization.cpp

@@ -102,8 +102,11 @@ PestppOptions::ARG_STATUS DynamicRegularization::assign_value_by_key(const std::
 		pest_utils::convert_ip(value, wf_min);
 	else if (key == "WFMAX")
 		pest_utils::convert_ip(value, wf_max);
-	else if (key == "WFINIT")
-		pest_utils::convert_ip(value, wf_init);
+	else if (key == "WFINIT") {
+        pest_utils::convert_ip(value, wf_init);
+        tikhonov_weight = wf_init;
+    }
+
 	else if (key == "WFTOL")
 		pest_utils::convert_ip(value, wftol);
 	else if (key == "WFFAC")
@@ -140,6 +143,7 @@ void DynamicRegularization::set_defaults()
 	wffac = 0;
 	wftol = 1000;
 	wf_init = 1.0;
+	tikhonov_weight = 1.0;
 	max_reg_iter = 20;
 }