From ad8c8b69d5bd07aa26a137beaa53f4261dba6853 Mon Sep 17 00:00:00 2001
From: jwhite <jtwhite1000@gmail.com>
Date: Fri, 5 Apr 2024 10:11:06 -0600
Subject: [PATCH] temp attempt to speed up IDW chance mappingh

---
 src/libs/pestpp_common/Jacobian_1to1.cpp |  3 +-
 src/libs/pestpp_common/MOEA.cpp          |  3 +-
 src/libs/pestpp_common/constraints.cpp   | 58 ++++++++++++++++--------
 3 files changed, 44 insertions(+), 20 deletions(-)

diff --git a/src/libs/pestpp_common/Jacobian_1to1.cpp b/src/libs/pestpp_common/Jacobian_1to1.cpp
index a488d57a..0b7a2d88 100644
--- a/src/libs/pestpp_common/Jacobian_1to1.cpp
+++ b/src/libs/pestpp_common/Jacobian_1to1.cpp
@@ -394,7 +394,8 @@ bool Jacobian_1to1::forward_diff(const string &par_name, double base_derivative_
 
 	// perturb derivative parameters
 	double incr = derivative_inc(par_name, group_info, base_derivative_val, false);
-	if (incr == 0.0) return false;
+	if (incr == 0.0)
+        return false;
 	new_par_val = new_par[par_name] = base_derivative_val + incr;
 	// try forward derivative
 	out_of_bound_forward = out_of_bounds(new_par, par_info_ptr);
diff --git a/src/libs/pestpp_common/MOEA.cpp b/src/libs/pestpp_common/MOEA.cpp
index e29c5558..9f825f16 100644
--- a/src/libs/pestpp_common/MOEA.cpp
+++ b/src/libs/pestpp_common/MOEA.cpp
@@ -2294,7 +2294,7 @@ void MOEA::initialize()
 		{
 			//this can be done, but we need to make sure the appropriate chance restart
 			//args were supplied: base_jacobian or obs_stack
-			throw_moea_error("chance constraints not yet supported with restart");
+			throw_moea_error("chance constraints/objectives not yet supported with restart");
 		}
 	
 		//since mou reqs strict linking of realization names, let's see if we can find an intersection set 
@@ -3078,6 +3078,7 @@ void MOEA::initialize_population_schedule()
         }
         in.close();
     }
+
     ofstream& frec = file_manager.rec_ofstream();
     frec << "...population schedule: generation,population size:" << endl;
     for (int i=0;i<pest_scenario.get_control_info().noptmax;i++)
diff --git a/src/libs/pestpp_common/constraints.cpp b/src/libs/pestpp_common/constraints.cpp
index e13fd580..9f5d854c 100644
--- a/src/libs/pestpp_common/constraints.cpp
+++ b/src/libs/pestpp_common/constraints.cpp
@@ -943,15 +943,15 @@ void Constraints::initial_report()
 		}
 		cout << "...opt_risk = " << risk;
 		if (std_weights)
-			cout << ", using FOSM-based chance constraints with the weights of non-zero-weighted constraints as standard deviation" << endl;
+			cout << ", using FOSM-based chance constraints/objectives with the weights of non-zero-weighted constraints as standard deviation" << endl;
 		else if (use_fosm)
-			cout << ", using FOSM-based chance constraints with " << adj_par_names.size() << " adjustable parameters" << endl;
+			cout << ", using FOSM-based chance constraints/objectives with " << adj_par_names.size() << " adjustable parameters" << endl;
 		else
-			cout << ", using stack-based chance constraints with " << stack_pe.shape().first << " realizations" << endl;
+			cout << ", using stack-based chance constraints/objectives with " << stack_pe.shape().first << " realizations" << endl;
 	}
 	else
 	{
-		cout << "...not using chance constraints" << endl;
+		cout << "...not using chance constraints/objectives" << endl;
 	}
 }
 
@@ -1209,6 +1209,16 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl
 			vector<double> factors,temp;
 			vector<string> dreal_names;
 
+            Eigen::MatrixXd stack_pe_mat(stack_pe_map.size(),dvnames.size());
+            vector<string> stack_pe_map_rnames;
+            int i = 0;
+            for (auto& p : stack_pe_map) {
+                stack_dv_vec = p.second.get_data_eigen_vec(dvnames);
+                stack_pe_mat.row(i) = stack_dv_vec;
+                stack_pe_map_rnames.push_back(p.first);
+                i++;
+            }
+
             Eigen::MatrixXd shifts;
             double factor_sum, factor, factor_sum2;
 			for (int i=0;i<missing.size();i++)
@@ -1220,24 +1230,36 @@ ObservationEnsemble Constraints::get_chance_shifted_constraints(ParameterEnsembl
 				distances.clear();
                 double _risk = risk;
 
-                if (risk_obj.size() > 0)
-                {
+                if (risk_obj.size() > 0) {
                     //_risk = stack_pe_map[min_real_name][risk_obj];
                     _risk = risk_map[missing_real_name];
                 }
 
-				for (auto& p : stack_pe_map)
-				{
-					stack_dv_vec = p.second.get_data_eigen_vec(dvnames);
-					dist = (stack_dv_vec - missing_dv_vec).squaredNorm();
-					distances.push_back(dist);
-					dreal_names.push_back(p.first);
-					if (dist < min_dist)
-					{
-						min_dist = dist;
-						min_real_name = p.first;
-					}
-				}
+//				for (auto& p : stack_pe_map)
+//				{
+//					stack_dv_vec = p.second.get_data_eigen_vec(dvnames);
+//					dist = (stack_dv_vec - missing_dv_vec).squaredNorm();
+//					distances.push_back(dist);
+//					dreal_names.push_back(p.first);
+//					if (dist < min_dist)
+//					{
+//						min_dist = dist;
+//						min_real_name = p.first;
+//					}
+//				}
+                for (i=0;i<stack_pe_map_rnames.size();i++)
+                {
+                    dist = (stack_pe_mat.row(i) - missing_dv_vec).squaredNorm();
+                    distances.push_back(dist);
+                    if (dist < min_dist)
+                    {
+                        min_dist = dist;
+                        min_real_name = stack_pe_map_rnames[i];
+                    }
+                }
+                dreal_names = stack_pe_map_rnames;
+
+
                 temp.clear();
 				max_dist = *max_element(distances.begin(),distances.end());
                 for (auto& d : distances)