// Copyright 2010-2025 Google LLC // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. // You may obtain a copy of the License at // // http://www.apache.org/licenses/LICENSE-2.0 // // Unless required by applicable law or agreed to in writing, software // distributed under the License is distributed on an "AS IS" BASIS, // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. #include "ortools/glop/variables_info.h" #include #include #include "absl/log/check.h" #include "ortools/lp_data/lp_types.h" #include "ortools/lp_data/sparse.h" namespace operations_research { namespace glop { VariablesInfo::VariablesInfo(const CompactSparseMatrix& matrix) : matrix_(matrix) {} bool VariablesInfo::LoadBoundsAndReturnTrueIfUnchanged( const DenseRow& new_lower_bounds, const DenseRow& new_upper_bounds) { const ColIndex num_cols = matrix_.num_cols(); DCHECK_EQ(num_cols, new_lower_bounds.size()); DCHECK_EQ(num_cols, new_upper_bounds.size()); // Optim if nothing changed. if (lower_bounds_ == new_lower_bounds && upper_bounds_ == new_upper_bounds) { return true; } lower_bounds_ = new_lower_bounds; upper_bounds_ = new_upper_bounds; variable_type_.resize(num_cols, VariableType::UNCONSTRAINED); for (ColIndex col(0); col < num_cols; ++col) { variable_type_[col] = ComputeVariableType(col); } return false; } void VariablesInfo::InitializeFromMutatedState() { const ColIndex num_cols = matrix_.num_cols(); DCHECK_EQ(num_cols, lower_bounds_.size()); DCHECK_EQ(num_cols, upper_bounds_.size()); variable_type_.resize(num_cols, VariableType::UNCONSTRAINED); for (ColIndex col(0); col < num_cols; ++col) { variable_type_[col] = ComputeVariableType(col); } } bool VariablesInfo::LoadBoundsAndReturnTrueIfUnchanged( const DenseRow& variable_lower_bounds, const DenseRow& variable_upper_bounds, const DenseColumn& constraint_lower_bounds, const DenseColumn& constraint_upper_bounds) { const ColIndex num_cols = matrix_.num_cols(); const ColIndex num_variables = variable_upper_bounds.size(); const RowIndex num_rows = constraint_lower_bounds.size(); bool is_unchanged = (num_cols == lower_bounds_.size()); DCHECK_EQ(num_cols, num_variables + RowToColIndex(num_rows)); lower_bounds_.resize(num_cols, 0.0); upper_bounds_.resize(num_cols, 0.0); variable_type_.resize(num_cols, VariableType::FIXED_VARIABLE); // Copy bounds of the variables. for (ColIndex col(0); col < num_variables; ++col) { if (lower_bounds_[col] != variable_lower_bounds[col] || upper_bounds_[col] != variable_upper_bounds[col]) { lower_bounds_[col] = variable_lower_bounds[col]; upper_bounds_[col] = variable_upper_bounds[col]; is_unchanged = false; variable_type_[col] = ComputeVariableType(col); } } // Copy bounds of the slack. for (RowIndex row(0); row < num_rows; ++row) { const ColIndex col = num_variables + RowToColIndex(row); if (lower_bounds_[col] != -constraint_upper_bounds[row] || upper_bounds_[col] != -constraint_lower_bounds[row]) { lower_bounds_[col] = -constraint_upper_bounds[row]; upper_bounds_[col] = -constraint_lower_bounds[row]; is_unchanged = false; variable_type_[col] = ComputeVariableType(col); } } return is_unchanged; } void VariablesInfo::ResetStatusInfo() { const ColIndex num_cols = matrix_.num_cols(); DCHECK_EQ(num_cols, lower_bounds_.size()); DCHECK_EQ(num_cols, upper_bounds_.size()); // TODO(user): These could just be Resized() but there is a bug with the // iteration and resize it seems. Investigate. I suspect the last bucket // is not cleared so you can still iterate on the ones there even if it all // positions before num_cols are set to zero. variable_status_.resize(num_cols, VariableStatus::FREE); can_increase_.ClearAndResize(num_cols); can_decrease_.ClearAndResize(num_cols); is_basic_.ClearAndResize(num_cols); not_basic_.ClearAndResize(num_cols); non_basic_boxed_variables_.ClearAndResize(num_cols); // This one cannot just be resized. boxed_variables_are_relevant_ = true; num_entries_in_relevant_columns_ = 0; relevance_.ClearAndResize(num_cols); } void VariablesInfo::InitializeFromBasisState(ColIndex first_slack_col, ColIndex num_new_cols, const BasisState& state) { ResetStatusInfo(); const ColIndex num_cols = lower_bounds_.size(); DCHECK_LE(num_new_cols, first_slack_col); const ColIndex first_new_col(first_slack_col - num_new_cols); // Compute the status for all the columns (note that the slack variables are // already added at the end of the matrix at this stage). const int state_size = state.statuses.size().value(); for (ColIndex col(0); col < num_cols; ++col) { // Start with the given "warm" status from the BasisState if it exists. VariableStatus status; if (col < first_new_col && col < state_size) { status = state.statuses[col]; } else if (col >= first_slack_col && col - num_new_cols < state_size) { status = state.statuses[col - num_new_cols]; } else { UpdateToNonBasicStatus(col, DefaultVariableStatus(col)); continue; } // Remove incompatibilities between the warm status and the current state. switch (status) { case VariableStatus::BASIC: // Because we just called ResetStatusInfo(), we optimize the call to // UpdateToNonBasicStatus(col) here. In an incremental setting with // almost no work per call, the update of all the DenseBitRow are // visible. variable_status_[col] = VariableStatus::BASIC; is_basic_.Set(col, true); break; case VariableStatus::AT_LOWER_BOUND: if (lower_bounds_[col] == upper_bounds_[col]) { UpdateToNonBasicStatus(col, VariableStatus::FIXED_VALUE); } else { UpdateToNonBasicStatus(col, lower_bounds_[col] == -kInfinity ? DefaultVariableStatus(col) : status); } break; case VariableStatus::AT_UPPER_BOUND: if (lower_bounds_[col] == upper_bounds_[col]) { UpdateToNonBasicStatus(col, VariableStatus::FIXED_VALUE); } else { UpdateToNonBasicStatus(col, upper_bounds_[col] == kInfinity ? DefaultVariableStatus(col) : status); } break; default: UpdateToNonBasicStatus(col, DefaultVariableStatus(col)); } } } int VariablesInfo::ChangeUnusedBasicVariablesToFree( const RowToColMapping& basis) { const ColIndex num_cols = lower_bounds_.size(); is_basic_.ClearAndResize(num_cols); for (const ColIndex col : basis) { UpdateToBasicStatus(col); } int num_no_longer_in_basis = 0; for (ColIndex col(0); col < num_cols; ++col) { if (!is_basic_[col] && variable_status_[col] == VariableStatus::BASIC) { ++num_no_longer_in_basis; if (variable_type_[col] == VariableType::FIXED_VARIABLE) { UpdateToNonBasicStatus(col, VariableStatus::FIXED_VALUE); } else { UpdateToNonBasicStatus(col, VariableStatus::FREE); } } } return num_no_longer_in_basis; } int VariablesInfo::SnapFreeVariablesToBound(Fractional distance, const DenseRow& starting_values) { int num_changes = 0; const ColIndex num_cols = lower_bounds_.size(); for (ColIndex col(0); col < num_cols; ++col) { if (variable_status_[col] != VariableStatus::FREE) continue; if (variable_type_[col] == VariableType::UNCONSTRAINED) continue; const Fractional value = col < starting_values.size() ? starting_values[col] : 0.0; const Fractional diff_ub = upper_bounds_[col] - value; const Fractional diff_lb = value - lower_bounds_[col]; if (diff_lb <= diff_ub) { if (diff_lb <= distance) { ++num_changes; UpdateToNonBasicStatus(col, VariableStatus::AT_LOWER_BOUND); } } else { if (diff_ub <= distance) { ++num_changes; UpdateToNonBasicStatus(col, VariableStatus::AT_UPPER_BOUND); } } } return num_changes; } void VariablesInfo::InitializeToDefaultStatus() { ResetStatusInfo(); const ColIndex num_cols = lower_bounds_.size(); for (ColIndex col(0); col < num_cols; ++col) { UpdateToNonBasicStatus(col, DefaultVariableStatus(col)); } } VariableStatus VariablesInfo::DefaultVariableStatus(ColIndex col) const { DCHECK_GE(col, 0); DCHECK_LT(col, lower_bounds_.size()); if (lower_bounds_[col] == upper_bounds_[col]) { return VariableStatus::FIXED_VALUE; } if (lower_bounds_[col] == -kInfinity && upper_bounds_[col] == kInfinity) { return VariableStatus::FREE; } // Returns the bound with the lowest magnitude. Note that it must be finite // because the VariableStatus::FREE case was tested earlier. DCHECK(IsFinite(lower_bounds_[col]) || IsFinite(upper_bounds_[col])); return std::abs(lower_bounds_[col]) <= std::abs(upper_bounds_[col]) ? VariableStatus::AT_LOWER_BOUND : VariableStatus::AT_UPPER_BOUND; } void VariablesInfo::MakeBoxedVariableRelevant(bool value) { if (value == boxed_variables_are_relevant_) return; boxed_variables_are_relevant_ = value; if (value) { for (const ColIndex col : non_basic_boxed_variables_) { SetRelevance(col, variable_type_[col] != VariableType::FIXED_VARIABLE); } } else { for (const ColIndex col : non_basic_boxed_variables_) { SetRelevance(col, false); } } } void VariablesInfo::UpdateToBasicStatus(ColIndex col) { if (in_dual_phase_one_) { // TODO(user): A bit annoying that we need to test this even if we // don't use the dual. But the cost is minimal. if (lower_bounds_[col] != 0.0) lower_bounds_[col] = -kInfinity; if (upper_bounds_[col] != 0.0) upper_bounds_[col] = +kInfinity; variable_type_[col] = ComputeVariableType(col); } variable_status_[col] = VariableStatus::BASIC; is_basic_.Set(col, true); not_basic_.Set(col, false); can_increase_.Set(col, false); can_decrease_.Set(col, false); non_basic_boxed_variables_.Set(col, false); SetRelevance(col, false); } void VariablesInfo::UpdateToNonBasicStatus(ColIndex col, VariableStatus status) { DCHECK_NE(status, VariableStatus::BASIC); variable_status_[col] = status; is_basic_.Set(col, false); not_basic_.Set(col, true); can_increase_.Set(col, status == VariableStatus::AT_LOWER_BOUND || status == VariableStatus::FREE); can_decrease_.Set(col, status == VariableStatus::AT_UPPER_BOUND || status == VariableStatus::FREE); const bool boxed = variable_type_[col] == VariableType::UPPER_AND_LOWER_BOUNDED; non_basic_boxed_variables_.Set(col, boxed); const bool relevance = status != VariableStatus::FIXED_VALUE && (boxed_variables_are_relevant_ || !boxed); SetRelevance(col, relevance); } const VariableTypeRow& VariablesInfo::GetTypeRow() const { return variable_type_; } const VariableStatusRow& VariablesInfo::GetStatusRow() const { return variable_status_; } const DenseBitRow& VariablesInfo::GetCanIncreaseBitRow() const { return can_increase_; } const DenseBitRow& VariablesInfo::GetCanDecreaseBitRow() const { return can_decrease_; } const DenseBitRow& VariablesInfo::GetIsRelevantBitRow() const { return relevance_; } const DenseBitRow& VariablesInfo::GetIsBasicBitRow() const { return is_basic_; } const DenseBitRow& VariablesInfo::GetNotBasicBitRow() const { return not_basic_; } const DenseBitRow& VariablesInfo::GetNonBasicBoxedVariables() const { return non_basic_boxed_variables_; } EntryIndex VariablesInfo::GetNumEntriesInRelevantColumns() const { return num_entries_in_relevant_columns_; } VariableType VariablesInfo::ComputeVariableType(ColIndex col) const { DCHECK_LE(lower_bounds_[col], upper_bounds_[col]); if (lower_bounds_[col] == -kInfinity) { if (upper_bounds_[col] == kInfinity) { return VariableType::UNCONSTRAINED; } return VariableType::UPPER_BOUNDED; } else if (upper_bounds_[col] == kInfinity) { return VariableType::LOWER_BOUNDED; } else if (lower_bounds_[col] == upper_bounds_[col]) { return VariableType::FIXED_VARIABLE; } else { return VariableType::UPPER_AND_LOWER_BOUNDED; } } void VariablesInfo::SetRelevance(ColIndex col, bool relevance) { if (relevance_.IsSet(col) == relevance) return; if (relevance) { relevance_.Set(col); num_entries_in_relevant_columns_ += matrix_.ColumnNumEntries(col); } else { relevance_.Clear(col); num_entries_in_relevant_columns_ -= matrix_.ColumnNumEntries(col); } } // This is really similar to InitializeFromBasisState() but there is less // cases to consider for TransformToDualPhaseIProblem()/EndDualPhaseI(). void VariablesInfo::UpdateStatusForNewType(ColIndex col) { switch (variable_status_[col]) { case VariableStatus::BASIC: UpdateToBasicStatus(col); break; case VariableStatus::AT_LOWER_BOUND: if (lower_bounds_[col] == upper_bounds_[col]) { UpdateToNonBasicStatus(col, VariableStatus::FIXED_VALUE); } else if (lower_bounds_[col] == -kInfinity) { UpdateToNonBasicStatus(col, DefaultVariableStatus(col)); } else { // TODO(user): This is only needed for boxed variable to update their // relevance. It should probably be done with the type and not the // status update. UpdateToNonBasicStatus(col, variable_status_[col]); } break; case VariableStatus::AT_UPPER_BOUND: if (lower_bounds_[col] == upper_bounds_[col]) { UpdateToNonBasicStatus(col, VariableStatus::FIXED_VALUE); } else if (upper_bounds_[col] == kInfinity) { UpdateToNonBasicStatus(col, DefaultVariableStatus(col)); } else { // TODO(user): Same as in the AT_LOWER_BOUND branch above. UpdateToNonBasicStatus(col, variable_status_[col]); } break; default: // TODO(user): boxed variable that become fixed in // TransformToDualPhaseIProblem() will be changed status twice. Once here, // and once when we make them dual feasible according to their reduced // cost. We should probably just do all at once. UpdateToNonBasicStatus(col, DefaultVariableStatus(col)); } } void VariablesInfo::TransformToDualPhaseIProblem( Fractional dual_feasibility_tolerance, DenseRow::ConstView reduced_costs) { DCHECK(!in_dual_phase_one_); in_dual_phase_one_ = true; saved_lower_bounds_ = lower_bounds_; saved_upper_bounds_ = upper_bounds_; // Transform the bound and type to get a new problem. If this problem has an // optimal value of 0.0, then the problem is dual feasible. And more // importantly, by keeping the same basis, we have a feasible solution of the // original problem. const ColIndex num_cols = matrix_.num_cols(); for (ColIndex col(0); col < num_cols; ++col) { switch (variable_type_[col]) { case VariableType::FIXED_VARIABLE: // ABSL_FALLTHROUGH_INTENDED case VariableType::UPPER_AND_LOWER_BOUNDED: lower_bounds_[col] = 0.0; upper_bounds_[col] = 0.0; variable_type_[col] = VariableType::FIXED_VARIABLE; break; case VariableType::LOWER_BOUNDED: lower_bounds_[col] = 0.0; upper_bounds_[col] = 1.0; variable_type_[col] = VariableType::UPPER_AND_LOWER_BOUNDED; break; case VariableType::UPPER_BOUNDED: lower_bounds_[col] = -1.0; upper_bounds_[col] = 0.0; variable_type_[col] = VariableType::UPPER_AND_LOWER_BOUNDED; break; case VariableType::UNCONSTRAINED: lower_bounds_[col] = -1000.0; upper_bounds_[col] = 1000.0; variable_type_[col] = VariableType::UPPER_AND_LOWER_BOUNDED; break; } // Make sure we start with a feasible dual solution. // If the reduced cost is close to zero, we keep the "default" status. if (variable_type_[col] == VariableType::UPPER_AND_LOWER_BOUNDED) { if (reduced_costs[col] > dual_feasibility_tolerance) { variable_status_[col] = VariableStatus::AT_LOWER_BOUND; } else if (reduced_costs[col] < -dual_feasibility_tolerance) { variable_status_[col] = VariableStatus::AT_UPPER_BOUND; } } UpdateStatusForNewType(col); } } void VariablesInfo::EndDualPhaseI(Fractional dual_feasibility_tolerance, DenseRow::ConstView reduced_costs) { DCHECK(in_dual_phase_one_); in_dual_phase_one_ = false; std::swap(saved_lower_bounds_, lower_bounds_); std::swap(saved_upper_bounds_, upper_bounds_); // This is to clear the memory of the saved bounds since it is no longer // needed. DenseRow empty1, empty2; std::swap(empty1, saved_lower_bounds_); std::swap(empty1, saved_upper_bounds_); // Restore the type and update all other fields. const ColIndex num_cols = matrix_.num_cols(); for (ColIndex col(0); col < num_cols; ++col) { variable_type_[col] = ComputeVariableType(col); // We make sure that the old fixed variables that are now boxed are dual // feasible. // // TODO(user): When there is a choice, use the previous status that might // have been warm-started ? but then this is not high priority since // warm-starting with a non-dual feasible basis seems unfrequent. if (variable_type_[col] == VariableType::UPPER_AND_LOWER_BOUNDED) { if (reduced_costs[col] > dual_feasibility_tolerance) { variable_status_[col] = VariableStatus::AT_LOWER_BOUND; } else if (reduced_costs[col] < -dual_feasibility_tolerance) { variable_status_[col] = VariableStatus::AT_UPPER_BOUND; } } UpdateStatusForNewType(col); } } } // namespace glop } // namespace operations_research