Merge pull request #137 from jtwhite79/hotfix_opt

Hotfix opt

src/libs/common/config_os.h

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

src/libs/pestpp_common/Jacobian.cpp

@@ -170,13 +170,51 @@ unordered_map<string, int> Jacobian::get_obs2row_map() const
 }
 
 
-Eigen::SparseMatrix<double> Jacobian::get_matrix(const vector<string> &obs_names, const vector<string> & par_names) const
+Eigen::SparseMatrix<double> Jacobian::get_matrix(const vector<string> &obs_names, const vector<string> & par_names, bool forgive_missing, int n_cols) const
 {
+    /* the n_cols arg is so you can reserve a sparse matrix with columns that are all zeros - this is for the LP solver in pestpp-opt when you have
+     * external dec vars - they wont be in the Jacobian instance but they need to be in the LP solution matrix.
+     */
+    stringstream ss;
 	int n_rows = obs_names.size();
-	int n_cols = par_names.size();
+	if (n_cols == 0)
+	    n_cols = par_names.size();
 	int irow_new;
 	int icol_new;
 
+	if (!forgive_missing) {
+
+        set<string> sobs_names(base_sim_obs_names.begin(), base_sim_obs_names.end());
+        set<string> spar_names(base_numeric_par_names.begin(), base_numeric_par_names.end());
+
+        set<string>::iterator send = sobs_names.end();
+        vector<string> missing;
+        for (auto &obs_name : obs_names)
+            if (sobs_names.find(obs_name) == send)
+                missing.push_back(obs_name);
+        if (missing.size() > 0) {
+            ss.str("");
+
+            ss << "Jco::get_matrix(): the following obs names are not in the matrix:";
+            for (auto &m : missing)
+                ss << " " << m;
+            throw runtime_error(ss.str());
+        }
+
+        send = spar_names.end();
+        for (auto &par_name : par_names)
+            if (spar_names.find(par_name) == send)
+                missing.push_back(par_name);
+        if (missing.size() > 0){
+            ss.str("");
+
+            ss << "Jco::get_matrix(): the following par names are not in the matrix:";
+            for (auto &m : missing)
+                ss << " " << m;
+            throw runtime_error(ss.str());
+        }
+    }
+
 	unordered_map<string, int> obs_name2newindex_map;
 	unordered_map<string, int> par_name2new_index_map;
 

src/libs/pestpp_common/Jacobian.h

@@ -67,7 +67,7 @@ class Jacobian {
 	virtual const vector<string>& observation_list() const {return  base_sim_obs_names;}
 	virtual const vector<string>& obs_and_reg_list() const;
 	virtual const Parameters &get_base_numeric_parameters() const{return base_numeric_parameters;};
-	Eigen::SparseMatrix<double> get_matrix(const vector<string> &obs_names, const vector<string> & par_name_vec) const;
+	Eigen::SparseMatrix<double> get_matrix(const vector<string> &obs_names, const vector<string> & par_name_vec, bool forgive_missing=false, int n_cols=0) const;
 
 	virtual bool build_runs(ModelRun &model_run, vector<string> numeric_par_names, ParamTransformSeq &par_transform,
 		const ParameterGroupInfo &group_info, const ParameterInfo &ctl_par_info,

src/libs/pestpp_common/MOEA.cpp

@@ -1081,19 +1081,21 @@ map<string, map<string, double>> MOEA::obj_func_report(ParameterEnsemble& _dp, O
 
 	int max_len = get_max_len_obj_name();
 	string dir;
-	frec << left << setw(max_len) << "objective function" << right << setw(10) << "direction" << setw(10) << "mean" << setw(20) << "standard devation" << setw(12) << "min" << setw(12) << "max" << endl;
+	frec << left << setw(max_len) << "objective function" << right << setw(12) << "direction" << setw(15) << "mean" << setw(15) << "std dev" << setw(15) << "min" << setw(15) << "max" << endl;
+    frec << setprecision(7);
 	for (auto obs_obj : obs_obj_names)
 	{
 
 		frec << left << setw(max_len) << obs_obj;
 		dir = "minimize";
 		if (obj_dir_mult[obs_obj] == -1)
 			dir = "maximize";
-		frec << right << setw(10) << dir;
-		frec << right << setw(10) << summary[obs_obj]["mean"];
-		frec << setw(20) << summary[obs_obj]["std"];
-		frec << setw(12) << summary[obs_obj]["min"];
-		frec << setw(12) << summary[obs_obj]["max"] << endl;
+
+		frec << right << setw(12) << dir;
+		frec << right << setw(15) << summary[obs_obj]["mean"];
+		frec << setw(15) << summary[obs_obj]["std"];
+		frec << setw(15) << summary[obs_obj]["min"];
+		frec << setw(15) << summary[obs_obj]["max"] << endl;
 	}
 
 
@@ -1103,11 +1105,11 @@ map<string, map<string, double>> MOEA::obj_func_report(ParameterEnsemble& _dp, O
 		dir = "minimize";
 		if (obj_dir_mult[pi_obj] == -1)
 			dir = "maximize";
-		frec << right << setw(10) << dir;
-		frec << right << setw(10) << summary[pi_obj]["mean"];
-		frec << setw(20) << summary[pi_obj]["std"];
-		frec << setw(12) << summary[pi_obj]["min"];
-		frec << setw(12) << summary[pi_obj]["max"] << endl;
+		frec << right << setw(12) << dir;
+		frec << right << setw(15) << summary[pi_obj]["mean"];
+		frec << setw(15) << summary[pi_obj]["std"];
+		frec << setw(15) << summary[pi_obj]["min"];
+		frec << setw(15) << summary[pi_obj]["max"] << endl;
 	}
 
 	frec << endl;
@@ -1912,7 +1914,7 @@ void MOEA::initialize()
 		pest_scenario.get_ctl_parameter_info_ptr_4_mod()->get_parameter_rec_ptr_4_mod(RISK_NAME)->lbnd = max(b,0.01);
 		b = pest_scenario.get_ctl_parameter_info().get_parameter_rec_ptr(RISK_NAME)->ubnd;
 		pest_scenario.get_ctl_parameter_info_ptr_4_mod()->get_parameter_rec_ptr_4_mod(RISK_NAME)->ubnd = min(b, 0.99);
-		//set this just to make sure everuything gets initialized right
+		//set this just to make sure everything gets initialized right
 		pest_scenario.get_pestpp_options_ptr()->set_opt_risk(0.95);
 
 
@@ -2266,8 +2268,14 @@ void MOEA::initialize()
 	message(0, "initial population objective function summary:");
 	previous_obj_summary = obj_func_report(dp, op);
 
+
+
+
 	if (constraints.get_use_chance())
 	{
+        string sum = constraints.mou_population_observation_constraint_summary(0,op,"pre-shift",obs_obj_names);
+        frec << sum << endl;
+        cout << sum << endl;
 	    string opt_member;
 		ObservationEnsemble shifted_op = get_chance_shifted_op(dp, op, opt_member);
 		ss.str("");
@@ -2578,6 +2586,9 @@ void MOEA::iterate_to_solution()
 		save_populations(new_dp, new_op);
 		if (constraints.get_use_chance())
 		{
+		    string csum = constraints.mou_population_observation_constraint_summary(iter,new_op,"pre-shift",obs_obj_names);
+		    cout << csum;
+		    file_manager.rec_ofstream() << csum;
 		    string opt_member;
 			pair<Parameters,Observations> po = get_optimal_solution(dp, op, opt_member);
 			constraints.presolve_chance_report(iter, po.second,true, "chance constraint summary (calculated at optimal/mean decision variable point)");

src/libs/pestpp_common/Pest.cpp

@@ -1048,6 +1048,12 @@ int Pest::process_ctl_file(ifstream& fin, string _pst_filename, ofstream& f_rec)
 					throw_control_file_error(f_rec, "'* control data keyword' cant be used with '++' args");
 				sections_found.insert("PLUSPLUS");
 				pestpp_input.push_back(line);
+				section = "PLUSPLUS";
+				if (!prior_info_string.empty())
+                {
+                    tokens_to_pi_rec(f_rec,prior_info_string);
+                    prior_info_string.clear();
+                }
 			}
 
 			else if (line_upper[0] == '*')

src/libs/pestpp_common/SQP.cpp

@@ -2601,7 +2601,12 @@ bool SeqQuadProgram::seek_feasible()
 		snames.erase(obj_obs);
 	}
 	ObservationInfo* oi = ies_pest_scenario.get_observation_info_ptr();
-	Observations shifted = constraints.get_chance_shifted_constraints(current_obs);
+
+	Observations shifted = pest_scenario.get_ctl_observations();
+	if (constraints.get_use_chance())
+    {
+	    shifted = constraints.get_chance_shifted_constraints(current_obs);
+    }
 	Observations ctl_obs = ies_pest_scenario.get_ctl_observations_4_mod();
 	for (auto& name : ies_pest_scenario.get_ctl_ordered_obs_names())
 	{