1package solver23// a pbSet is a set representation of a PB constraint.4// It indicates, for each variable in the problem, what its weight is in the constraint.5// For instance, in a problem with 5 vars, the lits in constraint 5 x1 +3 ~x2 +2 x4+ x5 >= 6 will be encoded as6// []int{5, -3, 0, 2, 5}.7type pbSet struct {8 weights []int // The weight of each variable in the constraint, or 0 if the var isn't in the constraint.9 card int10}1112// pbSet converts c to the psSet structure.13// buffer is a buffer to store the weights. This is a parameter so as to avoid too frequent allocations.14func (s *Solver) pbSet(c *Clause, buffer []int) *pbSet {15 res := &pbSet{weights: buffer, card: c.Cardinality()}16 for i := range buffer { // Buffer has to be cleaned first17 buffer[i] = 018 }19 for i := 0; i < c.Len(); i++ {20 lit := c.Get(i)21 v := lit.Var()22 w := c.Weight(i)23 if !lit.IsPositive() {24 w = -w25 }26 res.weights[v] = w27 }28 return res29}3031func (pb *pbSet) clause() *Clause {32 lits := make([]Lit, 0, len(pb.weights))33 weights := make([]int, 0, len(pb.weights))34 for i, w := range pb.weights {35 if w == 0 {36 continue37 }38 idx := int32(i + 1)39 absW := w40 if w < 0 {41 absW = -w42 idx = -idx43 }44 lit := IntToLit(idx)45 lits = append(lits, lit)46 weights = append(weights, absW)47 }48 return NewPBClause(lits, weights, pb.card)49}5051// clash will make pb1 and pb2 clash, and update the values in pb1.52// pb2 will be unmodified.53// There should be at least one variable whose weight becomes 0 in the process.54func (pb1 *pbSet) clash(s *Solver, pb2 *pbSet) {55 pb1.card += pb2.card56 for i, w1 := range pb1.weights {57 w2 := pb2.weights[i]58 pb1.weights[i] += w259 if w1*w2 < 0 { // vars don't have the same polarity in both constraints60 pb1.card -= min(abs(w1), abs(w2))61 }62 }63}6465// slack returns the slack of pb1 for the given decision level.66// It is defined as card - sum(falsified lits at decLevel <= lvl).67// If slack = 0, all unassigned lits shuld be propagated.68// If slack < 0, the constraint is falsified under current assignments.69func (s *Solver) slack(pb *pbSet, lvl decLevel) int {70 res := -pb.card71 for i, w := range pb.weights {72 if w == 0 {73 continue74 }75 l := IntToLit(int32(i + 1))76 absW := w77 if w < 0 {78 absW = -w79 l = l.Negation()80 }81 v := Var(i)82 litLvl := abs(s.model[v])83 if s.litStatus(l) != Unsat || litLvl > lvl {84 res += absW85 }86 }87 return res88}8990// falsifies returns true iff lit's negation appears in pb.91// Only then should clashes happen.92func (pb *pbSet) falsifies(lit Lit) bool {93 v := lit.Var()94 w := pb.weights[v]95 if w == 0 {96 return false97 }98 return (w < 0) == lit.IsPositive()99}100101// cuttingPlanes learns a new PB constraint using the cutting plane resolution system.102// This is usually less efficient than calling learnClause, but will dramatically improve103// efficiency in corner cases, such as the pigeonhole problem.104func (s *Solver) cuttingPlanes(confl *Clause, lvl decLevel) (learned *Clause, propagated []Lit, newLvl decLevel) {105 // fmt.Printf("conflict at level %d! Constraint is: %s, trail is %s\n", lvl, confl.PBString(), s.trailString())106 seen := make([]bool, s.nbVars) // Was the var seen in the resolution process, making it a candidate for bumping?107 s.clauseBumpActivity(confl)108 for _, lit := range confl.lits {109 seen[lit.Var()] = true110 }111 pb := s.pbSet(confl, s.pbSetBuf)112 ptr := len(s.trail) - 1113 for pb.onlyFalsified(s, ptr, lvl) < 0 {114 if lvl == 1 { // Top-level conflict: UNSAT115 return nil, nil, -1116 }117 lit := s.trail[ptr]118 for !pb.falsifies(lit) {119 if s.reason[lit.Var()] == nil {120 lvl--121 }122 s.model[lit.Var()] = 0123 // s.trail = s.trail[:ptr]124 ptr--125 lit = s.trail[ptr]126 }127 v := lit.Var()128 s.varBumpActivity(v) // RoundingSAT's strategy: eliminated variables are bumped twice129 pb.roundToOne(s, v, lvl)130 reason := s.reason[v]131 if reason == nil {132 lvl--133 continue134 }135 for _, lit := range reason.lits {136 seen[lit.Var()] = true137 }138 s.clauseBumpActivity(reason)139 pb2 := s.pbSet(reason, s.pbSetBuf2)140 pb2.roundToOne(s, v, lvl)141 pb.clash(s, pb2)142 }143 unit := pb.onlyFalsified(s, ptr, lvl).Negation()144 btLvl := pb.backtrackLevel(s, unit)145 pb.roundToOne(s, unit.Var(), lvl)146 for i := range seen {147 if seen[i] {148 s.varBumpActivity(Var(i))149 }150 }151 // s.varDecayActivity()152 // s.clauseDecayActivity()153 if propagated, learned, ok := pb.clause().SimplifyPB(); !ok {154 return nil, nil, -1155 } else if len(propagated) > 0 {156 return nil, propagated, 1157 } else {158 return learned, []Lit{unit}, btLvl159 }160}161162func (pb *pbSet) backtrackLevel(s *Solver, falsified Lit) decLevel {163 v := falsified.Var()164 lvl := abs(s.model[v])165 maxLvl := decLevel(1)166 for i, w := range pb.weights {167 if w == 0 || Var(i) == v {168 continue169 }170 if lvlI := abs(s.model[i]); lvlI > maxLvl && lvlI != lvl {171 maxLvl = lvlI172 }173 }174 return maxLvl175}176177// returns the only false literal at level lvl in pb, or -1 if no or several lits are false.178func (pb *pbSet) onlyFalsified(s *Solver, ptr int, lvl decLevel) Lit {179 var res Lit = -1180 for ptr >= 0 {181 lit := s.trail[ptr]182 if abs(s.model[lit.Var()]) != lvl { // We're out of lvl now: we're done183 return res184 }185 if pb.falsifies(lit) {186 if res != -1 {187 // We already had one: there are several non-false literals!188 return -1189 }190 res = lit191 }192 ptr--193 }194 return res195}196197// roundToOne weakens pb by rounding the falsified literal ('locked'), as described in the RoundingSAT paper.198func (pb *pbSet) roundToOne(s *Solver, locked Var, lvl decLevel) {199 wi := abs(pb.weights[locked])200 if wi == 1 {201 return202 }203 for j, wj := range pb.weights {204 if wj == 0 {205 continue206 }207 assign := s.model[j]208 if wj%wi != 0 && (assign == 0 || ((assign > 0) == (wj > 0))) {209 // lit j isn't falsified: weaken constraint by removing it210 pb.weights[j] = 0211 pb.card -= abs(wj)212 }213 }214 pb.divideBy(wi)215}216217// divideBy performs a division on the conflict var by applying rounded division on its weight.218// locked should be present in pb or a panic will ensue.219// Ideally, this function shouldn't be called when coeff == 1.220func (pb *pbSet) divideBy(coeff int) {221 for j, wj := range pb.weights {222 if wj == 0 {223 continue224 }225 if wj%coeff == 0 {226 pb.weights[j] = wj / coeff227 } else if wj > 0 {228 pb.weights[j] = (wj / coeff) + 1229 } else {230 pb.weights[j] = (wj / coeff) - 1231 }232 }233 if pb.card%coeff == 0 {234 pb.card = pb.card / coeff235 } else {236 pb.card = (pb.card / coeff) + 1237 }238}