1 /*
2 * Revision Control Information
3 *
4 * $Source$
5 * $Author$
6 * $Revision$
7 * $Date$
8 *
9 */
10 #include "mincov_int.h"
11
12 ABC_NAMESPACE_IMPL_START
13
14
15 /*
16 * mincov.c
17 */
18
19 #define USE_GIMPEL
20 #define USE_INDEP_SET
21
22 static int select_column();
23 static void select_essential();
24 static int verify_cover();
25
26 #define fail(why) {\
27 (void) fprintf(stderr, "Fatal error: file %s, line %d\n%s\n",\
28 __FILE__, __LINE__, why);\
29 (void) fflush(stdout);\
30 abort();\
31 }
32
33 sm_row *
sm_minimum_cover(A,weight,heuristic,debug_level)34 sm_minimum_cover(A, weight, heuristic, debug_level)
35 sm_matrix *A;
36 int *weight;
37 int heuristic; /* set to 1 for a heuristic covering */
38 int debug_level; /* how deep in the recursion to provide info */
39 {
40 stats_t stats;
41 solution_t *best, *select;
42 sm_row *prow, *sol;
43 sm_col *pcol;
44 sm_matrix *dup_A;
45 int nelem, bound;
46 double sparsity;
47
48 /* Avoid sillyness */
49 if (A->nrows <= 0) {
50 return sm_row_alloc(); /* easy to cover */
51 }
52
53 /* Initialize debugging structure */
54 stats.start_time = util_cpu_time();
55 stats.debug = debug_level > 0;
56 stats.max_print_depth = debug_level;
57 stats.max_depth = -1;
58 stats.nodes = 0;
59 stats.component = stats.comp_count = 0;
60 stats.gimpel = stats.gimpel_count = 0;
61 stats.no_branching = heuristic != 0;
62 stats.lower_bound = -1;
63
64 /* Check the matrix sparsity */
65 nelem = 0;
66 sm_foreach_row(A, prow) {
67 nelem += prow->length;
68 }
69 sparsity = (double) nelem / (double) (A->nrows * A->ncols);
70
71 /* Determine an upper bound on the solution */
72 bound = 1;
73 sm_foreach_col(A, pcol) {
74 bound += WEIGHT(weight, pcol->col_num);
75 }
76
77 /* Perform the covering */
78 select = solution_alloc();
79 dup_A = sm_dup(A);
80 best = sm_mincov(dup_A, select, weight, 0, bound, 0, &stats);
81 sm_free(dup_A);
82 solution_free(select);
83
84 if (stats.debug) {
85 if (stats.no_branching) {
86 (void) printf("**** heuristic covering ...\n");
87 (void) printf("lower bound = %d\n", stats.lower_bound);
88 }
89 (void) printf("matrix = %d by %d with %d elements (%4.3f%%)\n",
90 A->nrows, A->ncols, nelem, sparsity * 100.0);
91 (void) printf("cover size = %d elements\n", best->row->length);
92 (void) printf("cover cost = %d\n", best->cost);
93 (void) printf("time = %s\n",
94 util_print_time(util_cpu_time() - stats.start_time));
95 (void) printf("components = %d\n", stats.comp_count);
96 (void) printf("gimpel = %d\n", stats.gimpel_count);
97 (void) printf("nodes = %d\n", stats.nodes);
98 (void) printf("max_depth = %d\n", stats.max_depth);
99 }
100
101 sol = sm_row_dup(best->row);
102 if (! verify_cover(A, sol)) {
103 fail("mincov: internal error -- cover verification failed\n");
104 }
105 solution_free(best);
106 return sol;
107 }
108
109 /*
110 * Find the best cover for 'A' (given that 'select' already selected);
111 *
112 * - abort search if a solution cannot be found which beats 'bound'
113 *
114 * - if any solution meets 'lower_bound', then it is the optimum solution
115 * and can be returned without further work.
116 */
117
118 solution_t *
sm_mincov(A,select,weight,lb,bound,depth,stats)119 sm_mincov(A, select, weight, lb, bound, depth, stats)
120 sm_matrix *A;
121 solution_t *select;
122 int *weight;
123 int lb;
124 int bound;
125 int depth;
126 stats_t *stats;
127 {
128 sm_matrix *A1, *A2, *L, *R;
129 sm_element *p;
130 solution_t *select1, *select2, *best, *best1, *best2, *indep;
131 int pick, lb_new, debug;
132
133 /* Start out with some debugging information */
134 stats->nodes++;
135 if (depth > stats->max_depth) stats->max_depth = depth;
136 debug = stats->debug && (depth <= stats->max_print_depth);
137
138 /* Apply row dominance, column dominance, and select essentials */
139 select_essential(A, select, weight, bound);
140 if (select->cost >= bound) {
141 return NIL(solution_t);
142 }
143
144 /* See if gimpel's reduction technique applies ... */
145 #ifdef USE_GIMPEL
146 if ( weight == NIL(int)) { /* hack until we fix it */
147 if (gimpel_reduce(A, select, weight, lb, bound, depth, stats, &best)) {
148 return best;
149 }
150 }
151 #endif
152
153 #ifdef USE_INDEP_SET
154 /* Determine bound from here to final solution using independent-set */
155 indep = sm_maximal_independent_set(A, weight);
156
157 /* make sure the lower bound is monotonically increasing */
158 lb_new = MAX(select->cost + indep->cost, lb);
159 pick = select_column(A, weight, indep);
160 solution_free(indep);
161 #else
162 lb_new = select->cost + (A->nrows > 0);
163 pick = select_column(A, weight, NIL(solution_t));
164 #endif
165
166 if (depth == 0) {
167 stats->lower_bound = lb_new + stats->gimpel;
168 }
169
170 if (debug) {
171 (void) printf("ABSMIN[%2d]%s", depth, stats->component ? "*" : " ");
172 (void) printf(" %3dx%3d sel=%3d bnd=%3d lb=%3d %12s ",
173 A->nrows, A->ncols, select->cost + stats->gimpel,
174 bound + stats->gimpel, lb_new + stats->gimpel,
175 util_print_time(util_cpu_time()-stats->start_time));
176 }
177
178 /* Check for bounding based on no better solution possible */
179 if (lb_new >= bound) {
180 if (debug) (void) printf("bounded\n");
181 best = NIL(solution_t);
182
183
184 /* Check for new best solution */
185 } else if (A->nrows == 0) {
186 best = solution_dup(select);
187 if (debug) (void) printf("BEST\n");
188 if (stats->debug && stats->component == 0) {
189 (void) printf("new 'best' solution %d at level %d (time is %s)\n",
190 best->cost + stats->gimpel, depth,
191 util_print_time(util_cpu_time() - stats->start_time));
192 }
193
194
195 /* Check for a partition of the problem */
196 } else if (sm_block_partition(A, &L, &R)) {
197 /* Make L the smaller problem */
198 if (L->ncols > R->ncols) {
199 A1 = L;
200 L = R;
201 R = A1;
202 }
203 if (debug) (void) printf("comp %d %d\n", L->nrows, R->nrows);
204 stats->comp_count++;
205
206 /* Solve problem for L */
207 select1 = solution_alloc();
208 stats->component++;
209 best1 = sm_mincov(L, select1, weight, 0,
210 bound-select->cost, depth+1, stats);
211 stats->component--;
212 solution_free(select1);
213 sm_free(L);
214
215 /* Add best solution to the selected set */
216 if (best1 == NIL(solution_t)) {
217 best = NIL(solution_t);
218 } else {
219 for(p = best1->row->first_col; p != 0; p = p->next_col) {
220 solution_add(select, weight, p->col_num);
221 }
222 solution_free(best1);
223
224 /* recur for the remaining block */
225 best = sm_mincov(R, select, weight, lb_new, bound, depth+1, stats);
226 }
227 sm_free(R);
228
229 /* We've tried as hard as possible, but now we must split and recur */
230 } else {
231 if (debug) (void) printf("pick=%d\n", pick);
232
233 /* Assume we choose this column to be in the covering set */
234 A1 = sm_dup(A);
235 select1 = solution_dup(select);
236 solution_accept(select1, A1, weight, pick);
237 best1 = sm_mincov(A1, select1, weight, lb_new, bound, depth+1, stats);
238 solution_free(select1);
239 sm_free(A1);
240
241 /* Update the upper bound if we found a better solution */
242 if (best1 != NIL(solution_t) && bound > best1->cost) {
243 bound = best1->cost;
244 }
245
246 /* See if this is a heuristic covering (no branching) */
247 if (stats->no_branching) {
248 return best1;
249 }
250
251 /* Check for reaching lower bound -- if so, don't actually branch */
252 if (best1 != NIL(solution_t) && best1->cost == lb_new) {
253 return best1;
254 }
255
256 /* Now assume we cannot have that column */
257 A2 = sm_dup(A);
258 select2 = solution_dup(select);
259 solution_reject(select2, A2, weight, pick);
260 best2 = sm_mincov(A2, select2, weight, lb_new, bound, depth+1, stats);
261 solution_free(select2);
262 sm_free(A2);
263
264 best = solution_choose_best(best1, best2);
265 }
266
267 return best;
268 }
269
270 static int
select_column(A,weight,indep)271 select_column(A, weight, indep)
272 sm_matrix *A;
273 int *weight;
274 solution_t *indep;
275 {
276 register sm_col *pcol;
277 register sm_row *prow, *indep_cols;
278 register sm_element *p, *p1;
279 double w, best;
280 int best_col;
281
282 indep_cols = sm_row_alloc();
283 if (indep != NIL(solution_t)) {
284 /* Find which columns are in the independent sets */
285 for(p = indep->row->first_col; p != 0; p = p->next_col) {
286 prow = sm_get_row(A, p->col_num);
287 for(p1 = prow->first_col; p1 != 0; p1 = p1->next_col) {
288 (void) sm_row_insert(indep_cols, p1->col_num);
289 }
290 }
291 } else {
292 /* select out of all columns */
293 sm_foreach_col(A, pcol) {
294 (void) sm_row_insert(indep_cols, pcol->col_num);
295 }
296 }
297
298 /* Find the best column */
299 best_col = -1;
300 best = -1;
301
302 /* Consider only columns which are in some independent row */
303 sm_foreach_row_element(indep_cols, p1) {
304 pcol = sm_get_col(A, p1->col_num);
305
306 /* Compute the total 'value' of all things covered by the column */
307 w = 0.0;
308 for(p = pcol->first_row; p != 0; p = p->next_row) {
309 prow = sm_get_row(A, p->row_num);
310 w += 1.0 / ((double) prow->length - 1.0);
311 }
312
313 /* divide this by the relative cost of choosing this column */
314 w = w / (double) WEIGHT(weight, pcol->col_num);
315
316 /* maximize this ratio */
317 if (w > best) {
318 best_col = pcol->col_num;
319 best = w;
320 }
321 }
322
323 sm_row_free(indep_cols);
324 return best_col;
325 }
326
327 static void
select_essential(A,select,weight,bound)328 select_essential(A, select, weight, bound)
329 sm_matrix *A;
330 solution_t *select;
331 int *weight;
332 int bound; /* must beat this solution */
333 {
334 register sm_element *p;
335 register sm_row *prow, *essen;
336 int delcols, delrows, essen_count;
337
338 do {
339 /* Check for dominated columns */
340 delcols = sm_col_dominance(A, weight);
341
342 /* Find the rows with only 1 element (the essentials) */
343 essen = sm_row_alloc();
344 sm_foreach_row(A, prow) {
345 if (prow->length == 1) {
346 (void) sm_row_insert(essen, prow->first_col->col_num);
347 }
348 }
349
350 /* Select all of the elements */
351 sm_foreach_row_element(essen, p) {
352 solution_accept(select, A, weight, p->col_num);
353 /* Make sure solution still looks good */
354 if (select->cost >= bound) {
355 sm_row_free(essen);
356 return;
357 }
358 }
359 essen_count = essen->length;
360 sm_row_free(essen);
361
362 /* Check for dominated rows */
363 delrows = sm_row_dominance(A);
364
365 } while (delcols > 0 || delrows > 0 || essen_count > 0);
366 }
367
368 static int
verify_cover(A,cover)369 verify_cover(A, cover)
370 sm_matrix *A;
371 sm_row *cover;
372 {
373 sm_row *prow;
374
375 sm_foreach_row(A, prow) {
376 if (! sm_row_intersects(prow, cover)) {
377 return 0;
378 }
379 }
380 return 1;
381 }
382 ABC_NAMESPACE_IMPL_END
383
384