1 // Copyright 2010-2021 Google LLC 2 // Licensed under the Apache License, Version 2.0 (the "License"); 3 // you may not use this file except in compliance with the License. 4 // You may obtain a copy of the License at 5 // 6 // http://www.apache.org/licenses/LICENSE-2.0 7 // 8 // Unless required by applicable law or agreed to in writing, software 9 // distributed under the License is distributed on an "AS IS" BASIS, 10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 11 // See the License for the specific language governing permissions and 12 // limitations under the License. 13 14 #ifndef OR_TOOLS_SAT_UTIL_H_ 15 #define OR_TOOLS_SAT_UTIL_H_ 16 17 #include <cstdint> 18 #include <deque> 19 20 #include "absl/random/bit_gen_ref.h" 21 #include "absl/random/random.h" 22 #include "ortools/sat/model.h" 23 #include "ortools/sat/sat_base.h" 24 #include "ortools/sat/sat_parameters.pb.h" 25 #include "ortools/util/random_engine.h" 26 #include "ortools/util/time_limit.h" 27 28 #if !defined(__PORTABLE_PLATFORM__) 29 #include "google/protobuf/descriptor.h" 30 #endif // __PORTABLE_PLATFORM__ 31 32 namespace operations_research { 33 namespace sat { 34 35 // Returns a in [0, m) such that a * x = 1 modulo m. 36 // If gcd(x, m) != 1, there is no inverse, and it returns 0. 37 // 38 // This DCHECK that x is in [0, m). 39 // This is integer overflow safe. 40 // 41 // Note(user): I didn't find this in a easily usable standard library. 42 int64_t ModularInverse(int64_t x, int64_t m); 43 44 // Just returns x % m but with a result always in [0, m). 45 int64_t PositiveMod(int64_t x, int64_t m); 46 47 // If we know that X * coeff % mod = rhs % mod, this returns c such that 48 // PositiveMod(X, mod) = c. 49 // 50 // This requires coeff != 0, mod !=0 and gcd(coeff, mod) == 1. 51 // The result will be in [0, mod) but there is no other condition on the sign or 52 // magnitude of a and b. 53 // 54 // This is overflow safe, and when rhs == 0 or abs(mod) == 1, it returns 0. 55 int64_t ProductWithModularInverse(int64_t coeff, int64_t mod, int64_t rhs); 56 57 // Returns true if the equation a * X + b * Y = cte has some integer solutions. 58 // For now, we check that a and b are different from 0 and from int64_t min. 59 // 60 // There is actually always a solution if cte % gcd(|a|, |b|) == 0. And because 61 // a, b and cte fit on an int64_t, if there is a solution, there is one with X 62 // and Y fitting on an int64_t. 63 // 64 // We will divide everything by gcd(a, b) first, so it is why we take reference 65 // and the equation can change. 66 // 67 // If there are solutions, we return one of them (x0, y0). 68 // From any such solution, the set of all solutions is given for Z integer by: 69 // X = x0 + b * Z; 70 // Y = y0 - a * Z; 71 // 72 // Given a domain for X and Y, it is possible to compute the "exact" domain of Z 73 // with our Domain functions. Note however that this will only compute solution 74 // where both x-x0 and y-y0 do fit on an int64_t: 75 // DomainOf(x).SubtractionWith(x0).InverseMultiplicationBy(b).IntersectionWith( 76 // DomainOf(y).SubtractionWith(y0).InverseMultiplicationBy(-a)) 77 bool SolveDiophantineEquationOfSizeTwo(int64_t& a, int64_t& b, int64_t& cte, 78 int64_t& x0, int64_t& y0); 79 80 // The model "singleton" random engine used in the solver. 81 // 82 // In test, we usually set use_absl_random() so that the sequence is changed at 83 // each invocation. This way, clients do not relly on the wrong assumption that 84 // a particular optimal solution will be returned if they are many equivalent 85 // ones. 86 class ModelRandomGenerator : public absl::BitGenRef { 87 public: 88 // We seed the strategy at creation only. This should be enough for our use 89 // case since the SatParameters is set first before the solver is created. We 90 // also never really need to change the seed afterwards, it is just used to 91 // diversify solves with identical parameters on different Model objects. ModelRandomGenerator(Model * model)92 explicit ModelRandomGenerator(Model* model) 93 : absl::BitGenRef(deterministic_random_) { 94 const auto& params = *model->GetOrCreate<SatParameters>(); 95 deterministic_random_.seed(params.random_seed()); 96 if (params.use_absl_random()) { 97 absl_random_ = absl::BitGen(absl::SeedSeq({params.random_seed()})); 98 absl::BitGenRef::operator=(absl::BitGenRef(absl_random_)); 99 } 100 } 101 102 // This is just used to display ABSL_RANDOM_SALT_OVERRIDE in the log so that 103 // it is possible to reproduce a failure more easily while looking at a solver 104 // log. 105 // 106 // TODO(user): I didn't find a cleaner way to log this. LogSalt()107 void LogSalt() const {} 108 109 private: 110 random_engine_t deterministic_random_; 111 absl::BitGen absl_random_; 112 }; 113 114 // The model "singleton" shared time limit. 115 class ModelSharedTimeLimit : public SharedTimeLimit { 116 public: ModelSharedTimeLimit(Model * model)117 explicit ModelSharedTimeLimit(Model* model) 118 : SharedTimeLimit(model->GetOrCreate<TimeLimit>()) {} 119 }; 120 121 // Randomizes the decision heuristic of the given SatParameters. 122 void RandomizeDecisionHeuristic(absl::BitGenRef random, 123 SatParameters* parameters); 124 125 // Context: this function is not really generic, but required to be unit-tested. 126 // It is used in a clause minimization algorithm when we try to detect if any of 127 // the clause literals can be propagated by a subset of the other literal being 128 // false. For that, we want to enqueue in the solver all the subset of size n-1. 129 // 130 // This moves one of the unprocessed literal from literals to the last position. 131 // The function tries to do that while preserving the longest possible prefix of 132 // literals "amortized" through the calls assuming that we want to move each 133 // literal to the last position once. 134 // 135 // For a vector of size n, if we want to call this n times so that each literal 136 // is last at least once, the sum of the size of the changed suffixes will be 137 // O(n log n). If we were to use a simpler algorithm (like moving the last 138 // unprocessed literal to the last position), this sum would be O(n^2). 139 // 140 // Returns the size of the common prefix of literals before and after the move, 141 // or -1 if all the literals are already processed. The argument 142 // relevant_prefix_size is used as a hint when keeping more that this prefix 143 // size do not matter. The returned value will always be lower or equal to 144 // relevant_prefix_size. 145 int MoveOneUnprocessedLiteralLast(const std::set<LiteralIndex>& processed, 146 int relevant_prefix_size, 147 std::vector<Literal>* literals); 148 149 // ============================================================================ 150 // Implementation. 151 // ============================================================================ 152 153 // Manages incremental averages. 154 class IncrementalAverage { 155 public: 156 // Initializes the average with 'initial_average' and number of records to 0. IncrementalAverage(double initial_average)157 explicit IncrementalAverage(double initial_average) 158 : average_(initial_average) {} IncrementalAverage()159 IncrementalAverage() {} 160 161 // Sets the number of records to 0 and average to 'reset_value'. 162 void Reset(double reset_value); 163 CurrentAverage()164 double CurrentAverage() const { return average_; } NumRecords()165 int64_t NumRecords() const { return num_records_; } 166 167 void AddData(double new_record); 168 169 private: 170 double average_ = 0.0; 171 int64_t num_records_ = 0; 172 }; 173 174 // Manages exponential moving averages defined as 175 // new_average = decaying_factor * old_average 176 // + (1 - decaying_factor) * new_record. 177 // where 0 < decaying_factor < 1. 178 class ExponentialMovingAverage { 179 public: ExponentialMovingAverage(double decaying_factor)180 explicit ExponentialMovingAverage(double decaying_factor) 181 : decaying_factor_(decaying_factor) { 182 DCHECK_GE(decaying_factor, 0.0); 183 DCHECK_LE(decaying_factor, 1.0); 184 } 185 186 // Returns exponential moving average for all the added data so far. CurrentAverage()187 double CurrentAverage() const { return average_; } 188 189 // Returns the total number of added records so far. NumRecords()190 int64_t NumRecords() const { return num_records_; } 191 192 void AddData(double new_record); 193 194 private: 195 double average_ = 0.0; 196 int64_t num_records_ = 0; 197 const double decaying_factor_; 198 }; 199 200 // Utility to calculate percentile (First variant) for limited number of 201 // records. Reference: https://en.wikipedia.org/wiki/Percentile 202 // 203 // After the vector is sorted, we assume that the element with index i 204 // correspond to the percentile 100*(i+0.5)/size. For percentiles before the 205 // first element (resp. after the last one) we return the first element (resp. 206 // the last). And otherwise we do a linear interpolation between the two element 207 // around the asked percentile. 208 class Percentile { 209 public: Percentile(int record_limit)210 explicit Percentile(int record_limit) : record_limit_(record_limit) {} 211 212 void AddRecord(double record); 213 214 // Returns number of stored records. NumRecords()215 int64_t NumRecords() const { return records_.size(); } 216 217 // Note that this is not fast and runs in O(n log n) for n records. 218 double GetPercentile(double percent); 219 220 private: 221 std::deque<double> records_; 222 const int record_limit_; 223 }; 224 225 // This method tries to compress a list of tuples by merging complementary 226 // tuples, that is a set of tuples that only differ on one variable, and that 227 // cover the domain of the variable. In that case, it will keep only one tuple, 228 // and replace the value for variable by any_value, the equivalent of '*' in 229 // regexps. 230 // 231 // This method is exposed for testing purposes. 232 void CompressTuples(absl::Span<const int64_t> domain_sizes, int64_t any_value, 233 std::vector<std::vector<int64_t>>* tuples); 234 235 } // namespace sat 236 } // namespace operations_research 237 238 #endif // OR_TOOLS_SAT_UTIL_H_ 239