Merge pull request usgs#252 from jtwhite79/feat_emsolutionreporting

Feat emsolutionreporting
jtwhite79 · May 22, 2023 · a8a829d · a8a829d
2 parents 934413c + f29d2ba
commit a8a829d
Show file tree

Hide file tree

Showing 10 changed files with 1,438 additions and 1,378 deletions.
diff --git a/benchmarks/basic_tests.py b/benchmarks/basic_tests.py
@@ -1351,7 +1351,7 @@ def build_and_draw_prior(t_d="ends",num_reals=500):
     #glm_long_name_test()
     #sen_plusplus_test()
     #parchglim_test()
-    #unc_file_test()
+    unc_file_test()
     #cmdline_test()
     #secondary_marker_test()
     #basic_test("ies_10par_xsec")
@@ -1382,7 +1382,7 @@ def build_and_draw_prior(t_d="ends",num_reals=500):
     #mf6_v5_sen_test()
 
     #shutil.copy2(os.path.join("..","exe","windows","x64","Debug","pestpp-opt.exe"),os.path.join("..","bin","win","pestpp-opt.exe"))
-    mf6_v5_opt_stack_test()
+    #mf6_v5_opt_stack_test()
     #mf6_v5_glm_test()
     #mf6_v5_ies_test()
     #cmdline_test()

diff --git a/documentation/pestpp_users_guide_v5.2.3.docx → documentation/pestpp_users_guide_v5.2.4.docx b/documentation/pestpp_users_guide_v5.2.3.docx → documentation/pestpp_users_guide_v5.2.4.docx
diff --git a/documentation/pestpp_users_manual.md b/documentation/pestpp_users_manual.md
diff --git a/src/libs/common/config_os.h b/src/libs/common/config_os.h
@@ -2,7 +2,7 @@
 #define CONFIG_OS_H_
 
 
-#define PESTPP_VERSION "5.2.3";
+#define PESTPP_VERSION "5.2.4";
 
 #if defined(_WIN32) || defined(_WIN64)
 #define OS_WIN

diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.cpp b/src/libs/pestpp_common/EnsembleMethodUtils.cpp
diff --git a/src/libs/pestpp_common/EnsembleMethodUtils.h b/src/libs/pestpp_common/EnsembleMethodUtils.h
@@ -111,6 +111,8 @@ class L2PhiHandler
     map<string,double> get_actual_swr_map(ObservationEnsemble& oe, string real_name="");
 	map<string,map<string,double>> get_meas_phi_weight_ensemble(ObservationEnsemble& oe, ObservationEnsemble& weights);
 
+    vector<string> get_violating_realizations(ObservationEnsemble& oe, const vector<string>& viol_obs_names);
+
 private:
 	string tag;
 	map<string, double> get_summary_stats(phiType pt);
@@ -287,6 +289,14 @@ class UpgradeThread
 
 	virtual void work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector<string> par_names, vector<string> obs_names) { ; }
 
+
+    static void ensemble_solution(const int iter, const int verbose_level,const int maxsing,  const int thread_id,
+                           const int t_count, const bool use_prior_scaling,const bool use_approx, const bool use_glm,
+                           const double cur_lam,const double eigthresh, Eigen::MatrixXd& par_resid, Eigen::MatrixXd& par_diff,
+                           const Eigen::MatrixXd& Am, Eigen::MatrixXd& obs_resid,Eigen::MatrixXd& obs_diff, Eigen::MatrixXd& upgrade_1,
+                           Eigen::MatrixXd& obs_err, const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& weights,
+                           const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& parcov_inv,
+                           const vector<string>& act_obs_names,const vector<string>& act_par_names);
 protected:
 	PerformanceLog* performance_log;
 	Localizer::How how;
@@ -308,55 +318,35 @@ class UpgradeThread
 	mutex par_diff_lock, am_lock, put_lock, obs_err_lock;
 	mutex next_lock;
 
