Merge remote-tracking branch 'origin/develop'

benchmarks/basic_tests.py

@@ -354,6 +354,31 @@ def unc_file_test():
     print(diff.sum())
     assert diff.sum().max() < 1.0e-10
 
+    cov.to_uncfile(os.path.join(m_d, "pest.unc"), covmat_file=None)
+    pst.control_data.noptmax = -2
+    pst.pestpp_options["ies_num_reals"] = 100000
+    pst.pestpp_options["ies_enforce_bounds"] = False
+    pst.write(os.path.join(m_d, "pest_unc.pst"))
+    pyemu.os_utils.run("{0} {1}".format(exe_path, "pest_unc.pst"), cwd=m_d)
+    pe_3 = pd.read_csv(os.path.join(m_d, "pest_unc.0.par.csv"), index_col=0).apply(np.log10)
+    print(pe_3.std(ddof=0))
+    pe_std = pe_3.std(ddof=0)
+    for r,v  in zip(cov.row_names,cov.x):
+        d = np.abs(pe_std.loc[r] - np.sqrt(v))
+
+        print(r,v,np.sqrt(v),d)
+        assert d < 0.01
+    pst.control_data.noptmax = -1
+    pst.write(os.path.join(m_d, "pest_unc.pst"))
+    pyemu.os_utils.run("{0} {1}".format(exe_path.replace("-ies","-glm"), "pest_unc.pst"), cwd=m_d)
+    fosm_df = pd.read_csv(os.path.join(m_d,"pest_unc.par.usum.csv"),index_col=0)
+    cov_df = cov.to_dataframe()
+    for pname,prior_std in zip(fosm_df.index,fosm_df.prior_stdev):
+        d = np.abs(prior_std - np.sqrt(cov_df.loc[pname,pname]))
+        print(pname,d)
+        assert d < 1.0e-4
+
+
 def parchglim_test():
     model_d = "ies_10par_xsec"
 
@@ -1191,13 +1216,13 @@ def ins_missing_e_test():
     #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-glm.exe"),os.path.join("..","bin","win","pestpp-glm.exe"))
     #shutil.copy2(os.path.join("..", "exe", "windows", "x64", "Debug", "pestpp-ies.exe"),
     #             os.path.join("..", "bin", "win", "pestpp-ies.exe"))
-    ins_missing_e_test()
+    #ins_missing_e_test()
     #basic_test()
     #agnostic_path_test()
     #glm_long_name_test()
     #sen_plusplus_test()
     #parchglim_test()
-    #unc_file_test()
+    unc_file_test()
     #secondary_marker_test()
     #basic_test("ies_10par_xsec")
     #glm_save_binary_test()
@@ -1223,7 +1248,7 @@ def ins_missing_e_test():
     #da_mf6_freyberg_test_2()
     #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-ies.exe"),os.path.join("..","bin","win","pestpp-ies.exe"))
     #tplins1_test()
-    mf6_v5_ies_test()
+    #mf6_v5_ies_test()
     #mf6_v5_sen_test()
 
     #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-opt.exe"),os.path.join("..","bin","win","pestpp-opt.exe"))

documentation/tocgen.py

@@ -1,3 +1,4 @@
+import sys
 import os
 
 def processFile(inFile, outFile):
@@ -156,7 +157,7 @@ def clean_4_toc(docx_file,inFile,outFile,run_pandoc=True):
 
 if __name__ == "__main__":
 
-    clean_4_toc("pestpp_users_guide_v5.1.6.docx","file.md","temp.md",True)
+    clean_4_toc(sys.argv[1],"file.md","temp.md",True)
     processFile("temp.md","pestpp_users_manual.md")
 
 

src/libs/common/config_os.h

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

src/libs/pestpp_common/Pest.cpp

@@ -38,6 +38,7 @@
 #include "QSqrtMatrix.h"
 #include <limits>
 #include "network_package.h"
+#include <cmath>
 
 
 using namespace::std;
@@ -1329,11 +1330,16 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec)
 
 				}
 
-				string tcol;
+				string tcol,pcol;
 				if (cnames.find("PARTRANS") != cnames.end())
 					tcol = "PARTRANS";
 				else
 					tcol = "TRANSFORM";
+                if (cnames.find("PARNME") != cnames.end())
+                    pcol = "PARNME";
+                else
+                    pcol = "NAME";
+
 				vector<string> partrans = efile.get_col_string_vector(tcol);
 				set<string> s_partrans(partrans.begin(), partrans.end());
 				if (s_partrans.find("TIED") != s_partrans.end())
