Merge pull request PetroleumCyberneticsGroup#143 from PetroleumCybern…

…eticsGroup/feature/Penalty_function_improvements

Feature/penalty function and hybridization improvements

FieldOpt/Optimization/constraints/constraint_handler.cpp

@@ -219,8 +219,16 @@ long double ConstraintHandler::GetWeightedNormalizedPenalties(Case *c) {
     long double wnp = 0.0L;
     for (auto con : constraints_) {
         long double pen = con->PenaltyNormalized(c);
+        if (VERB_OPT >= 3 && pen > 0) {
+            Printer::ext_info("Penalty from constraint " + con->name()
+                    + ": " + Printer::num2str(pen), "Optimization", "ConstraintHandler");
+        }
         wnp += pen * con->GetPenaltyWeight();
     }
+    if (VERB_OPT >= 2) {
+        Printer::ext_info("Weighted, normalized penalty for case " + c->id().toString().toStdString()
+                + ": " + Printer::num2str(wnp), "Optimization", "ConstraintHandler");
+    }
     return wnp;
 }
 

FieldOpt/Optimization/constraints/interwell_distance.cpp

@@ -21,6 +21,8 @@
 #include "ConstraintMath/well_constraint_projections/well_constraint_projections.h"
 #include <boost/lexical_cast.hpp>
 #include <cmath>
