third_party/gophersat/solver is gophersat v1.4.0 (MIT) with bufLits moved from a package-level variable onto the Solver, so solvers run concurrently; the sat mutex is gone; go test -race ./sat passes. Measurement, idle G=60 scatter seed a: 32.8 s wall (baseline 34.4 s); time inside solve 32.7 s (was 167.7 s incl. lock wait); output identical to the baseline. Wall barely moved: the removal loops are CPU-bound on 108k solves each re-parsing the CNF, which step 3 addresses.
- commit
- 389b12137658cb53962f26fb8769c3f970cdfbe9
- parent
- 7a4c6ccdfe55395bf9a01b8cb8f79770d87a9ab0
- author
- dami <gazey55@gmail.com>
@@ -2,6 +2,4 @@ module mosaic-engine
22 33 go 1.24.544 5-require github.com/crillab/gophersat v1.4.06-75 replace github.com/FelineStateMachine/puzzletea => ../../puzzletea@@ -33,7 +33,7 @@ import (
3333 "sync/atomic"3434 "time"3535 36- "github.com/crillab/gophersat/solver"36+ "mosaic-engine/third_party/gophersat/solver"3737 3838 "mosaic-engine/rules"3939 )@@ -212,11 +212,6 @@ func keyFor(sh rules.Shape, mask rules.Mask) key {
212212 var (213213 cacheMu sync.Mutex214214 cache = map[key]*problem{}215- // gophersat keeps a package-level buffer for clause learning216- // (solver/learn.go: bufLits), so two solvers must never run at once.217- // Calls are milliseconds; serialising them costs less than the218- // generator's parallel waves would gain.219- solveMu sync.Mutex220215 )221216 222217 // call counters for the generator's timing table@@ -230,10 +225,10 @@ func Counters() (n int64, d time.Duration) {
230225 return calls.Load(), time.Duration(nanos.Load())231226 }232227 228+// solve runs one solver. The vendored gophersat keeps its clause-learning229+// buffer per solver (upstream: package-level), so solvers run concurrently.233230 func solve(s *solver.Solver) solver.Status {234231 t0 := time.Now()235- solveMu.Lock()236- defer solveMu.Unlock()237232 defer func() {238233 calls.Add(1)239234 nanos.Add(int64(time.Since(t0)))@@ -0,0 +1,7 @@
1+Copyright (c) 2017 Centre de Recherche en Informatique de Lens and contributors2+3+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:4+5+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.6+7+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.@@ -0,0 +1,28 @@
1+package solver2+3+// A CardConstr is a cardinality constraint, i.e a set of literals (represented with integer variables) associated with a minimal number of literals that must be true.4+// A propositional clause (i.e a disjunction of literals) is a cardinality constraint with a minimal cardinality of 1.5+type CardConstr struct {6+ Lits []int7+ AtLeast int8+}9+10+// AtLeast1 returns a cardinality constraint stating that at least one of the given lits must be true.11+// This is the equivalent of a propositional clause.12+func AtLeast1(lits ...int) CardConstr {13+ return CardConstr{Lits: lits, AtLeast: 1}14+}15+16+// AtMost1 returns a cardinality constraint stating that at most one of the given lits can be true.17+func AtMost1(lits ...int) CardConstr {18+ negated := make([]int, len(lits))19+ for i, lit := range lits {20+ negated[i] = -lit21+ }22+ return CardConstr{Lits: negated, AtLeast: len(lits) - 1}23+}24+25+// Exactly1 returns two cardinality constraints stating that exactly one of the given lits must be true.26+func Exactly1(lits ...int) []CardConstr {27+ return []CardConstr{AtLeast1(lits...), AtMost1(lits...)}28+}@@ -0,0 +1,280 @@
1+package solver2+3+import (4+ "fmt"5+ "sort"6+ "strings"7+)8+9+// data used in PB constraints.10+type pbData struct {11+ weights []int // weight of each literal. If nil, weights are all 1.12+ watched []bool // indices of watched literals.13+}14+15+// A Clause is a list of Lit, associated with possible data (for learned clauses).16+type Clause struct {17+ lits []Lit18+ // lbdValue's bits are as follow:19+ // leftmost bit: learned flag.20+ // second bit: locked flag (if learned).21+ // last 30 bits: LBD value (if learned) or minimal cardinality - 1 (if !learned).22+ // NOTE: actual cardinality is value + 1, since this is the default value and go defaults to 0.23+ lbdValue uint3224+ activity float3225+ pbData *pbData26+}27+28+const (29+ learnedMask uint32 = 1 << 3130+ lockedMask uint32 = 1 << 3031+ bothMasks uint32 = learnedMask | lockedMask32+)33+34+// NewClause returns a clause whose lits are given as an argument.35+func NewClause(lits []Lit) *Clause {36+ return &Clause{lits: lits}37+}38+39+// NewCardClause returns a clause whose lits are given as an argument and40+// for which at least 'card' literals must be true.41+// Note tha NewClause(lits) is equivalent to NewCardClause(lits, 1).42+func NewCardClause(lits []Lit, card int) *Clause {43+ if card < 1 || card > len(lits) {44+ panic("Invalid cardinality value")45+ }46+ return &Clause{lits: lits, lbdValue: uint32(card - 1)}47+}48+49+// Used to sort literals when constructing PB clause.50+type weightedLits struct {51+ lits []Lit52+ weights []int53+}54+55+func (wl weightedLits) Less(i, j int) bool { return wl.weights != nil && wl.weights[i] > wl.weights[j] }56+func (wl weightedLits) Len() int { return len(wl.lits) }57+func (wl *weightedLits) Swap(i, j int) {58+ wl.lits[i], wl.lits[j] = wl.lits[j], wl.lits[i]59+ if wl.weights != nil {60+ wl.weights[i], wl.weights[j] = wl.weights[j], wl.weights[i]61+ }62+}63+64+// NewPBClause returns a pseudo-boolean clause with the given lits, weights and minimal cardinality.65+func NewPBClause(lits []Lit, weights []int, card int) *Clause {66+ if card < 1 {67+ panic("Invalid cardinality value")68+ }69+ wl := &weightedLits{lits: lits, weights: weights}70+ sort.Sort(wl)71+ pbd := pbData{weights: weights, watched: make([]bool, len(lits))}72+ if pbd.weights == nil {73+ pbd.weights = make([]int, len(lits))74+ for i := range pbd.weights {75+ pbd.weights[i] = 176+ }77+ }78+ return &Clause{lits: lits, lbdValue: uint32(card - 1), pbData: &pbd}79+}80+81+// NewLearnedClause returns a new clause marked as learned.82+func NewLearnedClause(lits []Lit) *Clause {83+ return &Clause{lits: lits, lbdValue: learnedMask}84+}85+86+// Cardinality returns the minimum number of literals that must be true to satisfy the clause.87+func (c *Clause) Cardinality() int {88+ if c.Learned() {89+ return 190+ }91+ return int(c.lbdValue & ^bothMasks) + 192+}93+94+// Learned returns true iff c was a learned clause.95+func (c *Clause) Learned() bool {96+ return c.lbdValue&learnedMask == learnedMask97+}98+99+// PseudoBoolean returns true iff c is a pseudo boolean constraint, and not100+// just a propositional clause or cardinality constraint.101+func (c *Clause) PseudoBoolean() bool {102+ return c.pbData != nil103+}104+105+func (c *Clause) lock() {106+ c.lbdValue = c.lbdValue | lockedMask107+}108+109+func (c *Clause) unlock() {110+ c.lbdValue = c.lbdValue & ^lockedMask111+}112+113+func (c *Clause) lbd() int {114+ return int(c.lbdValue & ^bothMasks)115+}116+117+func (c *Clause) setLbd(lbd int) {118+ c.lbdValue = (c.lbdValue & bothMasks) | uint32(lbd)119+}120+121+func (c *Clause) incLbd() {122+ c.lbdValue++123+}124+125+func (c *Clause) isLocked() bool {126+ return c.lbdValue&bothMasks == bothMasks127+}128+129+// Len returns the nb of lits in the clause.130+func (c *Clause) Len() int {131+ return len(c.lits)132+}133+134+// First returns the first lit from the clause.135+func (c *Clause) First() Lit {136+ return c.lits[0]137+}138+139+// Second returns the second lit from the clause.140+func (c *Clause) Second() Lit {141+ return c.lits[1]142+}143+144+// Get returns the ith literal from the clause.145+func (c *Clause) Get(i int) Lit {146+ return c.lits[i]147+}148+149+// Set sets the ith literal of the clause.150+func (c *Clause) Set(i int, l Lit) {151+ c.lits[i] = l152+}153+154+// Weight returns the weight of the ith literal.155+// In a propositional clause or a cardinality constraint, that value will always be 1.156+func (c *Clause) Weight(i int) int {157+ if c.pbData == nil {158+ return 1159+ }160+ return c.pbData.weights[i]161+}162+163+// WeightSum returns the sum of the PB weights.164+// If c is a propositional clause, the function will return the length of the clause.165+func (c *Clause) WeightSum() int {166+ if c.pbData == nil {167+ return len(c.lits)168+ }169+ res := 0170+ for _, w := range c.pbData.weights {171+ res += w172+ }173+ return res174+}175+176+// swap swaps the ith and jth lits from the clause.177+func (c *Clause) swap(i, j int) {178+ c.lits[i], c.lits[j] = c.lits[j], c.lits[i]179+ if c.pbData != nil {180+ c.pbData.weights[i], c.pbData.weights[j] = c.pbData.weights[j], c.pbData.weights[i]181+ }182+}183+184+// updateCardinality adds "add" to c's cardinality.185+// Must not be called on learned clauses!186+func (c *Clause) updateCardinality(add int) {187+ if add < 0 && uint32(-add) > c.lbdValue {188+ c.lbdValue = 0189+ } else {190+ c.lbdValue += uint32(add)191+ }192+}193+194+// removeLit remove the idx'th lit from c.195+// The order of literals might be updated, so this function196+// should not be called once the whole solver was created,197+// only during simplification of the problem.198+func (c *Clause) removeLit(idx int) {199+ c.lits[idx] = c.lits[len(c.lits)-1]200+ c.lits = c.lits[:len(c.lits)-1]201+ if c.pbData != nil {202+ c.pbData.weights[idx] = c.pbData.weights[len(c.pbData.weights)-1]203+ c.pbData.weights = c.pbData.weights[:len(c.pbData.weights)-1]204+ }205+}206+207+// Shrink reduces the length of the clauses, by removing all lits208+// starting from position newLen.209+func (c *Clause) Shrink(newLen int) {210+ c.lits = c.lits[:newLen]211+ if c.pbData != nil {212+ c.pbData.weights = c.pbData.weights[:newLen]213+ c.pbData.watched = c.pbData.watched[:newLen]214+ }215+}216+217+// CNF returns a DIMACS CNF representation of the clause.218+func (c *Clause) CNF() string {219+ res := ""220+ for _, lit := range c.lits {221+ res += fmt.Sprintf("%d ", lit.Int())222+ }223+ return fmt.Sprintf("%s0", res)224+}225+226+// PBString returns a string representation of c as a pseudo-boolean expression.227+func (c *Clause) PBString() string {228+ terms := make([]string, c.Len())229+ for i, lit := range c.lits {230+ weight := 1231+ if c.pbData != nil {232+ weight = c.pbData.weights[i]233+ }234+ val := lit.Int()235+ sign := ""236+ if val < 0 {237+ val = -val238+ sign = "~"239+ }240+ terms[i] = fmt.Sprintf("%d %sx%d", weight, sign, val)241+ }242+ return fmt.Sprintf("%s >= %d ;", strings.Join(terms, " +"), c.Cardinality())243+}244+245+// SimplifyPB tries to simplify a pseudo boolean constraint by propagating all lits that can be propagated246+// not matter whet the assignment.247+// A PB constraint can also be false if all its coeficients are smaller than the cardinality.248+// It will return a list of unit lits (if any), a new, simplified clause (or nil if all lits can be propagated and/or ignored)249+// and a boolean ok indicating whether the clause can be satisfied or not. If ok is false, other return parameters are irrelevant.250+func (c *Clause) SimplifyPB() (units []Lit, c2 *Clause, ok bool) {251+ card := c.Cardinality()252+ thresh := c.WeightSum() - card253+ if thresh < 0 {254+ // Even if all lits are true, the clause is still unsatisfied255+ return nil, nil, false256+ }257+ // all lits whose weights are > thresh must be satisfied258+ i := 0259+ for i < c.Len() && c.Weight(i) > thresh {260+ units = append(units, c.Get(i))261+ card -= c.Weight(i)262+ i++263+ }264+ if card <= 0 {265+ // clause is now satisfied, other lits are irrelevant266+ return units, nil, true267+ }268+ // All lits starting from i must be kept269+ newLen := c.Len() - i270+ newLits := make([]Lit, newLen)271+ newWeights := make([]int, newLen)272+ copy(newLits, c.lits[i:])273+ copy(newWeights, c.pbData.weights[i:])274+ // saturate weights so that none is higher than card275+ i = 0276+ for newWeights[i] > card {277+ newWeights[i] = card278+ }279+ return units, NewPBClause(newLits, newWeights, card), true280+}@@ -0,0 +1,101 @@
1+/*2+Package solver gives access to a simple SAT and pseudo-boolean solver.3+Its input can be either a DIMACS CNF file, a PBS file or a solver.Problem object,4+containing the set of clauses to be solved. In the last case,5+the problem can be either a set of propositional clauses,6+or a set of pseudo-boolean constraints.7+8+No matter the input format,9+the solver.Solver will then solve the problem and indicate whether the problem is10+satisfiable or not. In the former case, it will be able to provide a model, i.e a set of bindings11+for all variables that makes the problem true.12+13+# Describing a problem14+15+A problem can be described in several ways:16+17+1. parse a DIMACS stream (io.Reader). If the io.Reader produces the following content:18+19+ p cnf 6 720+ 1 2 3 021+ 4 5 6 022+ -1 -4 023+ -2 -5 024+ -3 -6 025+ -1 -3 026+ -4 -6 027+28+the programmer can create the Problem by doing:29+30+ pb, err := solver.ParseCNF(f)31+32+2. create the equivalent list of list of literals. The problem above can be created programatically this way:33+34+ clauses := [][]int{35+ []int{1, 2, 3},36+ []int{4, 5, 6},37+ []int{-1, -4},38+ []int{-2, -5},39+ []int{-3, -6},40+ []int{-1, -3},41+ []int{-4, -6},42+ }43+ pb := solver.ParseSlice(clauses)44+45+3. create a list of cardinality constraints (CardConstr), if the problem to be solved is better represented this way.46+For instance, the problem stating that at least two literals must be true among the literals 1, 2, 3 and 4 could be described as a set of clauses:47+48+ clauses := [][]int{49+ []int{1, 2, 3},50+ []int{2, 3, 4},51+ []int{1, 2, 4},52+ []int{1, 3, 4},53+ }54+ pb := solver.ParseSlice(clauses)55+56+The number of clauses necessary to describe such a constrain can grow exponentially. Alternatively, it is possible to describe the same this way:57+58+ constrs := []CardConstr{59+ {Lits: []int{1, 2, 3, 4}, AtLeast: 2},60+ }61+ pb := solver.ParseCardConstrs(clauses)62+63+Note that a propositional clause has an implicit cardinality constraint of 1, since at least one of its literals must be true.64+65+4. parse a PBS stream (io.Reader). If the io.Reader contains the following problem:66+67+ 2 ~x1 +1 x2 +1 x3 >= 3 ;68+69+the programmer can create the Problem by doing:70+71+ pb, err := solver.ParsePBS(f)72+73+5. create a list of PBConstr. For instance, the following set of one PBConstrs will generate the same problem as above:74+75+ constrs := []PBConstr{GtEq([]int{1, 2, 3}, []int{2, 1, 1}, 3)}76+ pb := solver.ParsePBConstrs(constrs)77+78+# Solving a problem79+80+To solve a problem, one simply creates a solver with said problem.81+The solve() method then solves the problem and returns the corresponding status: Sat or Unsat.82+83+ s := solver.New(pb)84+ status := s.Solve()85+86+If the status was Sat, the programmer can ask for a model, i.e an assignment that makes all the clauses of the problem true:87+88+ m := s.Model()89+90+For the above problem, the status will be Sat and the model can be {false, true, false, true, false, false}.91+92+Alternatively, one can display on stdout the result and model (if any):93+94+ s.OutputModel()95+96+For the above problem described in the DIMACS format, the output can be:97+98+ SATISFIABLE99+ -1 2 -3 4 -5 -6100+*/101+package solver@@ -0,0 +1,35 @@
1+package solver2+3+// A Result is a status, either Sat, Unsat or Indet.4+// If the status is Sat, the Result also associates a ModelMap with an integer value.5+// This value is typically used in optimization processes.6+// If the weight is 0, that means all constraints could be solved.7+// By definition, in decision problems, the cost will always be 0.8+type Result struct {9+ Status Status10+ Model []bool11+ Weight int12+}13+14+// Interface is any type implementing a solver.15+// The basic Solver defined in this package implements it.16+// Any solver that uses the basic solver to solve more complex problems17+// (MAXSAT, MUS extraction, etc.) can implement it, too.18+type Interface interface {19+ // Optimal solves or optimizes the problem and returns the best result.20+ // If the results chan is non nil, it will write the associated model each time one is found.21+ // It will stop as soon as a model of cost 0 is found, or the problem is not satisfiable anymore.22+ // The last satisfying model, if any, will be returned with the Sat status.23+ // If no model at all could be found, the Unsat status will be returned.24+ // If the solver prematurely stopped, the Indet status will be returned.25+ // If data is sent to stop, the method may stop prematurely.26+ // In any case, results will be closed before the function returns.27+ // NOTE: data sent on stop may be ignored by an implementation.28+ Optimal(results chan Result, stop chan struct{}) Result29+ // Enumerate returns the number of models for the problem.30+ // If the models chan is non nil, it will write the associated model each time one is found.31+ // If data is sent to stop, the method may stop prematurely.32+ // In any case, models will be closed before the function returns.33+ // NOTE: data sent on stop may be ignored by an implementation.34+ Enumerate(models chan []bool, stop chan struct{}) int35+}@@ -0,0 +1,89 @@
1+package solver2+3+const (4+ nbMaxRecent = 50 // How many recent LBD values we consider; "X" in papers about LBD.5+ triggerRestartK = 0.86+ nbMaxTrail = 5_000 // How many elements in queueTrail we consider; "Y" in papers about LBD.7+ postponeRestartT = 1.48+)9+10+type queueData struct {11+ totalNb int // Current total nb of values considered12+ totalSum int // Sum of all values so far13+ nbRecent int // NB of values used in the array14+ ptr int // current index of oldest value in the array15+ recentAvg float64 // Average value16+}17+18+// lbdStats is a structure dealing with recent LBD evolutions.19+type lbdStats struct {20+ lbdData queueData21+ trailData queueData22+ recentVals [nbMaxRecent]int // Last LBD values23+ recentTrails [nbMaxTrail]int // Last trail lengths24+}25+26+// mustRestart is true iff recent LBDs are much smaller on average than average of all LBDs.27+func (l *lbdStats) mustRestart() bool {28+ if l.lbdData.nbRecent < nbMaxRecent {29+ return false30+ }31+ return l.lbdData.recentAvg*triggerRestartK > float64(l.lbdData.totalSum)/float64(l.lbdData.totalNb)32+}33+34+// addConflict adds information about a conflict that just happened.35+func (l *lbdStats) addConflict(trailSz int) {36+ td := &l.trailData37+ td.totalNb++38+ td.totalSum += trailSz39+ if td.nbRecent < nbMaxTrail {40+ l.recentTrails[td.nbRecent] = trailSz41+ old := float64(td.nbRecent)42+ new := old + 143+ td.recentAvg = (td.recentAvg*old)/new + float64(trailSz)/new44+ td.nbRecent++45+ } else {46+ old := l.recentTrails[td.ptr]47+ l.recentTrails[td.ptr] = trailSz48+ td.ptr++49+ if td.ptr == nbMaxTrail {50+ td.ptr = 051+ }52+ td.recentAvg = td.recentAvg - float64(old)/nbMaxTrail + float64(trailSz)/nbMaxTrail53+ }54+ if td.nbRecent == nbMaxTrail && l.lbdData.nbRecent == nbMaxRecent && trailSz > int(postponeRestartT*td.recentAvg) {55+ // Too many good assignments: postpone restart56+ l.clear()57+ }58+59+}60+61+// addLbd adds information about a recent learned clause's LBD.62+// TODO: this is very close to addConflicts's code, this should probably be rewritten/merged.63+func (l *lbdStats) addLbd(lbd int) {64+ ld := &l.lbdData65+ ld.totalNb++66+ ld.totalSum += lbd67+ if ld.nbRecent < nbMaxRecent {68+ l.recentVals[ld.nbRecent] = lbd69+ old := float64(ld.nbRecent)70+ new := old + 171+ ld.recentAvg = (ld.recentAvg*old)/new + float64(lbd)/new72+ ld.nbRecent++73+ } else {74+ old := l.recentVals[ld.ptr]75+ l.recentVals[ld.ptr] = lbd76+ ld.ptr++77+ if ld.ptr == nbMaxRecent {78+ ld.ptr = 079+ }80+ ld.recentAvg = ld.recentAvg - float64(old)/nbMaxRecent + float64(lbd)/nbMaxRecent81+ }82+}83+84+// clear clears last values. It should be called after a restart.85+func (l *lbdStats) clear() {86+ l.lbdData.ptr = 087+ l.lbdData.nbRecent = 088+ l.lbdData.recentAvg = 0.089+}@@ -0,0 +1,138 @@
1+package solver2+3+// computeLbd computes and sets c's LBD (Literal Block Distance).4+func (c *Clause) computeLbd(model Model) {5+ c.setLbd(1)6+ curLvl := abs(model[c.Get(0).Var()])7+ for i := 0; i < c.Len(); i++ {8+ lit := c.Get(i)9+ if lvl := abs(model[lit.Var()]); lvl != curLvl {10+ curLvl = lvl11+ c.incLbd()12+ }13+ }14+}15+16+// addClauseLits is a helper function for learnClause.17+// It deals with lits from the conflict clause.18+func (s *Solver) addClauseLits(confl *Clause, lvl decLevel, met, metLvl []bool, lits *[]Lit) int {19+ nbLvl := 020+ for i := 0; i < confl.Len(); i++ {21+ l := confl.Get(i)22+ v := l.Var()23+ if s.litStatus(l) != Unsat {24+ // In clauses where cardinality > 1, some lits might be true in the conflict clause: ignore them25+ continue26+ }27+ met[v] = true28+ s.varBumpActivity(v)29+ if abs(s.model[v]) == lvl {30+ metLvl[v] = true31+ nbLvl++32+ } else {33+ *lits = append(*lits, l)34+ }35+ }36+ return nbLvl37+}38+39+// bufLits (per Solver) is the buffer for lits in learnClause, kept to40+// reduce allocations; a package-level buffer made concurrent solvers41+// unsafe (mosaic: moved onto the Solver).42+43+// learnClause creates a conflict clause and returns either:44+// - the clause itself, if its len is at least 2,45+// - a nil clause and a unit literal, if its len is exactly 1,46+// - a nil clause and -1, if the empty clause was learned.47+func (s *Solver) learnClause(confl *Clause, lvl decLevel) (learned *Clause, unit Lit) {48+ s.clauseBumpActivity(confl)49+ if s.bufLits == nil {50+ s.bufLits = make([]Lit, 10_000)51+ }52+ lits := s.bufLits[:1] // Not 0: make room for asserting literal53+ buf := make([]bool, s.nbVars*2) // Buffer for met and metLvl; reduces allocs/deallocs54+ met := buf[:s.nbVars] // List of all vars already met55+ metLvl := buf[s.nbVars:] // List of all vars from current level to deal with56+ // nbLvl is the nb of vars in lvl currently used57+ nbLvl := s.addClauseLits(confl, lvl, met, metLvl, &lits)58+ ptr := len(s.trail) - 1 // Pointer in propagation trail59+ for nbLvl > 1 { // We will stop once we only have one lit from current level.60+ for !metLvl[s.trail[ptr].Var()] {61+ if abs(s.model[s.trail[ptr].Var()]) == lvl { // This var was deduced afterwards and was not a reason for the conflict62+ met[s.trail[ptr].Var()] = true63+ }64+ ptr--65+ }66+ v := s.trail[ptr].Var()67+ if s.assumptions[v] {68+ // Now we only have assumed lits: this is a top-level conflict69+ return nil, -170+ }71+ ptr--72+ nbLvl--73+ if reason := s.reason[v]; reason != nil {74+ s.clauseBumpActivity(reason)75+ for i := 0; i < reason.Len(); i++ {76+ lit := reason.Get(i)77+ if v2 := lit.Var(); !met[v2] {78+ if s.litStatus(lit) != Unsat { // In clauses where cardinality > 1, some lits might be true in the conflict clause: ignore them79+ continue80+ }81+ met[v2] = true82+ s.varBumpActivity(v2)83+ if abs(s.model[v2]) == lvl {84+ metLvl[v2] = true85+ nbLvl++86+ } else {87+ lits = append(lits, lit)88+ }89+ }90+ }91+ }92+ }93+ for _, l := range s.trail { // Look for last lit from lvl and use it as asserting lit94+ if metLvl[l.Var()] {95+ lits[0] = l.Negation()96+ break97+ }98+ }99+ s.varDecayActivity()100+ s.clauseDecayActivity()101+ sortLiterals(lits, s.model)102+ sz := s.minimizeLearned(met, lits)103+ if sz == 1 {104+ // fmt.Printf("learned unit %d, trail is %s\n", lits[0].Int(), s.trailString())105+ return nil, lits[0]106+ }107+ if cap(lits) > cap(s.bufLits) {108+ s.bufLits = lits[:cap(lits)] // the buffer grew: keep it109+ }110+ lits2 := make([]Lit, sz)111+ copy(lits2, lits[:sz])112+ learned = NewLearnedClause(lits2)113+ learned.computeLbd(s.model)114+ // fmt.Printf("learned clause %s, trail is %s\n", learned.CNF(), s.trailString())115+ return learned, -1116+}117+118+// minimizeLearned reduces (if possible) the length of the learned clause and returns the size119+// of the new list of lits.120+func (s *Solver) minimizeLearned(met []bool, learned []Lit) int {121+ sz := 1122+ for i := 1; i < len(learned); i++ {123+ if reason := s.reason[learned[i].Var()]; reason == nil {124+ learned[sz] = learned[i]125+ sz++126+ } else {127+ for k := 0; k < reason.Len(); k++ {128+ lit := reason.Get(k)129+ if !met[lit.Var()] /*&& abs(s.model[lit.Var()]) > 1*/ {130+ learned[sz] = learned[i]131+ sz++132+ break133+ }134+ }135+ }136+ }137+ return sz138+}@@ -0,0 +1,238 @@
1+package solver2+3+// 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}.7+type 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+}11+12+// 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.14+func (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+}30+31+func (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+}50+51+// 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.54+func (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+}64+65+// 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.69+func (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+}89+90+// falsifies returns true iff lit's negation appears in pb.91+// Only then should clashes happen.92+func (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+}100+101+// 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.104+func (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+}161+162+func (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+}176+177+// returns the only false literal at level lvl in pb, or -1 if no or several lits are false.178+func (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+}196+197+// roundToOne weakens pb by rounding the falsified literal ('locked'), as described in the RoundingSAT paper.198+func (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+}216+217+// 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.220+func (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+}@@ -0,0 +1,188 @@
1+package solver2+3+import (4+ "bufio"5+ "fmt"6+ "io"7+ "strconv"8+ "strings"9+)10+11+// ParseSlice parse a slice of slice of lits and returns the equivalent problem.12+// The argument is supposed to be a well-formed CNF.13+func ParseSlice(cnf [][]int) *Problem {14+ var pb Problem15+ pb.parseSlice(cnf)16+ return &pb17+}18+19+// ParseSliceNb parse a slice of slice of lits and returns the equivalent problem.20+// The argument is supposed to be a well-formed CNF.21+// The number of vars is provided because clauses might be added to it later.22+func ParseSliceNb(cnf [][]int, nbVars int) *Problem {23+ pb := Problem{NbVars: nbVars}24+ pb.parseSlice(cnf)25+ return &pb26+}27+28+func (pb *Problem) parseSlice(cnf [][]int) {29+ for _, line := range cnf {30+ switch len(line) {31+ case 0:32+ pb.Status = Unsat33+ return34+ case 1:35+ if line[0] == 0 {36+ panic("null unit clause")37+ }38+ lit := IntToLit(int32(line[0]))39+ v := lit.Var()40+ if int(v) >= pb.NbVars {41+ pb.NbVars = int(v) + 142+ }43+ pb.Units = append(pb.Units, lit)44+ default:45+ lits := make([]Lit, len(line))46+ for j, val := range line {47+ if val == 0 {48+ panic(fmt.Sprintf("null literal in clause %v", lits))49+ }50+ lits[j] = IntToLit(int32(val))51+ if v := int(lits[j].Var()); v >= pb.NbVars {52+ pb.NbVars = v + 153+ }54+ }55+ pb.Clauses = append(pb.Clauses, NewClause(lits))56+ }57+ }58+ pb.Model = make([]decLevel, pb.NbVars)59+ for _, unit := range pb.Units {60+ v := unit.Var()61+ if pb.Model[v] == 0 {62+ if unit.IsPositive() {63+ pb.Model[v] = 164+ } else {65+ pb.Model[v] = -166+ }67+ } else if pb.Model[v] > 0 != unit.IsPositive() {68+ pb.Status = Unsat69+ return70+ }71+ }72+ pb.simplify2()73+}74+75+func isSpace(b byte) bool {76+ return b == ' ' || b == '\t' || b == '\n' || b == '\r'77+}78+79+// readInt reads an int from r.80+// 'b' is the last read byte. It can be a space, a '-' or a digit.81+// The int can be negated.82+// All spaces before the int value are ignored.83+// Can return EOF.84+func readInt(b *byte, r *bufio.Reader) (res int, err error) {85+ for err == nil && isSpace(*b) {86+ *b, err = r.ReadByte()87+ }88+ if err == io.EOF {89+ return res, io.EOF90+ }91+ if err != nil {92+ return res, fmt.Errorf("could not read digit: %v", err)93+ }94+ neg := 195+ if *b == '-' {96+ neg = -197+ *b, err = r.ReadByte()98+ if err != nil {99+ return 0, fmt.Errorf("cannot read int: %v", err)100+ }101+ }102+ for err == nil {103+ if *b < '0' || *b > '9' {104+ return 0, fmt.Errorf("cannot read int: %q is not a digit", *b)105+ }106+ res = 10*res + int(*b-'0')107+ *b, err = r.ReadByte()108+ if isSpace(*b) {109+ break110+ }111+ }112+ res *= neg113+ return res, err114+}115+116+func parseHeader(r *bufio.Reader) (nbVars, nbClauses int, err error) {117+ line, err := r.ReadString('\n')118+ if err != nil {119+ return 0, 0, fmt.Errorf("cannot read header: %v", err)120+ }121+ fields := strings.Fields(line)122+ if len(fields) < 3 {123+ return 0, 0, fmt.Errorf("invalid syntax %q in header", line)124+ }125+ nbVars, err = strconv.Atoi(fields[1])126+ if err != nil {127+ return 0, 0, fmt.Errorf("nbvars not an int : %q", fields[1])128+ }129+ nbClauses, err = strconv.Atoi(fields[2])130+ if err != nil {131+ return 0, 0, fmt.Errorf("nbClauses not an int : '%s'", fields[2])132+ }133+ return nbVars, nbClauses, nil134+}135+136+// ParseCNF parses a CNF file and returns the corresponding Problem.137+func ParseCNF(f io.Reader) (*Problem, error) {138+ r := bufio.NewReader(f)139+ var (140+ nbClauses int141+ pb Problem142+ )143+ b, err := r.ReadByte()144+ for err == nil {145+ if b == 'c' { // Ignore comment146+ b, err = r.ReadByte()147+ for err == nil && b != '\n' {148+ b, err = r.ReadByte()149+ }150+ } else if b == 'p' { // Parse header151+ pb.NbVars, nbClauses, err = parseHeader(r)152+ if err != nil {153+ return nil, fmt.Errorf("cannot parse CNF header: %v", err)154+ }155+ pb.Model = make([]decLevel, pb.NbVars)156+ pb.Clauses = make([]*Clause, 0, nbClauses)157+ } else {158+ lits := make([]Lit, 0, 3) // Make room for some lits to improve performance159+ for {160+ val, err := readInt(&b, r)161+ if err == io.EOF {162+ if len(lits) != 0 { // This is not a trailing space at the end...163+ pb.Clauses = append(pb.Clauses, NewClause(lits))164+ }165+ break // When there are only several useless spaces at the end of the file, that is ok166+ }167+ if err != nil {168+ return nil, fmt.Errorf("cannot parse clause: %v", err)169+ }170+ if val == 0 {171+ pb.Clauses = append(pb.Clauses, NewClause(lits))172+ break173+ } else {174+ if val > pb.NbVars || -val > pb.NbVars {175+ return nil, fmt.Errorf("invalid literal %d for problem with %d vars only", val, pb.NbVars)176+ }177+ lits = append(lits, IntToLit(int32(val)))178+ }179+ }180+ }181+ b, err = r.ReadByte()182+ }183+ if err != io.EOF {184+ return nil, err185+ }186+ pb.simplify2()187+ return &pb, nil188+}@@ -0,0 +1,266 @@
1+package solver2+3+import (4+ "bufio"5+ "fmt"6+ "io"7+ "strconv"8+ "strings"9+)10+11+// ParseCardConstrs parses the given cardinality constraints.12+// Will panic if a zero value appears in the literals.13+func ParseCardConstrs(constrs []CardConstr) *Problem {14+ var pb Problem15+ for _, constr := range constrs {16+ card := constr.AtLeast17+ if card <= 0 { // Clause is trivially SAT, ignore18+ continue19+ }20+ if len(constr.Lits) < card { // Clause cannot be satsfied21+ pb.Status = Unsat22+ return &pb23+ }24+ if len(constr.Lits) == card { // All lits must be true25+ for i := range constr.Lits {26+ if constr.Lits[i] == 0 {27+ panic("literal 0 found in clause")28+ }29+ lit := IntToLit(int32(constr.Lits[i]))30+ v := lit.Var()31+ if int(v) >= pb.NbVars {32+ pb.NbVars = int(v) + 133+ }34+ pb.Units = append(pb.Units, lit)35+ }36+ } else {37+ lits := make([]Lit, len(constr.Lits))38+ for j, val := range constr.Lits {39+ if val == 0 {40+ panic("literal 0 found in clause")41+ }42+ lits[j] = IntToLit(int32(val))43+ if v := int(lits[j].Var()); v >= pb.NbVars {44+ pb.NbVars = v + 145+ }46+ }47+ pb.Clauses = append(pb.Clauses, NewCardClause(lits, card))48+ }49+ }50+ pb.Model = make([]decLevel, pb.NbVars)51+ for _, unit := range pb.Units {52+ v := unit.Var()53+ if pb.Model[v] == 0 {54+ if unit.IsPositive() {55+ pb.Model[v] = 156+ } else {57+ pb.Model[v] = -158+ }59+ } else if pb.Model[v] > 0 != unit.IsPositive() {60+ pb.Status = Unsat61+ return &pb62+ }63+ }64+ pb.simplifyCard()65+ return &pb66+}67+68+func (pb *Problem) appendClause(constr PBConstr) {69+ lits := make([]Lit, len(constr.Lits))70+ for j, val := range constr.Lits {71+ lits[j] = IntToLit(int32(val))72+ }73+ pb.Clauses = append(pb.Clauses, NewPBClause(lits, constr.Weights, constr.AtLeast))74+}75+76+// ParsePBConstrs parses and returns a PB problem from PBConstr values.77+func ParsePBConstrs(constrs []PBConstr) *Problem {78+ var pb Problem79+ for _, constr := range constrs {80+ for i := range constr.Lits {81+ lit := IntToLit(int32(constr.Lits[i]))82+ v := lit.Var()83+ if int(v) >= pb.NbVars {84+ pb.NbVars = int(v) + 185+ }86+ }87+ card := constr.AtLeast88+ if card <= 0 { // Clause is trivially SAT, ignore89+ continue90+ }91+ sumW := constr.WeightSum()92+ if sumW < card { // Clause cannot be satsfied93+ pb.Status = Unsat94+ return &pb95+ }96+ if sumW == card { // All lits must be true97+ for i := range constr.Lits {98+ lit := IntToLit(int32(constr.Lits[i]))99+ found := false100+ for _, u := range pb.Units {101+ if u == lit {102+ found = true103+ break104+ }105+ }106+ if !found {107+ pb.Units = append(pb.Units, lit)108+ }109+ }110+ } else {111+ pb.appendClause(constr)112+ }113+ }114+ pb.Model = make([]decLevel, pb.NbVars)115+ for _, unit := range pb.Units {116+ v := unit.Var()117+ if pb.Model[v] == 0 {118+ if unit.IsPositive() {119+ pb.Model[v] = 1120+ } else {121+ pb.Model[v] = -1122+ }123+ } else if pb.Model[v] > 0 != unit.IsPositive() {124+ pb.Status = Unsat125+ return &pb126+ }127+ }128+ pb.simplifyPB()129+ return &pb130+}131+132+// parsePBOptim parses the "min:" instruction.133+func (pb *Problem) parsePBOptim(fields []string, line string) error {134+ weights, lits, err := pb.parseTerms(fields[1:], line)135+ if err != nil {136+ return err137+ }138+ pb.minLits = make([]Lit, len(lits))139+ for i, lit := range lits {140+ pb.minLits[i] = IntToLit(int32(lit))141+ }142+ pb.minWeights = weights143+ return nil144+}145+146+func (pb *Problem) parsePBLine(line string) error {147+ if line[len(line)-1] != ';' {148+ return fmt.Errorf("line %q does not end with semicolon", line)149+ }150+ fields := strings.Fields(line[:len(line)-1])151+ if len(fields) == 0 {152+ return fmt.Errorf("empty line in file")153+ }154+ if fields[0] == "min:" { // Optimization constraint155+ return pb.parsePBOptim(fields, line)156+ }157+ return pb.parsePBConstrLine(fields, line)158+}159+160+func (pb *Problem) parsePBConstrLine(fields []string, line string) error {161+ if len(fields) < 3 {162+ return fmt.Errorf("invalid syntax %q", line)163+ }164+ operator := fields[len(fields)-2]165+ if operator != ">=" && operator != "=" {166+ return fmt.Errorf("invalid operator %q in %q: expected \">=\" or \"=\"", operator, line)167+ }168+ rhs, err := strconv.Atoi(fields[len(fields)-1])169+ if err != nil {170+ return fmt.Errorf("invalid value %q in %q: %v", fields[len(fields)-1], line, err)171+ }172+ weights, lits, err := pb.parseTerms(fields[:len(fields)-2], line)173+ if err != nil {174+ return err175+ }176+ var constrs []PBConstr177+ if operator == ">=" {178+ constrs = []PBConstr{GtEq(lits, weights, rhs)}179+ } else {180+ constrs = Eq(lits, weights, rhs)181+ }182+ for _, constr := range constrs {183+ card := constr.AtLeast184+ sumW := constr.WeightSum()185+ if sumW < card { // Clause cannot be satsfied186+ pb.Status = Unsat187+ return nil188+ }189+ if sumW == card { // All lits must be true190+ for i := range constr.Lits {191+ lit := IntToLit(int32(constr.Lits[i]))192+ pb.Units = append(pb.Units, lit)193+ }194+ } else {195+ lits := make([]Lit, len(constr.Lits))196+ for j, val := range constr.Lits {197+ lits[j] = IntToLit(int32(val))198+ }199+ pb.Clauses = append(pb.Clauses, NewPBClause(lits, constr.Weights, card))200+ }201+ }202+ return nil203+}204+205+func (pb *Problem) parseTerms(terms []string, line string) (weights []int, lits []int, err error) {206+ weights = make([]int, 0, len(terms)/2)207+ lits = make([]int, 0, len(terms)/2)208+ i := 0209+ for i < len(terms) {210+ var l string211+ w, err := strconv.Atoi(terms[i])212+ if err != nil {213+ l = terms[i]214+ if !strings.HasPrefix(l, "x") && !strings.HasPrefix(l, "~x") {215+ return nil, nil, fmt.Errorf("invalid weight %q in %q: %v", terms[i*2], line, err)216+ }217+ // This is a weightless lit, i.e a lit with weight 1.218+ weights = append(weights, 1)219+ } else {220+ weights = append(weights, w)221+ i++222+ l = terms[i]223+ if !strings.HasPrefix(l, "x") && !strings.HasPrefix(l, "~x") || len(l) < 2 {224+ return nil, nil, fmt.Errorf("invalid variable name %q in %q", l, line)225+ }226+ }227+ var lit int228+ if l[0] == '~' {229+ lit, err = strconv.Atoi(l[2:])230+ lits = append(lits, -lit)231+ } else {232+ lit, err = strconv.Atoi(l[1:])233+ lits = append(lits, lit)234+ }235+ if err != nil {236+ return nil, nil, fmt.Errorf("invalid variable %q in %q: %v", l, line, err)237+ }238+ if lit > pb.NbVars {239+ pb.NbVars = lit240+ }241+ i++242+ }243+ return weights, lits, nil244+}245+246+// ParseOPB parses a file corresponding to the OPB syntax.247+// See http://www.cril.univ-artois.fr/PB16/format.pdf for more details.248+func ParseOPB(f io.Reader) (*Problem, error) {249+ scanner := bufio.NewScanner(f)250+ var pb Problem251+ for scanner.Scan() {252+ line := scanner.Text()253+ if line == "" || line[0] == '*' {254+ continue255+ }256+ if err := pb.parsePBLine(line); err != nil {257+ return nil, err258+ }259+ }260+ if err := scanner.Err(); err != nil {261+ return nil, fmt.Errorf("could not parse OPB: %v", err)262+ }263+ pb.Model = make([]decLevel, pb.NbVars)264+ pb.simplifyPB()265+ return &pb, nil266+}@@ -0,0 +1,113 @@
1+package solver2+3+// A PBConstr is a Pseudo-Boolean constraint.4+type PBConstr struct {5+ Lits []int // List of literals, designed with integer values. A positive value means the literal is true, a negative one it is false.6+ Weights []int // Weight of each lit from Lits. If nil, all lits == 17+ AtLeast int // Sum of all lits must be at least this value8+}9+10+// WeightSum returns the sum of the weight of all terms.11+func (c PBConstr) WeightSum() int {12+ if c.Weights == nil { // All weights = 113+ return len(c.Lits)14+ }15+ res := 016+ for _, w := range c.Weights {17+ res += w18+ }19+ return res20+}21+22+// Clause returns the clause (in fact, a constraint but the type is named Clause) associated with the given constraint.23+func (c PBConstr) Clause() *Clause {24+ lits := make([]Lit, len(c.Lits))25+ for i, val := range c.Lits {26+ lits[i] = IntToLit(int32(val))27+ }28+ // Saturate weights29+ for i := range c.Weights {30+ if c.Weights[i] > c.AtLeast {31+ c.Weights[i] = c.AtLeast32+ }33+ }34+ return NewPBClause(lits, c.Weights, c.AtLeast)35+}36+37+// PropClause returns a PB constraint equivalent to a propositional clause: at least one of the given38+// literals must be true.39+// It takes ownership of lits.40+func PropClause(lits ...int) PBConstr {41+ return PBConstr{Lits: lits, AtLeast: 1}42+}43+44+// AtLeast returns a PB constraint stating that at least n literals must be true.45+// It takes ownership of lits.46+func AtLeast(lits []int, n int) PBConstr {47+ return PBConstr{Lits: lits, AtLeast: n}48+}49+50+// AtMost returns a PB constraint stating that at most n literals can be true.51+// It takes ownership of lits.52+func AtMost(lits []int, n int) PBConstr {53+ lits2 := make([]int, len(lits))54+ for i := range lits {55+ lits2[i] = -lits[i]56+ }57+ return PBConstr{Lits: lits2, AtLeast: len(lits2) - n}58+}59+60+// GtEq returns a PB constraint stating that the sum of all literals multiplied by their weight61+// must be at least n.62+// Will panic if len(weights) != len(lits).63+func GtEq(lits []int, weights []int, n int) PBConstr {64+ if len(weights) != 0 && len(lits) != len(weights) {65+ panic("not as many lits as weights")66+ }67+ for i := 0; i < len(weights); i++ {68+ if weights[i] < 0 {69+ weights[i] = -weights[i]70+ n += weights[i]71+ lits[i] = -lits[i]72+ }73+ if weights[i] == 0 {74+ weights = append(weights[:i], weights[i+1:]...)75+ lits = append(lits[:i], lits[i+1:]...)76+ i--77+ }78+ }79+ return PBConstr{Lits: lits, Weights: weights, AtLeast: n}80+}81+82+// LtEq returns a PB constraint stating that the sum of all literals multiplied by their weight83+// must be at most n.84+// Will panic if len(weights) != len(lits).85+func LtEq(lits []int, weights []int, n int) PBConstr {86+ sum := 087+ for i := range lits {88+ lits[i] = -lits[i]89+ sum += weights[i]90+ }91+ n = sum - n92+ return GtEq(lits, weights, n)93+}94+95+// Eq returns a set of PB constraints stating that the sum of all literals multiplied by their weight96+// must be exactly n.97+// Will panic if len(weights) != len(lits).98+func Eq(lits []int, weights []int, n int) []PBConstr {99+ lits2 := make([]int, len(lits))100+ weights2 := make([]int, len(weights))101+ copy(lits2, lits)102+ copy(weights2, weights)103+ ge := GtEq(lits2, weights2, n)104+ le := LtEq(lits, weights, n)105+ var res []PBConstr106+ if ge.AtLeast > 0 {107+ res = append(res, ge)108+ }109+ if le.AtLeast > 0 {110+ res = append(res, le)111+ }112+ return res113+}@@ -0,0 +1,183 @@
1+package solver2+3+//4+// import "log"5+//6+// // Subsumes returns true iff c subsumes c2.7+// func (c *Clause) Subsumes(c2 *Clause) bool {8+// if c.Len() > c2.Len() {9+// return false10+// }11+// for _, lit := range c.lits {12+// found := false13+// for _, lit2 := range c2.lits {14+// if lit == lit2 {15+// found = true16+// break17+// }18+// if lit2 > lit {19+// return false20+// }21+// }22+// if !found {23+// return false24+// }25+// }26+// return true27+// }28+//29+// // SelfSubsumes returns true iff c self-subsumes c2.30+// func (c *Clause) SelfSubsumes(c2 *Clause) bool {31+// oneNeg := false32+// for _, lit := range c.lits {33+// found := false34+// for _, lit2 := range c2.lits {35+// if lit == lit2 {36+// found = true37+// break38+// }39+// if lit == lit2.Negation() {40+// if oneNeg { // We want exactly one, but this is the second41+// return false42+// }43+// oneNeg = true44+// found = true45+// break46+// }47+// if lit2 > lit { // We won't find it anymore48+// return false49+// }50+// }51+// if !found {52+// return false53+// }54+// }55+// return oneNeg56+// }57+//58+// // Simplify simplifies the given clause by removing redundant lits.59+// // If the clause is trivially satisfied (i.e contains both a lit and its negation),60+// // true is returned. Otherwise, false is returned.61+// func (c *Clause) Simplify() (isSat bool) {62+// c.Sort()63+// lits := make([]Lit, 0, len(c.lits))64+// i := 065+// for i < len(c.lits) {66+// if i < len(c.lits)-1 && c.lits[i] == c.lits[i+1].Negation() {67+// return true68+// }69+// lit := c.lits[i]70+// lits = append(lits, lit)71+// i++72+// for i < len(c.lits) && c.lits[i] == lit {73+// i++74+// }75+// }76+// if len(lits) < len(c.lits) {77+// c.lits = lits78+// }79+// return false80+// }81+//82+// // Generate returns a subsumed clause from c and c2, by removing v.83+// func (c *Clause) Generate(c2 *Clause, v Var) *Clause {84+// c3 := &Clause{lits: make([]Lit, 0, len(c.lits)+len(c2.lits)-2)}85+// for _, lit := range c.lits {86+// if lit.Var() != v {87+// c3.lits = append(c3.lits, lit)88+// }89+// }90+// for _, lit2 := range c2.lits {91+// if lit2.Var() != v {92+// c3.lits = append(c3.lits, lit2)93+// }94+// }95+// return c396+// }97+//98+// func (pb *Problem) preprocess() {99+// log.Printf("Preprocessing... %d clauses currently", len(pb.Clauses))100+// occurs := make([][]int, pb.NbVars*2)101+// for i, c := range pb.Clauses {102+// for j := 0; j < c.Len(); j++ {103+// occurs[c.Get(j)] = append(occurs[c.Get(j)], i)104+// }105+// }106+// modified := true107+// neverModified := true108+// for modified {109+// modified = false110+// for i := 0; i < pb.NbVars; i++ {111+// if pb.Model[i] != 0 {112+// continue113+// }114+// v := Var(i)115+// lit := v.Lit()116+// nbLit := len(occurs[lit])117+// nbLit2 := len(occurs[lit.Negation()])118+// if (nbLit < 10 || nbLit2 < 10) && (nbLit != 0 || nbLit2 != 0) {119+// modified = true120+// neverModified = false121+// // pb.deleted[v] = true122+// log.Printf("%d can be removed: %d and %d", lit.Int(), len(occurs[lit]), len(occurs[lit.Negation()]))123+// for _, idx1 := range occurs[lit] {124+// for _, idx2 := range occurs[lit.Negation()] {125+// c1 := pb.Clauses[idx1]126+// c2 := pb.Clauses[idx2]127+// newC := c1.Generate(c2, v)128+// if !newC.Simplify() {129+// switch newC.Len() {130+// case 0:131+// log.Printf("Inferred UNSAT")132+// pb.Status = Unsat133+// return134+// case 1:135+// log.Printf("Unit %d", newC.First().Int())136+// lit2 := newC.First()137+// if lit2.IsPositive() {138+// if pb.Model[lit2.Var()] == -1 {139+// pb.Status = Unsat140+// return141+// }142+// pb.Model[lit2.Var()] = 1143+// } else {144+// if pb.Model[lit2.Var()] == 1 {145+// pb.Status = Unsat146+// return147+// }148+// pb.Model[lit2.Var()] = -1149+// }150+// pb.Units = append(pb.Units, lit2)151+// default:152+// pb.Clauses = append(pb.Clauses, newC)153+// }154+// }155+// }156+// }157+// nbRemoved := 0158+// for _, idx := range occurs[lit] {159+// pb.Clauses[idx] = pb.Clauses[len(pb.Clauses)-nbRemoved-1]160+// nbRemoved++161+// }162+// for _, idx := range occurs[lit.Negation()] {163+// pb.Clauses[idx] = pb.Clauses[len(pb.Clauses)-nbRemoved-1]164+// nbRemoved++165+// }166+// pb.Clauses = pb.Clauses[:len(pb.Clauses)-nbRemoved]167+// log.Printf("clauses=%s", pb.CNF())168+// // Redo occurs169+// occurs = make([][]int, pb.NbVars*2)170+// for i, c := range pb.Clauses {171+// for j := 0; j < c.Len(); j++ {172+// occurs[c.Get(j)] = append(occurs[c.Get(j)], i)173+// }174+// }175+// continue176+// }177+// }178+// }179+// if !neverModified {180+// pb.simplify()181+// }182+// log.Printf("Done. %d clauses now", len(pb.Clauses))183+// }@@ -0,0 +1,396 @@
1+package solver2+3+import (4+ "fmt"5+)6+7+// A Problem is a list of clauses & a nb of vars.8+type Problem struct {9+ NbVars int // Total nb of vars10+ Clauses []*Clause // List of non-empty, non-unit clauses11+ Status Status // Status of the problem. Can be trivially UNSAT (if empty clause was met or inferred by UP) or Indet.12+ Units []Lit // List of unit literal found in the problem.13+ Model []decLevel // For each var, its inferred binding. 0 means unbound, 1 means bound to true, -1 means bound to false.14+ minLits []Lit // For an optimisation problem, the list of lits whose sum must be minimized15+ minWeights []int // For an optimisation problem, the weight of each lit.16+}17+18+// Optim returns true iff pb is an optimisation problem, ie19+// a problem for which we not only want to find a model, but also20+// the best possible model according to an optimization constraint.21+func (pb *Problem) Optim() bool {22+ return pb.minLits != nil23+}24+25+// CNF returns a DIMACS CNF representation of the problem.26+func (pb *Problem) CNF() string {27+ res := fmt.Sprintf("p cnf %d %d\n", pb.NbVars, len(pb.Clauses)+len(pb.Units))28+ for _, unit := range pb.Units {29+ res += fmt.Sprintf("%d 0\n", unit.Int())30+ }31+ for _, clause := range pb.Clauses {32+ res += fmt.Sprintf("%s\n", clause.CNF())33+ }34+ return res35+}36+37+// PBString returns a representation of the problem as a pseudo-boolean problem.38+func (pb *Problem) PBString() string {39+ res := pb.costFuncString()40+ for _, unit := range pb.Units {41+ sign := ""42+ if !unit.IsPositive() {43+ sign = "~"44+ unit = unit.Negation()45+ }46+ res += fmt.Sprintf("1 %sx%d = 1 ;\n", sign, unit.Int())47+ }48+ for _, clause := range pb.Clauses {49+ res += fmt.Sprintf("%s\n", clause.PBString())50+ }51+ return res52+}53+54+// SetCostFunc sets the function to minimize when optimizing the problem.55+// If all weights are 1, weights can be nil.56+// In all other cases, len(lits) must be the same as len(weights).57+func (pb *Problem) SetCostFunc(lits []Lit, weights []int) {58+ if weights != nil && len(lits) != len(weights) {59+ panic("length of lits and of weights don't match")60+ }61+ pb.minLits = lits62+ pb.minWeights = weights63+}64+65+// costFuncString returns a string representation of the cost function of the problem, if any, followed by a \n.66+// If there is no cost function, the empty string will be returned.67+func (pb *Problem) costFuncString() string {68+ if pb.minLits == nil {69+ return ""70+ }71+ res := "min: "72+ for i, lit := range pb.minLits {73+ w := 174+ if pb.minWeights != nil {75+ w = pb.minWeights[i]76+ }77+ sign := ""78+ if w >= 0 && i != 0 { // No plus sign for the first term or for negative terms.79+ sign = "+"80+ }81+ val := lit.Int()82+ neg := ""83+ if val < 0 {84+ val = -val85+ neg = "~"86+ }87+ res += fmt.Sprintf("%s%d %sx%d", sign, w, neg, val)88+ }89+ res += " ;\n"90+ return res91+}92+93+func (pb *Problem) updateStatus(nbClauses int) {94+ pb.Clauses = pb.Clauses[:nbClauses]95+ if pb.Status == Indet && nbClauses == 0 {96+ pb.Status = Sat97+ }98+}99+100+// simplify simplifies the pure SAT problem, i.e runs unit propagation if possible.101+func (pb *Problem) simplify2() {102+ nbClauses := len(pb.Clauses)103+ restart := true104+ for restart {105+ restart = false106+ i := 0107+ for i < nbClauses {108+ c := pb.Clauses[i]109+ nbLits := c.Len()110+ clauseSat := false111+ j := 0112+ for j < nbLits {113+ lit := c.Get(j)114+ k := j + 1115+ for k < nbLits {116+ lit2 := c.Get(k)117+ if lit2 == lit.Negation() {118+ clauseSat = true119+ break120+ }121+ if lit2 == lit { // duplicate lit122+ nbLits--123+ c.Set(k, c.Get(nbLits))124+ } else {125+ k++126+ }127+ }128+ if clauseSat {129+ clauseSat = true130+ break131+ }132+ if pb.Model[lit.Var()] == 0 {133+ j++134+ } else if (pb.Model[lit.Var()] == 1) == lit.IsPositive() {135+ clauseSat = true136+ break137+ } else {138+ nbLits--139+ c.Set(j, c.Get(nbLits))140+ }141+ }142+ if clauseSat {143+ nbClauses--144+ pb.Clauses[i] = pb.Clauses[nbClauses]145+ } else if nbLits == 0 {146+ pb.Status = Unsat147+ return148+ } else if nbLits == 1 { // UP149+ pb.addUnit(c.First())150+ if pb.Status == Unsat {151+ return152+ }153+ nbClauses--154+ pb.Clauses[i] = pb.Clauses[nbClauses]155+ restart = true // Must restart, since this lit might have made one more clause Unit or SAT.156+ } else { // nb lits unbound > cardinality157+ if c.Len() != nbLits {158+ c.Shrink(nbLits)159+ }160+ i++161+ }162+ }163+ }164+ pb.updateStatus(nbClauses)165+}166+167+// DetectAtMostOne tries to detect AtMostOne constraints encoded using the pairwise encoding.168+// It replaces those binary clauses by a single cardinality constraint.169+// This should mostly be called using the CuttingPlanes option, as it can dramatically improve the resolution process in some cases.170+func (pb *Problem) DetectAtMostOne() {171+ considered := make([]bool, pb.NbVars*2) // Has lit 1 already been detected in a clique?172+ propagates := make([][]Lit, pb.NbVars*2) // For each lit, the literals it propagates in a binary clause173+ indexes := make([][]int, len(propagates)) // Indexes of binary clauses, to remove them efficiently174+ toRemove := make([]int, 0, 1_000) // Indexes of clauses that have to be removed afterwards175+ for i, c := range pb.Clauses {176+ if c.Len() == 2 {177+ lit1 := c.First()178+ neg1 := lit1.Negation()179+ lit2 := c.Second()180+ neg2 := lit2.Negation()181+ propagates[neg1] = append(propagates[neg1], lit2)182+ propagates[neg2] = append(propagates[neg2], lit1)183+ indexes[neg1] = append(indexes[neg1], i)184+ indexes[neg2] = append(indexes[neg2], i)185+ }186+ }187+ for i := range propagates {188+ if considered[i] {189+ continue190+ }191+ lit := Lit(i)192+ others := propagates[lit]193+ if len(others) < 2 { // We won't find a cardinality constraint here194+ continue195+ }196+ constr := []Lit{lit.Negation()}197+ for j, other := range others {198+ if considered[other] {199+ continue200+ }201+ ok := true202+ for j := 1; j < len(constr); j++ {203+ lit2 := constr[j].Negation()204+ found := false205+ for _, lit3 := range propagates[lit2] {206+ if lit3 == other {207+ found = true208+ break209+ }210+ }211+ if !found {212+ ok = false213+ break214+ }215+ }216+ if ok { // other was found in a binary clause with each literal in constr217+ constr = append(constr, other)218+ toRemove = append(toRemove, indexes[lit][j])219+ }220+ }221+ if len(constr) > 2 { // We detected a stronger cardinality constraint222+ for _, lit := range constr {223+ considered[lit.Negation()] = true224+ }225+ pb.Clauses = append(pb.Clauses, NewCardClause(constr, len(constr)-1))226+ }227+ }228+ pb.removeBinaries(toRemove)229+}230+231+// removeBinaries removes the binary clauses that were used to build the232+// hidden cardinality constraint whose lits are given as a parameter.233+func (pb *Problem) removeBinaries(toRemove []int) {234+ if len(toRemove) == 0 {235+ return236+ }237+ newClauses := make([]*Clause, 0, len(pb.Clauses)-len(toRemove))238+ i := 0239+ for j, c := range pb.Clauses {240+ if j == toRemove[i] {241+ i++242+ if i == len(toRemove) {243+ break244+ }245+ } else {246+ newClauses = append(newClauses, c)247+ }248+ }249+ pb.Clauses = newClauses250+}251+252+// simplifyCard simplifies the problem, i.e runs unit propagation if possible.253+func (pb *Problem) simplifyCard() {254+ nbClauses := len(pb.Clauses)255+ restart := true256+ for restart {257+ restart = false258+ i := 0259+ for i < nbClauses {260+ c := pb.Clauses[i]261+ nbLits := c.Len()262+ card := c.Cardinality()263+ clauseSat := false264+ nbSat := 0265+ j := 0266+ for j < nbLits {267+ lit := c.Get(j)268+ if pb.Model[lit.Var()] == 0 {269+ j++270+ } else if (pb.Model[lit.Var()] == 1) == lit.IsPositive() {271+ nbSat++272+ if nbSat == card {273+ clauseSat = true274+ break275+ }276+ } else {277+ nbLits--278+ c.Set(j, c.Get(nbLits))279+ }280+ }281+ if clauseSat {282+ nbClauses--283+ pb.Clauses[i] = pb.Clauses[nbClauses]284+ } else if nbLits < card {285+ pb.Status = Unsat286+ return287+ } else if nbLits == card { // UP288+ pb.addUnits(c, nbLits)289+ if pb.Status == Unsat {290+ return291+ }292+ nbClauses--293+ pb.Clauses[i] = pb.Clauses[nbClauses]294+ restart = true // Must restart, since this lit might have made one more clause Unit or SAT.295+ } else { // nb lits unbound > cardinality296+ if c.Len() != nbLits {297+ c.Shrink(nbLits)298+ }299+ i++300+ }301+ }302+ }303+ pb.updateStatus(nbClauses)304+}305+306+func (pb *Problem) simplifyPB() {307+ pb.replicateUnits()308+ modified := true309+ for modified {310+ modified = false311+ i := 0312+ for i < len(pb.Clauses) {313+ c := pb.Clauses[i]314+ j := 0315+ card := c.Cardinality()316+ wSum := c.WeightSum()317+ for j < c.Len() {318+ lit := c.Get(j)319+ v := lit.Var()320+ w := c.Weight(j)321+ if pb.Model[v] == 0 { // Literal not assigned: is it unit?322+ if wSum-w < card { // Lit must be true for the clause to be satisfiable323+ pb.addUnit(lit)324+ if pb.Status == Unsat {325+ return326+ }327+ c.removeLit(j)328+ card -= w329+ c.updateCardinality(-w)330+ wSum -= w331+ modified = true332+ } else {333+ j++334+ }335+ } else { // Bound literal: remove it and update, if needed, cardinality336+ wSum -= w337+ if (pb.Model[v] == 1) == lit.IsPositive() {338+ card -= w339+ c.updateCardinality(-w)340+ }341+ c.removeLit(j)342+ modified = true343+ }344+ }345+ if card <= 0 { // Clause is Sat346+ pb.Clauses[i] = pb.Clauses[len(pb.Clauses)-1]347+ pb.Clauses = pb.Clauses[:len(pb.Clauses)-1]348+ modified = true349+ } else if wSum < card {350+ pb.Clauses = nil351+ pb.Status = Unsat352+ return353+ } else {354+ i++355+ }356+ }357+ }358+ if pb.Status == Indet && len(pb.Clauses) == 0 {359+ pb.Status = Sat360+ }361+}362+363+func (pb *Problem) addUnit(lit Lit) {364+ if lit.IsPositive() {365+ if pb.Model[lit.Var()] == -1 {366+ pb.Status = Unsat367+ return368+ }369+ pb.Model[lit.Var()] = 1370+ } else {371+ if pb.Model[lit.Var()] == 1 {372+ pb.Status = Unsat373+ return374+ }375+ pb.Model[lit.Var()] = -1376+ }377+ pb.Units = append(pb.Units, lit)378+}379+380+func (pb *Problem) addUnits(c *Clause, nbLits int) {381+ for i := 0; i < nbLits; i++ {382+ lit := c.Get(i)383+ pb.addUnit(lit)384+ }385+}386+387+func (pb *Problem) replicateUnits() {388+ for _, unit := range pb.Units {389+ v := unit.Var()390+ if unit.IsPositive() {391+ pb.Model[v] = 1392+ } else {393+ pb.Model[v] = -1394+ }395+ }396+}@@ -0,0 +1,128 @@
1+/******************************************************************************************[Heap.h]2+Copyright (c) 2003-2006, Niklas Een, Niklas Sorensson3+Copyright (c) 2007-2010, Niklas Sorensson4+5+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and6+associated documentation files (the "Software"), to deal in the Software without restriction,7+including without limitation the rights to use, copy, modify, merge, publish, distribute,8+sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is9+furnished to do so, subject to the following conditions:10+11+The above copyright notice and this permission notice shall be included in all copies or12+substantial portions of the Software.13+14+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT15+NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND16+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,17+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT18+OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.19+**************************************************************************************************/20+21+package solver22+23+// A heap implementation with support for decrease/increase key. This is24+// strongly inspired from Minisat's mtl/Heap.h.25+26+type queue struct {27+ activity []float64 // Activity of each variable. This should be the solver's slice, not a copy.28+ content []int // Actual content.29+ indices []int // Reverse queue, i.e position of each item in content; -1 means absence.30+}31+32+func newQueue(activity []float64) queue {33+ q := queue{34+ activity: activity,35+ }36+ for i := range q.activity {37+ q.insert(i)38+ }39+ return q40+}41+42+func (q *queue) lt(i, j int) bool {43+ return q.activity[i] > q.activity[j]44+}45+46+// Traversal functions.47+func left(i int) int { return i*2 + 1 }48+func right(i int) int { return (i + 1) * 2 }49+func parent(i int) int { return (i - 1) >> 1 }50+51+func (q *queue) percolateUp(i int) {52+ x := q.content[i]53+ p := parent(i)54+ for i != 0 && q.lt(x, q.content[p]) {55+ q.content[i] = q.content[p]56+ q.indices[q.content[p]] = i57+ i = p58+ p = parent(p)59+ }60+ q.content[i] = x61+ q.indices[x] = i62+}63+64+func (q *queue) percolateDown(i int) {65+ x := q.content[i]66+ for left(i) < len(q.content) {67+ var child int68+ if right(i) < len(q.content) && q.lt(q.content[right(i)], q.content[left(i)]) {69+ child = right(i)70+ } else {71+ child = left(i)72+ }73+ if !q.lt(q.content[child], x) {74+ break75+ }76+ q.content[i] = q.content[child]77+ q.indices[q.content[i]] = i78+ i = child79+ }80+ q.content[i] = x81+ q.indices[x] = i82+}83+84+func (q *queue) empty() bool { return len(q.content) == 0 }85+86+func (q *queue) contains(n int) bool {87+ return n < len(q.indices) && q.indices[n] >= 088+}89+90+func (q *queue) decrease(n int) {91+ q.percolateUp(q.indices[n])92+}93+94+func (q *queue) insert(n int) {95+ for i := len(q.indices); i <= n; i++ {96+ q.indices = append(q.indices, -1)97+ }98+ q.indices[n] = len(q.content)99+ q.content = append(q.content, n)100+ q.percolateUp(q.indices[n])101+}102+103+func (q *queue) removeMin() int {104+ x := q.content[0]105+ q.content[0] = q.content[len(q.content)-1]106+ q.indices[q.content[0]] = 0107+ q.indices[x] = -1108+ q.content = q.content[:len(q.content)-1]109+ if len(q.content) > 1 {110+ q.percolateDown(0)111+ }112+ return x113+}114+115+// Rebuild the heap from scratch, using the elements in 'ns':116+func (q *queue) build(ns []int) {117+ for i := range q.content {118+ q.indices[q.content[i]] = -1119+ }120+ q.content = q.content[:0]121+ for i, val := range ns {122+ q.indices[val] = i123+ q.content = append(q.content, val)124+ }125+ for i := len(q.content)/2 - 1; i >= 0; i-- {126+ q.percolateDown(i)127+ }128+}@@ -0,0 +1,1061 @@
1+package solver2+3+import (4+ "fmt"5+ "sort"6+ "strings"7+ "time"8+)9+10+const (11+ initNbMaxClauses = 2_000 // Maximum # of learned clauses, at first.12+ incrNbMaxClauses = 300 // By how much # of learned clauses is incremented at each conflict.13+ incrPostponeNbMax = 1_000 // By how much # of learned is increased when lots of good clauses are currently learned.14+ clauseDecay = 0.999 // By how much clauses bumping decays over time.15+ defaultVarDecay = 0.8 // On each var decay, how much the varInc should be decayed at startup16+)17+18+// Stats are statistics about the resolution of the problem.19+// They are provided for information purpose only.20+type Stats struct {21+ NbRestarts int22+ NbConflicts int23+ NbDecisions int24+ NbUnitLearned int // How many unit clauses were learned25+ NbBinaryLearned int // How many binary clauses were learned26+ NbLearned int // How many clauses were learned27+ NbDeleted int // How many clauses were deleted28+}29+30+// The level a decision was made.31+// A negative value means "negative assignement at that level".32+// A positive value means "positive assignment at that level".33+type decLevel int34+35+// A Model is a binding for several variables.36+// It can be totally bound (i.e all vars have a true or false binding)37+// or only partially (i.e some vars have no binding yet or their binding has no impact).38+// Each var, in order, is associated with a binding. Binding are implemented as39+// decision levels:40+// - a 0 value means the variable is free,41+// - a positive value means the variable was set to true at the given decLevel,42+// - a negative value means the variable was set to false at the given decLevel.43+type Model []decLevel44+45+func (m Model) String() string {46+ bound := make(map[int]decLevel)47+ for i := range m {48+ if m[i] != 0 {49+ bound[i+1] = m[i]50+ }51+ }52+ return fmt.Sprintf("%v", bound)53+}54+55+// A Solver solves a given problem. It is the main data structure.56+type Solver struct {57+ Verbose bool // Indicates whether the solver should display information during solving or not. False by default58+ Certified bool // Indicates whether a certificate should be generated during solving or not, using the RUP notation. This is useful to prove UNSAT instances. False by default.59+ CertChan chan string // Indicates where to write the certificate. If Certified is true but CertChan is nil, the certificate will be written on stdout.60+ CuttingPlanes bool // Indicates that the cutting planes resolution method should be used. Note that this is only efficient on PB problems.61+ nbVars int62+ status Status63+ wl watcherList64+ trail []Lit // Current assignment stack65+ model Model // 0 means unbound, other value is a binding66+ lastModel Model // Placeholder for last model found, useful when looking for several models67+ activity []float64 // How often each var is involved in conflicts68+ polarity []bool // Preferred sign for each var69+ assumptions []bool // True iff the var's binding is assumed70+ // For each var, clause considered when it was unified71+ // If the var is not bound yet, or if it was bound by a decision, value is nil.72+ reason []*Clause73+ varQueue queue74+ varInc float64 // On each var bump, how big the increment should be75+ clauseInc float32 // On each var bump, how big the increment should be76+ lbdStats lbdStats77+ lubyNextRestart int // When will the next restart happen when using Luby's strategy?78+ Stats Stats // Statistics about the solving process.79+ minLits []Lit // Lits to minimize if the problem was an optimization problem.80+ minWeights []int // Weight of each lit to minimize if the problem was an optimization problem.81+ hypothesis []Lit // Literals that are, ideally, true. Useful when trying to minimize a function.82+ localNbRestarts int // How many restarts since Solve() was called?83+ varDecay float64 // On each var decay, how much the varInc should be decayed84+ trailBuf []int // A buffer while cleaning bindings85+ pbSetBuf []int // A buffer to reduce allocation when performing cutting planes86+ pbSetBuf2 []int // A buffer to reduce allocation when performing cutting planes87+ bufLits []Lit // Buffer for lits in learnClause (per solver; see learn.go)88+}89+90+// New makes a solver, given a number of variables and a set of clauses.91+// nbVars should be consistent with the content of clauses, i.e.92+// the biggest variable in clauses should be >= nbVars.93+func New(problem *Problem) *Solver {94+ if problem.Status == Unsat {95+ return &Solver{status: Unsat}96+ }97+ nbVars := problem.NbVars98+99+ trailCap := nbVars100+ if len(problem.Units) > trailCap {101+ trailCap = len(problem.Units)102+ }103+104+ s := &Solver{105+ nbVars: nbVars,106+ status: problem.Status,107+ trail: make([]Lit, len(problem.Units), trailCap),108+ model: problem.Model,109+ activity: make([]float64, nbVars),110+ polarity: make([]bool, nbVars),111+ assumptions: make([]bool, nbVars),112+ reason: make([]*Clause, nbVars),113+ varInc: 1.0,114+ clauseInc: 1.0,115+ lubyNextRestart: int(lubyConstant * luby(1)),116+ minLits: problem.minLits,117+ minWeights: problem.minWeights,118+ varDecay: defaultVarDecay,119+ trailBuf: make([]int, nbVars),120+ pbSetBuf: make([]int, nbVars),121+ pbSetBuf2: make([]int, nbVars),122+ }123+ s.resetOptimPolarity()124+ s.initOptimActivity()125+ s.initWatcherList(problem.Clauses)126+ s.varQueue = newQueue(s.activity)127+ for i, lit := range problem.Units {128+ if lit.IsPositive() {129+ s.model[lit.Var()] = 1130+ } else {131+ s.model[lit.Var()] = -1132+ }133+ s.trail[i] = lit134+ }135+ return s136+}137+138+// newVar is used to indicate a new variable must be added to the solver.139+// This can be used when new clauses are appended and these clauses contain vars that were unseen so far.140+// If the var already existed, nothing will happen.141+func (s *Solver) newVar(v Var) {142+ if cnfVar := int(v.Int()); cnfVar > s.nbVars {143+ // If the var already existed, do nothing144+ for i := s.nbVars; i < cnfVar; i++ {145+ s.model = append(s.model, 0)146+ s.activity = append(s.activity, 0.)147+ s.polarity = append(s.polarity, false)148+ s.reason = append(s.reason, nil)149+ s.trailBuf = append(s.trailBuf, 0)150+ }151+ s.varQueue = newQueue(s.activity)152+ s.addVarWatcherList(v)153+ s.nbVars = cnfVar154+ }155+}156+157+// sets initial activity for optimization variables, if any.158+func (s *Solver) initOptimActivity() {159+ for i, lit := range s.minLits {160+ w := 1161+ if s.minWeights != nil {162+ w = s.minWeights[i]163+ }164+ s.activity[lit.Var()] += float64(w)165+ }166+}167+168+// resets polarity of optimization lits so that they are negated by default.169+func (s *Solver) resetOptimPolarity() {170+ if s.minLits != nil {171+ for _, lit := range s.minLits {172+ s.polarity[lit.Var()] = !lit.IsPositive() // Try to make lits from the optimization clause false173+ }174+ }175+}176+177+// Optim returns true iff the underlying problem is an optimization problem (rather than a satisfaction one).178+func (s *Solver) Optim() bool {179+ return s.minLits != nil180+}181+182+// OutputModel outputs the model for the problem on stdout.183+func (s *Solver) OutputModel() {184+ if s.status == Sat || s.lastModel != nil {185+ fmt.Printf("s SATISFIABLE\nv ")186+ model := s.model187+ if s.lastModel != nil {188+ model = s.lastModel189+ }190+ for i, val := range model {191+ if val < 0 {192+ fmt.Printf("%d ", -i-1)193+ } else {194+ fmt.Printf("%d ", i+1)195+ }196+ }197+ fmt.Printf("\n")198+ } else if s.status == Unsat {199+ fmt.Printf("s UNSATISFIABLE\n")200+ } else {201+ fmt.Printf("s INDETERMINATE\n")202+ }203+}204+205+// litStatus returns whether the literal is made true (Sat) or false (Unsat) by the206+// current bindings, or if it is unbounded (Indet).207+func (s *Solver) litStatus(l Lit) Status {208+ assign := s.model[l.Var()]209+ if assign == 0 {210+ return Indet211+ }212+ if assign > 0 == l.IsPositive() {213+ return Sat214+ }215+ return Unsat216+}217+218+func (s *Solver) varDecayActivity() {219+ s.varInc *= 1 / s.varDecay220+}221+222+func (s *Solver) varBumpActivity(v Var) {223+ // fmt.Printf("bumping var %d\n", v.Int())224+ s.activity[v] += s.varInc225+ if s.activity[v] > 1e100 { // Rescaling is needed to avoid overflowing226+ for i := range s.activity {227+ s.activity[i] *= 1e-100228+ }229+ s.varInc *= 1e-100230+ }231+ if s.varQueue.contains(int(v)) {232+ s.varQueue.decrease(int(v))233+ }234+}235+236+// Decays each clause's activity237+func (s *Solver) clauseDecayActivity() {238+ s.clauseInc *= 1 / clauseDecay239+}240+241+// Bumps the given clause's activity.242+func (s *Solver) clauseBumpActivity(c *Clause) {243+ if c.Learned() {244+ c.activity += s.clauseInc245+ if c.activity > 1e30 { // Rescale to avoid overflow246+ for _, c2 := range s.wl.learned {247+ c2.activity *= 1e-30248+ }249+ s.clauseInc *= 1e-30250+ }251+ }252+}253+254+// Chooses an unbound literal to be tested, or -1255+// if all the variables are already bound.256+func (s *Solver) chooseLit() Lit {257+ v := Var(-1)258+ for v == -1 && !s.varQueue.empty() {259+ if v2 := Var(s.varQueue.removeMin()); s.model[v2] == 0 { // Ignore already bound vars260+ v = v2261+ }262+ }263+ if v == -1 {264+ return Lit(-1)265+ }266+ s.Stats.NbDecisions++267+ return v.SignedLit(!s.polarity[v])268+}269+270+type number interface {271+ int | int32 | decLevel272+}273+274+func abs[T number](val T) T {275+ if val < 0 {276+ return -val277+ }278+ return val279+}280+281+func min[T number](a, b T) T {282+ if a < b {283+ return a284+ }285+ return b286+}287+288+// Reinitializes bindings (both model & reason) for all variables bound at a decLevel >= lvl.289+// TODO: check this method as it has a weird behavior regarding performance.290+// TODO: clean-up commented-out code and understand underlying performance pattern.291+func (s *Solver) cleanupBindings(lvl decLevel) {292+ i := 0293+ for i < len(s.trail) && abs(s.model[s.trail[i].Var()]) <= lvl {294+ i++295+ }296+ /*297+ for j := len(s.trail) - 1; j >= i; j-- {298+ lit2 := s.trail[j]299+ v := lit2.Var()300+ s.model[v] = 0301+ if s.reason[v] != nil {302+ s.reason[v].unlock()303+ s.reason[v] = nil304+ }305+ s.polarity[v] = lit2.IsPositive()306+ if !s.varQueue.contains(int(v)) {307+ s.varQueue.insert(int(v))308+ }309+ }310+ s.trail = s.trail[:i]311+ */312+ toInsert := s.trailBuf[:0] // make([]int, 0, len(s.trail)-i)313+ for j := i; j < len(s.trail); j++ {314+ lit2 := s.trail[j]315+ v := lit2.Var()316+ s.model[v] = 0317+ if s.reason[v] != nil {318+ s.reason[v].unlock()319+ s.reason[v] = nil320+ }321+ s.polarity[v] = lit2.IsPositive()322+ if !s.varQueue.contains(int(v)) {323+ toInsert = append(toInsert, int(v))324+ s.varQueue.insert(int(v))325+ }326+ }327+ s.trail = s.trail[:i]328+ for i := len(toInsert) - 1; i >= 0; i-- {329+ s.varQueue.insert(toInsert[i])330+ }331+ /*for i := len(s.trail) - 1; i >= 0; i-- {332+ lit := s.trail[i]333+ v := lit.Var()334+ if abs(s.model[v]) <= lvl { // All lits in trail before here must keep their status.335+ s.trail = s.trail[:i+1]336+ break337+ }338+ s.model[v] = 0339+ if s.reason[v] != nil {340+ s.reason[v].unlock()341+ s.reason[v] = nil342+ }343+ s.polarity[v] = lit.IsPositive()344+ if !s.varQueue.contains(int(v)) {345+ s.varQueue.insert(int(v))346+ }347+ }*/348+ s.resetOptimPolarity()349+}350+351+func (s *Solver) trailString() string {352+ res := ""353+ for i, lit := range s.trail {354+ v := lit.Var()355+ assign := s.model[v]356+ if assign == 0 {357+ panic(fmt.Sprintf("error: literal in trail at position %d but not assigned: %d", i, lit.Int()))358+ }359+ res += fmt.Sprintf("%d@%d ", lit.Int(), abs(assign))360+ }361+ return res362+}363+364+// Given the last learnt clause and the levels at which vars were bound,365+// Returns the level to bt to and the literal to bind366+func backtrackData(c *Clause, model []decLevel) (btLevel decLevel, lit Lit) {367+ btLevel = abs(model[c.Get(1).Var()])368+ return btLevel, c.Get(0)369+}370+371+func (s *Solver) rebuildOrderHeap() {372+ ints := make([]int, s.nbVars)373+ for v := 0; v < s.nbVars; v++ {374+ if s.model[v] == 0 {375+ ints = append(ints, int(v))376+ }377+ }378+ s.varQueue.build(ints)379+}380+381+// propagate binds the given lit, propagates it and searches for a solution,382+// until it is found or a restart is needed.383+func (s *Solver) propagateAndSearch(lit Lit, lvl decLevel) Status {384+ if s.CuttingPlanes {385+ return s.propagateAndSearchPB(lit, lvl)386+ }387+ for lit != -1 {388+ // log.Printf("picked %d at lvl %d", lit.Int(), lvl)389+ if conflict := s.unifyLiteral(lit, lvl); conflict == nil { // Pick new branch or restart390+ if s.lbdStats.mustRestart() {391+ s.lbdStats.clear()392+ s.cleanupBindings(1)393+ return Indet394+ }395+ if s.Stats.NbConflicts >= s.wl.idxReduce*s.wl.nbMax {396+ s.wl.idxReduce = s.Stats.NbConflicts/s.wl.nbMax + 1397+ s.reduceLearned()398+ s.bumpNbMax()399+ }400+ lvl++401+ lit = s.chooseLit()402+ } else { // Deal with conflict403+ s.Stats.NbConflicts++404+ if s.Stats.NbConflicts%5_000 == 0 && s.varDecay < 0.95 {405+ s.varDecay += 0.01406+ }407+ s.lbdStats.addConflict(len(s.trail))408+ learnt, unit := s.learnClause(conflict, lvl)409+ if learnt == nil { // Unit clause was learned: this lit is known for sure410+ if unit == -1 || (abs(s.model[unit.Var()]) == 1 && s.litStatus(unit) == Unsat) { // Top-level conflict411+ return s.setUnsat()412+ }413+ s.Stats.NbUnitLearned++414+ s.lbdStats.addLbd(1)415+ s.cleanupBindings(1)416+ s.addLearnedUnit(unit)417+ s.model[unit.Var()] = lvlToSignedLvl(unit, 1)418+ if conflict = s.unifyLiteral(unit, 1); conflict != nil { // top-level conflict419+ return s.setUnsat()420+ }421+ s.rebuildOrderHeap()422+ lit = s.chooseLit()423+ lvl = 2424+ } else {425+ if learnt.Len() == 2 {426+ s.Stats.NbBinaryLearned++427+ }428+ s.Stats.NbLearned++429+ s.lbdStats.addLbd(learnt.lbd())430+ s.addLearned(learnt)431+ lvl, lit = backtrackData(learnt, s.model)432+ s.cleanupBindings(lvl)433+ s.reason[lit.Var()] = learnt434+ learnt.lock()435+ }436+ }437+ }438+ return Sat439+}440+441+// propagateAndSearchPB performs pseudo-boolean constraint learning whenever a conflcit arises.442+func (s *Solver) propagateAndSearchPB(lit Lit, lvl decLevel) Status {443+ for lit != -1 {444+ // log.Printf("picked %d at lvl %d", lit.Int(), lvl)445+ if conflict := s.unifyLiteral(lit, lvl); conflict == nil { // Pick new branch or restart446+ if s.Stats.NbConflicts >= s.lubyNextRestart {447+ s.lubyNextRestart += int(lubyConstant * luby(uint(s.Stats.NbRestarts)+2))448+ s.cleanupBindings(1)449+ return Indet450+ }451+ if s.Stats.NbConflicts >= s.wl.idxReduce*s.wl.nbMax {452+ s.wl.idxReduce = s.Stats.NbConflicts/s.wl.nbMax + 1453+ s.reduceLearnedPB()454+ s.bumpNbMax()455+ }456+ lvl++457+ lit = s.chooseLit()458+ } else { // Deal with conflict459+ for conflict != nil {460+ // log.Printf("conflict: %s", conflict.PBString())461+ s.Stats.NbConflicts++462+ if s.Stats.NbConflicts%5_000 == 0 && s.varDecay < 0.95 {463+ s.varDecay += 0.01464+ }465+ s.lbdStats.addConflict(len(s.trail))466+ learnt, propagated, newLvl := s.cuttingPlanes(conflict, lvl)467+ // log.Printf("learnt=%v, propagated=%v, newLvl=%d", learnt, propagated, newLvl)468+ if newLvl == -1 { // Generated constraint is false469+ return s.setUnsat()470+ }471+ if newLvl == 1 {472+ for _, unit := range propagated {473+ if unit == -1 || (abs(s.model[unit.Var()]) == 1 && s.litStatus(unit) == Unsat) { // Top-level conflict474+ return s.setUnsat()475+ }476+ s.Stats.NbUnitLearned++477+ s.lbdStats.addLbd(1)478+ s.cleanupBindings(1)479+ s.addLearnedUnit(unit)480+ s.model[unit.Var()] = lvlToSignedLvl(unit, 1)481+ if conflict = s.unifyLiteral(unit, 1); conflict != nil { // top-level conflict482+ return s.setUnsat()483+ }484+ }485+ s.rebuildOrderHeap()486+ lit = s.chooseLit()487+ lvl = 2488+ } else {489+ lvl = newLvl490+ // A constraint was learned and lits have to be propagated at lvl > 1491+ // if learnt != nil {492+ // log.Printf("propagated the following lits at lvl %d because of %s:", lvl, learnt.PBString())493+ // } else {494+ // log.Printf("propagated the following lits at lvl %d with empty constr:", lvl)495+ // }496+ // for _, lit := range propagated {497+ // log.Printf("%d ", lit.Int())498+ // }499+ s.Stats.NbLearned++500+ s.addLearned(learnt)501+ learnt.lock()502+ lvl = newLvl503+ s.cleanupBindings(lvl)504+ for _, lit := range propagated {505+ s.reason[lit.Var()] = learnt506+ }507+ conflict = s.unifyLiterals(propagated, lvl)508+ if conflict == nil {509+ lit = s.chooseLit()510+ lvl++511+ }512+ }513+ }514+ }515+ }516+ return Sat517+}518+519+// Sets the status to unsat and do cleanup tasks.520+func (s *Solver) setUnsat() Status {521+ if s.Certified {522+ if s.CertChan == nil {523+ fmt.Printf("0\n")524+ } else {525+ s.CertChan <- "0"526+ }527+ }528+ s.status = Unsat529+ return Unsat530+}531+532+// Searches until a restart is needed.533+func (s *Solver) search() Status {534+ s.localNbRestarts++535+ lvl := decLevel(2) // Level starts at 2, for implementation reasons : 1 is for top-level bindings; 0 means "no level assigned yet"536+ // s.status = s.propagateAndSearch(s.chooseLit(), lvl)537+ s.status = s.propagateAndSearch(s.chooseLit(), lvl)538+ return s.status539+}540+541+// Solve solves the problem associated with the solver and returns the appropriate status.542+func (s *Solver) Solve() Status {543+ if s.status == Unsat {544+ return s.status545+ }546+ s.status = Indet547+ //s.lbdStats.clear()548+ s.localNbRestarts = 0549+ var end chan struct{}550+ if s.Verbose {551+ end = make(chan struct{})552+ defer close(end)553+ go func() { // Function displaying stats during resolution554+ fmt.Printf("c ======================================================================================\n")555+ fmt.Printf("c | Restarts | Conflicts | Learned | Deleted | Del%% | Reduce | Units learned |\n")556+ fmt.Printf("c ======================================================================================\n")557+ ticker := time.NewTicker(3 * time.Second)558+ defer ticker.Stop()559+ for { // There might be concurrent access in a few places but this is okay since we are very conservative and don't modify state.560+ select {561+ case <-ticker.C:562+ case <-end:563+ return564+ }565+ if s.status == Indet {566+ iter := s.Stats.NbRestarts + 1567+ nbConfl := s.Stats.NbConflicts568+ nbReduce := s.wl.idxReduce - 1569+ nbLearned := len(s.wl.learned)570+ nbDel := s.Stats.NbDeleted571+ pctDel := int(100 * float64(nbDel) / float64(s.Stats.NbLearned))572+ nbUnit := s.Stats.NbUnitLearned573+ fmt.Printf("c | %8d | %11d | %9d | %9d | %3d%% | %6d | %8d/%8d |\n", iter, nbConfl, nbLearned, nbDel, pctDel, nbReduce, nbUnit, s.nbVars)574+ }575+ }576+ }()577+ }578+ for s.status == Indet {579+ s.search()580+ if s.status == Indet {581+ s.Stats.NbRestarts++582+ s.rebuildOrderHeap()583+ }584+ }585+ if s.status == Sat {586+ s.lastModel = make(Model, len(s.model))587+ copy(s.lastModel, s.model)588+ }589+ if s.Verbose {590+ end <- struct{}{}591+ fmt.Printf("c ======================================================================================\n")592+ }593+ return s.status594+}595+596+// Assume adds unit literals to the solver.597+// This is useful when calling the solver several times, e.g to keep it "hot" while removing clauses.598+func (s *Solver) Assume(lits []Lit) Status {599+ s.cleanupBindings(0)600+ s.trail = s.trail[:0]601+ s.assumptions = make([]bool, s.nbVars)602+603+ for _, lit := range lits {604+605+ s.addLearnedUnit(lit)606+ s.assumptions[lit.Var()] = true607+ s.trail = append(s.trail, lit)608+ }609+ s.status = Indet610+ if confl := s.propagate(0, 1); confl != nil {611+ // Conflict after unit propagation612+ s.status = Unsat613+ return s.status614+ }615+ return s.status616+}617+618+// Enumerate returns the total number of models for the given problems.619+// if "models" is non-nil, it will write models on it as soon as it discovers them.620+// models will be closed at the end of the method.621+func (s *Solver) Enumerate(models chan []bool, stop chan struct{}) int {622+ if models != nil {623+ defer close(models)624+ }625+ s.lastModel = make(Model, len(s.model))626+ nb := 0627+ var lit Lit628+ var lvl decLevel629+ for s.status != Unsat {630+ for s.status == Indet {631+ s.search()632+ if s.status == Indet {633+ s.Stats.NbRestarts++634+ }635+ }636+ if s.status == Sat {637+ copy(s.lastModel, s.model)638+ if models != nil {639+ nb += s.addCurrentModels(models)640+ } else {641+ nb += s.countCurrentModels()642+ }643+ s.status = Indet644+ lits := s.decisionLits()645+ switch len(lits) {646+ case 0:647+ s.status = Unsat648+ case 1:649+ s.propagateUnits(lits)650+ default:651+ c := NewClause(lits)652+ s.appendClause(c)653+ lit = lits[len(lits)-1]654+ v := lit.Var()655+ lvl = abs(s.model[v]) - 1656+ s.cleanupBindings(lvl)657+ s.reason[v] = c // Must do it here because it won't be made by propagateAndSearch658+ s.propagateAndSearch(lit, lvl)659+ }660+ }661+ }662+ return nb663+}664+665+// CountModels returns the total number of models for the given problem.666+func (s *Solver) CountModels() int {667+ var end chan struct{}668+ if s.Verbose {669+ end = make(chan struct{})670+ defer close(end)671+ go func() { // Function displaying stats during resolution672+ fmt.Printf("c ======================================================================================\n")673+ fmt.Printf("c | Restarts | Conflicts | Learned | Deleted | Del%% | Reduce | Units learned |\n")674+ fmt.Printf("c ======================================================================================\n")675+ ticker := time.NewTicker(3 * time.Second)676+ defer ticker.Stop()677+ for { // There might be concurrent access in a few places but this is okay since we are very conservative and don't modify state.678+ select {679+ case <-ticker.C:680+ case <-end:681+ return682+ }683+ if s.status == Indet {684+ iter := s.Stats.NbRestarts + 1685+ nbConfl := s.Stats.NbConflicts686+ nbReduce := s.wl.idxReduce - 1687+ nbLearned := len(s.wl.learned)688+ nbDel := s.Stats.NbDeleted689+ pctDel := int(100 * float64(nbDel) / float64(s.Stats.NbLearned))690+ nbUnit := s.Stats.NbUnitLearned691+ fmt.Printf("c | %8d | %11d | %9d | %9d | %3d%% | %6d | %8d/%8d |\n", iter, nbConfl, nbLearned, nbDel, pctDel, nbReduce, nbUnit, s.nbVars)692+ }693+ }694+ }()695+ }696+ nb := 0697+ var lit Lit698+ var lvl decLevel699+ for s.status != Unsat {700+ for s.status == Indet {701+ s.search()702+ if s.status == Indet {703+ s.Stats.NbRestarts++704+ }705+ }706+ if s.status == Sat {707+ s.lastModel = s.model708+ nb += s.countCurrentModels()709+ if s.Verbose {710+ fmt.Printf("c found %d model(s)\n", nb)711+ }712+ s.status = Indet713+ lits := s.decisionLits()714+ switch len(lits) {715+ case 0:716+ s.status = Unsat717+ case 1:718+ s.propagateUnits(lits)719+ default:720+ c := NewClause(lits)721+ s.appendClause(c)722+ lit = lits[len(lits)-1]723+ v := lit.Var()724+ lvl = abs(s.model[v]) - 1725+ s.cleanupBindings(lvl)726+ s.reason[v] = c // Must do it here because it won't be made by propagateAndSearch727+ s.propagateAndSearch(lit, lvl)728+ }729+ }730+ }731+ if s.Verbose {732+ end <- struct{}{}733+ fmt.Printf("c ======================================================================================\n")734+ }735+ return nb736+}737+738+// decisionLits returns the negation of all decision values once a model was found, ordered by decision levels.739+// This will allow for searching other models.740+func (s *Solver) decisionLits() []Lit {741+ lastLit := s.trail[len(s.trail)-1]742+ lvls := abs(s.model[lastLit.Var()])743+ lits := make([]Lit, lvls-1)744+ for i, r := range s.reason {745+ if lvl := abs(s.model[i]); r == nil && lvl > 1 {746+ if s.model[i] < 0 {747+ // lvl-2 : levels beside unit clauses start at 2, not 0 or 1!748+ lits[lvl-2] = IntToLit(int32(i + 1))749+ } else {750+ lits[lvl-2] = IntToLit(int32(-i - 1))751+ }752+ }753+ }754+ return lits755+}756+757+func (s *Solver) propagateUnits(units []Lit) {758+ for _, unit := range units {759+ s.lbdStats.addLbd(1)760+ s.Stats.NbUnitLearned++761+ s.cleanupBindings(1)762+ s.model[unit.Var()] = lvlToSignedLvl(unit, 1)763+ if s.unifyLiteral(unit, 1) != nil {764+ s.status = Unsat765+ return766+ }767+ s.rebuildOrderHeap()768+ }769+}770+771+// PBString returns a representation of the solver's state as a pseudo-boolean problem.772+func (s *Solver) PBString() string {773+ meta := fmt.Sprintf("* #variable= %d #constraint= %d #learned= %d\n", s.nbVars, len(s.wl.origClauses), len(s.wl.learned))774+ minLine := ""775+ if s.minLits != nil {776+ terms := make([]string, len(s.minLits))777+ for i, lit := range s.minLits {778+ weight := 1779+ if s.minWeights != nil {780+ weight = s.minWeights[i]781+ }782+ val := lit.Int()783+ sign := ""784+ if val < 0 {785+ val = -val786+ sign = "~"787+ }788+ terms[i] = fmt.Sprintf("%d %sx%d", weight, sign, val)789+ }790+ minLine = fmt.Sprintf("min: %s ;\n", strings.Join(terms, " +"))791+ }792+ clauses := make([]string, len(s.wl.origClauses)+len(s.wl.learned))793+ for i, c := range s.wl.origClauses {794+ clauses[i] = c.PBString()795+ }796+ for i, c := range s.wl.learned {797+ clauses[i+len(s.wl.origClauses)] = c.PBString()798+ }799+ for i := 0; i < len(s.model); i++ {800+ if s.model[i] == 1 {801+ clauses = append(clauses, fmt.Sprintf("1 x%d = 1 ;", i+1))802+ } else if s.model[i] == -1 {803+ clauses = append(clauses, fmt.Sprintf("1 x%d = 0 ;", i+1))804+ }805+ }806+ return meta + minLine + strings.Join(clauses, "\n")807+}808+809+// AppendClause appends a new clause to the set of clauses.810+// This is not a learned clause, but a clause that is part of the problem added afterwards (during model counting, for instance).811+func (s *Solver) AppendClause(clause *Clause) {812+ s.cleanupBindings(1)813+ card := clause.Cardinality()814+ minW := 0815+ maxW := 0816+ i := 0817+ for i < clause.Len() {818+ lit := clause.Get(i)819+ s.newVar(lit.Var())820+ switch s.litStatus(lit) {821+ case Sat:822+ w := clause.Weight(i)823+ minW += w824+ maxW += w825+ clause.removeLit(i)826+ clause.updateCardinality(-w)827+ case Unsat:828+ clause.removeLit(i)829+ default:830+ maxW += clause.Weight(i)831+ i++832+ }833+ }834+ if minW >= card { // clause is already sat835+ return836+ }837+ if maxW < card { // clause cannot be satisfied838+ s.status = Unsat839+ return840+ }841+ if maxW == card { // Unit842+ s.propagateUnits(clause.lits)843+ } else {844+ s.appendClause(clause)845+ }846+}847+848+// Model returns a slice that associates, to each variable, its binding.849+// If s's status is not Sat, the method will panic.850+func (s *Solver) Model() []bool {851+ if s.lastModel == nil {852+ panic("cannot call Model() from a non-Sat solver")853+ }854+ res := make([]bool, s.nbVars)855+ for i, lvl := range s.lastModel {856+ res[i] = lvl > 0857+ }858+ return res859+}860+861+// addCurrentModels is called when a model was found.862+// It returns the total number of models from this point, and sends all models on ch.863+// The number can be different of 1 if there are unbound variables.864+// For instance, if there are 4 variables in the problem and only 1, 3 and 4 are bound,865+// there are actually 2 models currently: one with 2 set to true, the other with 2 set to false.866+func (s *Solver) addCurrentModels(ch chan []bool) int {867+ unbound := make([]int, 0, s.nbVars) // indices of unbound variables868+ var nb uint64 = 1 // total number of models found869+ model := make([]bool, s.nbVars) // partial model870+ for i, lvl := range s.lastModel {871+ if lvl == 0 {872+ unbound = append(unbound, i)873+ nb *= 2874+ } else {875+ model[i] = lvl > 0876+ }877+ }878+ for i := uint64(0); i < nb; i++ {879+ for j := range unbound {880+ mask := uint64(1 << j)881+ cur := i & mask882+ idx := unbound[j]883+ model[idx] = cur != 0884+ }885+ model2 := make([]bool, len(model))886+ copy(model2, model)887+ ch <- model2888+ }889+ return int(nb)890+}891+892+// countCurrentModels is called when a model was found.893+// It returns the total number of models from this point.894+// The number can be different of 1 if there are unbound variables.895+// For instance, if there are 4 variables in the problem and only 1, 3 and 4 are bound,896+// there are actually 2 models currently: one with 2 set to true, the other with 2 set to false.897+func (s *Solver) countCurrentModels() int {898+ var nb uint64 = 1 // total number of models found899+ for _, lvl := range s.lastModel {900+ if lvl == 0 {901+ nb *= 2902+ }903+ }904+ return int(nb)905+}906+907+// Optimal returns the optimal solution, if any.908+// If results is non-nil, all solutions will be written to it.909+// In any case, results will be closed at the end of the call.910+func (s *Solver) Optimal(results chan Result, stop chan struct{}) (res Result) {911+ if results != nil {912+ defer close(results)913+ }914+ status := s.Solve()915+ if status == Unsat { // Problem cannot be satisfied at all916+ res.Status = Unsat917+ if results != nil {918+ results <- res919+ }920+ return res921+ }922+ if s.minLits == nil { // No optimization clause: this is a decision problem, solution is optimal923+ s.lastModel = make(Model, len(s.model))924+ copy(s.lastModel, s.model)925+ res := Result{926+ Status: Sat,927+ Model: s.Model(),928+ Weight: 0,929+ }930+ if results != nil {931+ results <- res932+ }933+ return res934+ }935+ // log.Printf("found a solution, now minimizing...")936+ maxCost := 0937+ if s.minWeights == nil {938+ maxCost = len(s.minLits)939+ } else {940+ for _, w := range s.minWeights {941+ maxCost += w942+ }943+ }944+ s.hypothesis = make([]Lit, len(s.minLits))945+ for i, lit := range s.minLits {946+ s.hypothesis[i] = lit.Negation()947+ }948+ weights := make([]int, len(s.minWeights))949+ copy(weights, s.minWeights)950+ sort.Sort(wLits{lits: s.hypothesis, weights: weights})951+ s.lastModel = make(Model, len(s.model))952+ var cost int953+ for status == Sat {954+ copy(s.lastModel, s.model) // Save this model: it might be the last one955+ cost = 0956+ for i, lit := range s.minLits {957+ if s.model[lit.Var()] > 0 == lit.IsPositive() {958+ if s.minWeights == nil {959+ cost++960+ } else {961+ cost += s.minWeights[i]962+ }963+ }964+ }965+ res = Result{966+ Status: Sat,967+ Model: s.Model(),968+ Weight: cost,969+ }970+ // log.Printf("result=%v", res)971+ if results != nil {972+ results <- res973+ }974+ if cost == 0 {975+ break976+ }977+ // Add a constraint incrementing current best cost978+ lits2 := make([]Lit, len(s.minLits))979+ weights2 := make([]int, len(s.minWeights))980+ copy(lits2, s.hypothesis)981+ copy(weights2, weights)982+ s.AppendClause(NewPBClause(lits2, weights2, maxCost-cost+1))983+ s.rebuildOrderHeap()984+ status = s.Solve()985+ }986+ return res987+}988+989+// Minimize tries to find a model that minimizes the weight of the clause defined as the optimisation clause in the problem.990+// If no model can be found, it will return a cost of -1.991+// Otherwise, calling s.Model() afterwards will return the model that satisfy the formula, such that no other model with a smaller cost exists.992+// If this function is called on a non-optimization problem, it will either return -1, or a cost of 0 associated with a993+// satisfying model (ie any model is an optimal model).994+func (s *Solver) Minimize() int {995+ status := s.Solve()996+ if status == Unsat { // Problem cannot be satisfied at all997+ return -1998+ }999+ if s.minLits == nil { // No optimization clause: this is a decision problem, solution is optimal1000+ return 01001+ }1002+ maxCost := 01003+ if s.minWeights == nil {1004+ maxCost = len(s.minLits)1005+ } else {1006+ for _, w := range s.minWeights {1007+ maxCost += w1008+ }1009+ }1010+ s.hypothesis = make([]Lit, len(s.minLits))1011+ for i, lit := range s.minLits {1012+ s.hypothesis[i] = lit.Negation()1013+ }1014+ weights := make([]int, len(s.minWeights))1015+ copy(weights, s.minWeights)1016+ sort.Sort(wLits{lits: s.hypothesis, weights: weights})1017+ s.lastModel = make(Model, len(s.model))1018+ var cost int1019+ for status == Sat {1020+ copy(s.lastModel, s.model) // Save this model: it might be the last one1021+ cost = 01022+ for i, lit := range s.minLits {1023+ if s.model[lit.Var()] > 0 == lit.IsPositive() {1024+ if s.minWeights == nil {1025+ cost++1026+ } else {1027+ cost += s.minWeights[i]1028+ }1029+ }1030+ }1031+ if cost == 0 {1032+ return 01033+ }1034+ if s.Verbose {1035+ fmt.Printf("o %d\n", cost)1036+ }1037+ // Add a constraint incrementing current best cost1038+ lits2 := make([]Lit, len(s.minLits))1039+ weights2 := make([]int, len(s.minWeights))1040+ copy(lits2, s.hypothesis)1041+ copy(weights2, weights)1042+ s.AppendClause(NewPBClause(lits2, weights2, maxCost-cost+1))1043+ s.rebuildOrderHeap()1044+ status = s.Solve()1045+ }1046+ return cost1047+}1048+1049+// functions to sort hypothesis for pseudo-boolean minimization clause.1050+type wLits struct {1051+ lits []Lit1052+ weights []int1053+}1054+1055+func (wl wLits) Len() int { return len(wl.lits) }1056+func (wl wLits) Less(i, j int) bool { return wl.weights[i] > wl.weights[j] }1057+1058+func (wl wLits) Swap(i, j int) {1059+ wl.lits[i], wl.lits[j] = wl.lits[j], wl.lits[i]1060+ wl.weights[i], wl.weights[j] = wl.weights[j], wl.weights[i]1061+}@@ -0,0 +1,23 @@
1+package solver2+3+import "sort"4+5+// clauseSorter is a structure to facilitate the sorting of lits in a learned clause6+// according to their respective decision levels.7+type clauseSorter struct {8+ lits []Lit9+ model Model10+}11+12+func (cs *clauseSorter) Len() int { return len(cs.lits) }13+func (cs *clauseSorter) Less(i, j int) bool {14+ return abs(cs.model[cs.lits[i].Var()]) > abs(cs.model[cs.lits[j].Var()])15+}16+func (cs *clauseSorter) Swap(i, j int) { cs.lits[i], cs.lits[j] = cs.lits[j], cs.lits[i] }17+18+// sortLiterals sorts the literals depending on the decision level they were bound19+// i.e. abs(model[lits[i]]) <= abs(model[lits[i+1]]).20+func sortLiterals(lits []Lit, model []decLevel) {21+ cs := &clauseSorter{lits, model}22+ sort.Sort(cs)23+}@@ -0,0 +1,111 @@
1+package solver2+3+// Describes basic types and constants that are used in the solver4+5+// Status is the status of a given problem or clause at a given moment.6+type Status byte7+8+const (9+ // Indet means the problem is not proven sat or unsat yet.10+ Indet = Status(iota)11+ // Sat means the problem or clause is satisfied.12+ Sat13+ // Unsat means the problem or clause is unsatisfied.14+ Unsat15+ // Unit is a constant meaning the clause contains only one unassigned literal.16+ Unit17+ // Many is a constant meaning the clause contains at least 2 unassigned literals.18+ Many19+)20+21+func (s Status) String() string {22+ switch s {23+ case Indet:24+ return "INDETERMINATE"25+ case Sat:26+ return "SAT"27+ case Unsat:28+ return "UNSAT"29+ case Unit:30+ return "UNIT"31+ case Many:32+ return "MANY"33+ default:34+ panic("invalid status")35+ }36+}37+38+// Var start at 0 ; thus the CNF variable 1 is encoded as the Var 0.39+type Var int3240+41+// Lit start at 0 and are positive ; the sign is the last bit.42+// Thus the CNF literal -3 is encoded as 2 * (3-1) + 1 = 5.43+type Lit int3244+45+// IntsToLits converts a list of CNF literals to a []Lit.46+// This is a helper function as the same result can be achieved by calling47+// IntToLit several times.48+func IntsToLits(vals ...int32) []Lit {49+ res := make([]Lit, len(vals))50+ for i, val := range vals {51+ res[i] = IntToLit(val)52+ }53+ return res54+}55+56+// IntToLit converts a CNF literal to a Lit.57+func IntToLit(i int32) Lit {58+ if i < 0 {59+ return Lit(2*(-i-1) + 1)60+ }61+ return Lit(2 * (i - 1))62+}63+64+// IntToVar converts a CNF variable to a Var.65+func IntToVar(i int32) Var {66+ return Var(i - 1)67+}68+69+// Lit returns the positive Lit associated to v.70+func (v Var) Lit() Lit {71+ return Lit(v * 2)72+}73+74+// Int converts a Var to a CNF variable.75+func (v Var) Int() int32 {76+ return int32(v + 1)77+}78+79+// SignedLit returns the Lit associated to v, negated if 'signed', positive else.80+func (v Var) SignedLit(signed bool) Lit {81+ if signed {82+ return Lit(v*2) + 183+ }84+ return Lit(v * 2)85+}86+87+// Var returns the variable of l.88+func (l Lit) Var() Var {89+ return Var(l / 2)90+}91+92+// Int returns the equivalent CNF literal.93+func (l Lit) Int() int32 {94+ sign := l&1 == 195+ res := int32((l / 2) + 1)96+ if sign {97+ return -res98+ }99+ return res100+}101+102+// IsPositive is true iff l is > 0103+func (l Lit) IsPositive() bool {104+ return l%2 == 0105+}106+107+// Negation returns -l, i.e the positive version of l if it is negative,108+// or the negative version otherwise.109+func (l Lit) Negation() Lit {110+ return l ^ 1111+}@@ -0,0 +1,597 @@
1+package solver2+3+import (4+ "fmt"5+ "sort"6+)7+8+type watcher struct {9+ other Lit // Another lit from the clause10+ clause *Clause11+}12+13+// A watcherList is a structure used to store clauses and propagate unit literals efficiently.14+type watcherList struct {15+ nbMax int // Max # of learned clauses at current moment16+ idxReduce int // # of calls to reduce + 117+ wlistBin [][]watcher // For each literal, a list of binary clauses where its negation appears18+ wlist [][]watcher // For each literal, a list of non-binary clauses where its negation appears at position 1 or 219+ wlistPb [][]*Clause // For each literal, a list of PB or cardinality constraints.20+ wlistCardAMO [][]*Clause // For each literal, a list of cardinality constraints where card = length - 1, meaning any false literal propagates all others.21+ origClauses []*Clause // All the problem clauses.22+ learned []*Clause23+}24+25+// initWatcherList makes a new watcherList for the solver.26+func (s *Solver) initWatcherList(clauses []*Clause) {27+ nbMax := initNbMaxClauses28+ newClauses := make([]*Clause, len(clauses))29+ copy(newClauses, clauses)30+ s.wl = watcherList{31+ nbMax: nbMax,32+ idxReduce: 1,33+ wlistBin: make([][]watcher, s.nbVars*2),34+ wlist: make([][]watcher, s.nbVars*2),35+ wlistPb: make([][]*Clause, s.nbVars*2),36+ wlistCardAMO: make([][]*Clause, s.nbVars*2),37+ origClauses: newClauses,38+ }39+ for _, c := range clauses {40+ s.watchClause(c)41+ }42+}43+44+// Should be called when new vars are added to the problem (see Solver.newVar)45+func (s *Solver) addVarWatcherList(v Var) {46+ cnfVar := int(v.Int())47+ for i := s.nbVars; i < cnfVar; i++ {48+ s.wl.wlistBin = append(s.wl.wlistBin, nil, nil)49+ s.wl.wlist = append(s.wl.wlist, nil, nil)50+ s.wl.wlistPb = append(s.wl.wlistPb, nil, nil)51+ s.wl.wlistCardAMO = append(s.wl.wlistCardAMO, nil, nil)52+ }53+}54+55+// appendClause appends the clause without checking whether the clause is already satisfiable, unit, or unsatisfiable.56+// To perform those checks, call s.AppendClause.57+// clause is supposed to be a problem clause, not a learned one.58+func (s *Solver) appendClause(clause *Clause) {59+ s.wl.origClauses = append(s.wl.origClauses, clause)60+ // log.Printf("appending (and watching) %s", clause.PBString())61+ s.watchClause(clause)62+}63+64+// bumpNbMax increases the max nb of clauses used.65+// It is typically called after a restart.66+func (s *Solver) bumpNbMax() {67+ s.wl.nbMax += incrNbMaxClauses68+}69+70+// postponeNbMax increases the max nb of clauses used.71+// It is typically called when too many good clauses were learned and a cleaning was expected.72+func (s *Solver) postponeNbMax() {73+ s.wl.nbMax += incrPostponeNbMax74+}75+76+// Utilities for sorting according to clauses' LBD and activities.77+func (wl *watcherList) Len() int { return len(wl.learned) }78+func (wl *watcherList) Swap(i, j int) { wl.learned[i], wl.learned[j] = wl.learned[j], wl.learned[i] }79+80+func (wl *watcherList) Less(i, j int) bool {81+ ci := wl.learned[i]82+ cj := wl.learned[j]83+ lbdI := ci.lbd()84+ lbdJ := cj.lbd()85+ // Sort by lbd, break ties by activity86+ return lbdI > lbdJ || (lbdI == lbdJ && wl.learned[i].activity < wl.learned[j].activity)87+}88+89+// Watches the provided clause.90+func (s *Solver) watchClause(c *Clause) {91+ if c.PseudoBoolean() {92+ s.watchPB(c)93+ } else if card := c.Cardinality(); card > 1 {94+ if card == c.Len()+1 {95+ s.watchCardAMO(c, card)96+ } else {97+ // log.Printf("watching cardinality %s", c.PBString())98+ for i := 0; i < c.Cardinality()+1; i++ {99+ lit := c.Get(i)100+ neg := lit.Negation()101+ s.wl.wlistPb[neg] = append(s.wl.wlistPb[neg], c)102+ }103+ }104+ } else if c.Len() == 2 {105+ // log.Printf("watching binary %s", c.PBString())106+ first := c.First()107+ second := c.Second()108+ neg0 := first.Negation()109+ neg1 := second.Negation()110+ s.wl.wlistBin[neg0] = append(s.wl.wlistBin[neg0], watcher{clause: c, other: second})111+ s.wl.wlistBin[neg1] = append(s.wl.wlistBin[neg1], watcher{clause: c, other: first})112+ } else { // Regular, propositional clause113+ // log.Printf("watching regular %s", c.PBString())114+ first := c.First()115+ second := c.Second()116+ neg0 := first.Negation()117+ neg1 := second.Negation()118+ s.wl.wlist[neg0] = append(s.wl.wlist[neg0], watcher{clause: c, other: second})119+ s.wl.wlist[neg1] = append(s.wl.wlist[neg1], watcher{clause: c, other: first})120+ }121+}122+123+func (s *Solver) watchPB(c *Clause) {124+ // log.Printf("watching PB %s", c.PBString())125+ goal := c.Weight(0) + c.Cardinality() // We'll keep watching vars until the max weight at least reaches this value126+ sum := 0127+ i := 0128+ // log.Printf("goal is %d", goal)129+ for sum < goal && i < c.Len() {130+ lit := c.Get(i)131+ neg := lit.Negation()132+ s.wl.wlistPb[neg] = append(s.wl.wlistPb[neg], c)133+ c.pbData.watched[i] = true134+ sum += c.Weight(i)135+ i++136+ }137+}138+139+func (s *Solver) watchCardAMO(c *Clause, card int) {140+ // This is an AtMostOne constraint. At most one of the literals is false.141+ // Any falsified literal propagates all other lits.142+ // log.Printf("watching AMO %s", c.PBString())143+ for i := 0; i < card+1; i++ {144+ lit := c.Get(i)145+ neg := lit.Negation()146+ s.wl.wlistCardAMO[neg] = append(s.wl.wlistCardAMO[neg], c)147+ }148+}149+150+// unwatch the given learned clause.151+// NOTE: since it is only called when c.lbd() > 2, we know for sure152+// that c is not a binary clause.153+// We also know for sure this is a propositional clause, since only those are learned.154+func (s *Solver) unwatchClause(c *Clause) {155+ for i := 0; i < 2; i++ {156+ neg := c.Get(i).Negation()157+ j := 0158+ length := len(s.wl.wlist[neg])159+ // We're looking for the index of the clause.160+ // This will panic if c is not in wlist[neg], but this shouldn't happen.161+ for s.wl.wlist[neg][j].clause != c {162+ j++163+ }164+ s.wl.wlist[neg][j] = s.wl.wlist[neg][length-1]165+ s.wl.wlist[neg] = s.wl.wlist[neg][:length-1]166+ }167+}168+169+// unwatch the given learned PB constraint.170+// Note: this should only be called when c.PseudoBoolean() is true.171+func (s *Solver) unwatchPB(c *Clause) {172+ for i := 0; i < c.Len(); i++ {173+ if !c.pbData.watched[i] {174+ continue175+ }176+ neg := c.Get(i).Negation()177+ j := 0178+ length := len(s.wl.wlistPb[neg])179+ // We're looking for the index of the clause.180+ // This will panic if c is not in wlist[neg], but this shouldn't happen.181+ for s.wl.wlistPb[neg][j] != c {182+ j++183+ }184+ s.wl.wlistPb[neg][j] = s.wl.wlistPb[neg][length-1]185+ s.wl.wlistPb[neg] = s.wl.wlistPb[neg][:length-1]186+ }187+}188+189+// reduceLearned removes a few learned clauses that are deemed useless.190+func (s *Solver) reduceLearned() {191+ sort.Sort(&s.wl)192+ nbLearned := len(s.wl.learned)193+ length := nbLearned / 2194+ if s.wl.learned[length].lbd() <= 3 { // Lots of good clauses, postpone reduction195+ s.postponeNbMax()196+ }197+ nbRemoved := 0198+ for i := 0; i < length; i++ {199+ c := s.wl.learned[i]200+ if c.lbd() <= 2 || c.isLocked() {201+ continue202+ }203+ nbRemoved++204+ s.Stats.NbDeleted++205+ s.wl.learned[i] = s.wl.learned[nbLearned-nbRemoved]206+ s.unwatchClause(c)207+ }208+ nbLearned -= nbRemoved209+ s.wl.learned = s.wl.learned[:nbLearned]210+}211+212+type watcherListPB watcherList // A type synonymous to sort PB constraints a little more efficiently.213+214+func (wl *watcherListPB) Len() int { return len(wl.learned) }215+func (wl *watcherListPB) Swap(i, j int) { wl.learned[i], wl.learned[j] = wl.learned[j], wl.learned[i] }216+217+func (wl *watcherListPB) Less(i, j int) bool {218+ return wl.learned[i].activity < wl.learned[j].activity219+}220+221+func (s *Solver) reduceLearnedPB() {222+ wlpb := watcherListPB(s.wl)223+ sort.Sort(&wlpb)224+ nbLearned := len(s.wl.learned)225+ length := nbLearned / 2226+ nbRemoved := 0227+ for i := 0; i < length; i++ {228+ c := s.wl.learned[i]229+ if c.isLocked() {230+ continue231+ }232+ nbRemoved++233+ s.Stats.NbDeleted++234+ s.wl.learned[i] = s.wl.learned[nbLearned-nbRemoved]235+ s.unwatchPB(c)236+ }237+ nbLearned -= nbRemoved238+ s.wl.learned = s.wl.learned[:nbLearned]239+}240+241+// Adds the given learned clause and updates watchers.242+// If too many clauses have been learned yet, one will be removed.243+func (s *Solver) addLearned(c *Clause) {244+ s.wl.learned = append(s.wl.learned, c)245+ s.watchClause(c)246+ s.clauseBumpActivity(c)247+ if s.Certified {248+ if s.CertChan == nil {249+ fmt.Printf("%s\n", c.CNF())250+ } else {251+ s.CertChan <- c.CNF()252+ }253+ }254+}255+256+// Adds the given unit literal to the model at the top level.257+func (s *Solver) addLearnedUnit(unit Lit) {258+ s.model[unit.Var()] = lvlToSignedLvl(unit, 1)259+ if s.Certified {260+ if s.CertChan == nil {261+ fmt.Printf("%d 0\n", unit.Int())262+ } else {263+ s.CertChan <- fmt.Sprintf("%d 0", unit.Int())264+ }265+ }266+}267+268+// If l is negative, -lvl is returned. Else, lvl is returned.269+func lvlToSignedLvl(l Lit, lvl decLevel) decLevel {270+ if l.IsPositive() {271+ return lvl272+ }273+ return -lvl274+}275+276+// Removes the first occurrence of c from lst.277+// The element *must* be present into lst.278+func removeFrom(lst []*Clause, c *Clause) []*Clause {279+ i := 0280+ for lst[i] != c {281+ i++282+ }283+ last := len(lst) - 1284+ lst[i] = lst[last]285+ return lst[:last]286+}287+288+// Propagates literals in the trail starting from the ptrth, and returns a conflict clause, or nil if none arose.289+func (s *Solver) propagate(ptr int, lvl decLevel) *Clause {290+ for ptr < len(s.trail) {291+ lit := s.trail[ptr]292+ // log.Printf("propagating %d", lit.Int())293+ for _, w := range s.wl.wlistBin[lit] {294+ v2 := w.other.Var()295+ if assign := s.model[v2]; assign == 0 { // Other was unbounded: propagate296+ s.reason[v2] = w.clause297+ s.model[v2] = lvlToSignedLvl(w.other, lvl)298+ s.trail = append(s.trail, w.other)299+ } else if (assign > 0) != w.other.IsPositive() { // Conflict here300+ return w.clause301+ }302+ }303+ if confl := s.simplifyPropClauses(lit, lvl); confl != nil {304+ return confl305+ }306+ for _, c := range s.wl.wlistPb[lit] {307+ if c.PseudoBoolean() {308+ if !s.simplifyPseudoBool(c, lvl) {309+ return c310+ }311+ } else {312+ if !s.simplifyCardConstr(c, lvl) {313+ return c314+ }315+ }316+ }317+ for _, c := range s.wl.wlistCardAMO[lit] {318+ if !s.simplifyCardAMOConstr(c, lvl) {319+ return c320+ }321+ }322+ ptr++323+ }324+ // No unsat clause was met325+ return nil326+}327+328+// Unifies the given literal and returns a conflict clause, or nil if no conflict arose.329+func (s *Solver) unifyLiteral(lit Lit, lvl decLevel) *Clause {330+ s.model[lit.Var()] = lvlToSignedLvl(lit, lvl)331+ s.trail = append(s.trail, lit)332+ return s.propagate(len(s.trail)-1, lvl)333+}334+335+func (s *Solver) unifyLiterals(lits []Lit, lvl decLevel) *Clause {336+ for _, lit := range lits {337+ s.model[lit.Var()] = lvlToSignedLvl(lit, lvl)338+ s.trail = append(s.trail, lit)339+ }340+ for i := 0; i < len(lits); i++ {341+ if confl := s.propagate(len(s.trail)-len(lits)-i, lvl); confl != nil {342+ return confl343+ }344+ }345+ return nil346+}347+348+func (s *Solver) propagateUnit(c *Clause, lvl decLevel, unit Lit) {349+ // log.Printf("propagating unit %d", unit.Int())350+ v := unit.Var()351+ s.reason[v] = c352+ c.lock()353+ s.model[v] = lvlToSignedLvl(unit, lvl)354+ s.trail = append(s.trail, unit)355+}356+357+func (s *Solver) simplifyPropClauses(lit Lit, lvl decLevel) *Clause {358+ wl := s.wl.wlist[lit]359+ j := 0360+ for i, w := range wl {361+ if s.litStatus(w.other) == Sat { // blocking literal is SAT? Don't explore clause!362+ wl[j] = w363+ j++364+ continue365+ }366+ c := w.clause367+ // make sure c.Second() is lit368+ if c.First() == lit.Negation() {369+ c.swap(0, 1)370+ }371+ w2 := watcher{clause: c, other: c.First()}372+ firstStatus := s.litStatus(c.First())373+ if firstStatus == Sat { // Clause is already sat374+ wl[j] = w2375+ j++376+ } else {377+ found := false378+ for k := 2; k < c.Len(); k++ {379+ if litK := c.Get(k); s.litStatus(litK) != Unsat {380+ c.swap(1, k)381+ neg := litK.Negation()382+ s.wl.wlist[neg] = append(s.wl.wlist[neg], w2)383+ found = true384+ break385+ }386+ }387+ if !found { // No free lit found: unit propagation or UNSAT388+ wl[j] = w2389+ j++390+ if firstStatus == Unsat {391+ copy(wl[j:], wl[i+1:]) // Keep remaining clauses392+ s.wl.wlist[lit] = wl[:len(wl)-((i+1)-j)]393+ return c394+ }395+ s.propagateUnit(c, lvl, c.First())396+ }397+ }398+ }399+ s.wl.wlist[lit] = wl[:j]400+ return nil401+}402+403+// simplifyCardConstr simplifies a constraint of cardinality > 1, but with all weights = 1.404+// returns false iff the clause cannot be satisfied.405+func (s *Solver) simplifyCardConstr(clause *Clause, lvl decLevel) bool {406+ length := clause.Len()407+ card := clause.Cardinality()408+ nbTrue := 0409+ nbFalse := 0410+ nbUnb := 0411+ for i := 0; i < length; i++ {412+ lit := clause.Get(i)413+ switch s.litStatus(lit) {414+ case Indet:415+ nbUnb++416+ case Sat:417+ nbTrue++418+ if nbTrue == card {419+ return true420+ }421+ case Unsat:422+ nbFalse++423+ if length-nbFalse < card {424+ return false425+ }426+ }427+ if nbUnb+nbTrue > card {428+ break429+ }430+ }431+ if nbUnb+nbTrue == card {432+ // All unbounded lits must be bound to make the clause true433+ i := 0434+ for nbUnb > 0 {435+ lit := clause.Get(i)436+ if s.model[lit.Var()] == 0 {437+ s.propagateUnit(clause, lvl, lit)438+ nbUnb--439+ } else {440+ i++441+ }442+ }443+ return true444+ }445+ s.swapFalse(clause)446+ return true447+}448+449+// simplifyCardAMOConstr simplifies the special cardinality constraints where card == length -1, and returns false iff the constraint is UNSAT.450+// Whenever a literal is false, all other literals must be true.451+// This is a special case, which can be dealt with slightly more efficiently than more general cases.452+func (s *Solver) simplifyCardAMOConstr(clause *Clause, lvl decLevel) bool {453+ card := clause.Cardinality()454+ length := card + 1455+ foundFalse := false456+ for i := 0; i < length; i++ {457+ lit := clause.Get(i)458+ if s.litStatus(lit) == Unsat {459+ if foundFalse { // A second false lit460+ return false461+ }462+ foundFalse = true463+ }464+ }465+ // All unbounded lits must be bound to make the clause true466+ for i := 0; i < length; i++ {467+ lit := clause.Get(i)468+ if s.model[lit.Var()] == 0 {469+ s.propagateUnit(clause, lvl, lit)470+ }471+ }472+ return true473+}474+475+// swapFalse swaps enough literals from the clause so that all watching literals are either true or unbounded lits.476+// Must only be called when there a at least cardinality + 1 true and unbounded lits.477+func (s *Solver) swapFalse(clause *Clause) {478+ card := clause.Cardinality()479+ i := 0480+ j := card + 1481+ for i < card+1 {482+ lit := clause.Get(i)483+ for s.litStatus(lit) != Unsat {484+ i++485+ if i == card+1 {486+ return487+ }488+ lit = clause.Get(i)489+ }490+ lit = clause.Get(j)491+ for s.litStatus(lit) == Unsat {492+ j++493+ lit = clause.Get(j)494+ }495+ ni := &s.wl.wlistPb[clause.Get(i).Negation()]496+ nj := &s.wl.wlistPb[clause.Get(j).Negation()]497+ clause.swap(i, j)498+ *ni = removeFrom(*ni, clause)499+ *nj = append(*nj, clause)500+ i++501+ j++502+ }503+}504+505+// slackSum returns slack value for c and whether the clause is already sat or not.506+// The slack value is defined as sum of weights - cardinality - sum of weights of falsified lits.507+// It can be negative, meaning the whole constraint is falsified.508+// If it's 0 or above, it means all literals with a weight >= slack must be propagated.509+// If the clause is already satisfied, the slack value shall not be used.510+// This is mostly useful for PB constraints.511+func (s *Solver) slackSum(c *Clause) (slack int, sat bool) {512+ card := c.Cardinality()513+ slack = -card514+ sum := 0515+ for i, w := range c.pbData.weights {516+ status := s.litStatus(c.Get(i))517+ switch status {518+ case Indet:519+ slack += w520+ case Sat:521+ slack += w522+ sum += w523+ if sum >= card {524+ return slack, true525+ }526+ }527+ }528+ return slack, false529+}530+531+// propagateAll propagates all unbounded literals from c as unit literals532+func (s *Solver) propagateAll(c *Clause, lvl decLevel) {533+ for i := 0; i < c.Len(); i++ {534+ if lit := c.Get(i); s.litStatus(lit) == Indet {535+ s.propagateUnit(c, lvl, lit)536+ }537+ }538+}539+540+func (s *Solver) simplifyPseudoBool(clause *Clause, lvl decLevel) bool {541+ foundUnit := true542+ for foundUnit {543+ slack, sat := s.slackSum(clause)544+ if sat {545+ return true546+ }547+ if slack < 0 {548+ return false549+ }550+ if slack == 0 {551+ s.propagateAll(clause, lvl)552+ return true553+ }554+ foundUnit = false555+ for i := 0; i < clause.Len(); i++ {556+ lit := clause.Get(i)557+ if s.litStatus(lit) == Indet && clause.Weight(i) > slack { // lit will be propagated558+ s.propagateUnit(clause, lvl, lit)559+ foundUnit = true560+ }561+ }562+ }563+ s.updateWatchPB(clause)564+ return true565+}566+567+func (s *Solver) updateWatchPB(clause *Clause) {568+ weightWatched := 0569+ i := 0570+ card := clause.Cardinality()571+ for weightWatched <= card && i < clause.Len() {572+ lit := clause.Get(i)573+ if s.litStatus(lit) == Unsat {574+ if clause.pbData.watched[i] {575+ ni := &s.wl.wlistPb[lit.Negation()]576+ *ni = removeFrom(*ni, clause)577+ clause.pbData.watched[i] = false578+ }579+ } else {580+ weightWatched += clause.Weight(i)581+ if !clause.pbData.watched[i] {582+ ni := &s.wl.wlistPb[lit.Negation()]583+ *ni = append(*ni, clause)584+ clause.pbData.watched[i] = true585+ }586+ }587+ i++588+ }589+ // If there are some more watched literals, they are now useless590+ for i := i; i < clause.Len(); i++ {591+ if clause.pbData.watched[i] {592+ ni := &s.wl.wlistPb[clause.Get(i).Negation()]593+ *ni = removeFrom(*ni, clause)594+ clause.pbData.watched[i] = false595+ }596+ }597+}