@@ -1353,8 +1359,8 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec)
 					//save any tied pars for processing later bc the par its tied to
 					//might not have been processed yet.
 					row_map = efile.get_row_map(ro);
-					if (row_map["PARTRANS"] == "TIED")
-						temp_tied_map[row_map["PARNME"]] = row_map["PARTIED"];
+					if (row_map.at(tcol) == "TIED")
+						temp_tied_map[row_map.at(pcol)] = row_map.at("PARTIED");
 				}
 				efile.keep_cols(efile_keep_cols);
 				if (efiles_map.find(section) == efiles_map.end())
@@ -1689,11 +1695,17 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec)
 
 	// handle any tied pars found in external files
 	double numer, demon, ratio;
+	vector<string> missing;
 	for (auto p: temp_tied_map)
 	{
 		name = p.first;
 		string name_tied = p.second;
 		numer = ctl_parameters[name];
+//		if (ctl_parameters.find(name_tied) == ctl_parameters.end())
+//        {
+//		    missing.push_back(name_tied);
+//		    continue;
+//        }
 		demon = ctl_parameters[name_tied];
 		if (demon == 0.0)
 		{
@@ -1713,6 +1725,16 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec)
 		tied_names.insert(name_tied);
 	}
 
+	if (missing.size() > 0)
+    {
+	    ss.str("");
+	    ss << "Error: the following `partied` parameters were not found in the control file:";
+	    for (auto& m: missing)
+	        ss << m << ",";
+	    f_rec << ss.str() << endl;
+	    throw runtime_error(ss.str());
+    }
+
 	//process pestpp options
 	map<string, PestppOptions::ARG_STATUS> arg_map, line_arg_map;
 	vector<string> dup;
