rm step

Author: Yanyu000

highs/pdlp/hipdlp/Meld

@@ -9,6 +9,4 @@ meld restart.hpp $CU_PDLP_CPP_ROOT/include/restart.hpp
 meld restart.cc $CU_PDLP_CPP_ROOT/src/restart.cc
 meld scaling.hpp $CU_PDLP_CPP_ROOT/include/scaling.hpp
 meld scaling.cc $CU_PDLP_CPP_ROOT/src/scaling.cc
-meld solver_results.hpp $CU_PDLP_CPP_ROOT/include/solver_results.hpp
-meld step.hpp $CU_PDLP_CPP_ROOT/include/step.hpp
-meld step.cc $CU_PDLP_CPP_ROOT/src/step.cc
+meld solver_results.hpp $CU_PDLP_CPP_ROOT/include/solver_results.hpp

highs/pdlp/hipdlp/pdhg.cc

@@ -392,6 +392,10 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
   // --- 2. Main PDHG Loop ---
   debugPdlpIterHeaderLog(debug_pdlp_log_file_);
   debugPdlpDataInitialise(&debug_pdlp_data_);
+  debug_pdlp_data_.ax_average_norm = 0.0;
+  debug_pdlp_data_.aty_average_norm = 0.0;
+  debug_pdlp_data_.x_average_norm = 0.0;
+  debug_pdlp_data_.ax_norm = 0.0;
 
   for (int iter = 0; iter < params_.max_iterations; ++iter) {
     debugPdlpIterLog(debug_pdlp_log_file_, iter, &debug_pdlp_data_,
@@ -412,7 +416,6 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
                          (iter > 0 && (iter % PDHG_CHECK_INTERVAL == 0));
 
     if (bool_checking) {
-      std::cout << "sum_weights_ = " << sum_weights_ << "\n";
       ComputeAverageIterate(Ax_avg, ATy_avg);
       // Compute matrix-vector products for current and average iterates
       linalg::Ax(lp, x_current_, Ax_cache_);
@@ -468,8 +471,6 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
           restart_scheme_.Check(iter, current_results, average_results);
 
       if (restart_info.should_restart) {
-        logger_.info("Restarting at iteration " + std::to_string(iter));
-
         std::vector<double> restart_x, restart_y;
         if (restart_info.restart_to_average) {
           // Restart from the average iterate
@@ -501,19 +502,19 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
 
         // Recompute Ax and ATy for the restarted iterates
         linalg::Ax(lp, x_current_, Ax_cache_);
+        debug_pdlp_data_.ax_norm = linalg::vector_norm(Ax_cache_);
         linalg::ATy(lp, y_current_, ATy_cache_);
-
+        debug_pdlp_data_.aty_norm = linalg::vector_norm(ATy_cache_);
         restart_scheme_.SetLastRestartIter(iter);
       }
     } else {
       // If not checking, still need Ax and ATy for the update
       linalg::Ax(lp, x_current_, Ax_cache_);
+      debug_pdlp_data_.ax_norm = linalg::vector_norm(Ax_cache_);
       linalg::ATy(lp, y_current_, ATy_cache_);
+      debug_pdlp_data_.aty_norm = linalg::vector_norm(ATy_cache_);
     }
 
-    // Record Ax norm for debugging
-    debug_pdlp_data_.ax_norm = linalg::vector_norm(Ax_cache_);
-
     // --- 5. Core PDHG Update Step ---
     bool step_success = true;
 
@@ -654,8 +655,6 @@ void PDLPSolver::UpdateAverageIterates(const std::vector<double>& x,
   }
 
   sum_weights_ += dMeanStepSize;
-  printf("x_sum norm before %g after %g \n", x_norm_before,
-         linalg::vector_norm(x_sum_));
 }
 
 void PDLPSolver::ComputeAverageIterate(std::vector<double>& ax_avg,
@@ -670,10 +669,11 @@ void PDLPSolver::ComputeAverageIterate(std::vector<double>& ax_avg,
     y_avg_[i] = y_sum_[i] * dDualScale; 
   }
 
-  debug_pdlp_data_.x_average_norm = linalg::vector_norm(x_avg_);
+  debug_pdlp_data_.x_average_norm = linalg::vector_norm_squared(x_avg_);
   linalg::Ax(lp_, x_avg_, ax_avg);
   linalg::ATy(lp_, y_avg_, aty_avg);
-  debug_pdlp_data_.ax_average_norm = linalg::vector_norm(ax_avg);
+  debug_pdlp_data_.ax_average_norm = linalg::vector_norm_squared(ax_avg);
+  debug_pdlp_data_.aty_average_norm = linalg::vector_norm_squared(aty_avg);
 }
 
 // lambda = c - proj_{\Lambda}(c - K^T y)

highs/pdlp/hipdlp/restart.cc

@@ -17,6 +17,13 @@
 #include "io/HighsIO.h"          // For pdlpLogging
 #include "pdlp/cupdlp/cupdlp.h"  // For pdlpLogging
 
+// Helper function to compute the score
+static double compute_score(double beta, double primal_feas, double dual_feas, double duality_gap) {
+    return std::sqrt(beta * primal_feas * primal_feas +
+                     dual_feas * dual_feas / beta +
+                     duality_gap * duality_gap);
+}
+
 // Initializes the restart scheme with parameters and initial results
 void RestartScheme::Initialize(const SolverResults& results) {
   strategy_ = params_->restart_strategy;
@@ -28,34 +35,16 @@ void RestartScheme::Initialize(const SolverResults& results) {
   //    in the paper, primal weight is the ω:
   //    ω = √β
   beta_ = params_->omega * params_->omega;
-
-  last_restart_iter_ = 0;
-  last_restart_score_ = std::numeric_limits<double>::infinity();
-  last_candidate_score_ = std::numeric_limits<double>::infinity();
-}
-
-// Computes a weighted score to evaluate solution quality
-double RestartScheme::ComputeRestartScore(const SolverResults& results) {
-  double weight_squared = beta_;
-  double primal_feas_sq =
-      results.primal_feasibility * results.primal_feasibility;
-  double dual_feas_sq = results.dual_feasibility * results.dual_feasibility;
-  double gap_sq = results.duality_gap * results.duality_gap;
-
-  // debugPdlpRestarScoretLog(debug_pdlp_log_file_, weight_squared,
-  // results.primal_feasibility, results.dual_feasibility, results.duality_gap);
-
-  return std::sqrt(weight_squared * primal_feas_sq +
-                   dual_feas_sq / weight_squared + gap_sq);
 }
 
+/*
 // Main logic to check if a restart is needed
 RestartInfo RestartScheme::Check(int current_iter,
                                  const SolverResults& current_results,
                                  const SolverResults& average_results) {
   RestartInfo info;
 
-  if (current_iter == last_restart_iter_) {
+  if (current_iter <= last_restart_iter_) {
     double current_score = ComputeRestartScore(current_results);
     last_restart_score_ = current_score;
     last_candidate_score_ = current_score;
@@ -77,8 +66,8 @@ RestartInfo RestartScheme::Check(int current_iter,
     case RestartStrategy::ADAPTIVE_RESTART: {
       double current_score = ComputeRestartScore(current_results);
       double average_score = ComputeRestartScore(average_results);
-      debugPdlpRestartLog(debug_pdlp_log_file_, current_iter, current_score,
-                          average_score);
+      //debugPdlpRestartLog(debug_pdlp_log_file_, current_iter, current_score,
+      //                    average_score);
 
       // Choose the best candidate (current vs. average) based on the score
       double candidate_score = std::min(current_score, average_score);
@@ -133,3 +122,88 @@ RestartInfo RestartScheme::Check(int current_iter,
 
   return info;
 }
+*/
+
+RestartInfo RestartScheme::Check(int current_iter,
+                                 const SolverResults& current_results,
+                                 const SolverResults& average_results) {
+  if (strategy_ != RestartStrategy::ADAPTIVE_RESTART) {
+    if (strategy_ == RestartStrategy::FIXED_RESTART && 
+        current_iter - last_restart_iter_ >= fixed_restart_interval_) {
+        return RestartInfo(true, true);
+    }
+    return RestartInfo(false, false); // ✨ Corrected return
+  }
+
+  // Base Case: First check after a restart. Initialize state and exit.
+  if (current_iter == last_restart_iter_) {
+    primal_feas_last_restart_ = current_results.primal_feasibility;
+    dual_feas_last_restart_ = current_results.dual_feasibility;
+    duality_gap_last_restart_ = current_results.duality_gap;
+
+    primal_feas_last_candidate_ = current_results.primal_feasibility;
+    dual_feas_last_candidate_ = current_results.dual_feasibility;
+    duality_gap_last_candidate_ = current_results.duality_gap;
+    
+    return {false, false};
+  }
+
+  // 1. Calculate scores and determine the best candidate
+  double mu_current = compute_score(beta_, current_results.primal_feasibility,
+                                    current_results.dual_feasibility, current_results.duality_gap);
+  double mu_average = compute_score(beta_, average_results.primal_feasibility,
+                                    average_results.dual_feasibility, average_results.duality_gap);
+  debugPdlpRestartLog(debug_pdlp_log_file_, current_iter, mu_current,
+                          mu_average);
+  double candidate_score;
+  bool restart_to_average;
+  if (mu_current < mu_average) {
+    candidate_score = mu_current;
+    restart_to_average = false;
+  } else {
+    candidate_score = mu_average;
+    restart_to_average = true;
+  }
+
+  // 2. Check the three restart conditions in order
+  bool should_restart = false;
+  PDHG_restart_choice restart_choice_for_logic = PDHG_NO_RESTART;
+
+  if ((current_iter - last_restart_iter_) >= 0.36 * current_iter) {
+    // Condition 1: Artificial Restart
+    should_restart = true;
+  } else {
+    double mu_last_restart = compute_score(beta_, primal_feas_last_restart_,
+                                           dual_feas_last_restart_, duality_gap_last_restart_);
+    
+    if (candidate_score < sufficient_decay_factor_ * mu_last_restart) {
+      // Condition 2: Sufficient Decay
+      should_restart = true;
+    } else {
+      double mu_last_candidate = compute_score(beta_, primal_feas_last_candidate_,
+                                                 dual_feas_last_candidate_, duality_gap_last_candidate_);
+      
+      if (candidate_score < necessary_decay_factor_ * mu_last_restart && candidate_score > mu_last_candidate) {
+        // Condition 3: Necessary Decay
+        should_restart = true;
+      }
+    }
+  }
+  
+  // 3. ALWAYS update the "last candidate" metrics for the next check
+  if (restart_to_average) {
+    primal_feas_last_candidate_ = average_results.primal_feasibility;
+    dual_feas_last_candidate_ = average_results.dual_feasibility;
+    duality_gap_last_candidate_ = average_results.duality_gap;
+  } else {
+    primal_feas_last_candidate_ = current_results.primal_feasibility;
+    dual_feas_last_candidate_ = current_results.dual_feasibility;
+    duality_gap_last_candidate_ = current_results.duality_gap;
+  }
+
+  // NOTE: Your main solver loop is responsible for calling `SetLastRestartIter(current_iter)`
+  // when it receives `should_restart = true`. This function does not modify that state directly,
+  // it only determines if it *should* happen.
+  
+  return RestartInfo(should_restart, restart_to_average);
+}

highs/pdlp/hipdlp/restart.hpp

@@ -25,12 +25,16 @@ struct RestartInfo {
   bool should_restart = false;
   bool restart_to_average =
       true;  // If true, restart to average; otherwise, to current
+
+  RestartInfo(bool should_restart_, bool restart_to_average_)
+      : should_restart(should_restart_), restart_to_average(restart_to_average_){};
+
+  RestartInfo() = default;
 };
 
 class RestartScheme {
  public:
   RestartScheme() = default;
-
   void Initialize(const SolverResults& results);
 
   // Checks if a restart should be performed based on the chosen strategy
@@ -61,17 +65,20 @@ class RestartScheme {
   FILE* debug_pdlp_log_file_ = nullptr;
   DebugPdlpData* debug_pdlp_data_;
 
-  // Computes a merit score for a given set of residuals
-  double ComputeRestartScore(const SolverResults& results);
-
   RestartStrategy strategy_ = RestartStrategy::NO_RESTART;
   int fixed_restart_interval_ = 100;
   int last_restart_iter_ = 0;
   double beta_;
 
   // State for adaptive restart
-  double last_restart_score_ = 1;
-  double last_candidate_score_ = 1;
+  double primal_feas_last_restart_ = 0.0;
+  double dual_feas_last_restart_ = 0.0;
+  double duality_gap_last_restart_ = 0.0;
+  
+  double primal_feas_last_candidate_ = 0.0;
+  double dual_feas_last_candidate_ = 0.0;
+  double duality_gap_last_candidate_ = 0.0;
+
   double sufficient_decay_factor_ = 0.2;
   double necessary_decay_factor_ = 0.8;
 };