-};
-
-void ensemble_solution(const int iter, const int verbose_level,const int maxsing,  const int thread_id,
-                       const int t_count, const bool
-                  use_prior_scaling,const bool use_approx, const bool use_glm, const double cur_lam,
-                  const double eigthresh, Eigen::MatrixXd& par_resid, Eigen::MatrixXd& par_diff,
-                  const Eigen::MatrixXd& Am, Eigen::MatrixXd& obs_resid,Eigen::MatrixXd& obs_diff, Eigen::MatrixXd& upgrade_1,
-                  Eigen::MatrixXd& obs_err,
-                  const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& weights,
-                  const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& parcov_inv);
-
 
-class CovLocalizationUpgradeThread : public UpgradeThread
-{
-public:
-	using UpgradeThread::UpgradeThread;
 
-	void work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector<string> par_names, vector<string> obs_names);
 };
 
+//void ensemble_solution(const int iter, const int verbose_level,const int maxsing,  const int thread_id,
+//                       const int t_count, const bool
+//                  use_prior_scaling,const bool use_approx, const bool use_glm, const double cur_lam,
+//                  const double eigthresh, Eigen::MatrixXd& par_resid, Eigen::MatrixXd& par_diff,
+//                  const Eigen::MatrixXd& Am, Eigen::MatrixXd& obs_resid,Eigen::MatrixXd& obs_diff, Eigen::MatrixXd& upgrade_1,
+//                  Eigen::MatrixXd& obs_err,
+//                  const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& weights,
+//                  const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& parcov_inv,
+//                  const vector<string>& act_obs_names,const vector<string>& act_par_names);
+
+
+//class CovLocalizationUpgradeThread : public UpgradeThread
+//{
+//public:
+//	using UpgradeThread::UpgradeThread;
+//
+//	void work(int thread_id, int iter, double cur_lam, bool use_glm_form, vector<string> par_names, vector<string> obs_names);
+//};
+
 class LocalAnalysisUpgradeThread: public UpgradeThread
 {
 public:
 	using UpgradeThread::UpgradeThread;
 
 	void work(int thread_id, int iter, double cur_lam,bool use_glm_form, vector<string> par_names, vector<string> obs_names);
-
-
-//private:
-//	PerformanceLog* performance_log;
-//	Localizer::How how;
-//	vector<string> keys;
-//	int count, total;
-//
-//	unordered_map<string, pair<vector<string>, vector<string>>>& cases;
-//
-//	ParameterEnsemble& pe_upgrade;
-//	Localizer& localizer;
-//	unordered_map<string, double>& parcov_inv_map;
-//	unordered_map<string, double>& weight_map;
-//
-//	unordered_map<string, Eigen::VectorXd>& par_resid_map, & par_diff_map, & Am_map;
-//	unordered_map<string, Eigen::VectorXd>& obs_resid_map, & obs_diff_map, obs_err_map;
-//
-//	mutex ctrl_lock, weight_lock, loc_lock, parcov_lock;
-//	mutex obs_resid_lock, obs_diff_lock, par_resid_lock;
-//	mutex par_diff_lock, am_lock, put_lock, obs_err_lock;
-//	mutex next_lock;
-
 };
 
 class EnsembleMethod
@@ -456,11 +446,13 @@ class EnsembleMethod
 	vector<double> lam_mults;
 	vector<string> oe_org_real_names, pe_org_real_names;
 	vector<string> act_obs_names, act_par_names;
+	vector<string> violation_obs;
 	ParameterEnsemble pe, pe_base;
 	ObservationEnsemble oe, oe_base, weights;
 	Eigen::DiagonalMatrix<double, Eigen::Dynamic> obscov_inv_sqrt, parcov_inv_sqrt;
 	bool oe_drawn, pe_drawn;
 
+
 	bool solve_glm(int cycle = NetPackage::NULL_DA_CYCLE);
 
 	bool solve_mda(bool last_iter, int cycle = NetPackage::NULL_DA_CYCLE);
@@ -477,7 +469,7 @@ class EnsembleMethod
 
 	void initialize_restart();
 
-	void drop_bad_phi(ParameterEnsemble& _pe, ObservationEnsemble& _oe, vector<int> subset_idxs = vector<int>());
+	void drop_bad_reals(ParameterEnsemble& _pe, ObservationEnsemble& _oe, vector<int> subset_idxs = vector<int>());
 
 	void add_bases();
 