@@ -3155,7 +3177,7 @@ void Pest::tokens_to_par_rec(ofstream &f_rec, const vector<string>& tokens, Tran
         {
 	        float f;
 	        convert_ip(tokens[9],f);
-	        pi.dercom = std::round(f);
+	        pi.dercom = (int)std::floor(f + 0.5);
         }
 	else
 		pi.dercom = 1;

src/libs/pestpp_common/SVDSolver.cpp

@@ -238,9 +238,19 @@ ModelRun SVDSolver::solve(RunManagerAbstract &run_manager, TerminationController
 		ModelRun prev_run(best_upgrade_run);
 		bool upgrade_start = (restart_controller.get_restart_option() == RestartController::RestartOption::RESUME_UPGRADE_RUNS);
 		best_upgrade_run = iteration_upgrd(run_manager, termination_ctl, prev_run, upgrade_start);
+        if ((global_iter_num == 2) && (pest_scenario.get_pestpp_options().get_glm_debug_high_2nd_iter_phi()))
+        {
+            Observations fake_obs(optimum_run.get_obs());
+            for (auto& oname : fake_obs.get_keys())
+            {
+                fake_obs.update_rec(oname,1.0e+10);
+            }
+            best_upgrade_run.set_observations(fake_obs);
+        }
 		// reload best parameters and set flag to switch to central derivatives next iteration
 		double prev_phi = prev_run.get_phi(*regul_scheme_ptr);
 		double best_new_phi = best_upgrade_run.get_phi(*regul_scheme_ptr);
+
 		double phi_ratio = best_new_phi / prev_phi;
 
 		cout << endl << "  ...Lambda testing complete for iteration " << termination_ctl.get_iteration_number() + 1 << endl;
@@ -288,11 +298,6 @@ ModelRun SVDSolver::solve(RunManagerAbstract &run_manager, TerminationController
 			best_upgrade_run.get_obs(), *(best_upgrade_run.get_obj_func_ptr()),
 			best_upgrade_run.get_ctl_pars());
 		file_manager.close_file(filename.str());
-		// par file for this iteration
-		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("par");
-
 		filename.str(""); // reset the stringstream
 		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()),
@@ -328,6 +333,10 @@ ModelRun SVDSolver::solve(RunManagerAbstract &run_manager, TerminationController
 			// jacobian calculated next iteration will be at the current parameters and
 			// will be more accurate than the one caluculated at the begining of this iteration
 			save_nextjac = true;
+            // par file for this iteration
+            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("par");
 		}
 		os << endl;
 		iteration_update_and_report(os, prev_run, best_upgrade_run, termination_ctl, run_manager);
@@ -1140,21 +1149,24 @@ ModelRun SVDSolver::iteration_upgrd(RunManagerAbstract &run_manager, Termination
 		RestartController::write_upgrade_runs_built(fout_restart);
 	}
 	//instance of a Mat for the jco
+    Mat j;
+    LinearAnalysis la(j, pest_scenario, file_manager, *performance_log, parcov, rand_gen_ptr);
+
 	pair<ParameterEnsemble, map<int, int>> fosm_real_info;
-	Mat j(jacobian.get_sim_obs_names(), jacobian.get_base_numeric_par_names(),
-		jacobian.get_matrix_ptr());
-	if (pest_scenario.get_prior_info_ptr()->get_nnz_pi() > 0)
-	{
-		vector<string> pi_names = pest_scenario.get_ctl_ordered_pi_names();
-		j.drop_rows(pi_names);
-	}
-	LinearAnalysis la(j, pest_scenario, file_manager, *performance_log, parcov, rand_gen_ptr);
+
 	if (pest_scenario.get_pestpp_options().get_uncert_flag())
 	{
 		cout << "-->starting iteration FOSM process..." << endl;
 
 		performance_log->log_event("LinearAnalysis::glm_iter_fosm");
-
+        Mat j(jacobian.get_sim_obs_names(), jacobian.get_base_numeric_par_names(),
+              jacobian.get_matrix_ptr());
+        if (pest_scenario.get_prior_info_ptr()->get_nnz_pi() > 0)
+        {
+            vector<string> pi_names = pest_scenario.get_ctl_ordered_pi_names();
+            j.drop_rows(pi_names);
+        }
+        LinearAnalysis la(j, pest_scenario, file_manager, *performance_log, parcov, rand_gen_ptr);
 		try
 		{
 			la.glm_iter_fosm(base_run, output_file_writer, termination_ctl.get_iteration_number(), &run_manager);

src/libs/pestpp_common/covariance.cpp

@@ -1334,7 +1334,7 @@ void Covariance::from_uncertainty_file(const string &filename, vector<string> &o
 					if (find(names.begin(), names.end(), name) != names.end())
 						throw runtime_error(name + " listed more than once in uncertainty file:" + filename);
 					names.push_back(tokens[0]);
-					triplet_list.push_back(Eigen::Triplet<double>(irow, jcol, val * std_mlt));
+					triplet_list.push_back(Eigen::Triplet<double>(irow, jcol, (val * std_mlt)*(val * std_mlt)));
 					irow++, jcol++;
 				}
 

src/libs/pestpp_common/linear_analysis.cpp

@@ -1140,7 +1140,8 @@ void LinearAnalysis::write_par_credible_range(ofstream &fout, string sum_filenam
 	pair<double, double> range;
 	for (auto &pname : ordered_names)
 	{
-		if (find(jacobian.cn_ptr()->begin(), jacobian.cn_ptr()->end(), pname) == jacobian.cn_ptr()->end())
+		//if (find(jacobian.cn_ptr()->begin(), jacobian.cn_ptr()->end(), pname) == jacobian.cn_ptr()->end())
+		if (prior_vars.find(pname) == prior_vars.end())
 			missing.push_back(pname);
 		else
 		{
@@ -1149,7 +1150,7 @@ void LinearAnalysis::write_par_credible_range(ofstream &fout, string sum_filenam
 			//if (parinfo.get_parameter_rec_ptr(pname)->tranform_type == ParameterRec::TRAN_TYPE::LOG)
 			//	value = log10(value);
 			//range = get_range(value, prior_vars[pname], parinfo.get_parameter_rec_ptr(pname)->tranform_type);
-			stdev = sqrt(prior_vars[pname]);
+			stdev = sqrt(prior_vars.at(pname));
 
 			fout << setw(20) << pest_utils::lower_cp(pname) << setw(20) << value << setw(20) << stdev << setw(20) <<
 				value - (2.0*stdev) << setw(20) << value + (2.0*stdev);
@@ -1158,7 +1159,7 @@ void LinearAnalysis::write_par_credible_range(ofstream &fout, string sum_filenam
 
 			//posterior
 			value = opt_pars.get_rec(pname);
-			stdev = sqrt(post_vars[pname]);
+			stdev = sqrt(post_vars.at(pname));
 			//range = get_range(value, post_vars[pname], parinfo.get_parameter_rec_ptr(pname)->tranform_type);
 
 			fout << setw(20) << value << setw(20) << stdev << setw(20) <<

src/libs/pestpp_common/pest_data_structs.cpp

@@ -477,7 +477,7 @@ PestppOptions::ARG_STATUS PestppOptions::assign_value_by_key(string key, const s
 		passed_args.insert("BASE_JACOBIAN");
 		passed_args.insert("BASE_JACOBIAN_FILENAME");
 
-		//convert_ip(org_value, basejac_filename);
+		//convert_ip(org_value, basejac_filename);f
 		basejac_filename = org_value;
 	}
 
@@ -579,6 +579,10 @@ PestppOptions::ARG_STATUS PestppOptions::assign_value_by_key(string key, const s
 	{
 		glm_debug_real_fail = pest_utils::parse_string_arg_to_bool(value);
 	}
+    else if (key == "GLM_DEBUG_HIGH_2ND_ITER_PHI")
+    {
+        glm_debug_high_2nd_iter_phi = pest_utils::parse_string_arg_to_bool(value);
+    }
 	else if (key == "UPGRADE_AUGMENT")
 	{
 		cout << "++UPGRADE_AUGMENT is deprecated and no longer supported...ignoring" << endl;
@@ -1565,6 +1569,7 @@ void PestppOptions::summary(ostream& os) const
 	os << "glm_accept_mc_phi: " << glm_accept_mc_phi << endl;
 	os << "glm_rebase_super: " << glm_rebase_super << endl;
 	os << "glm_iter_mc: " << glm_iter_mc << endl;
+	os << "glm_high_2nd_iter_phi: " << glm_debug_high_2nd_iter_phi << endl;
 
 //	if (global_opt == OPT_DE)
 //	{
@@ -1774,6 +1779,7 @@ void PestppOptions::set_defaults()
 	set_glm_accept_mc_phi(false);
 	set_glm_rebase_super(false);
 	set_glm_iter_mc(false);
+    set_glm_debug_high_2nd_iter_phi(false);
 	set_prediction_names(vector<string>());
 	set_parcov_filename(string());
 	set_obscov_filename(string());

src/libs/pestpp_common/pest_data_structs.h

@@ -295,6 +295,9 @@ class PestppOptions {
 	void set_glm_rebase_super(bool _flag) { glm_rebase_super = _flag; }
     bool get_glm_iter_mc() const { return glm_iter_mc; }
     void set_glm_iter_mc(bool _flag) { glm_iter_mc = _flag; }
+    bool get_glm_debug_high_2nd_iter_phi() const {return glm_debug_high_2nd_iter_phi;}
+    void set_glm_debug_high_2nd_iter_phi(bool _flag) {glm_debug_high_2nd_iter_phi = _flag;}
+
 
 
 
@@ -638,6 +641,7 @@ class PestppOptions {
 	bool glm_accept_mc_phi;
 	bool glm_rebase_super;
 	bool glm_iter_mc;
+	bool glm_debug_high_2nd_iter_phi;
 
 	vector<double> base_lambda_vec;
 	vector<double> lambda_scale_vec;

src/programs/pestpp-ies/pestpp-ies.cpp

@@ -35,9 +35,9 @@ using namespace pest_utils;
 
 int main(int argc, char* argv[])
 {
-#ifndef _DEBUG
-	try {
-#endif
+//#ifndef _DEBUG
+//	try {
+//#endif
         string version = PESTPP_VERSION;
         cout << endl << endl;
         cout << "             pestpp-ies: a GLM iterative ensemble smoother" << endl << endl;
@@ -301,20 +301,20 @@ int main(int argc, char* argv[])
         fout_rec.close();
         return 0;
 
-#ifndef _DEBUG
-	}
-	catch (exception &e)
-	{
-		cout << "Error condition prevents further execution: " << endl << e.what() << endl;
-		//cout << "press enter to continue" << endl;
-		//char buf[256];
-		//OperSys::gets_s(buf, sizeof(buf));
-		return 1;
-	}
-	catch (...)
-	{
-		cout << "Error condition prevents further execution" << endl;
-		return 1;
-	}
-#endif
+//#ifndef _DEBUG
+//	}
+//	catch (exception &e)
+//	{
+//		cout << "Error condition prevents further execution: " << endl << e.what() << endl;
+//		//cout << "press enter to continue" << endl;
+//		//char buf[256];
+//		//OperSys::gets_s(buf, sizeof(buf));
+//		return 1;
+//	}
+//	catch (...)
+//	{
+//		cout << "Error condition prevents further execution" << endl;
+//		return 1;
+//	}
+//#endif
 }