Skip to content

Commit

Permalink
Merge pull request #252 from jtwhite79/feat_emsolutionreporting
Browse files Browse the repository at this point in the history
Feat emsolutionreporting
  • Loading branch information
jtwhite79 authored May 22, 2023
2 parents 934413c + f29d2ba commit a8a829d
Show file tree
Hide file tree
Showing 10 changed files with 1,438 additions and 1,378 deletions.
4 changes: 2 additions & 2 deletions benchmarks/basic_tests.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down Expand Up @@ -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()
Expand Down
Binary file not shown.
16 changes: 11 additions & 5 deletions documentation/pestpp_users_manual.md

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion src/libs/common/config_os.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2,659 changes: 1,350 additions & 1,309 deletions src/libs/pestpp_common/EnsembleMethodUtils.cpp

Large diffs are not rendered by default.

74 changes: 34 additions & 40 deletions src/libs/pestpp_common/EnsembleMethodUtils.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand All @@ -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();

Expand Down Expand Up @@ -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
51 changes: 35 additions & 16 deletions src/libs/pestpp_common/Localizer.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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')
{
Expand All @@ -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();
Expand Down
2 changes: 1 addition & 1 deletion src/libs/pestpp_common/Pest.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
6 changes: 3 additions & 3 deletions src/libs/pestpp_common/SVDPackage.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);
Expand All @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion src/programs/pestpp-ies/pestpp-ies.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit a8a829d

Please sign in to comment.