1package solver23import (4 "fmt"5 "sort"6 "strings"7 "time"8)910const (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)1718// Stats are statistics about the resolution of the problem.19// They are provided for information purpose only.20type 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}2930// 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".33type decLevel int3435// 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.43type Model []decLevel4445func (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}5455// A Solver solves a given problem. It is the main data structure.56type 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}8990// 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.93func New(problem *Problem) *Solver {94 if problem.Status == Unsat {95 return &Solver{status: Unsat}96 }97 nbVars := problem.NbVars9899 trailCap := nbVars100 if len(problem.Units) > trailCap {101 trailCap = len(problem.Units)102 }103104 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}137138// 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.141func (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}156157// sets initial activity for optimization variables, if any.158func (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}167168// resets polarity of optimization lits so that they are negated by default.169func (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}176177// Optim returns true iff the underlying problem is an optimization problem (rather than a satisfaction one).178func (s *Solver) Optim() bool {179 return s.minLits != nil180}181182// OutputModel outputs the model for the problem on stdout.183func (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}204205// litStatus returns whether the literal is made true (Sat) or false (Unsat) by the206// current bindings, or if it is unbounded (Indet).207func (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}217218func (s *Solver) varDecayActivity() {219 s.varInc *= 1 / s.varDecay220}221222func (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}235236// Decays each clause's activity237func (s *Solver) clauseDecayActivity() {238 s.clauseInc *= 1 / clauseDecay239}240241// Bumps the given clause's activity.242func (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}253254// Chooses an unbound literal to be tested, or -1255// if all the variables are already bound.256func (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}269270type number interface {271 int | int32 | decLevel272}273274func abs[T number](val T) T {275 if val < 0 {276 return -val277 }278 return val279}280281func min[T number](a, b T) T {282 if a < b {283 return a284 }285 return b286}287288// 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.291func (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}350351func (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}363364// Given the last learnt clause and the levels at which vars were bound,365// Returns the level to bt to and the literal to bind366func backtrackData(c *Clause, model []decLevel) (btLevel decLevel, lit Lit) {367 btLevel = abs(model[c.Get(1).Var()])368 return btLevel, c.Get(0)369}370371func (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}380381// propagate binds the given lit, propagates it and searches for a solution,382// until it is found or a restart is needed.383func (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}440441// propagateAndSearchPB performs pseudo-boolean constraint learning whenever a conflcit arises.442func (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}518519// Sets the status to unsat and do cleanup tasks.520func (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}531532// Searches until a restart is needed.533func (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}540541// Solve solves the problem associated with the solver and returns the appropriate status.542func (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}595596// 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.598func (s *Solver) Assume(lits []Lit) Status {599 s.cleanupBindings(0)600 s.trail = s.trail[:0]601 s.assumptions = make([]bool, s.nbVars)602603 for _, lit := range lits {604605 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}617618// 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.621func (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}664665// CountModels returns the total number of models for the given problem.666func (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}737738// 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.740func (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}756757func (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}770771// PBString returns a representation of the solver's state as a pseudo-boolean problem.772func (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}808809// 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).811func (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}847848// Model returns a slice that associates, to each variable, its binding.849// If s's status is not Sat, the method will panic.850func (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}860861// 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.866func (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}891892// 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.897func (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}906907// 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.910func (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}988989// 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).994func (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}10481049// functions to sort hypothesis for pseudo-boolean minimization clause.1050type wLits struct {1051 lits []Lit1052 weights []int1053}10541055func (wl wLits) Len() int { return len(wl.lits) }1056func (wl wLits) Less(i, j int) bool { return wl.weights[i] > wl.weights[j] }10571058func (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}