1 #include <chuffed/mdd/weighted_dfa.h>
2 #include <chuffed/mdd/wmdd_prop.h>
3
4
5 //#define FULL_PROP
6
7 // #define DEBUG_INFER
8 //#define DEBUG_EXPLN
9 #define INC_EXPLAIN
10 //#define WEAKNOGOOD
11
12 #define RELAX_UBC_LATER
13 #define EARLY_CUTOFF
14
mkDisjRef(vec<EdgeID> & es)15 DisjRef mkDisjRef(vec<EdgeID>& es)
16 {
17 int mem_size = sizeof(Disj) + (es.size()-1)*sizeof(EdgeID);
18 void* mem = malloc(mem_size);
19 Disj* disj = new (mem) Disj;
20
21 for(int ei = 0; ei < es.size(); ei++)
22 disj->edges[ei] = es[ei];
23 disj->sz = es.size();
24 disj->curr_sz = disj->sz;
25 return disj;
26 }
27
28 enum WatchFlag { W_ABOVE = 1, W_BELOW = 2, W_VAL = 4, W_ANY = 7 };
29 #define SET_WATCH(eid, flag) edges[(eid)].watch_flags |= (flag)
30 #define CLEAR_WATCH(eid, flag) edges[(eid)].watch_flags &= (~(flag))
31 #define CHECK_WATCH(eid, flag) (edges[(eid)].watch_flags&(flag))
32
33 #ifdef DEBUG_GRAPH
34 static int __count = 0;
35 #endif
36
37 // Flag used to indicate that an edge
38 // needs to be processed.
39 #define EDGE_PROC 8
40
41 #ifdef DEBUG_EXPLN
42 #include <iostream>
operator <<(std::ostream & out,Lit l)43 std::ostream& operator<<(std::ostream& out, Lit l)
44 {
45 if(!sign(l))
46 out << "~";
47 out << var(l);
48 return out;
49 }
50
51 template<class T>
operator <<(std::ostream & out,vec<T> & v)52 std::ostream& operator<<(std::ostream& out, vec<T>& v)
53 {
54 out << "[";
55 if(v.size() > 0)
56 out << v[0];
57 for(int ii = 1; ii < v.size(); ii++)
58 out << ", " << v[ii];
59 out << "]";
60 return out;
61 }
62 #endif
63
WMDDProp(vec<IntView<>> & _vs,IntView<> _c,vec<int> & _levels,vec<Edge> & _edges,const MDDOpts & _opts)64 WMDDProp::WMDDProp(vec< IntView<> >& _vs, IntView<> _c, vec<int>& _levels, vec<Edge>& _edges, const MDDOpts& _opts)
65 : intvars(_vs), cost(_c),
66 nodes(),
67 root(1), T(0),
68 edges(_edges),
69 dead_edges(_edges.size()),
70 fixedvars(0), // Do we want to use this, or backtrack timestamps?
71 expl_care(0),
72 opts(_opts)
73 {
74 // Attach variables
75 // Literals for each assignment [| x != i |]
76 vec< vec<int> > val_edges;
77 for( int i = 0; i < intvars.size(); i++ )
78 {
79 int min = intvars[i].getMin();
80 int max = intvars[i].getMax();
81 int dom = max - min + 1;
82 int offset = boolvars.size();
83 VarInfo info(min, offset, dom);
84 varinfo.push(info);
85
86 for( int j = min; j <= max; j++ )
87 {
88 boolvars.push(intvars[i].getLit(j,1)); // v[i] \eq j
89 boolvars.last().attach(this, boolvars.size()-1, EVENT_U);
90 val_edges.push();
91 }
92 }
93 fixedvars.growTo(val_edges.size());
94
95 // Assumptions:
96 // ============
97 // Node 0 is T.
98 // Node 1 is the root.
99 // There are no long edges.
100 vec< vec<int> > in_vec;
101 vec< vec<int> > out_vec;
102 for(int ni = 0; ni < _levels.size(); ni++)
103 {
104 in_vec.push();
105 out_vec.push();
106 }
107
108 for(int ei = 0; ei < edges.size(); ei++)
109 {
110 Edge& e(edges[ei]);
111 e.kill_flags = 0;
112 e.watch_flags = 0;
113 in_vec[e.end].push(ei);
114 out_vec[e.begin].push(ei);
115
116 VarInfo vinfo(varinfo[_levels[e.begin]]);
117 int vidx = vinfo.offset + (e.val - vinfo.min);
118 val_edges[vidx].push(ei);
119 // Rename the edge value to be the unique (var,val) index.
120 e.val = vidx;
121 }
122
123 for(int ni = 0; ni < _levels.size(); ni++)
124 {
125 // Construct each node.
126 DisjRef in_edges(mkDisjRef(in_vec[ni]));
127 if(in_edges->sz > 0)
128 {
129 SET_WATCH(in_edges->edges[0], W_BELOW);
130 } else {
131 // Only the root node has no parents.
132 assert(ni == 1);
133 }
134
135 DisjRef out_edges(mkDisjRef(out_vec[ni]));
136 if(out_edges->sz > 0)
137 {
138 SET_WATCH(out_edges->edges[0], W_ABOVE);
139 } else {
140 // Similarly, only the T node has no children.
141 assert(ni == 0);
142 }
143
144 Node curr = {
145 // Basic properties
146 _levels[ni],
147 in_edges,
148 out_edges
149 };
150 nodes.push(curr);
151
152 in_base.push(0);
153 out_base.push(0);
154 }
155 expl_care.growTo(nodes.size());
156
157 for(int vv = 0; vv < intvars.size(); vv++)
158 {
159 VarInfo vinfo(varinfo[vv]);
160 int vidx = vinfo.offset;
161 int val = vinfo.min;
162 int end = vinfo.min + vinfo.dom;
163 for(; val < end; val++, vidx++)
164 {
165 // Create the array of supports.
166 DisjRef disj(mkDisjRef(val_edges[vidx]));
167 Val vref = {
168 vv,
169 val,
170 disj
171 };
172 vals.push(vref);
173
174 // Set up the watches.
175 if(disj->sz > 0)
176 {
177 SET_WATCH(disj->edges[0], W_VAL);
178 } else {
179 if(intvars[vv].remValNotR(val))
180 intvars[vv].remVal(val);
181 }
182 }
183 }
184
185 // Also need to wake up if the maximum cost drops.
186 cost.attach(this, boolvars.size(), EVENT_U);
187
188 // Ensure all the paths are initialized.
189 bool okay = fullProp();
190 if(!okay)
191 TL_FAIL();
192
193 for(int nID = 0; nID < nodes.size(); nID++)
194 {
195 int inV = nodes[nID].in_value;
196 int outV = nodes[nID].out_value;
197 nodes[nID].in_pathC = inV;
198 nodes[nID].out_pathC = outV;
199
200 in_base[nID] = inV;
201 out_base[nID] = outV;
202 }
203
204 #ifdef DEBUG_GRAPH
205 if(__count == 0)
206 debugStateDot();
207 __count++;
208 #endif
209 }
210
211 // Shrink the constraint by eliminating any infeasible
212 // edges/nodes.
213 // Makes the assumption that in_value and out_value are up
214 // to date for all reachable nodes.
215 // FIXME: Currently doesn't update the watches.
compact()216 void WMDDProp::compact()
217 {
218 #ifdef FULL_PROP
219 int maxC = cost.getMax();
220 #endif
221 for(int nID = 0; nID < nodes.size(); nID++)
222 {
223 Node& node(nodes[nID]);
224 #ifndef FULL_PROP
225 in_base[nID] = node.in_pathC;
226 out_base[nID] = node.out_pathC;
227 #else
228 in_base[nID] = node.in_value;
229 out_base[nID] = node.out_value;
230 #endif
231
232 int jj = 0;
233 for(int ei = 0; ei < node.out->sz; ei++)
234 {
235 int eid = node.out->edges[ei];
236 #ifdef FULL_PROP
237 Edge& edge(edges[eid]);
238 if(boolvars[edge.val].isFalse()
239 || nodes[edge.end].out_value == INT_MAX
240 || node.in_value + edge.weight + nodes[edge.end].out_value > maxC)
241 continue;
242 #else
243 if(dead_edges.elem(eid))
244 continue;
245 #endif
246 node.out->edges[jj++] = eid;
247 }
248 node.out->sz = jj;
249 node.out->curr_sz = jj;
250
251 jj = 0;
252 for(int ei = 0; ei < node.in->sz; ei++)
253 {
254 int eid = node.in->edges[ei];
255 #ifdef FULL_PROP
256 Edge& edge(edges[eid]);
257 if(boolvars[edge.val].isFalse()
258 || node.out_value == INT_MAX
259 || nodes[edge.begin].in_value + edge.weight + node.out_value > maxC)
260 continue;
261 #else
262 if(dead_edges.elem(eid))
263 continue;
264 #endif
265 node.in->edges[jj++] = eid;
266 }
267 node.in->sz = jj;
268 node.in->curr_sz = jj;
269 }
270 #ifndef FULL_PROP
271 for(int vv = 0; vv < vals.size(); vv++)
272 {
273 Val& vinfo(vals[vv]);
274 int jj = 0;
275 for(int ii = 0; ii < vinfo.edges->sz; ii++)
276 {
277 if(dead_edges.elem(vinfo.edges->edges[ii]))
278 continue;
279 vinfo.edges->edges[jj++] = vinfo.edges->edges[ii];
280 }
281 vinfo.edges->sz = jj;
282 vinfo.edges->curr_sz = jj;
283 if(jj > 0)
284 {
285 int eid = vinfo.edges->edges[0];
286 if(!CHECK_WATCH(eid, W_VAL))
287 SET_WATCH(eid, W_VAL);
288 }
289 }
290 #endif
291 }
292
propagate()293 bool WMDDProp::propagate()
294 {
295 #ifdef FULL_PROP
296 bool ret = fullProp();
297 #else
298 nodes[T].status = 0;
299 nodes[root].status = 0;
300 bool ret = incProp();
301 #endif
302
303 if(sat.decisionLevel() == 0 && ret)
304 compact();
305
306 return ret;
307 }
308
309 // Full propagation algorithm.
fullProp(void)310 bool WMDDProp::fullProp(void)
311 {
312 // for(int nID = 1; nID < nodes.size(); nID++)
313 // assert(nodes[nID].status == 0);
314 for(int vv = 0; vv < vals.size(); vv++)
315 {
316 if(!boolvars[vv].isFalse())
317 vals[vv].status = VAL_UNSUPP;
318 }
319
320 clear_queue.clear();
321
322 // Ensure T is never enqueued.
323 nodes[T].status = 1;
324 nodes[T].in_value = INT_MAX;
325 nodes[T].out_value = 0;
326
327 // Walk forward, propagating the shortest path from r.
328 vec<int> stateQ;
329 stateQ.push(root);
330 int qidx = 0;
331
332 while(qidx < stateQ.size())
333 {
334 int nID = stateQ[qidx];
335 Node& node(nodes[nID]);
336
337 // for(int ei = 0; ei < node.out->sz; ei++)
338 for(int ei = 0; ei < node.out->curr_sz; ei++)
339 {
340 int eid = node.out->edges[ei];
341 Edge& e(edges[eid]);
342
343 // if(fixedvars.elem(e.val))
344 if(boolvars[e.val].isFalse())
345 continue;
346
347 Node& dest(nodes[e.end]);
348 if(!dest.status)
349 {
350 // Not yet touched.
351 dest.status = 1;
352 dest.in_value = node.in_value + e.weight;
353 stateQ.push(e.end);
354 } else {
355 dest.in_value = std::min(dest.in_value, node.in_value + e.weight);
356 }
357 }
358 qidx++;
359 }
360
361 int minC = nodes[T].in_value;
362 int maxC = cost.getMax();
363 // Check feasibility.
364 if(minC > maxC)
365 {
366 // Clear the status flags.
367 for(qidx = 0; qidx < stateQ.size(); qidx++)
368 nodes[stateQ[qidx]].status = 0;
369
370 if( so.lazy )
371 {
372 // Generate an explanation.
373 for(int vv = 0; vv < vals.size(); vv++)
374 vals[vv].status = 0;
375
376 Clause* r = explainConflict();
377 sat.confl = r;
378 }
379 return false;
380 }
381 if(cost.setMinNotR(minC))
382 {
383 Reason r = createReason((minC<<1)|1);
384 if(!cost.setMin(minC,r))
385 return false;
386 }
387
388 // Walk backwards, propagating the shortest path to T.
389 for(qidx = stateQ.size()-1; qidx >= 0; qidx--)
390 {
391 int nID = stateQ[qidx];
392 Node& node(nodes[nID]);
393
394 // Reset the status flags on the way back.
395 node.status = 0;
396 int dist_from = INT_MAX;
397
398 int sz = node.out->curr_sz;
399 // for(int ei = 0; ei < node.out->sz; ei++)
400 for(int ei = sz-1; ei >= 0; ei--)
401 {
402 int eid = node.out->edges[ei];
403
404 Edge& e(edges[eid]);
405 // if(fixedvars.elem(e.val))
406 if(boolvars[e.val].isFalse())
407 {
408 std::swap(node.out->edges[ei], node.out->edges[--sz]);
409 continue;
410 }
411
412 Node& dest(nodes[e.end]);
413 // Need to be careful adding anything to INT_MAX when it's to be
414 // stored as an integer.
415 if(dest.out_value == INT_MAX
416 || node.in_value + e.weight + dest.out_value > maxC)
417 {
418 // No point updating distances if they're infeasible.
419 std::swap(node.out->edges[ei], node.out->edges[--sz]);
420 continue;
421 }
422 vals[e.val].status &= (~VAL_UNSUPP);
423 dist_from = std::min(dist_from, dest.out_value + e.weight);
424 }
425 node.out_value = dist_from;
426 if(sz < node.out->curr_sz)
427 trailChange(node.out->curr_sz, sz);
428 }
429
430 #ifdef DEBUG_INFER
431 bool changes = false;
432 #endif
433 for(int vv = 0; vv < vals.size(); vv++)
434 {
435 if(vals[vv].status)
436 {
437 int var(vals[vv].var);
438 int val(vals[vv].val);
439 if(intvars[var].remValNotR(val))
440 {
441 #ifdef DEBUG_INFER
442 if(!changes)
443 {
444 changes = true;
445 printf("[");
446 } else {
447 printf(",");
448 }
449 printf("%d", vv);
450 #endif
451 Reason r = createReason(vv<<1);
452 if(!intvars[var].remVal(val,r))
453 return false;
454 // fixedvars.insert(vv);
455 }
456 vals[vv].status = 0;
457 }
458 }
459 #ifdef DEBUG_INFER
460 if(changes)
461 printf("]\n");
462 #endif
463
464 return true;
465 }
466
incPropDown(vec<int> & clear_queue,int maxC,vec<int> & valQ)467 void WMDDProp::incPropDown(vec<int>& clear_queue, int maxC, vec<int>& valQ)
468 {
469 // for(int ni = 0; ni < nodes.size(); ni++)
470 // nodes[ni].status = 0;
471 vec<int> downQ;
472 int clear_idx = 0;
473 int down_idx = 0;
474 int level = -1;
475 while(clear_idx < clear_queue.size() || down_idx < downQ.size())
476 {
477 level = (down_idx == downQ.size()) ? vals[clear_queue[clear_idx]].var : level+1;
478
479 // Remove any additional edges corresponding to values removed from the current variable, and
480 // enqueue the affected nodes.
481 for(; clear_idx < clear_queue.size(); clear_idx++)
482 {
483 int vv = clear_queue[clear_idx];
484 if(vals[vv].var != level)
485 break;
486
487 DisjRef es(vals[vv].edges);
488 for(int ei = 0; ei < es->sz; ei++)
489 {
490 int eid = es->edges[ei];
491 assert(eid < edges.size());
492 if(dead_edges.elem(eid))
493 continue;
494
495 dead_edges.insert(eid);
496 Edge& e(edges[eid]);
497 e.kill_flags |= EDGE_PROC;
498
499 // Enqueue e.end if the loss of e might change the shortest path.
500 if(!nodes[e.end].status)
501 {
502 nodes[e.end].status = 1;
503 downQ.push(e.end);
504 }
505 }
506 }
507
508 int down_next = downQ.size();
509 for(; down_idx < down_next; down_idx++)
510 {
511 int nID = downQ[down_idx];
512 Node& node(nodes[nID]);
513 node.status = 0;
514 #ifdef EARLY_CUTOFF
515 if(node.in_pathC + node.out_pathC > maxC)
516 continue;
517 #endif
518
519 // Compute the updated cost c(r -> n)
520 int oldC = node.in_pathC;
521 int bestC = INT_MAX;
522 for(int ei = 0; ei < node.in->sz; ei++)
523 {
524 int eid = node.in->edges[ei];
525 assert(eid < edges.size());
526 // if(dead_edges.elem(eid))
527 if(fixedvars.elem(edges[eid].val))
528 continue;
529 Edge& e(edges[eid]);
530 int eC = nodes[e.begin].in_pathC == INT_MAX ? INT_MAX : nodes[e.begin].in_pathC + e.weight;
531 if(eC < bestC)
532 {
533 bestC = eC;
534 if(bestC == oldC)
535 break;
536 }
537 }
538 // If the path length remains the same,
539 // we can stop.
540 if(bestC == oldC)
541 continue;
542
543 // Separated out this case so we don't have to check if
544 // the node is reachable every time we update an edge.
545 trailChange(node.in_pathC, bestC);
546 if(bestC == INT_MAX)
547 {
548 for(int ei = 0; ei < node.out->sz; ei++)
549 {
550 int eid = node.out->edges[ei];
551 assert(eid < edges.size());
552 if(dead_edges.elem(eid))
553 continue;
554 dead_edges.insert(eid);
555 if(CHECK_WATCH(eid, W_VAL))
556 valQ.push(edges[eid].val);
557
558 Edge& e(edges[eid]);
559
560 if(!nodes[e.end].status)
561 {
562 nodes[e.end].status = 1;
563 downQ.push(e.end);
564 }
565 }
566 } else {
567 // Propagate the affected nodes, updating
568 for(int ei = 0; ei < node.out->sz; ei++)
569 {
570 int eid = node.out->edges[ei];
571 assert(eid < edges.size());
572 if(dead_edges.elem(eid))
573 continue;
574 // If nodes[e.end].out_pathC was INT_MAX, e would already
575 // be dead; so we don't need to worry about INT_MAX-y nodes.
576 Edge& e(edges[eid]);
577 if(bestC + e.weight + nodes[e.end].out_pathC > maxC)
578 {
579 dead_edges.insert(eid);
580 if(CHECK_WATCH(eid, W_VAL))
581 valQ.push(edges[eid].val);
582 }
583
584 if(!nodes[e.end].status)
585 {
586 nodes[e.end].status = 1;
587 downQ.push(e.end);
588 }
589 }
590 }
591 }
592 }
593 }
594
incPropUp(vec<int> & clear_queue,int maxC,vec<int> & valQ)595 void WMDDProp::incPropUp(vec<int>& clear_queue, int maxC, vec<int>& valQ)
596 {
597 // for(int ni = 0; ni < nodes.size(); ni++)
598 // nodes[ni].status = 0;
599
600 vec<int> upQ;
601 int clear_idx = clear_queue.size()-1;
602 int up_idx = 0;
603 int level = nodes[T].var;
604 while(clear_idx >= 0 || up_idx < upQ.size())
605 {
606 level = (up_idx == upQ.size()) ? vals[clear_queue[clear_idx]].var : level-1;
607
608 // Remove any additional edges corresponding to values removed from the current variable, and
609 // enqueue the affected nodes.
610 for(; clear_idx >= 0; clear_idx--)
611 {
612 int vv = clear_queue[clear_idx];
613 if(vals[vv].var != level)
614 break;
615
616 DisjRef es(vals[vv].edges);
617 for(int ei = 0; ei < es->sz; ei++)
618 {
619 int eid = es->edges[ei];
620 assert(eid < edges.size());
621 // All these edges are already dead, since we've run incPropDown.
622 // We use kill-flags instead.
623 Edge& e(edges[eid]);
624 if(!(e.kill_flags&EDGE_PROC))
625 continue;
626 e.kill_flags &= (~EDGE_PROC);
627
628 // Enqueue e.begin if the loss of e might change the shortest path.
629 // FIXME: This actually needs to be the value of pathC
630 // before propagation started, not at the current time.
631 if(!nodes[e.begin].status)
632 {
633 assert(e.begin < nodes.size());
634 nodes[e.begin].status = 1;
635 upQ.push(e.begin);
636 }
637 }
638 }
639
640 int up_next = upQ.size();
641 for(; up_idx < up_next; up_idx++)
642 {
643 int nID = upQ[up_idx];
644 assert(nID < nodes.size());
645 Node& node(nodes[nID]);
646 node.status = 0;
647 #ifdef EARLY_CUTOFF
648 if(node.in_pathC + node.out_pathC > maxC)
649 continue;
650 #endif
651
652 // Compute the updated cost c(r -> n)
653 int oldC = node.out_pathC;
654 int bestC = INT_MAX;
655 for(int ei = 0; ei < node.out->sz; ei++)
656 {
657 int eid = node.out->edges[ei];
658 Edge& e(edges[eid]);
659 // if(dead_edges.elem(eid))
660 if(fixedvars.elem(e.val))
661 continue;
662 int eC = nodes[e.end].out_pathC == INT_MAX ? INT_MAX : nodes[e.end].out_pathC + e.weight;
663 // int eC = nodes[e.end].out_pathC + e.weight;
664 if(eC < bestC)
665 {
666 bestC = eC;
667 if(bestC == oldC)
668 break;
669 }
670 }
671 // If the path length remains the same,
672 // we can stop.
673 if(bestC == oldC)
674 continue;
675
676 // Separated out this case so we don't have to check if
677 // the node is reachable every time we update an edge.
678 trailChange(node.out_pathC, bestC);
679 if(bestC == INT_MAX)
680 {
681 for(int ei = 0; ei < node.in->sz; ei++)
682 {
683 int eid = node.in->edges[ei];
684 assert(eid < edges.size());
685 if(dead_edges.elem(eid))
686 continue;
687 dead_edges.insert(eid);
688 if(CHECK_WATCH(eid, W_VAL))
689 valQ.push(edges[eid].val);
690
691 Edge& e(edges[eid]);
692 if(!nodes[e.begin].status)
693 {
694 assert(e.begin < nodes.size());
695 nodes[e.begin].status = 1;
696 upQ.push(e.begin);
697 }
698 }
699 } else {
700 // Propagate the affected nodes, updating
701 for(int ei = 0; ei < node.in->sz; ei++)
702 {
703 int eid = node.in->edges[ei];
704 assert(eid < edges.size());
705 if(dead_edges.elem(eid))
706 continue;
707 // If nodes[e.end].in_pathC was INT_MAX, e would already
708 // be dead; so we don't need to worry about INT_MAX-y nodes.
709 Edge& e(edges[eid]);
710 if(bestC + e.weight + nodes[e.begin].in_pathC > maxC)
711 {
712 dead_edges.insert(eid);
713 if(CHECK_WATCH(eid, W_VAL))
714 valQ.push(edges[eid].val);
715 }
716
717 if(!nodes[e.begin].status)
718 {
719 assert(e.begin < nodes.size());
720 nodes[e.begin].status = 1;
721 upQ.push(e.begin);
722 }
723 }
724 }
725 }
726 }
727 }
728
729 // Incremental propagation algorithm
incProp(void)730 bool WMDDProp::incProp(void)
731 {
732 vec<int> valQ;
733 // Changing ub(cost) doesn't change the distance
734 // along any paths; it just eliminates some paths.
735 // So it should be sufficient to check the value watches
736 // for any edges such that
737 // c(r -> begin) + weight + c(dest -> T) > ub(cost)
738 int maxC = cost.getMax();
739 if(cost_changed)
740 {
741 // Pretty sure this should never fire, since
742 // we should have propagated lb(cost) >= c(r -> T)
743 // earlier
744 if(nodes[root].out_pathC > maxC)
745 {
746 if( so.lazy )
747 {
748 // Generate an explanation.
749 Clause* r = explainConflict();
750 sat.confl = r;
751 }
752 return false;
753 }
754
755 // Check if the watch for each value is still
756 // alive.
757 for(int vv = 0; vv < vals.size(); vv++)
758 {
759 if(boolvars[vv].isFalse())
760 continue;
761
762 // If the watch is no longer on a satisfying path, mark the
763 // value for further investigation.
764 if(edge_pathC(vals[vv].edges->edges[0]) > maxC)
765 {
766 valQ.push(vv);
767 }
768 }
769 }
770
771 if(clear_queue.size() == 0)
772 return true;
773
774 // We want to process the affected edges level-at-a-time.
775 // So, sort the values in ascending order.
776 std::sort((int*) clear_queue, ((int*) clear_queue) + clear_queue.size());
777
778 // Compute the downward phase first, computing shortest
779 // paths c(r -> n).
780 incPropDown(clear_queue, maxC, valQ);
781
782 // Check if a path c(r -> n) is still feasible.
783 int minC = nodes[T].in_pathC;
784 if(minC > maxC)
785 {
786 if( so.lazy )
787 {
788 // Generate an explanation.
789 Clause* r = explainConflict();
790 sat.confl = r;
791 }
792 return false;
793 }
794 if(cost.setMinNotR(minC))
795 {
796 Reason r = createReason((minC<<1)|1);
797 if(!cost.setMin(minC,r)) return false;
798 }
799
800 // Compute the upward phase, shortest paths c(n -> T).
801 incPropUp(clear_queue, maxC, valQ);
802
803 // Determine which values are now missing supports.
804 #ifdef DEBUG_INFER
805 bool changes = false;
806 #endif
807 std::sort((int*) valQ, ((int*) valQ) + valQ.size());
808
809 for(int vi = 0; vi < valQ.size(); vi++)
810 {
811 int vv = valQ[vi];
812 Val& vinfo(vals[vv]);
813 vinfo.status = 0;
814
815 // Search for a new watch.
816 for(int ei = 0; ei < vinfo.edges->sz; ei++)
817 {
818 int eid = vinfo.edges->edges[ei];
819 if(!dead_edges.elem(eid))
820 {
821 // Update the watches.
822 CLEAR_WATCH(vinfo.edges->edges[0], W_VAL);
823 SET_WATCH(eid, W_VAL);
824 vinfo.edges->edges[ei] = vinfo.edges->edges[0];
825 vinfo.edges->edges[0] = eid;
826 goto vwatch_found;
827 }
828 }
829
830 // No watch. Kill the val.
831 {
832 int var(vinfo.var);
833 int val(vinfo.val);
834 if(intvars[var].remValNotR(val))
835 {
836 #ifdef DEBUG_INFER
837 if(!changes)
838 {
839 changes = true;
840 printf("[");
841 } else {
842 printf(",");
843 }
844 printf("%d", vv);
845 #endif
846 Reason r = createReason(vv<<1);
847 if(!intvars[var].remVal(val,r)) return false;
848 }
849 }
850 vwatch_found:
851 continue;
852 }
853 #ifdef DEBUG_INFER
854 if(changes)
855 printf("]\n");
856 #endif
857
858 return true;
859 }
860
861 /*
862 void WMDDProp::checkIncProp(void)
863 {
864 for(int nID = 0; nID < nodes.size(); nID++)
865 {
866 if(nID != root)
867 {
868 // Check the up-cost has been updated.
869 Node& n(nodes[nID]);
870
871 int minC = INT_MAX;
872 for(int ei = 0; ei < n.in->sz; ei++)
873 {
874 if(dead_edges.elem( ))
875 }
876 }
877 if(nID != T)
878 {
879
880 }
881 }
882 }
883 */
884
885 // (Locally) minimal explanation algorithm
886 // =======================================
887 // (1) Walk the graph breadth-first (assuming no long edges),
888 // computing the shortest path to T through [|x = v|].
889 // This tells us the weakest cost which remains infeasible.
890 // (2) Walk backwards in the reverse order, annotating each
891 // node with the shortest distance from r to n
892 // which does not make r -> n -> T feasible.
893 // Note that this pass traverses *all* nodes, not just those
894 // reachable under the current assignment (cf. (1))
895 // (3) Walk forward, collecting nodes in the same manner as
896 // the MDD algorithm.
897 // GKG: May be faster to use an inline linked list
898 // instead of a queue.
899
900 // Compute the least upper bound necessary to restore
901 // feasibility.
compute_minC(int var,int val)902 int WMDDProp::compute_minC(int var, int val)
903 {
904 assert(0 && "No longer used.");
905 // for(int nID = 0; nID < nodes.size(); nID++)
906 // nodes[nID].status = 0;
907
908 vec<int> stateQ;
909 stateQ.push(root);
910 int qidx = 0;
911
912 nodes[root].in_value = 0;
913 nodes[T].in_value = INT_MAX;
914
915 while(qidx < stateQ.size())
916 {
917 int nID = stateQ[qidx];
918 Node& node(nodes[nID]);
919
920 // The path must pass through [|x = v|].
921 if(node.var == var)
922 {
923 for(int ei = 0; ei < node.out->sz; ei++)
924 {
925 int eid = node.out->edges[ei];
926 Edge& e(edges[eid]);
927 if(e.val == val)
928 {
929 Node& dest(nodes[e.end]);
930 if(!dest.status)
931 {
932 dest.status = 1;
933 dest.in_value = node.in_value + e.weight;
934 stateQ.push(e.end);
935 } else {
936 dest.in_value = std::min(dest.in_value, node.in_value + e.weight);
937 }
938 }
939 }
940 } else {
941 // For any other node, just compute the distance similarly.
942 for(int ei = 0; ei < node.out->sz; ei++)
943 {
944 int eid = node.out->edges[ei];
945 Edge& e(edges[eid]);
946
947 // if(boolvars[e.val].isFalse())
948 if(fixedvars.elem(e.val))
949 continue;
950
951 Node& dest(nodes[e.end]);
952 if(!dest.status)
953 {
954 // Not yet touched.
955 dest.status = 1;
956 dest.in_value = node.in_value + e.weight;
957 stateQ.push(e.end);
958 } else {
959 dest.in_value = std::min(dest.in_value, node.in_value + e.weight);
960 }
961 }
962 }
963 qidx++;
964 }
965
966 return nodes[T].in_value;
967 }
968
969 // Same as compute_minC, but with respect to a subset of
970 // the fixed variables.
late_minC(int var,int val)971 int WMDDProp::late_minC(int var, int val)
972 {
973 vec<int> stateQ;
974 stateQ.push(root);
975 int qidx = 0;
976
977 nodes[root].in_value = 0;
978 nodes[T].in_value = INT_MAX;
979
980 while(qidx < stateQ.size())
981 {
982 int nID = stateQ[qidx];
983 Node& node(nodes[nID]);
984 node.status = 0;
985
986 // The path must pass through [|x = v|].
987 if(node.var == var)
988 {
989 for(int ei = 0; ei < node.out->sz; ei++)
990 {
991 int eid = node.out->edges[ei];
992 Edge& e(edges[eid]);
993 if(e.val == val)
994 {
995 Node& dest(nodes[e.end]);
996 if(!dest.status)
997 {
998 dest.status = 1;
999 dest.in_value = node.in_value + e.weight;
1000 stateQ.push(e.end);
1001 } else {
1002 dest.in_value = std::min(dest.in_value, node.in_value + e.weight);
1003 }
1004 }
1005 }
1006 } else {
1007 // For any other node, just compute the distance similarly.
1008 for(int ei = 0; ei < node.out->sz; ei++)
1009 {
1010 int eid = node.out->edges[ei];
1011 Edge& e(edges[eid]);
1012
1013 // if(boolvars[e.val].isFalse())
1014 if(vals[e.val].status)
1015 continue;
1016
1017 Node& dest(nodes[e.end]);
1018 if(!dest.status)
1019 {
1020 // Not yet touched.
1021 dest.status = 1;
1022 dest.in_value = node.in_value + e.weight;
1023 stateQ.push(e.end);
1024 } else {
1025 dest.in_value = std::min(dest.in_value, node.in_value + e.weight);
1026 }
1027 }
1028 }
1029 qidx++;
1030 }
1031
1032 return nodes[T].in_value;
1033 }
1034
1035 // For each node, mark the shortest distance r -> n such that
1036 // c(r -> n -> T) > maxC under [|x = v|].
1037 // Assumption: Nodes are ordered topologically.
1038 // Also resets all node status flags.
mark_frontier(int var,int val)1039 int WMDDProp::mark_frontier(int var, int val)
1040 {
1041 assert(T == 0);
1042 assert(root == 1);
1043 nodes[T].out_value = 0;
1044
1045 for(int nID = nodes.size()-1; nID > 0; nID--)
1046 {
1047 Node& node(nodes[nID]);
1048
1049 // Reset status flags
1050 node.status = 0;
1051 int dist_from = INT_MAX;
1052
1053 if(node.var == var)
1054 {
1055 for(int ei = 0; ei < node.out->sz; ei++)
1056 {
1057 int eid = node.out->edges[ei];
1058 Edge& e(edges[eid]);
1059 if(e.val == val)
1060 {
1061 int dC = nodes[e.end].out_value;
1062 if(dC < INT_MAX)
1063 dist_from = dC + e.weight;
1064 }
1065 }
1066 } else {
1067 for(int ei = 0; ei < node.out->sz; ei++)
1068 {
1069 int eid = node.out->edges[ei];
1070
1071 Edge& e(edges[eid]);
1072 // if(boolvars[e.val].isFalse())
1073 if(fixedvars.elem(e.val))
1074 continue;
1075
1076 Node& dest(nodes[e.end]);
1077 // Need to be careful adding anything to INT_MAX when it's to be
1078 // stored as an integer.
1079 if(dest.out_value == INT_MAX)
1080 continue;
1081 dist_from = std::min(dist_from, dest.out_value + e.weight);
1082 }
1083 }
1084 node.out_value = dist_from;
1085 }
1086 nodes[T].status = 0;
1087
1088 return nodes[root].out_value;
1089 }
1090
1091 // In the presence of early-cutoff, this is not guaranteed
1092 // to be exact (and therefore minimal). However, it should
1093 // be correct.
1094 // We only need to explain nodes that have some path through
1095 // x = v; these are recorded in expl_care.
1096 #define NODE_QUEUED 1
1097 #define NODE_CHECK 2
1098 /*
1099 int WMDDProp::mark_frontier_phased(int var, int val)
1100 {
1101 assert(var >= 0);
1102 expl_care.clear();
1103
1104 // Check which set of edges we start with.
1105 Val& vinfo(vals[val]);
1106
1107 vec<int> upQ;
1108
1109 // Any nodes past the given edge have correct values
1110 for(int ei = 0; ei < vinfo.edges->sz; ei++)
1111 {
1112 EdgeID eid(vinfo.edges->edges[ei]);
1113 Edge& e(edges[eid]);
1114
1115 if(!nodes[e.begin].status)
1116 {
1117 nodes[e.begin].status = 1;
1118 nodes[e.begin].dist_from = nodes[e.end].out_pathC + e.weight;
1119 upQ.push(e.begin);
1120 } else {
1121 nodes[e.begin].dist_from = std::min(nodes[e.begin].dist_from, nodes[e.end].out_pathC + e.weight);
1122 }
1123 }
1124
1125 int qhead = 0;
1126 while(qhead < upQ.size())
1127 {
1128 int nID = upQ[qhead];
1129 Node& n(nodes[nID]);
1130
1131 // Reset the status flags.
1132 n.status = 0;
1133
1134 // We need to know which nodes to ignore.
1135 // We'll reset the status flags during collection.
1136 expl_care.insert(nID);
1137
1138 for(int ei = 0; ei < n.in->sz; ei++)
1139 {
1140 Edge& e(edges[eid]);
1141 if(!nodes[e.begin].status)
1142 {
1143 nodes[e.begin].status = 1;
1144 nodes[e.begin].dist_from = INT_MAX;
1145 upQ.push(e.begin);
1146 }
1147 if(!fixedvars.elem(e.val))
1148 {
1149 nodes[e.begin].dist_from = std::min(nodes[e.begin].dist_from, n.dist_from + e.weight);
1150 }
1151 }
1152 }
1153 }
1154 */
1155
1156 #define VAL_LOCKED 1
1157 #define VAL_WEAK 2
1158
1159 #define WEAKEN_THRESHOLD 2
1160 // Phased version of minimal explanation.
1161 // For layers [0, var), path cost is given by dest_from.
1162 // For layer [var, N), cost is given by out_pathC.
1163 //
1164 /*
1165 void WMDDProp::minimize_expln_phased(int var, int val, int maxC)
1166 {
1167 // A variable can be relaxed if its status is 0.
1168 for(int vi = 0; vi < vals.size(); vi++)
1169 vals[vi].status = 0;
1170
1171 vec<int> stateQ;
1172 stateQ.push(root);
1173 int qhead = 0;
1174
1175 nodes[root].in_value = 0;
1176
1177 int level = 0;
1178
1179 while(qhead < stateQ.size())
1180 {
1181 // Ensure that we're keeping track of the levels correctly.
1182 assert(nodes[stateQ[qhead]].var == level);
1183
1184 int qnext = stateQ.size();
1185
1186 if(level == var)
1187 {
1188 for(; qhead < qnext; qhead++)
1189 {
1190 int nID = stateQ[qhead];
1191 Node& node(nodes[nID]);
1192 for(int ei = 0; ei < node.out->sz; ei++)
1193 {
1194 int eid = node.out->edges[ei];
1195 Edge& e(edges[eid]);
1196 if(e.val == val)
1197 {
1198 if(!expl_care.elem(e.end))
1199 continue;
1200 if(!nodes[e.end].status)
1201 {
1202 nodes[e.end].status = 1;
1203 stateQ.push(e.end);
1204 nodes[e.end].in_value = node.in_value + e.weight;
1205 } else {
1206 nodes[e.end].in_value =
1207 std::min(nodes[e.end].in_value, node.in_value + e.weight);
1208 }
1209 }
1210 }
1211 node.status = 0;
1212 }
1213 } else if(level < var) {
1214 // Determine which values must be locked.
1215 for(int qidx = qhead; qidx < qnext; qidx++)
1216 {
1217 int nID = stateQ[qidx];
1218 Node& node(nodes[nID]);
1219 for(int ei = 0; ei < node.out->sz; ei++)
1220 {
1221 int eid = node.out->edges[ei];
1222 Edge& e(edges[eid]);
1223 int dC = nodes[e.end].out_value;
1224 if(dC < INT_MAX && node.in_value + e.weight + dC <= maxC)
1225 {
1226 assert(boolvars[e.val].isFalse());
1227 vals[e.val].status = VAL_LOCKED;
1228 }
1229 }
1230 }
1231
1232 // Add any nodes along unlocked edges to the next level.
1233 for(; qhead < qnext; qhead++)
1234 {
1235 int nID = stateQ[qhead];
1236 Node& node(nodes[nID]);
1237 for(int ei = 0; ei < node.out->sz; ei++)
1238 {
1239 int eid = node.out->edges[ei];
1240 Edge& e(edges[eid]);
1241 if(!vals[e.val].status)
1242 {
1243 if(!nodes[e.end].status)
1244 {
1245 nodes[e.end].status = 1;
1246 stateQ.push(e.end);
1247 nodes[e.end].in_value = node.in_value + e.weight;
1248 } else {
1249 nodes[e.end].in_value =
1250 std::min(nodes[e.end].in_value, node.in_value + e.weight);
1251 }
1252 }
1253 }
1254 // Clear the status flags.
1255 node.status = 0;
1256 }
1257 } else {
1258 // Level > var.
1259 }
1260 level++;
1261 }
1262 }
1263 */
1264
1265 // Walk the graph var-at-a-time, unfixing any values which will not
1266 // expose a path to T of length <= maxC.
1267 // Assumption:
1268 // - We've just run mark_frontier. As such, at each node:
1269 // + n.out_value is the length of the shortest path n -> T
1270 // + (n.out_value = INT_MAX and n.in_value = 0) or n.in_value + n.out_value = maxC
minimize_expln(int var,int val,int maxC)1271 void WMDDProp::minimize_expln(int var, int val, int maxC)
1272 {
1273 // A variable can be relaxed if its status is 0.
1274 for(int vi = 0; vi < vals.size(); vi++)
1275 vals[vi].status = 0;
1276
1277 vec<int> stateQ;
1278 stateQ.push(root);
1279 int qhead = 0;
1280
1281 nodes[root].in_value = 0;
1282
1283 int level = 0;
1284
1285 while(qhead < stateQ.size())
1286 {
1287 // Ensure that we're keeping track of the levels correctly.
1288 assert(nodes[stateQ[qhead]].var == level);
1289
1290 int qnext = stateQ.size();
1291
1292 if(level == var)
1293 {
1294 for(; qhead < qnext; qhead++)
1295 {
1296 int nID = stateQ[qhead];
1297 Node& node(nodes[nID]);
1298 for(int ei = 0; ei < node.out->sz; ei++)
1299 {
1300 int eid = node.out->edges[ei];
1301 Edge& e(edges[eid]);
1302 if(e.val == val)
1303 {
1304 if(!nodes[e.end].status)
1305 {
1306 nodes[e.end].status = 1;
1307 stateQ.push(e.end);
1308 nodes[e.end].in_value = node.in_value + e.weight;
1309 } else {
1310 nodes[e.end].in_value =
1311 std::min(nodes[e.end].in_value, node.in_value + e.weight);
1312 }
1313 }
1314 }
1315 node.status = 0;
1316 }
1317 } else {
1318 // Determine which values must be locked.
1319 #ifdef WEAKNOGOOD
1320 int val_count = 0;
1321 #endif
1322 for(int qidx = qhead; qidx < qnext; qidx++)
1323 {
1324 int nID = stateQ[qidx];
1325 Node& node(nodes[nID]);
1326 if(node.in_value + out_base[nID] > maxC)
1327 continue;
1328 for(int ei = 0; ei < node.out->sz; ei++)
1329 {
1330 int eid = node.out->edges[ei];
1331 Edge& e(edges[eid]);
1332 int dC = nodes[e.end].out_value;
1333 if(dC < INT_MAX && node.in_value + e.weight + dC <= maxC)
1334 {
1335 assert(boolvars[e.val].isFalse());
1336 #ifdef WEAKNOGOOD
1337 if(!vals[e.val].status)
1338 val_count++;
1339 #endif
1340 vals[e.val].status = VAL_LOCKED;
1341 }
1342 }
1343 }
1344
1345 #ifdef WEAKNOGOOD
1346 bool invert = (val_count > WEAKEN_THRESHOLD) && intvars[level].isFixed();
1347 if(invert)
1348 {
1349 VarInfo vinfo(varinfo[level]);
1350 int lockval = vinfo.offset + (intvars[level].getVal() - vinfo.min);
1351 int vidx = vinfo.offset;
1352 int end = vinfo.offset + vinfo.dom;
1353 for(; vidx < end; vidx++)
1354 {
1355 vals[vidx].status = (vidx == lockval) ? VAL_WEAK : (VAL_LOCKED|VAL_WEAK);
1356 }
1357 }
1358 #endif
1359 // Add any nodes along unlocked edges to the next level.
1360 for(; qhead < qnext; qhead++)
1361 {
1362 int nID = stateQ[qhead];
1363 Node& node(nodes[nID]);
1364 // Clear the status flags.
1365 node.status = 0;
1366
1367 if(node.in_value + out_base[nID] > maxC)
1368 continue;
1369 for(int ei = 0; ei < node.out->sz; ei++)
1370 {
1371 int eid = node.out->edges[ei];
1372 Edge& e(edges[eid]);
1373 if(!vals[e.val].status&VAL_LOCKED)
1374 {
1375 if(!nodes[e.end].status)
1376 {
1377 nodes[e.end].status = 1;
1378 stateQ.push(e.end);
1379 nodes[e.end].in_value = node.in_value + e.weight;
1380 } else {
1381 nodes[e.end].in_value =
1382 std::min(nodes[e.end].in_value, node.in_value + e.weight);
1383 }
1384 }
1385 }
1386 }
1387 }
1388 level++;
1389 }
1390 }
1391
collect_lits(vec<Lit> & expln)1392 void WMDDProp::collect_lits(vec<Lit>& expln)
1393 {
1394 // Determine the set of literals that must be included in the explanation.
1395 for(int vv = 0; vv < vals.size(); vv++)
1396 {
1397 Val& vinfo(vals[vv]);
1398 #ifndef WEAKNOGOOD
1399 if(vinfo.status)
1400 {
1401 expln.push(intvars[vinfo.var].getLit(vinfo.val, 1));
1402 vinfo.status = 0;
1403 }
1404 #else
1405 switch(vinfo.status)
1406 {
1407 case VAL_LOCKED:
1408 // vinfo.status is set to 2 if the literal is negated
1409 expln.push(intvars[vinfo.var].getLit(vinfo.val, 1));
1410 break;
1411 case VAL_WEAK:
1412 expln.push(intvars[vinfo.var].getLit(vinfo.val, 0));
1413 break;
1414 default:
1415 break;
1416 }
1417 vinfo.status = 0;
1418 #endif
1419 }
1420 }
1421
explainConflict(void)1422 Clause* WMDDProp::explainConflict(void)
1423 {
1424 vec<Lit> expln;
1425
1426 #ifndef RELAX_UBC_LATER
1427 int maxC = mark_frontier(-1, -1);
1428 #if 1
1429 if(maxC < INT_MAX)
1430 {
1431 expln.push(cost.getLit(maxC, 2));
1432 maxC--; // Force the assignment to be infeasible.
1433 }
1434 #else
1435 maxC = cost.getMax();
1436 expln.push(cost.getMaxLit());
1437 #endif
1438
1439 minimize_expln(-1, -1, maxC);
1440 collect_lits(expln);
1441 #else
1442 int currC = cost.getMax();
1443 mark_frontier(-1, -1);
1444 minimize_expln(-1, -1, currC);
1445
1446 int maxC = late_minC(-1, -1);
1447 if(maxC < INT_MAX)
1448 expln.push(cost.getLit(maxC, 2));
1449 collect_lits(expln);
1450 #endif
1451
1452 /*
1453 for(int vv = 0; vv < vals.size(); vv++)
1454 assert(vals[vv].status == 0);
1455 for(int nID = 1; nID < nodes.size(); nID++)
1456 assert(nodes[nID].status == 0);
1457 */
1458
1459 // if(++count == 2)
1460 // debugStateDot();
1461
1462 Clause* r = Reason_new(expln.size());
1463 for( int ii = 0; ii < expln.size(); ii++ )
1464 (*r)[ii] = expln[ii];
1465
1466 return r;
1467 }
1468
explain(Lit p,int inf)1469 Clause* WMDDProp::explain(Lit p, int inf)
1470 {
1471 // Backtrack to the correct time.
1472 /*
1473 Pinfo pi = p_info[inf];
1474 engine.btToPos(pi.btpos);
1475 int tag = pi.tag;
1476 */
1477 int tag = inf;
1478
1479 vec<Lit> expln;
1480 expln.push();
1481
1482 if(tag&1)
1483 {
1484 if(opts.expl_alg == MDDOpts::E_GREEDY)
1485 {
1486 vec<int> upQ;
1487 nodes[T].status = (tag>>1);
1488 assert(nodes[T].status > 0);
1489 assert(nodes[T].status <= nodes[T].in_pathC);
1490 upQ.push(T);
1491 incExplainUp(upQ, expln);
1492 } else {
1493 // Explaining lb(c) >= k
1494 int maxC = (tag>>1)-1;
1495 mark_frontier(-1, -1);
1496 minimize_expln(-1, -1, maxC);
1497 collect_lits(expln);
1498 }
1499 } else {
1500 int val = (tag>>1);
1501 int var = vals[val].var;
1502
1503 if(opts.expl_alg == MDDOpts::E_GREEDY)
1504 return incExplain(p, var, val);
1505
1506 #ifndef RELAX_UBC_LATER
1507 int maxC = mark_frontier(var, val);
1508 #if 1
1509 if(maxC < INT_MAX)
1510 {
1511 expln.push(cost.getLit(maxC, 2));
1512 maxC--;
1513 }
1514 #else
1515 maxC = cost.getMax();
1516 expln.push(cost.getMaxLit());
1517 #endif
1518
1519 minimize_expln(var, val, maxC);
1520 collect_lits(expln);
1521 #else
1522 // Relax ub(c) later.
1523 int currC = cost.getMax();
1524 mark_frontier(var, val);
1525 minimize_expln(var, val, currC);
1526
1527 int maxC = late_minC(var, val);
1528 if(maxC < INT_MAX)
1529 expln.push(cost.getLit(maxC, 2));
1530 collect_lits(expln);
1531 #endif
1532 }
1533
1534 #ifdef DEBUG_EXPLN
1535 expln[0] = p;
1536 std::cout << expln << std::endl;
1537 #endif
1538
1539 if(opts.expl_strat == MDDOpts::E_KEEP)
1540 {
1541 expln[0] = p;
1542 Clause* c = Clause_new(expln, true);
1543 c->learnt = true;
1544 sat.addClause(*c);
1545 return c;
1546 } else {
1547 Clause* r = Reason_new(expln.size());
1548 for( int ii = 1; ii < expln.size(); ii++ )
1549 (*r)[ii] = expln[ii];
1550 (*r)[0] = p;
1551 return r;
1552 }
1553 }
1554
1555 // Incremental explanation
1556 // Starts at the set of edges to be explained, then
1557 // walks only as far as is needed for the explanation.
1558
1559 // Precondition:
1560 // For each node, in_pathC and out_pathC is correct,
1561 // and status is 0.
1562 // FIXME: We need to store the unconstrained (level 0)
1563 // cost to/from each node.
1564 // Currently assuming this is stored in in_value
1565 // and out_value.
incExplain(Lit p,int var,int val)1566 Clause* WMDDProp::incExplain(Lit p, int var, int val)
1567 {
1568 // Check which set of edges we need to explain
1569 Val& vinfo(vals[val]);
1570
1571 vec<Lit> expln;
1572 expln.push(p);
1573 expln.push(cost.getMaxLit());
1574
1575 vec<int> upQ;
1576 vec<int> downQ;
1577
1578 int minC = cost.getMax()+1;
1579 for(int ei = 0; ei < vinfo.edges->sz; ei++)
1580 {
1581 EdgeID eid(vinfo.edges->edges[ei]);
1582 Edge& e(edges[eid]);
1583
1584 // Compute the amount by which the original minimum cost must be exceeded.
1585 // NOTE: Assumes that inBase and outBase are smallish (not INT_MAX).
1586 int inBase = in_base[e.begin];
1587 int outBase = out_base[e.end];
1588 int surplus = minC - (inBase + outBase + e.weight);
1589
1590 // If the edge was infeasible to begin with, ignore it.
1591 if(surplus <= 0)
1592 continue;
1593
1594 // In our generated explanation, we need to maintain
1595 // the invariant that c(r, s) + w + c(s, T) > maxC.
1596
1597 // Compute the minimum cost through the edge
1598 int inC = nodes[e.begin].in_pathC;
1599 int outC = nodes[e.end].out_pathC;
1600
1601 // How are we going to allocate the minC cost?
1602 int in_part = 0;
1603 int out_part = 0;
1604 if(inC == INT_MAX)
1605 {
1606 in_part = minC - (e.weight + outBase);
1607 } else if(outC == INT_MAX) {
1608 out_part = minC - (e.weight + inBase);
1609 } else {
1610 assert(inC + e.weight + outC >= minC);
1611 // Explain as much as possible in the upward pass.
1612 in_part = std::min(inC, minC - e.weight - outBase);
1613 // Then put the remaining slack into the downward pass.
1614 out_part = std::max(0, minC - e.weight - in_part);
1615 }
1616 assert(in_part >= 0);
1617 assert(out_part >= 0);
1618 if(in_part > inBase)
1619 {
1620 if(!nodes[e.begin].status)
1621 {
1622 upQ.push(e.begin);
1623 }
1624 nodes[e.begin].status = std::max(nodes[e.begin].status, in_part);
1625 assert(nodes[e.begin].status > 0);
1626 }
1627 if(out_part > outBase)
1628 {
1629 if(!nodes[e.end].status)
1630 {
1631 downQ.push(e.end);
1632 }
1633 nodes[e.end].status = std::max(nodes[e.end].status, out_part);
1634 assert(nodes[e.end].status > 0);
1635 }
1636 }
1637
1638 incExplainUp(upQ, expln);
1639 incExplainDown(downQ, expln);
1640
1641 #ifdef DEBUG_EXPLN
1642 std::cout << expln << std::endl;
1643 #endif
1644
1645 if(opts.expl_strat == MDDOpts::E_KEEP)
1646 {
1647 expln[0] = p;
1648 Clause* c = Clause_new(expln, true);
1649 c->learnt = true;
1650 sat.addClause(*c);
1651 return c;
1652 } else {
1653 Clause* r = Reason_new(expln.size());
1654 for( int ii = 0; ii < expln.size(); ii++ )
1655 (*r)[ii] = expln[ii];
1656 return r;
1657 }
1658 }
1659
incExplainUp(vec<int> & upQ,vec<Lit> & expln)1660 void WMDDProp::incExplainUp(vec<int>& upQ, vec<Lit>& expln)
1661 {
1662 vec<int> explVals;
1663
1664 int qhead = 0;
1665
1666 while(qhead < upQ.size())
1667 {
1668 int qnext = qhead;
1669
1670 // First scan, determine which values must remain
1671 // fixed.
1672 for(; qnext < upQ.size(); qnext++)
1673 {
1674 int nID(upQ[qnext]);
1675 Node& node(nodes[nID]);
1676
1677 int node_minC = node.status;
1678
1679 for(int ei = 0; ei < node.in->sz; ei++)
1680 {
1681 int eID(node.in->edges[ei]);
1682 Edge& e(edges[eID]);
1683 if(nodes[e.begin].in_pathC == INT_MAX)
1684 continue;
1685 if(nodes[e.begin].in_pathC + e.weight < node_minC)
1686 {
1687 // Including this edge could result in a feasible path
1688 // through some edge.
1689 if(!vals[e.val].status)
1690 {
1691 assert(boolvars[e.val].isFalse());
1692 explVals.push(e.val);
1693 vals[e.val].status = 1;
1694 }
1695 }
1696 }
1697 }
1698 // Second scan; check which edges
1699 for(; qhead < qnext; qhead++)
1700 {
1701 int nID(upQ[qhead]);
1702 Node& node(nodes[nID]);
1703
1704 int node_minC = node.status;
1705
1706 for(int ei = 0; ei < node.in->sz; ei++)
1707 {
1708 int eID(node.in->edges[ei]);
1709 Edge& e(edges[eID]);
1710
1711 // If the value is included in the explanation, ignore the edge.
1712 if(vals[e.val].status)
1713 continue;
1714 int dest_minC = node_minC - e.weight;
1715
1716 // If the required cost is explained at the top level, ignore the edge.
1717 if(dest_minC <= in_base[e.begin])
1718 continue;
1719
1720 // If the node already has a greater required cost, no need to update it.
1721 if(dest_minC <= nodes[e.begin].status)
1722 continue;
1723
1724 if(!nodes[e.begin].status)
1725 {
1726 upQ.push(e.begin);
1727 }
1728 nodes[e.begin].status = node_minC - e.weight;
1729 assert(nodes[e.begin].status > 0);
1730 }
1731
1732 // Clear the status flag on node.
1733 node.status = 0;
1734 }
1735 }
1736
1737 // Add the freshly fixed values to the explanation,
1738 // and reset the status flags.
1739 for(int vi = 0; vi < explVals.size(); vi++)
1740 {
1741 int vv = explVals[vi];
1742 Val& vinfo(vals[vv]);
1743 expln.push(intvars[vinfo.var].getLit(vinfo.val, 1));
1744
1745 vinfo.status = 0;
1746 }
1747 }
1748
incExplainDown(vec<int> & downQ,vec<Lit> & expln)1749 void WMDDProp::incExplainDown(vec<int>& downQ, vec<Lit>& expln)
1750 {
1751 vec<int> explVals;
1752
1753 int qhead = 0;
1754
1755 while(qhead < downQ.size())
1756 {
1757 int qnext = qhead;
1758
1759 // First scan, determine which values must remain
1760 // fixed.
1761 for(; qnext < downQ.size(); qnext++)
1762 {
1763 int nID(downQ[qnext]);
1764 Node& node(nodes[nID]);
1765
1766 int node_minC = node.status;
1767
1768 for(int ei = 0; ei < node.out->sz; ei++)
1769 {
1770 int eID(node.out->edges[ei]);
1771 Edge& e(edges[eID]);
1772 if(nodes[e.end].out_pathC == INT_MAX)
1773 continue;
1774 if(nodes[e.end].out_pathC + e.weight < node_minC)
1775 {
1776 // Including this edge could result in a feasible path
1777 // through some edge.
1778 if(!vals[e.val].status)
1779 {
1780 assert(boolvars[e.val].isFalse());
1781 explVals.push(e.val);
1782 vals[e.val].status = 1;
1783 }
1784 }
1785 }
1786 }
1787 // Second scan; check which edges
1788 for(; qhead < qnext; qhead++)
1789 {
1790 int nID(downQ[qhead]);
1791 Node& node(nodes[nID]);
1792
1793 int node_minC = node.status;
1794
1795 bool flow = false;
1796 for(int ei = 0; ei < node.out->sz; ei++)
1797 {
1798 int eID(node.out->edges[ei]);
1799 Edge& e(edges[eID]);
1800
1801 // If the value is included in the explanation, ignore the edge.
1802 if(vals[e.val].status)
1803 {
1804 flow = true;
1805 continue;
1806 }
1807 int dest_minC = node_minC - e.weight;
1808
1809 // If the required cost is explained at the top level, ignore the edge.
1810 if(dest_minC <= out_base[e.end])
1811 continue;
1812
1813 flow = true;
1814
1815 // If the node already has a greater required cost, no need to update it.
1816 if(dest_minC <= nodes[e.end].status)
1817 continue;
1818
1819 if(!nodes[e.end].status)
1820 {
1821 downQ.push(e.end);
1822 }
1823 nodes[e.end].status = node_minC - e.weight;
1824 assert(nodes[e.end].status > 0);
1825 }
1826 assert(flow);
1827
1828 // Clear the status flag on node.
1829 node.status = 0;
1830 }
1831 }
1832
1833 // Add the freshly fixed values to the explanation,
1834 // and reset the status flags.
1835 for(int vi = 0; vi < explVals.size(); vi++)
1836 {
1837 int vv = explVals[vi];
1838 Val& vinfo(vals[vv]);
1839 expln.push(intvars[vinfo.var].getLit(vinfo.val, 1));
1840
1841 vinfo.status = 0;
1842 }
1843 }
1844
1845 // Output the current state of the graph
1846 // in dot format.
debugStateDot(void)1847 void WMDDProp::debugStateDot(void)
1848 {
1849 printf("digraph ingraph { graph [ranksep=\"1.0 equally\"] \n");
1850
1851 int nID = 1;
1852
1853 while(nID < nodes.size())
1854 {
1855 // Scan through the nodes in each level.
1856 int level = nodes[nID].var;
1857 for(; nID < nodes.size() && nodes[nID].var == level; nID++)
1858 {
1859 Node& node(nodes[nID]);
1860
1861 // printf(" { node [shape=record label=\"{<prefix>%d: (%d, %d) | {", nID, node.in_value, node.out_value);
1862 printf(" { node [shape=record label=\"{<prefix>%d: (%d, %d) | {", nID, node.in_pathC, node.out_pathC);
1863 bool first = true;
1864 for(int ei = 0; ei < node.out->sz; ei++)
1865 {
1866 int eid = node.out->edges[ei];
1867 Edge& edge(edges[eid]);
1868
1869 if(first)
1870 first = false;
1871 else
1872 printf("|");
1873 printf("<p%d>%d(%d)", eid, edge.val, edge.weight);
1874 if(boolvars[edge.val].isFalse())
1875 printf("X");
1876 }
1877 printf("} }\"] %d };\n", nID);
1878 }
1879 }
1880 for(int eid = 0; eid < edges.size(); eid++)
1881 {
1882 printf("\t%d:p%d -> %d;\n", edges[eid].begin, eid, edges[eid].end);
1883 }
1884 printf("};\n");
1885 }
1886
1887 // Ergh. It feels like I'm sort of triple-handling. Because we traverse the graph
1888 // to extract the set of edges, and then traverse the edges to create the propagator.
1889 // WARNING: At no point have we actually checked to ensure there aren't any long edges.
1890 // Or, indeed, that there aren't any repetitions or inversions in variable order.
evgraph_to_wmdd(vec<IntVar * > _vs,IntVar * _cost,EVLayerGraph & g,EVLayerGraph::NodeID rootID,const MDDOpts & opts)1891 WMDDProp* evgraph_to_wmdd(vec<IntVar*> _vs, IntVar* _cost, EVLayerGraph& g, EVLayerGraph::NodeID rootID, const MDDOpts& opts)
1892 {
1893 int nNodes = g.traverse(rootID);
1894
1895 // Level for each node.
1896 vec<int> level;
1897
1898 for(int ni = 0; ni < nNodes; ni++)
1899 {
1900 // We can't initialize the levels until we traverse.
1901 level.push(0);
1902 }
1903
1904 // Traversal currently skips T.
1905 // Set up the level indicator.
1906 level[0] = _vs.size();
1907
1908 // T is 0, root is 1.
1909 vec<Edge> edges;
1910 for(EVNode curr = g.travBegin(); curr != g.travEnd(); ++curr)
1911 {
1912 level[curr.id()] = curr.var();
1913 for(int ei = 0; ei < curr.size(); ei++)
1914 {
1915 EVEdge edge(curr[ei]);
1916
1917 Edge e = {
1918 edge.val,
1919 edge.weight,
1920 curr.id(),
1921 edge.dest.id()
1922 };
1923 edges.push(e);
1924 }
1925 }
1926
1927 // Set up the views.
1928 vec< IntView<> > vs;
1929 for (int i = 0; i < _vs.size(); i++) _vs[i]->specialiseToEL();
1930 for (int i = 0; i < _vs.size(); i++) vs.push(IntView<>(_vs[i],1,0));
1931 IntView<> cost(_cost, 1, 0);
1932
1933 // Create the propagator.
1934 return new WMDDProp(vs, cost, level, edges, opts);
1935 }
1936