+#include "Utilities/verbosity.h"
+#include "Utilities/printer.hpp"
 
 namespace Optimization {
 namespace Constraints {
@@ -111,9 +113,9 @@ void InterwellDistance::InitializeNormalizer(QList<Case *> cases) {
     long double minimum_distance = 1e20;
     for (auto c : cases) {
         vector<double> endp_dist = endpointDistances(c);
-        double min_dist = *min_element(endp_dist.begin(), endp_dist.end());
-        if (min_dist < minimum_distance){
-            minimum_distance = min_dist;
+        for (double dist : endp_dist) {
+            if (abs(dist) < minimum_distance)
+                minimum_distance = abs(dist);
         }
     }
     normalizer_.set_max(1.0L);
@@ -125,9 +127,13 @@ double InterwellDistance::Penalty(Case *c) {
     double violation = 0.0;
     for (auto distance : endpoint_distances) {
         if (distance < distance_) {
+            if (VERB_OPT >= 2) Printer::ext_info("Interwell distance penalty for case " + c->id().toString().toStdString()
+                    + ": " + Printer::num2str(distance - distance_) + " m too close");
             violation += abs(distance - distance_);
         }
     }
+    if (VERB_OPT >= 2) Printer::ext_info("Interwell distance total violation for case "
+            + c->id().toString().toStdString() + ": " + Printer::num2str(violation));
     return violation;
 }
 
@@ -148,7 +154,11 @@ vector<double> InterwellDistance::endpointDistances(Case *c) {
     return endpoint_distances;
 }
 long double InterwellDistance::PenaltyNormalized(Case *c) {
-    return normalizer_.normalize(Penalty(c));
+    double penalty = Penalty(c);
+    if (penalty > 0.0)
+        return normalizer_.normalize(penalty);
+    else
+        return 0.0;
 }
 
 }

FieldOpt/Optimization/constraints/reservoir_boundary.cpp

@@ -20,6 +20,8 @@
 #include "ConstraintMath/well_constraint_projections/well_constraint_projections.h"
 #include <iomanip>
 #include "Utilities/math.hpp"
+#include "Utilities/verbosity.h"
+#include "Utilities/printer.hpp"
 
 namespace Optimization {
 namespace Constraints {
@@ -455,14 +457,20 @@ void ReservoirBoundary::InitializeNormalizer(QList<Case *> cases) {
 double ReservoirBoundary::Penalty(Case *c) {
     if (CaseSatisfiesConstraint(c))
         return 0.0;
-    else
+    else {
+        if (VERB_OPT >= 2) Printer::ext_info("Reservoir boundary penalty (infty) for case " + c->id().toString().toStdString(),
+                "Optimization", "ReservoirBoundary");
         return INFINITY;
+    }
 }
 long double ReservoirBoundary::PenaltyNormalized(Case *c) {
     if (CaseSatisfiesConstraint(c))
         return 0.0L;
-    else
+    else {
+        if (VERB_OPT >= 2) Printer::ext_info("Reservoir boundary penalty (infty) for case " + c->id().toString().toStdString(),
+                "Optimization", "ReservoirBoundary");
         return 1.0L;
+    }
 }
 
 }

FieldOpt/Optimization/constraints/well_spline_length.cpp

@@ -18,6 +18,8 @@
 ******************************************************************************/
 
 #include "Utilities/math.hpp"
+#include "Utilities/verbosity.h"
+#include "Utilities/printer.hpp"
 #include "well_spline_length.h"
 #include "ConstraintMath/well_constraint_projections/well_constraint_projections.h"
 
@@ -103,10 +105,16 @@ void WellSplineLength::InitializeNormalizer(QList<Case *> cases) {
 double WellSplineLength::Penalty(Case *c) {
     auto endpts = GetEndpointValueVectors(c, affected_well_);
     double well_length =  (endpts.first - endpts.second).norm();
-    if (well_length > max_length_)
+    if (well_length > max_length_) {
+        if (VERB_OPT >= 2) Printer::ext_info("Well length penalty for case " + c->id().toString().toStdString()
+                + ": " + Printer::num2str(well_length - max_length_) + " m too long", "Optimization", "WellSplineLength");
         return well_length - max_length_;
-    else if (well_length < min_length_)
+    }
+    else if (well_length < min_length_) {
+        if (VERB_OPT >= 2) Printer::ext_info("Well length penalty for case " + c->id().toString().toStdString()
+                + ":" + Printer::num2str(min_length_ - well_length) + " m too short", "Optimization", "WellSplineLength");
         return min_length_ - well_length;
+    }
     else
         return 0.0;
 }

FieldOpt/Optimization/optimizer.cpp

@@ -47,6 +47,7 @@ Optimizer::Optimizer(Settings::Optimizer *settings, Case *base_case,
     else {
         Printer::ext_info("Using shared CaseHandler.", "Optimizer", "Optimization");
         case_handler_ = case_handler;
+        is_hybrid_component_ = true;
     }
     if (constraint_handler == 0) {
         constraint_handler_ = new Constraints::ConstraintHandler(settings->constraints(), variables, grid);
@@ -64,6 +65,25 @@ Optimizer::Optimizer(Settings::Optimizer *settings, Case *base_case,
     enable_logging_ = true;
     verbosity_level_ = 0;
     penalize_ = settings->objective().use_penalty_function;
+
+    if (penalize_) {
+        if (!normalizer_ofv_.is_ready()) {
+            if (VERB_OPT >=1) {
+                Printer::ext_info("Initializing normalizers", "Optimization", "Optimizer");
+            }
+            initializeNormalizers();
+
+            // penalize the base case
+            double org_ofv = tentative_best_case_->objective_function_value();
+            double pen_ofv = PenalizedOFV(tentative_best_case_);
+            tentative_best_case_->set_objective_function_value(pen_ofv);
+            if (VERB_OPT >=1) {
+                Printer::ext_info("Penalized base case. " 
+                        "Original value: " + Printer::num2str(org_ofv) + "; "
+                        + "Penalized value: " + Printer::num2str(pen_ofv), "Optimization", "Optimizer");
+            }
+        }
+    }
 }
 
 Case *Optimizer::GetCaseForEvaluation()
@@ -81,7 +101,7 @@ Case *Optimizer::GetCaseForEvaluation()
 void Optimizer::SubmitEvaluatedCase(Case *c)
 {
     evaluated_cases_++;
-    if (penalize_ && iteration_ > 0) {
+    if (penalize_) {
         double penalized_ofv = PenalizedOFV(c);
         c->set_objective_function_value(penalized_ofv);
     }
@@ -243,13 +263,23 @@ double Optimizer::PenalizedOFV(Case *c) {
     long double norm_ofv = normalizer_ofv_.normalize(c->objective_function_value());
     long double penalty = constraint_handler_->GetWeightedNormalizedPenalties(c);
     long double norm_pen_ovf = norm_ofv - penalty;
+    double denormalized_ofv = normalizer_ofv_.denormalize(norm_pen_ovf);
+
+    if (VERB_OPT >= 2) {
+        Printer::ext_info("Penalized case "     + c->id().toString().toStdString()                + ". "
+                "Initial OFV: "                 + Printer::num2str(c->objective_function_value()) + "; "
+                "Normalized OFV :"              + Printer::num2str(norm_ofv)                      + "; "
+                "Normalized penalty: "          + Printer::num2str(penalty)                       + "; "
+                "Penalized, normalized OFV: "   + Printer::num2str(norm_pen_ovf)                  + "; "
+                "Denormalized, penalized OFV: " + Printer::num2str(denormalized_ofv), "Optimization", "Optimizer");
+    }
 
     if (norm_pen_ovf <= 0.0L) {
         cout << "RETURNING ZERO OFV" << endl;
         return 0.0;
     }
     else {
-        return normalizer_ofv_.denormalize(norm_pen_ovf);
+        return denormalized_ofv;
     }
 }
 

FieldOpt/Optimization/optimizer.h

@@ -177,6 +177,7 @@ class Optimizer : public Loggable
   bool penalize_; //!< Switch for whether or not to use penalty function to account for constraints.
   Case *tentative_best_case_; //!< The best case encountered thus far.
   int tentative_best_case_iteration_; //!< The iteration in which the current tentative best case was found.
+  bool is_hybrid_component_; //!< Indicates that this object is a hybrid optimization component.
 
   Normalizer normalizer_ofv_; //!< Normalizer for objective function values.
 

FieldOpt/Optimization/optimizers/APPS.cpp

@@ -44,13 +44,6 @@ APPS::APPS(Settings::Optimizer *settings,
     assert(settings->parameters().max_queue_size >= 1.0);
     max_queue_length_ = directions_.size() * settings->parameters().max_queue_size;
     is_async_ = true;
-    if (penalize_) {
-        if (!normalizer_ofv_.is_ready()) initializeNormalizers();
-
-        // penalize the base case
-        double pen_ofv = PenalizedOFV(tentative_best_case_);
-        tentative_best_case_->set_objective_function_value(pen_ofv);
-    }
     if (enable_logging_) {
         logger_->AddEntry(this);
     }
@@ -60,11 +53,16 @@ void APPS::iterate() {
     if (enable_logging_) {
         logger_->AddEntry(this);
     }
-    if (inactive().size() > 0) {
+    if (inactive().size() > 0 && evaluated_cases_ < max_evaluations_) {
         case_handler_->AddNewCases(generate_trial_points(inactive()));
         set_active(inactive());
     }
     iteration_++;
+    if (is_hybrid_component_) {
+        // Increment this here if this object is a hybrid optimization component,
+        // as it will not be incremented elsewhere.
+        evaluated_cases_ = iteration_;
+    }
 }
 
 void APPS::handleEvaluatedCase(Case *c) {
@@ -89,6 +87,7 @@ void APPS::unsuccessful_iteration(Case *c) {
         set_inactive(unsuccessful_direction);
         contract(unsuccessful_direction);
     }
+    if (VERB_OPT >= 3) print_state("Unsuccessful iteration");
     if (!is_converged()) iterate();
 }
 
@@ -117,10 +116,13 @@ vector<int> APPS::inactive() {
 }
 
 void APPS::prune_queue() {
-    if (case_handler_->QueuedCases().size() <= max_queue_length_ - directions_.size())
+    if (case_handler_->QueuedCases().size() <= max_queue_length_ - directions_.size()) {
         return;
+    }
     else {
-        while (case_handler_->QueuedCases().size() > max_queue_length_ - directions_.size()) {
+        int queue_size = max_queue_length_ - directions_.size();
+        if (evaluated_cases_ >= max_evaluations_) queue_size = 1;
+        while (case_handler_->QueuedCases().size() > queue_size) {
             auto dequeued_case = dequeue_case_with_worst_origin();
             if (dequeued_case->origin_case()->id() == GetTentativeBestCase()->id())
                 set_inactive(vector<int>{dequeued_case->origin_direction_index()});
@@ -132,10 +134,12 @@ void APPS::prune_queue() {
 void APPS::print_state(string header) {
     std::stringstream ss;
     ss << header << "|";
-    ss << "Step lengths  : " << vec_to_str(vector<double>(step_lengths_.data(), step_lengths_.data() + step_lengths_.size())) << "|";
-    ss << "Iteration: " << iteration_ << "|";
+    ss << "Iteration:         " << iteration_ << "|";
+    ss << "Evaluated cases:   " << evaluated_cases_ << "|";
+    ss << "Queued cases:      " << case_handler_->QueuedCases().size() << "|";
     ss << "Current best case: " << tentative_best_case_->id().toString().toStdString() << "|";
-    ss << "              OFV: " << tentative_best_case_->objective_function_value();
+    ss << "OFV:               " << tentative_best_case_->objective_function_value();
+    ss << "Step lengths  :    " << vec_to_str(vector<double>(step_lengths_.data(), step_lengths_.data() + step_lengths_.size())) << "|";
     Printer::ext_info(ss.str(), "Optimization", "APPS");
 }
 }

FieldOpt/Optimization/optimizers/GeneticAlgorithm.cpp

@@ -131,12 +131,6 @@ Case *GeneticAlgorithm::generateRandomCase() {
     new_case->SetRealVarValues(erands);
     return  new_case;
 }
-void GeneticAlgorithm::penalizeInitialGeneration() {
-    initializeNormalizers();
-    for (auto chrom : population_) {
-        double pen_ofv = PenalizedOFV(chrom.case_pointer);
-        chrom.case_pointer->set_objective_function_value(pen_ofv);
-    }
-}
+
 }
 }

FieldOpt/Optimization/optimizers/GeneticAlgorithm.h

@@ -120,11 +120,6 @@ class GeneticAlgorithm : public Optimizer {
    */
   Case *generateRandomCase();
 
-  /*!
-   * @brief Initialize normalizers and penalize the initial generation after
-   * the fact (should be start of the first iteration).
-   */
-  void penalizeInitialGeneration();
 };
 
 }

FieldOpt/Optimization/optimizers/RGARDD.cpp

@@ -39,6 +39,7 @@ RGARDD::RGARDD(Settings::Optimizer *settings,
     stagnation_limit_ = settings->parameters().stagnation_limit;
     mating_pool_ = population_;
     if (enable_logging_) {
+        logger_->AddEntry(this);
         logger_->AddEntry(new ConfigurationSummary(this));
     }
 }
@@ -51,9 +52,6 @@ void RGARDD::iterate() {
     if (enable_logging_) {
         logger_->AddEntry(this);
     }
-    if (iteration_ == 0 && penalize_) { // If we're done evaluating the initial population ...
-        penalizeInitialGeneration();
-    }
     population_ = sortPopulation(population_);
 
     if (is_stagnant()) {
@@ -110,6 +108,9 @@ void RGARDD::handleEvaluatedCase(Case *c) {
         population_[index] = mating_pool_[index];
         if (isImprovement(c)) {
             updateTentativeBestCase(c);
+            if (enable_logging_) {
+                logger_->AddEntry(this);
+            }
             if (VERB_OPT >= 1) {
                 Printer::ext_info("New best in generation " + Printer::num2str(iteration_) + ": "
                 + Printer::num2str(GetTentativeBestCase()->objective_function_value()), "RGARDD", "Optimization");

FieldOpt/Optimization/optimizers/bayesian_optimization/EGO.cpp

@@ -44,6 +44,20 @@ EGO::EGO(Settings::Optimizer *settings,
     time_af_opt_ = 0;
     settings_ = settings;
 
+    initializeNormalizers();
+
+    // penalize the base case
+    if (penalize_) {
+        double org_ofv = tentative_best_case_->objective_function_value();
+        double pen_ofv = PenalizedOFV(tentative_best_case_);
+        tentative_best_case_->set_objective_function_value(pen_ofv);
+        if (VERB_OPT >=1) {
+            Printer::ext_info("Penalized base case. " 
+                    "Original value: " + Printer::num2str(org_ofv) + "; "
+                    + "Penalized value: " + Printer::num2str(pen_ofv), "Optimization", "Optimizer");
+        }
+    }
+
     int n_cont_vars = variables->ContinousVariableSize();
 
     if (constraint_handler_->HasBoundaryConstraints()) {
@@ -149,22 +163,13 @@ Optimization::Optimizer::TerminationCondition EGO::IsFinished() {
     return tc;
 }
 void EGO::handleEvaluatedCase(Case *c) {
-    if (!normalizer_ofv_.is_ready())
-        initializeNormalizers();
-
     gp_->add_pattern(c->GetRealVarVector().data(), normalizer_ofv_.normalize(c->objective_function_value()));
     if (isImprovement(c)) {
         updateTentativeBestCase(c);
         Printer::ext_info("Found new tentative best case: " + Printer::num2str(c->objective_function_value()), "Optimization", "EGO");
     }
 }
 void EGO::iterate() {
-    if (!normalizer_ofv_.is_ready()) {
-        initializeNormalizers();
-        // Add base case to GP
-        gp_->add_pattern(tentative_best_case_->GetRealVarVector().data(), normalizer_ofv_.normalize(tentative_best_case_->objective_function_value()));
-    }
-
     if (enable_logging_) {
         logger_->AddEntry(this);
     }