@@ -505,5 +497,7 @@ class EnsembleMethod
                                     map<string,map<string,double>>& phi_fracs_by_real,
                                     vector<string>& index, bool check_reals);
 
+    void prep_drop_violations();
+
 };
 #endif
diff --git a/src/libs/pestpp_common/Localizer.cpp b/src/libs/pestpp_common/Localizer.cpp
@@ -29,7 +29,12 @@ bool Localizer::initialize(PerformanceLog *performance_log, ofstream& frec, bool
 	loc_typ = pest_scenario_ptr->get_pestpp_options().get_ies_loc_type();
 	if (loc_typ[0] == 'C')
 	{
-		loctyp = LocTyp::COVARIANCE;
+	    ss.str("");
+	    ss << "WARNING: covariance type localization deprecated, resetting to local-analysis";
+        frec << ss.str() << endl;
+        cout << ss.str() << endl;
+        performance_log->log_event(ss.str());
+		loctyp = LocTyp::LOCALANALYSIS;
 	}
 	else if (loc_typ[0] == 'L')
 	{
@@ -42,28 +47,42 @@ bool Localizer::initialize(PerformanceLog *performance_log, ofstream& frec, bool
 		throw runtime_error("Localizer error: 'ies_localization_type' must start with 'C' (covariance) or 'L' (local analysis) not " + loc_typ);
 	}
 
-	if (how_str[0] == 'P')
-	{
-		how = How::PARAMETERS;
-	}
-	else if (how_str[0] == 'O')
-	{
-		how = How::OBSERVATIONS;
-	}
+    how = How::PARAMETERS;
+    if (how_str[0] == 'P')
+    {
+
+    }
+    else if (how_str[0] == 'O')
+    {
+        ss.str("");
+        ss << "WARNING: localization 'how' by observations is deprecated, resetting to 'how' by parameters";
+        frec << ss.str() << endl;
+        cout << ss.str() << endl;
+        performance_log->log_event(ss.str());
+        how = How::PARAMETERS;
+    }
+//	if (how_str[0] == 'P')
+//	{
+//		how = How::PARAMETERS;
+//	}
+//	else if (how_str[0] == 'O')
+//	{
+//		how = How::OBSERVATIONS;
+//	}
 	else
 	{
 	    performance_log->log_event("Localizer error: 'ies_localize_how' or 'da_localize_how' must start with 'P' (pars) or 'O' (obs) not " + how_str);
 	    cout << "Localizer error: 'ies_localize_how' or 'da_localize_how' must start with 'P' (pars) or 'O' (obs) not \" + how_str" << endl;
 	    throw runtime_error("Localizer error: 'ies_localize_how' or 'da_localize_how' must start with 'P' (pars) or 'O' (obs) not " + how_str);
 	}
 
-	if ((loctyp == LocTyp::COVARIANCE) && (how == How::OBSERVATIONS)) {
-        performance_log->log_event("Localizer error: covariance localization can only be used with localization by parameters");
-        cout << "Localizer error: covariance localization can only be used with localization by parameters" << endl;
-
-        throw runtime_error(
-                "Localizer error: covariance localization can only be used with localization by parameters");
-    }
+//	if ((loctyp == LocTyp::COVARIANCE) && (how == How::OBSERVATIONS)) {
+//        performance_log->log_event("Localizer error: covariance localization can only be used with localization by parameters");
+//        cout << "Localizer error: covariance localization can only be used with localization by parameters" << endl;
+//
+//        throw runtime_error(
+//                "Localizer error: covariance localization can only be used with localization by parameters");
+//    }
 	string filename = pest_scenario_ptr->get_pestpp_options().get_ies_localizer();
 	autoadaloc = pest_scenario_ptr->get_pestpp_options().get_ies_autoadaloc();
 	sigma_dist = pest_scenario_ptr->get_pestpp_options().get_ies_autoadaloc_sigma_dist();

diff --git a/src/libs/pestpp_common/Pest.h b/src/libs/pestpp_common/Pest.h
@@ -121,7 +121,7 @@ class Pest {
 
 protected:
 	//this is the list of external file cols that have meaning...
-	set<string> efile_keep_cols{ "standard_deviation", "obsnme","parnme","name", "upper_bound","lower_bound", "cycle", "state_par_link" };
+	set<string> efile_keep_cols{ "standard_deviation", "obsnme","parnme","name", "upper_bound","lower_bound", "cycle", "state_par_link","drop_violations" };
 	int n_adj_par = 0;
 	string prior_info_string;
 	ControlInfo control_info;

diff --git a/src/libs/pestpp_common/SVDPackage.h b/src/libs/pestpp_common/SVDPackage.h
@@ -27,7 +27,7 @@
 class SVDPackage
 {
 public:
-	SVDPackage(std::string _descritpion="undefined", int _n_max_sing=1000, double _eign_thres=1.0e-7);
+	SVDPackage(std::string _descritpion="undefined", int _n_max_sing=1000000, double _eign_thres=1.0e-7);
 	SVDPackage(const SVDPackage &rhs) : description(rhs.description), n_max_sing(rhs.n_max_sing), eign_thres(rhs.eign_thres), performance_log(rhs.performance_log) {}
 	virtual void solve_ip(Eigen::SparseMatrix<double>& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix<double> & U, Eigen::SparseMatrix<double>& VT, Eigen::VectorXd &Sigma_trunc) = 0;
 	virtual void solve_ip(Eigen::SparseMatrix<double>& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix<double> & U, Eigen::SparseMatrix<double>& VT, Eigen::VectorXd &Sigma_trunc, double _eigen_thres) = 0;
@@ -48,7 +48,7 @@ class SVDPackage
 class SVD_EIGEN : public SVDPackage
 {
 public:
-	SVD_EIGEN(int _n_max_sing = 1000, double _eign_thres = 1.0e-7) : SVDPackage("Eigen JacobiSVD", _n_max_sing, _eign_thres)  {}
+	SVD_EIGEN(int _n_max_sing = 1000000, double _eign_thres = 1.0e-7) : SVDPackage("Eigen JacobiSVD", _n_max_sing, _eign_thres)  {}
 	SVD_EIGEN(const SVD_EIGEN &rhs) : SVDPackage(rhs) {}
 	virtual void solve_ip(Eigen::SparseMatrix<double>& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix<double>& U, Eigen::SparseMatrix<double>& VT, Eigen::VectorXd &Sigma_trunc);
 	virtual void solve_ip(Eigen::SparseMatrix<double>& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix<double>& U, Eigen::SparseMatrix<double>& VT, Eigen::VectorXd &Sigma_trunc, double _eigen_thres);
@@ -62,7 +62,7 @@ class SVD_EIGEN : public SVDPackage
 class SVD_REDSVD : public SVDPackage
 {
 public:
-	SVD_REDSVD(int _n_max_sing = 1000, double _eign_thres = 1.0e-7) : SVDPackage("RedSVD", _n_max_sing, _eign_thres) {}
+	SVD_REDSVD(int _n_max_sing = 1000000, double _eign_thres = 1.0e-7) : SVDPackage("RedSVD", _n_max_sing, _eign_thres) {}
 	SVD_REDSVD(const SVD_REDSVD &rhs) : SVDPackage(rhs) {}
 	virtual void solve_ip(Eigen::SparseMatrix<double>& A, Eigen::VectorXd &Sigma, Eigen::SparseMatrix<double>& U,
 		Eigen::SparseMatrix<double>& VT, Eigen::VectorXd &Sigma_trunc);

diff --git a/src/programs/pestpp-ies/pestpp-ies.cpp b/src/programs/pestpp-ies/pestpp-ies.cpp
@@ -133,7 +133,7 @@ int main(int argc, char* argv[])
 
         if (!restart_flag || save_restart_rec_header) {
             fout_rec << "             pestpp-ies - a GLM iterative Ensemble Smoother" << endl
-                     << "for PEST(++) datasets " << endl << endl;
+                     << "                      for PEST(++) datasets " << endl << endl;
             fout_rec << "                 by the PEST++ developement team" << endl << endl << endl;
             fout_rec << endl;
             fout_rec << endl << endl << "version: " << version << endl;