forked from google/or-tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
preprocessor.cc
3947 lines (3576 loc) · 165 KB
/
preprocessor.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// Copyright 2010-2021 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/preprocessor.h"
#include <cstdint>
#include <limits>
#include "absl/strings/str_format.h"
#include "ortools/base/iterator_adaptors.h"
#include "ortools/base/strong_vector.h"
#include "ortools/glop/revised_simplex.h"
#include "ortools/glop/status.h"
#include "ortools/lp_data/lp_data_utils.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/lp_data/lp_utils.h"
#include "ortools/lp_data/matrix_utils.h"
namespace operations_research {
namespace glop {
using ::util::Reverse;
namespace {
// Returns an interval as an human readable string for debugging.
std::string IntervalString(Fractional lb, Fractional ub) {
return absl::StrFormat("[%g, %g]", lb, ub);
}
#if defined(_MSC_VER)
double trunc(double d) { return d > 0 ? floor(d) : ceil(d); }
#endif
} // namespace
// --------------------------------------------------------
// Preprocessor
// --------------------------------------------------------
Preprocessor::Preprocessor(const GlopParameters* parameters)
: status_(ProblemStatus::INIT),
parameters_(*parameters),
in_mip_context_(false),
infinite_time_limit_(TimeLimit::Infinite()),
time_limit_(infinite_time_limit_.get()) {}
Preprocessor::~Preprocessor() {}
// --------------------------------------------------------
// MainLpPreprocessor
// --------------------------------------------------------
#define RUN_PREPROCESSOR(name) \
RunAndPushIfRelevant(std::unique_ptr<Preprocessor>(new name(¶meters_)), \
#name, time_limit_, lp)
bool MainLpPreprocessor::Run(LinearProgram* lp) {
RETURN_VALUE_IF_NULL(lp, false);
default_logger_.EnableLogging(parameters_.log_search_progress());
default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
SOLVER_LOG(logger_, "");
SOLVER_LOG(logger_, "Starting presolve...");
initial_num_rows_ = lp->num_constraints();
initial_num_cols_ = lp->num_variables();
initial_num_entries_ = lp->num_entries();
if (parameters_.use_preprocessing()) {
RUN_PREPROCESSOR(ShiftVariableBoundsPreprocessor);
// We run it a few times because running one preprocessor may allow another
// one to remove more stuff.
const int kMaxNumPasses = 20;
for (int i = 0; i < kMaxNumPasses; ++i) {
const int old_stack_size = preprocessors_.size();
RUN_PREPROCESSOR(FixedVariablePreprocessor);
RUN_PREPROCESSOR(SingletonPreprocessor);
RUN_PREPROCESSOR(ForcingAndImpliedFreeConstraintPreprocessor);
RUN_PREPROCESSOR(FreeConstraintPreprocessor);
RUN_PREPROCESSOR(ImpliedFreePreprocessor);
RUN_PREPROCESSOR(UnconstrainedVariablePreprocessor);
RUN_PREPROCESSOR(DoubletonFreeColumnPreprocessor);
RUN_PREPROCESSOR(DoubletonEqualityRowPreprocessor);
// Abort early if none of the preprocessors did something. Technically
// this is true if none of the preprocessors above needs postsolving,
// which has exactly the same meaning for these particular preprocessors.
if (preprocessors_.size() == old_stack_size) {
// We use i here because the last pass did nothing.
SOLVER_LOG(logger_, "Reached fixed point after presolve pass #", i);
break;
}
}
RUN_PREPROCESSOR(EmptyColumnPreprocessor);
RUN_PREPROCESSOR(EmptyConstraintPreprocessor);
// TODO(user): Run them in the loop above if the effect on the running time
// is good. This needs more investigation.
RUN_PREPROCESSOR(ProportionalColumnPreprocessor);
RUN_PREPROCESSOR(ProportionalRowPreprocessor);
// If DualizerPreprocessor was run, we need to do some extra preprocessing.
// This is because it currently adds a lot of zero-cost singleton columns.
const int old_stack_size = preprocessors_.size();
// TODO(user): We probably want to scale the costs before and after this
// preprocessor so that the rhs/objective of the dual are with a good
// magnitude.
RUN_PREPROCESSOR(DualizerPreprocessor);
if (old_stack_size != preprocessors_.size()) {
RUN_PREPROCESSOR(SingletonPreprocessor);
RUN_PREPROCESSOR(FreeConstraintPreprocessor);
RUN_PREPROCESSOR(UnconstrainedVariablePreprocessor);
RUN_PREPROCESSOR(EmptyColumnPreprocessor);
RUN_PREPROCESSOR(EmptyConstraintPreprocessor);
}
RUN_PREPROCESSOR(SingletonColumnSignPreprocessor);
}
// The scaling is controlled by use_scaling, not use_preprocessing.
RUN_PREPROCESSOR(ScalingPreprocessor);
return !preprocessors_.empty();
}
#undef RUN_PREPROCESSOR
void MainLpPreprocessor::RunAndPushIfRelevant(
std::unique_ptr<Preprocessor> preprocessor, const std::string& name,
TimeLimit* time_limit, LinearProgram* lp) {
RETURN_IF_NULL(preprocessor);
RETURN_IF_NULL(time_limit);
if (status_ != ProblemStatus::INIT || time_limit->LimitReached()) return;
const double start_time = time_limit->GetElapsedTime();
preprocessor->SetTimeLimit(time_limit);
// No need to run the preprocessor if the lp is empty.
// TODO(user): without this test, the code is failing as of 2013-03-18.
if (lp->num_variables() == 0 && lp->num_constraints() == 0) {
status_ = ProblemStatus::OPTIMAL;
return;
}
if (preprocessor->Run(lp)) {
const EntryIndex new_num_entries = lp->num_entries();
const double preprocess_time = time_limit->GetElapsedTime() - start_time;
SOLVER_LOG(logger_,
absl::StrFormat(
"%-45s: %d(%d) rows, %d(%d) columns, %d(%d) entries. (%fs)",
name, lp->num_constraints().value(),
(lp->num_constraints() - initial_num_rows_).value(),
lp->num_variables().value(),
(lp->num_variables() - initial_num_cols_).value(),
// static_cast<int64_t> is needed because the Android port
// uses int32_t.
static_cast<int64_t>(new_num_entries.value()),
static_cast<int64_t>(new_num_entries.value() -
initial_num_entries_.value()),
preprocess_time));
status_ = preprocessor->status();
preprocessors_.push_back(std::move(preprocessor));
return;
} else {
// Even if a preprocessor returns false (i.e. no need for postsolve), it
// can detect an issue with the problem.
status_ = preprocessor->status();
if (status_ != ProblemStatus::INIT) {
SOLVER_LOG(logger_, name, " detected that the problem is ",
GetProblemStatusString(status_));
}
}
}
void MainLpPreprocessor::RecoverSolution(ProblemSolution* solution) const {
SCOPED_INSTRUCTION_COUNT(time_limit_);
for (const auto& p : gtl::reversed_view(preprocessors_)) {
p->RecoverSolution(solution);
}
}
void MainLpPreprocessor::DestructiveRecoverSolution(ProblemSolution* solution) {
SCOPED_INSTRUCTION_COUNT(time_limit_);
while (!preprocessors_.empty()) {
preprocessors_.back()->RecoverSolution(solution);
preprocessors_.pop_back();
}
}
// --------------------------------------------------------
// ColumnDeletionHelper
// --------------------------------------------------------
void ColumnsSaver::SaveColumn(ColIndex col, const SparseColumn& column) {
const int index = saved_columns_.size();
CHECK(saved_columns_index_.insert({col, index}).second);
saved_columns_.push_back(column);
}
void ColumnsSaver::SaveColumnIfNotAlreadyDone(ColIndex col,
const SparseColumn& column) {
const int index = saved_columns_.size();
const bool inserted = saved_columns_index_.insert({col, index}).second;
if (inserted) saved_columns_.push_back(column);
}
const SparseColumn& ColumnsSaver::SavedColumn(ColIndex col) const {
const auto it = saved_columns_index_.find(col);
CHECK(it != saved_columns_index_.end());
return saved_columns_[it->second];
}
const SparseColumn& ColumnsSaver::SavedOrEmptyColumn(ColIndex col) const {
const auto it = saved_columns_index_.find(col);
return it == saved_columns_index_.end() ? empty_column_
: saved_columns_[it->second];
}
void ColumnDeletionHelper::Clear() {
is_column_deleted_.clear();
stored_value_.clear();
}
void ColumnDeletionHelper::MarkColumnForDeletion(ColIndex col) {
MarkColumnForDeletionWithState(col, 0.0, VariableStatus::FREE);
}
void ColumnDeletionHelper::MarkColumnForDeletionWithState(
ColIndex col, Fractional fixed_value, VariableStatus status) {
DCHECK_GE(col, 0);
if (col >= is_column_deleted_.size()) {
is_column_deleted_.resize(col + 1, false);
stored_value_.resize(col + 1, 0.0);
stored_status_.resize(col + 1, VariableStatus::FREE);
}
is_column_deleted_[col] = true;
stored_value_[col] = fixed_value;
stored_status_[col] = status;
}
void ColumnDeletionHelper::RestoreDeletedColumns(
ProblemSolution* solution) const {
DenseRow new_primal_values;
VariableStatusRow new_variable_statuses;
ColIndex old_index(0);
for (ColIndex col(0); col < is_column_deleted_.size(); ++col) {
if (is_column_deleted_[col]) {
new_primal_values.push_back(stored_value_[col]);
new_variable_statuses.push_back(stored_status_[col]);
} else {
new_primal_values.push_back(solution->primal_values[old_index]);
new_variable_statuses.push_back(solution->variable_statuses[old_index]);
++old_index;
}
}
// Copy the end of the vectors and swap them with the ones in solution.
const ColIndex num_cols = solution->primal_values.size();
DCHECK_EQ(num_cols, solution->variable_statuses.size());
for (; old_index < num_cols; ++old_index) {
new_primal_values.push_back(solution->primal_values[old_index]);
new_variable_statuses.push_back(solution->variable_statuses[old_index]);
}
new_primal_values.swap(solution->primal_values);
new_variable_statuses.swap(solution->variable_statuses);
}
// --------------------------------------------------------
// RowDeletionHelper
// --------------------------------------------------------
void RowDeletionHelper::Clear() { is_row_deleted_.clear(); }
void RowDeletionHelper::MarkRowForDeletion(RowIndex row) {
DCHECK_GE(row, 0);
if (row >= is_row_deleted_.size()) {
is_row_deleted_.resize(row + 1, false);
}
is_row_deleted_[row] = true;
}
void RowDeletionHelper::UnmarkRow(RowIndex row) {
if (row >= is_row_deleted_.size()) return;
is_row_deleted_[row] = false;
}
const DenseBooleanColumn& RowDeletionHelper::GetMarkedRows() const {
return is_row_deleted_;
}
void RowDeletionHelper::RestoreDeletedRows(ProblemSolution* solution) const {
DenseColumn new_dual_values;
ConstraintStatusColumn new_constraint_statuses;
RowIndex old_index(0);
const RowIndex end = is_row_deleted_.size();
for (RowIndex row(0); row < end; ++row) {
if (is_row_deleted_[row]) {
new_dual_values.push_back(0.0);
new_constraint_statuses.push_back(ConstraintStatus::BASIC);
} else {
new_dual_values.push_back(solution->dual_values[old_index]);
new_constraint_statuses.push_back(
solution->constraint_statuses[old_index]);
++old_index;
}
}
// Copy the end of the vectors and swap them with the ones in solution.
const RowIndex num_rows = solution->dual_values.size();
DCHECK_EQ(num_rows, solution->constraint_statuses.size());
for (; old_index < num_rows; ++old_index) {
new_dual_values.push_back(solution->dual_values[old_index]);
new_constraint_statuses.push_back(solution->constraint_statuses[old_index]);
}
new_dual_values.swap(solution->dual_values);
new_constraint_statuses.swap(solution->constraint_statuses);
}
// --------------------------------------------------------
// EmptyColumnPreprocessor
// --------------------------------------------------------
namespace {
// Computes the status of a variable given its value and bounds. This only works
// with a value exactly at one of the bounds, or a value of 0.0 for free
// variables.
VariableStatus ComputeVariableStatus(Fractional value, Fractional lower_bound,
Fractional upper_bound) {
if (lower_bound == upper_bound) {
DCHECK_EQ(value, lower_bound);
DCHECK(IsFinite(lower_bound));
return VariableStatus::FIXED_VALUE;
}
if (value == lower_bound) {
DCHECK_NE(lower_bound, -kInfinity);
return VariableStatus::AT_LOWER_BOUND;
}
if (value == upper_bound) {
DCHECK_NE(upper_bound, kInfinity);
return VariableStatus::AT_UPPER_BOUND;
}
// TODO(user): restrict this to unbounded variables with a value of zero.
// We can't do that when postsolving infeasible problem. Don't call postsolve
// on an infeasible problem?
return VariableStatus::FREE;
}
// Returns the input with the smallest magnitude or zero if both are infinite.
Fractional MinInMagnitudeOrZeroIfInfinite(Fractional a, Fractional b) {
const Fractional value = std::abs(a) < std::abs(b) ? a : b;
return IsFinite(value) ? value : 0.0;
}
Fractional MagnitudeOrZeroIfInfinite(Fractional value) {
return IsFinite(value) ? std::abs(value) : 0.0;
}
// Returns the maximum magnitude of the finite variable bounds of the given
// linear program.
Fractional ComputeMaxVariableBoundsMagnitude(const LinearProgram& lp) {
Fractional max_bounds_magnitude = 0.0;
const ColIndex num_cols = lp.num_variables();
for (ColIndex col(0); col < num_cols; ++col) {
max_bounds_magnitude = std::max(
max_bounds_magnitude,
std::max(MagnitudeOrZeroIfInfinite(lp.variable_lower_bounds()[col]),
MagnitudeOrZeroIfInfinite(lp.variable_upper_bounds()[col])));
}
return max_bounds_magnitude;
}
} // namespace
bool EmptyColumnPreprocessor::Run(LinearProgram* lp) {
SCOPED_INSTRUCTION_COUNT(time_limit_);
RETURN_VALUE_IF_NULL(lp, false);
column_deletion_helper_.Clear();
const ColIndex num_cols = lp->num_variables();
for (ColIndex col(0); col < num_cols; ++col) {
if (lp->GetSparseColumn(col).IsEmpty()) {
const Fractional lower_bound = lp->variable_lower_bounds()[col];
const Fractional upper_bound = lp->variable_upper_bounds()[col];
const Fractional objective_coefficient =
lp->GetObjectiveCoefficientForMinimizationVersion(col);
Fractional value;
if (objective_coefficient == 0) {
// Any feasible value will do.
if (upper_bound != kInfinity) {
value = upper_bound;
} else {
if (lower_bound != -kInfinity) {
value = lower_bound;
} else {
value = Fractional(0.0);
}
}
} else {
value = objective_coefficient > 0 ? lower_bound : upper_bound;
if (!IsFinite(value)) {
VLOG(1) << "Problem INFEASIBLE_OR_UNBOUNDED, empty column " << col
<< " has a minimization cost of " << objective_coefficient
<< " and bounds"
<< " [" << lower_bound << "," << upper_bound << "]";
status_ = ProblemStatus::INFEASIBLE_OR_UNBOUNDED;
return false;
}
lp->SetObjectiveOffset(lp->objective_offset() +
value * lp->objective_coefficients()[col]);
}
column_deletion_helper_.MarkColumnForDeletionWithState(
col, value, ComputeVariableStatus(value, lower_bound, upper_bound));
}
}
lp->DeleteColumns(column_deletion_helper_.GetMarkedColumns());
return !column_deletion_helper_.IsEmpty();
}
void EmptyColumnPreprocessor::RecoverSolution(ProblemSolution* solution) const {
SCOPED_INSTRUCTION_COUNT(time_limit_);
RETURN_IF_NULL(solution);
column_deletion_helper_.RestoreDeletedColumns(solution);
}
// --------------------------------------------------------
// ProportionalColumnPreprocessor
// --------------------------------------------------------
namespace {
// Subtracts 'multiple' times the column col of the given linear program from
// the constraint bounds. That is, for a non-zero entry of coefficient c,
// c * multiple is subtracted from both the constraint upper and lower bound.
void SubtractColumnMultipleFromConstraintBound(ColIndex col,
Fractional multiple,
LinearProgram* lp) {
DenseColumn* lbs = lp->mutable_constraint_lower_bounds();
DenseColumn* ubs = lp->mutable_constraint_upper_bounds();
for (const SparseColumn::Entry e : lp->GetSparseColumn(col)) {
const RowIndex row = e.row();
const Fractional delta = multiple * e.coefficient();
(*lbs)[row] -= delta;
(*ubs)[row] -= delta;
}
// While not needed for correctness, this allows the presolved problem to
// have the same objective value as the original one.
lp->SetObjectiveOffset(lp->objective_offset() +
lp->objective_coefficients()[col] * multiple);
}
// Struct used to detect proportional columns with the same cost. For that, a
// vector of such struct will be sorted, and only the columns that end up
// together need to be compared.
struct ColumnWithRepresentativeAndScaledCost {
ColumnWithRepresentativeAndScaledCost(ColIndex _col, ColIndex _representative,
Fractional _scaled_cost)
: col(_col), representative(_representative), scaled_cost(_scaled_cost) {}
ColIndex col;
ColIndex representative;
Fractional scaled_cost;
bool operator<(const ColumnWithRepresentativeAndScaledCost& other) const {
if (representative == other.representative) {
if (scaled_cost == other.scaled_cost) {
return col < other.col;
}
return scaled_cost < other.scaled_cost;
}
return representative < other.representative;
}
};
} // namespace
bool ProportionalColumnPreprocessor::Run(LinearProgram* lp) {
SCOPED_INSTRUCTION_COUNT(time_limit_);
RETURN_VALUE_IF_NULL(lp, false);
ColMapping mapping = FindProportionalColumns(
lp->GetSparseMatrix(), parameters_.preprocessor_zero_tolerance());
// Compute some statistics and make each class representative point to itself
// in the mapping. Also store the columns that are proportional to at least
// another column in proportional_columns to iterate on them more efficiently.
//
// TODO(user): Change FindProportionalColumns for this?
int num_proportionality_classes = 0;
std::vector<ColIndex> proportional_columns;
for (ColIndex col(0); col < mapping.size(); ++col) {
const ColIndex representative = mapping[col];
if (representative != kInvalidCol) {
if (mapping[representative] == kInvalidCol) {
proportional_columns.push_back(representative);
++num_proportionality_classes;
mapping[representative] = representative;
}
proportional_columns.push_back(col);
}
}
if (proportional_columns.empty()) return false;
VLOG(1) << "The problem contains " << proportional_columns.size()
<< " columns which belong to " << num_proportionality_classes
<< " proportionality classes.";
// Note(user): using the first coefficient may not give the best precision.
const ColIndex num_cols = lp->num_variables();
column_factors_.assign(num_cols, 0.0);
for (const ColIndex col : proportional_columns) {
const SparseColumn& column = lp->GetSparseColumn(col);
column_factors_[col] = column.GetFirstCoefficient();
}
// This is only meaningful for column representative.
//
// The reduced cost of a column is 'cost - dual_values.column' and we know
// that for all proportional columns, 'dual_values.column /
// column_factors_[col]' is the same. Here, we bound this quantity which is
// related to the cost 'slope' of a proportional column:
// cost / column_factors_[col].
DenseRow slope_lower_bound(num_cols, -kInfinity);
DenseRow slope_upper_bound(num_cols, +kInfinity);
for (const ColIndex col : proportional_columns) {
const ColIndex representative = mapping[col];
// We reason in terms of a minimization problem here.
const bool is_rc_positive_or_zero =
(lp->variable_upper_bounds()[col] == kInfinity);
const bool is_rc_negative_or_zero =
(lp->variable_lower_bounds()[col] == -kInfinity);
bool is_slope_upper_bounded = is_rc_positive_or_zero;
bool is_slope_lower_bounded = is_rc_negative_or_zero;
if (column_factors_[col] < 0.0) {
std::swap(is_slope_lower_bounded, is_slope_upper_bounded);
}
const Fractional slope =
lp->GetObjectiveCoefficientForMinimizationVersion(col) /
column_factors_[col];
if (is_slope_lower_bounded) {
slope_lower_bound[representative] =
std::max(slope_lower_bound[representative], slope);
}
if (is_slope_upper_bounded) {
slope_upper_bound[representative] =
std::min(slope_upper_bound[representative], slope);
}
}
// Deal with empty slope intervals.
for (const ColIndex col : proportional_columns) {
const ColIndex representative = mapping[col];
// This is only needed for class representative columns.
if (representative == col) {
if (!IsSmallerWithinFeasibilityTolerance(
slope_lower_bound[representative],
slope_upper_bound[representative])) {
VLOG(1) << "Problem INFEASIBLE_OR_UNBOUNDED, no feasible dual values"
<< " can satisfy the constraints of the proportional columns"
<< " with representative " << representative << "."
<< " the associated quantity must be in ["
<< slope_lower_bound[representative] << ","
<< slope_upper_bound[representative] << "].";
status_ = ProblemStatus::INFEASIBLE_OR_UNBOUNDED;
return false;
}
}
}
// Now, fix the columns that can be fixed to one of their bounds.
for (const ColIndex col : proportional_columns) {
const ColIndex representative = mapping[col];
const Fractional slope =
lp->GetObjectiveCoefficientForMinimizationVersion(col) /
column_factors_[col];
// The scaled reduced cost is slope - quantity.
bool variable_can_be_fixed = false;
Fractional target_bound = 0.0;
const Fractional lower_bound = lp->variable_lower_bounds()[col];
const Fractional upper_bound = lp->variable_upper_bounds()[col];
if (!IsSmallerWithinFeasibilityTolerance(slope_lower_bound[representative],
slope)) {
// The scaled reduced cost is < 0.
variable_can_be_fixed = true;
target_bound = (column_factors_[col] >= 0.0) ? upper_bound : lower_bound;
} else if (!IsSmallerWithinFeasibilityTolerance(
slope, slope_upper_bound[representative])) {
// The scaled reduced cost is > 0.
variable_can_be_fixed = true;
target_bound = (column_factors_[col] >= 0.0) ? lower_bound : upper_bound;
}
if (variable_can_be_fixed) {
// Clear mapping[col] so this column will not be considered for the next
// stage.
mapping[col] = kInvalidCol;
if (!IsFinite(target_bound)) {
VLOG(1) << "Problem INFEASIBLE_OR_UNBOUNDED.";
status_ = ProblemStatus::INFEASIBLE_OR_UNBOUNDED;
return false;
} else {
SubtractColumnMultipleFromConstraintBound(col, target_bound, lp);
column_deletion_helper_.MarkColumnForDeletionWithState(
col, target_bound,
ComputeVariableStatus(target_bound, lower_bound, upper_bound));
}
}
}
// Merge the variables with the same scaled cost.
std::vector<ColumnWithRepresentativeAndScaledCost> sorted_columns;
for (const ColIndex col : proportional_columns) {
const ColIndex representative = mapping[col];
// This test is needed because we already removed some columns.
if (mapping[col] != kInvalidCol) {
sorted_columns.push_back(ColumnWithRepresentativeAndScaledCost(
col, representative,
lp->objective_coefficients()[col] / column_factors_[col]));
}
}
std::sort(sorted_columns.begin(), sorted_columns.end());
// All this will be needed during postsolve.
merged_columns_.assign(num_cols, kInvalidCol);
lower_bounds_.assign(num_cols, -kInfinity);
upper_bounds_.assign(num_cols, kInfinity);
new_lower_bounds_.assign(num_cols, -kInfinity);
new_upper_bounds_.assign(num_cols, kInfinity);
for (int i = 0; i < sorted_columns.size();) {
const ColIndex target_col = sorted_columns[i].col;
const ColIndex target_representative = sorted_columns[i].representative;
const Fractional target_scaled_cost = sorted_columns[i].scaled_cost;
// Save the initial bounds before modifying them.
lower_bounds_[target_col] = lp->variable_lower_bounds()[target_col];
upper_bounds_[target_col] = lp->variable_upper_bounds()[target_col];
int num_merged = 0;
for (++i; i < sorted_columns.size(); ++i) {
if (sorted_columns[i].representative != target_representative) break;
if (std::abs(sorted_columns[i].scaled_cost - target_scaled_cost) >=
parameters_.preprocessor_zero_tolerance()) {
break;
}
++num_merged;
const ColIndex col = sorted_columns[i].col;
const Fractional lower_bound = lp->variable_lower_bounds()[col];
const Fractional upper_bound = lp->variable_upper_bounds()[col];
lower_bounds_[col] = lower_bound;
upper_bounds_[col] = upper_bound;
merged_columns_[col] = target_col;
// This is a bit counter intuitive, but when a column is divided by x,
// the corresponding bounds have to be multiplied by x.
const Fractional bound_factor =
column_factors_[col] / column_factors_[target_col];
// We need to shift the variable so that a basic solution of the new
// problem can easily be converted to a basic solution of the original
// problem.
// A feasible value for the variable must be chosen, and the variable must
// be shifted by this value. This is done to make sure that it will be
// possible to recreate a basic solution of the original problem from a
// basic solution of the pre-solved problem during post-solve.
const Fractional target_value =
MinInMagnitudeOrZeroIfInfinite(lower_bound, upper_bound);
Fractional lower_diff = (lower_bound - target_value) * bound_factor;
Fractional upper_diff = (upper_bound - target_value) * bound_factor;
if (bound_factor < 0.0) {
std::swap(lower_diff, upper_diff);
}
lp->SetVariableBounds(
target_col, lp->variable_lower_bounds()[target_col] + lower_diff,
lp->variable_upper_bounds()[target_col] + upper_diff);
SubtractColumnMultipleFromConstraintBound(col, target_value, lp);
column_deletion_helper_.MarkColumnForDeletionWithState(
col, target_value,
ComputeVariableStatus(target_value, lower_bound, upper_bound));
}
// If at least one column was merged, the target column must be shifted like
// the other columns in the same equivalence class for the same reason (see
// above).
if (num_merged > 0) {
merged_columns_[target_col] = target_col;
const Fractional target_value = MinInMagnitudeOrZeroIfInfinite(
lower_bounds_[target_col], upper_bounds_[target_col]);
lp->SetVariableBounds(
target_col, lp->variable_lower_bounds()[target_col] - target_value,
lp->variable_upper_bounds()[target_col] - target_value);
SubtractColumnMultipleFromConstraintBound(target_col, target_value, lp);
new_lower_bounds_[target_col] = lp->variable_lower_bounds()[target_col];
new_upper_bounds_[target_col] = lp->variable_upper_bounds()[target_col];
}
}
lp->DeleteColumns(column_deletion_helper_.GetMarkedColumns());
return !column_deletion_helper_.IsEmpty();
}
void ProportionalColumnPreprocessor::RecoverSolution(
ProblemSolution* solution) const {
SCOPED_INSTRUCTION_COUNT(time_limit_);
RETURN_IF_NULL(solution);
column_deletion_helper_.RestoreDeletedColumns(solution);
// The rest of this function is to unmerge the columns so that the solution be
// primal-feasible.
const ColIndex num_cols = merged_columns_.size();
DenseBooleanRow is_representative_basic(num_cols, false);
DenseBooleanRow is_distance_to_upper_bound(num_cols, false);
DenseRow distance_to_bound(num_cols, 0.0);
DenseRow wanted_value(num_cols, 0.0);
// The first pass is a loop over the representatives to compute the current
// distance to the new bounds.
for (ColIndex col(0); col < num_cols; ++col) {
if (merged_columns_[col] == col) {
const Fractional value = solution->primal_values[col];
const Fractional distance_to_upper_bound = new_upper_bounds_[col] - value;
const Fractional distance_to_lower_bound = value - new_lower_bounds_[col];
if (distance_to_upper_bound < distance_to_lower_bound) {
distance_to_bound[col] = distance_to_upper_bound;
is_distance_to_upper_bound[col] = true;
} else {
distance_to_bound[col] = distance_to_lower_bound;
is_distance_to_upper_bound[col] = false;
}
is_representative_basic[col] =
solution->variable_statuses[col] == VariableStatus::BASIC;
// Restore the representative value to a feasible value of the initial
// variable. Now all the merged variable are at a feasible value.
wanted_value[col] = value;
solution->primal_values[col] = MinInMagnitudeOrZeroIfInfinite(
lower_bounds_[col], upper_bounds_[col]);
solution->variable_statuses[col] = ComputeVariableStatus(
solution->primal_values[col], lower_bounds_[col], upper_bounds_[col]);
}
}
// Second pass to correct the values.
for (ColIndex col(0); col < num_cols; ++col) {
const ColIndex representative = merged_columns_[col];
if (representative != kInvalidCol) {
if (IsFinite(distance_to_bound[representative])) {
// If the distance is finite, then each variable is set to its
// corresponding bound (the one from which the distance is computed) and
// is then changed by as much as possible until the distance is zero.
const Fractional bound_factor =
column_factors_[col] / column_factors_[representative];
const Fractional scaled_distance =
distance_to_bound[representative] / std::abs(bound_factor);
const Fractional width = upper_bounds_[col] - lower_bounds_[col];
const bool to_upper_bound =
(bound_factor > 0.0) == is_distance_to_upper_bound[representative];
if (width <= scaled_distance) {
solution->primal_values[col] =
to_upper_bound ? lower_bounds_[col] : upper_bounds_[col];
solution->variable_statuses[col] =
ComputeVariableStatus(solution->primal_values[col],
lower_bounds_[col], upper_bounds_[col]);
distance_to_bound[representative] -= width * std::abs(bound_factor);
} else {
solution->primal_values[col] =
to_upper_bound ? upper_bounds_[col] - scaled_distance
: lower_bounds_[col] + scaled_distance;
solution->variable_statuses[col] =
is_representative_basic[representative]
? VariableStatus::BASIC
: ComputeVariableStatus(solution->primal_values[col],
lower_bounds_[col],
upper_bounds_[col]);
distance_to_bound[representative] = 0.0;
is_representative_basic[representative] = false;
}
} else {
// If the distance is not finite, then only one variable needs to be
// changed from its current feasible value in order to have a
// primal-feasible solution.
const Fractional error = wanted_value[representative];
if (error == 0.0) {
if (is_representative_basic[representative]) {
solution->variable_statuses[col] = VariableStatus::BASIC;
is_representative_basic[representative] = false;
}
} else {
const Fractional bound_factor =
column_factors_[col] / column_factors_[representative];
const bool use_this_variable =
(error * bound_factor > 0.0) ? (upper_bounds_[col] == kInfinity)
: (lower_bounds_[col] == -kInfinity);
if (use_this_variable) {
wanted_value[representative] = 0.0;
solution->primal_values[col] += error / bound_factor;
if (is_representative_basic[representative]) {
solution->variable_statuses[col] = VariableStatus::BASIC;
is_representative_basic[representative] = false;
} else {
// This should not happen on an OPTIMAL or FEASIBLE solution.
DCHECK(solution->status != ProblemStatus::OPTIMAL &&
solution->status != ProblemStatus::PRIMAL_FEASIBLE);
solution->variable_statuses[col] = VariableStatus::FREE;
}
}
}
}
}
}
}
// --------------------------------------------------------
// ProportionalRowPreprocessor
// --------------------------------------------------------
bool ProportionalRowPreprocessor::Run(LinearProgram* lp) {
SCOPED_INSTRUCTION_COUNT(time_limit_);
RETURN_VALUE_IF_NULL(lp, false);
const RowIndex num_rows = lp->num_constraints();
const SparseMatrix& transpose = lp->GetTransposeSparseMatrix();
// Use the first coefficient of each row to compute the proportionality
// factor. Note that the sign is important.
//
// Note(user): using the first coefficient may not give the best precision.
row_factors_.assign(num_rows, 0.0);
for (RowIndex row(0); row < num_rows; ++row) {
const SparseColumn& row_transpose = transpose.column(RowToColIndex(row));
if (!row_transpose.IsEmpty()) {
row_factors_[row] = row_transpose.GetFirstCoefficient();
}
}
// The new row bounds (only meaningful for the proportional rows).
DenseColumn lower_bounds(num_rows, -kInfinity);
DenseColumn upper_bounds(num_rows, +kInfinity);
// Where the new bounds are coming from. Only for the constraints that stay
// in the lp and are modified, kInvalidRow otherwise.
upper_bound_sources_.assign(num_rows, kInvalidRow);
lower_bound_sources_.assign(num_rows, kInvalidRow);
// Initialization.
// We need the first representative of each proportional row class to point to
// itself for the loop below. TODO(user): Already return such a mapping from
// FindProportionalColumns()?
ColMapping mapping = FindProportionalColumns(
transpose, parameters_.preprocessor_zero_tolerance());
DenseBooleanColumn is_a_representative(num_rows, false);
int num_proportional_rows = 0;
for (RowIndex row(0); row < num_rows; ++row) {
const ColIndex representative_row_as_col = mapping[RowToColIndex(row)];
if (representative_row_as_col != kInvalidCol) {
mapping[representative_row_as_col] = representative_row_as_col;
is_a_representative[ColToRowIndex(representative_row_as_col)] = true;
++num_proportional_rows;
}
}
// Compute the bound of each representative as implied by the rows
// which are proportional to it. Also keep the source row of each bound.
for (RowIndex row(0); row < num_rows; ++row) {
const ColIndex row_as_col = RowToColIndex(row);
if (mapping[row_as_col] != kInvalidCol) {
// For now, delete all the rows that are proportional to another one.
// Note that we will unmark the final representative of this class later.
row_deletion_helper_.MarkRowForDeletion(row);
const RowIndex representative_row = ColToRowIndex(mapping[row_as_col]);
const Fractional factor =
row_factors_[representative_row] / row_factors_[row];
Fractional implied_lb = factor * lp->constraint_lower_bounds()[row];
Fractional implied_ub = factor * lp->constraint_upper_bounds()[row];
if (factor < 0.0) {
std::swap(implied_lb, implied_ub);
}
// TODO(user): if the bounds are equal, use the largest row in magnitude?
if (implied_lb >= lower_bounds[representative_row]) {
lower_bounds[representative_row] = implied_lb;
lower_bound_sources_[representative_row] = row;
}
if (implied_ub <= upper_bounds[representative_row]) {
upper_bounds[representative_row] = implied_ub;
upper_bound_sources_[representative_row] = row;
}
}
}
// For maximum precision, and also to simplify the postsolve, we choose
// a representative for each class of proportional columns that has at least
// one of the two tightest bounds.
for (RowIndex row(0); row < num_rows; ++row) {
if (!is_a_representative[row]) continue;
const RowIndex lower_source = lower_bound_sources_[row];
const RowIndex upper_source = upper_bound_sources_[row];
lower_bound_sources_[row] = kInvalidRow;
upper_bound_sources_[row] = kInvalidRow;
DCHECK_NE(lower_source, kInvalidRow);
DCHECK_NE(upper_source, kInvalidRow);
if (lower_source == upper_source) {
// In this case, a simple change of representative is enough.
// The constraint bounds of the representative will not change.
DCHECK_NE(lower_source, kInvalidRow);
row_deletion_helper_.UnmarkRow(lower_source);
} else {
// Report ProblemStatus::PRIMAL_INFEASIBLE if the new lower bound is not
// lower than the new upper bound modulo the default tolerance.
if (!IsSmallerWithinFeasibilityTolerance(lower_bounds[row],
upper_bounds[row])) {
status_ = ProblemStatus::PRIMAL_INFEASIBLE;
return false;
}
// Special case for fixed rows.
if (lp->constraint_lower_bounds()[lower_source] ==
lp->constraint_upper_bounds()[lower_source]) {
row_deletion_helper_.UnmarkRow(lower_source);
continue;
}
if (lp->constraint_lower_bounds()[upper_source] ==
lp->constraint_upper_bounds()[upper_source]) {
row_deletion_helper_.UnmarkRow(upper_source);
continue;
}
// This is the only case where a more complex postsolve is needed.
// To maximize precision, the class representative is changed to either
// upper_source or lower_source depending of which row has the largest
// proportionality factor.
RowIndex new_representative = lower_source;
RowIndex other = upper_source;
if (std::abs(row_factors_[new_representative]) <
std::abs(row_factors_[other])) {
std::swap(new_representative, other);
}
// Initialize the new bounds with the implied ones.
const Fractional factor =
row_factors_[new_representative] / row_factors_[other];
Fractional new_lb = factor * lp->constraint_lower_bounds()[other];
Fractional new_ub = factor * lp->constraint_upper_bounds()[other];
if (factor < 0.0) {
std::swap(new_lb, new_ub);
}
lower_bound_sources_[new_representative] = new_representative;
upper_bound_sources_[new_representative] = new_representative;
if (new_lb > lp->constraint_lower_bounds()[new_representative]) {
lower_bound_sources_[new_representative] = other;
} else {
new_lb = lp->constraint_lower_bounds()[new_representative];
}
if (new_ub < lp->constraint_upper_bounds()[new_representative]) {
upper_bound_sources_[new_representative] = other;
} else {
new_ub = lp->constraint_upper_bounds()[new_representative];
}
const RowIndex new_lower_source =
lower_bound_sources_[new_representative];
if (new_lower_source == upper_bound_sources_[new_representative]) {
row_deletion_helper_.UnmarkRow(new_lower_source);
lower_bound_sources_[new_representative] = kInvalidRow;
upper_bound_sources_[new_representative] = kInvalidRow;
continue;
}
// Take care of small numerical imprecision by making sure that lb <= ub.
// Note that if the imprecision was greater than the tolerance, the code
// at the beginning of this block would have reported
// ProblemStatus::PRIMAL_INFEASIBLE.
DCHECK(IsSmallerWithinFeasibilityTolerance(new_lb, new_ub));
if (new_lb > new_ub) {
if (lower_bound_sources_[new_representative] == new_representative) {
new_ub = lp->constraint_lower_bounds()[new_representative];
} else {
new_lb = lp->constraint_upper_bounds()[new_representative];
}
}
row_deletion_helper_.UnmarkRow(new_representative);
lp->SetConstraintBounds(new_representative, new_lb, new_ub);
}
}
lp_is_maximization_problem_ = lp->IsMaximizationProblem();
lp->DeleteRows(row_deletion_helper_.GetMarkedRows());