/farmR/src/glpscg.c

  1/* glpscg.c */
  2
  3/***********************************************************************
  4*  This code is part of GLPK (GNU Linear Programming Kit).
  5*
  6*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
  7*  Department for Applied Informatics, Moscow Aviation Institute,
  8*  Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
  9*
 10*  GLPK is free software: you can redistribute it and/or modify it
 11*  under the terms of the GNU General Public License as published by
 12*  the Free Software Foundation, either version 3 of the License, or
 13*  (at your option) any later version.
 14*
 15*  GLPK is distributed in the hope that it will be useful, but WITHOUT
 16*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 17*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
 18*  License for more details.
 19*
 20*  You should have received a copy of the GNU General Public License
 21*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
 22***********************************************************************/
 23
 24#include "glpscg.h"
 25
 26#define _GLPSCG_DEBUG 0
 27
 28/**********************************************************************/
 29
 30SCG *scg_create_graph(int n)
 31{     /* create cliqued graph */
 32      SCG *g;
 33      xassert(n >= 0);
 34      g = xmalloc(sizeof(SCG));
 35      g->pool = dmp_create_pool();
 36      g->n_max = 50;
 37      g->nc_max = 10;
 38      g->n = g->nc = 0;
 39      g->i_ptr = xcalloc(1+g->n_max, sizeof(SCGRIB *));
 40      g->j_ptr = xcalloc(1+g->n_max, sizeof(SCGRIB *));
 41      g->c_ptr = xcalloc(1+g->nc_max, sizeof(SCGCQE *));
 42      g->v_ptr = xcalloc(1+g->n_max, sizeof(SCGCQE *));
 43      g->flag = xcalloc(1+g->n_max, sizeof(char));
 44      if (n > 0) scg_add_nodes(g, n);
 45      return g;
 46}
 47
 48/**********************************************************************/
 49
 50int scg_add_nodes(SCG *g, int num)
 51{     /* add new nodes to cliqued graph */
 52      int n_new, i;
 53      xassert(num > 0);
 54      /* determine new number of nodes */
 55      n_new = g->n + num;
 56      xassert(n_new > 0);
 57      /* enlarge the room, if necessary */
 58      if (g->n_max < n_new)
 59      {  void **save;
 60         while (g->n_max < n_new)
 61         {  g->n_max += g->n_max;
 62            xassert(g->n_max > 0);
 63         }
 64         save = (void **)g->i_ptr;
 65         g->i_ptr = xcalloc(1+g->n_max, sizeof(SCGRIB *));
 66         memcpy(&g->i_ptr[1], &save[1], g->n * sizeof(SCGRIB *));
 67         xfree(save);
 68         save = (void **)g->j_ptr;
 69         g->j_ptr = xcalloc(1+g->n_max, sizeof(SCGRIB *));
 70         memcpy(&g->j_ptr[1], &save[1], g->n * sizeof(SCGRIB *));
 71         xfree(save);
 72         save = (void **)g->v_ptr;
 73         g->v_ptr = xcalloc(1+g->n_max, sizeof(SCGCQE *));
 74         memcpy(&g->v_ptr[1], &save[1], g->n * sizeof(SCGCQE *));
 75         xfree(save);
 76         xfree(g->flag);
 77         g->flag = xcalloc(1+g->n_max, sizeof(char));
 78         memset(&g->flag[1], 0, g->n);
 79      }
 80      /* add new nodes to the end of the node list */
 81      for (i = g->n+1; i <= n_new; i++)
 82      {  g->i_ptr[i] = NULL;
 83         g->j_ptr[i] = NULL;
 84         g->v_ptr[i] = NULL;
 85         g->flag[i] = 0;
 86      }
 87      /* set new number of nodes */
 88      g->n = n_new;
 89      /* return number of first node added */
 90      return n_new - num + 1;
 91}
 92
 93/**********************************************************************/
 94
 95SCGRIB *scg_add_edge(SCG *g, int i, int j)
 96{     /* add new edge (i,j) to cliqued graph */
 97      SCGRIB *e;
 98      int t;
 99      xassert(1 <= i && i <= g->n);
100      xassert(1 <= j && j <= g->n);
101      if (i > j) t = i, i = j, j = t;
102      xassert(i < j);
103      e = dmp_get_atom(g->pool, sizeof(SCGRIB));
104      e->i = i;
105      e->j = j;
106      e->i_prev = NULL;
107      e->i_next = g->i_ptr[i];
108      e->j_prev = NULL;
109      e->j_next = g->j_ptr[j];
110      if (e->i_next != NULL) e->i_next->i_prev = e;
111      if (e->j_next != NULL) e->j_next->j_prev = e;
112      g->i_ptr[i] = g->j_ptr[j] = e;
113      return e;
114}
115
116/***********************************************************************
117*  NAME
118*
119*  scg_adj_list - get adjacency list for a node of cliqued graph
120*
121*  SYNOPSIS
122*
123*  #include "glpscg.h"
124*  int scg_adj_list(SCG *g, int i, int adj[]);
125*
126*  DESCRIPTION
127*
128*  For a given node i of the graph g the routine scg_adj_list stores
129*  numbers of adjacent nodes to locations adj[1], ..., adj[nadj], where
130*  0 <= nadj <= n-1 is the number of adjacent nodes, n is the number of
131*  nodes in the graph. Note that the list of adjacent nodes does not
132*  include node i.
133*
134*  RETURNS
135*
136*  The routine scg_adj_list returns nadj, which is the number of nodes
137*  adjacent to node i. */
138
139int scg_adj_list(SCG *g, int i, int adj[])
140{     int n = g->n;
141      char *flag = g->flag;
142      SCGRIB *e;
143      SCGCQE *p, *q;
144      int j, c, nadj = 0;
145      xassert(1 <= i && i <= n);
146#if _GLPSCG_DEBUG
147      for (j = 1; j <= n; j++) xassert(!flag[j]);
148#endif
149      /* go through the list of edges (i,j) and add node j to the
150         adjacency list */
151      for (e = g->i_ptr[i]; e != NULL; e = e->i_next)
152      {  j = e->j;
153#if _GLPSCG_DEBUG
154         xassert(i < j && j <= n);
155#endif
156         if (!flag[j]) adj[++nadj] = j, flag[j] = 1;
157      }
158      /* go through the list of edges (j,i) and add node j to the
159         adjacency list */
160      for (e = g->j_ptr[i]; e != NULL; e = e->j_next)
161      {  j = e->i;
162#if _GLPSCG_DEBUG
163         xassert(1 <= j && j < i);
164#endif
165         if (!flag[j]) adj[++nadj] = j, flag[j] = 1;
166      }
167#if 1
168      /* not tested yet */
169      xassert(g->v_ptr[i] == NULL);
170#endif
171      /* go through the list of cliques containing node i */
172      for (p = g->v_ptr[i]; p != NULL; p = p->v_next)
173      {  c = p->c;
174         /* clique c contains node i */
175#if _GLPSCG_DEBUG
176         xassert(1 <= c && c <= g->nc);
177#endif
178         /* go through the list of nodes of clique c and add them
179            (except node i) to the adjacency list */
180         for (q = g->c_ptr[c]; q != NULL; q = q->c_next)
181         {  j = q->i;
182            /* clique c contains node j */
183#if _GLPSCG_DEBUG
184            xassert(1 <= j && j <= n);
185#endif
186            if (j != i && !flag[j]) adj[++nadj] = j, flag[j] = 1;
187         }
188      }
189      /* reset the working array */
190      for (j = 1; j <= nadj; j++) flag[adj[j]] = 0;
191#if _GLPSCG_DEBUG
192      for (j = 1; j <= n; j++) xassert(!flag[j]);
193#endif
194      return nadj;
195}
196
197/***********************************************************************
198*  MAXIMUM WEIGHT CLIQUE
199*
200*  Two subroutines sub and wclique below are used to find a maximum
201*  weight clique in a given undirected graph. These subroutines are
202*  slightly modified version of the program WCLIQUE developed by Patric
203*  Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based on
204*  ideas from the article "P. R. J. Ostergard, A new algorithm for the
205*  maximum-weight clique problem, submitted for publication", which, in
206*  turn, is a generalization of the algorithm for unweighted graphs
207*  presented in "P. R. J. Ostergard, A fast algorithm for the maximum
208*  clique problem, submitted for publication".
209*
210*  USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. */
211
212struct dsa
213{     /* dynamic storage area */
214      SCG *g;
215      /* given (cliqued) graph */
216      int i;
217      /* node number, for which the adjacency list is built and stored
218         in the arrays adj and flag */
219      int nadj;
220      /* size of the adjacency list, 0 <= nadj <= n-1 */
221      int *adj; /* int list[1+n]; */
222      /* the adjacency list: adj[1], ..., adj[nadj] are nodes adjacent
223         to node i */
224      char *flag; /* int flag[1+n]; */
225      /* flag[j] means that there is edge (i,j) in the graph */
226      const int *wt; /* int wt[0:n-1]; */
227      /* weights */
228      int record;
229      /* weight of best clique */
230      int rec_level;
231      /* number of vertices in best clique */
232      int *rec; /* int rec[0:n-1]; */
233      /* best clique so far */
234      int *clique; /* int clique[0:n-1]; */
235      /* table for pruning */
236      int *set; /* int set[0:n-1]; */
237      /* current clique */
238};
239
240static int is_edge(struct dsa *dsa, int i, int j)
241{     /* returns non-zero if there is edge (i,j) in the graph */
242      SCG *g = dsa->g;
243      int n = g->n;
244      int *adj = dsa->adj;
245      char *flag = dsa->flag;
246      int k;
247      i++, j++;
248      xassert(1 <= i && i <= n);
249      xassert(1 <= j && j <= n);
250      /* build the adjacency list, if necessary */
251      if (dsa->i != i)
252      {  for (k = dsa->nadj; k >= 1; k--) flag[adj[k]] = 0;
253         dsa->i = i;
254         dsa->nadj = scg_adj_list(g, i, adj);
255         for (k = dsa->nadj; k >= 1; k--) flag[adj[k]] = 1;
256      }
257      return flag[j];
258}
259
260static void sub(struct dsa *dsa, int ct, int table[], int level,
261      int weight, int l_weight)
262{     int n = dsa->g->n;
263      const int *wt = dsa->wt;
264      int *rec = dsa->rec;
265      int *clique = dsa->clique;
266      int *set = dsa->set;
267      int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
268      newtable = xcalloc(n, sizeof(int));
269      if (ct <= 0)
270      {  /* 0 or 1 elements left; include these */
271         if (ct == 0)
272         {  set[level++] = table[0];
273            weight += l_weight;
274         }
275         if (weight > dsa->record)
276         {  dsa->record = weight;
277            dsa->rec_level = level;
278            for (i = 0; i < level; i++) rec[i] = set[i];
279         }
280         goto done;
281      }
282      for (i = ct; i >= 0; i--)
283      {  if ((level == 0) && (i < ct)) goto done;
284         k = table[i];
285         if ((level > 0) && (clique[k] <= (dsa->record - weight)))
286            goto done; /* prune */
287         set[level] = k;
288         curr_weight = weight + wt[k];
289         l_weight -= wt[k];
290         if (l_weight <= (dsa->record - curr_weight))
291            goto done; /* prune */
292         p1 = newtable;
293         p2 = table;
294         left_weight = 0;
295         while (p2 < table + i)
296         {  j = *p2++;
297#if 0
298            if (is_edge(dsa, j, k))
299#else
300            /* to minimize building adjacency lists (by mao) */
301            if (is_edge(dsa, k, j))
302#endif
303            {  *p1++ = j;
304               left_weight += wt[j];
305            }
306         }
307         if (left_weight <= (dsa->record - curr_weight)) continue;
308         sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight,
309            left_weight);
310      }
311done: xfree(newtable);
312      return;
313}
314
315static int wclique(SCG *g, const int w[], int sol[])
316{     int n = g->n;
317      const *wt = &w[1];
318      struct dsa _dsa, *dsa = &_dsa;
319      int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
320      xlong_t timer;
321      xassert(n > 0);
322      dsa->g = g;
323      dsa->i = 0;
324      dsa->nadj = 0;
325      dsa->adj = xcalloc(1+n, sizeof(int));
326      dsa->flag = xcalloc(1+n, sizeof(char));
327      memset(&dsa->flag[1], 0, n);
328      dsa->wt = wt;
329      dsa->record = 0;
330      dsa->rec_level = 0;
331      dsa->rec = &sol[1];
332      dsa->clique = xcalloc(n, sizeof(int));
333      dsa->set = xcalloc(n, sizeof(int));
334      used = xcalloc(n, sizeof(int));
335      nwt = xcalloc(n, sizeof(int));
336      pos = xcalloc(n, sizeof(int));
337      /* start timer */
338      timer = xtime();
339      /* order vertices */
340      for (i = 0; i < n; i++)
341      {  nwt[i] = 0;
342         for (j = 0; j < n; j++)
343            if (is_edge(dsa, i, j)) nwt[i] += wt[j];
344      }
345      for (i = 0; i < n; i++)
346         used[i] = 0;
347      for (i = n-1; i >= 0; i--)
348      {  max_wt = -1;
349         max_nwt = -1;
350         for (j = 0; j < n; j++)
351         {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
352               && nwt[j] > max_nwt)))
353            {  max_wt = wt[j];
354               max_nwt = nwt[j];
355               p = j;
356            }
357         }
358         pos[i] = p;
359         used[p] = 1;
360         for (j = 0; j < n; j++)
361            if ((!used[j]) && (j != p) && (is_edge(dsa, p, j)))
362               nwt[j] -= wt[p];
363      }
364      /* main routine */
365      wth = 0;
366      for (i = 0; i < n; i++)
367      {  wth += wt[pos[i]];
368         sub(dsa, i, pos, 0, 0, wth);
369         dsa->clique[pos[i]] = dsa->record;
370#if _GLPSCG_DEBUG
371         ;
372#else
373         if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
374#endif
375         {  /* print current record and reset timer */
376            xprintf("level = %d (%d); best = %d\n", i+1, n,
377               dsa->record);
378            timer = xtime();
379         }
380      }
381      xfree(dsa->adj);
382      xfree(dsa->flag);
383      xfree(dsa->clique);
384      xfree(dsa->set);
385      xfree(used);
386      xfree(nwt);
387      xfree(pos);
388      /* return the solution found */
389      for (i = 1; i <= dsa->rec_level; i++) sol[i]++;
390      return dsa->rec_level;
391}
392
393/***********************************************************************
394*  NAME
395*
396*  scg_max_clique - find maximum weight clique in given graph
397*
398*  SYNOPSIS
399*
400*  #include "glpscg.h"
401*  int scg_max_clique(SCG *g, const int w[], int list[]);
402*
403*  DESCRIPTION
404*
405*  The routine scg_max_clique finds an exact maximum weight clique in
406*  the given (cliqued) graph.
407*
408*  On entry the array w specifies node weights in locations w[1], ...
409*  w[n], where n is the number of nodes in the graph.
410*
411*  On exit the routine stores numbers of nodes included in the clique
412*  in locations list[1], ..., list[size], where 0 <= size <= n is the
413*  clique size.
414*
415*  RETURNS
416*
417*  The routine scg_max_clique returns the size of the clique found. */
418
419int scg_max_clique(SCG *g, const int w[], int list[])
420{     int size;
421      if (g->n == 0)
422         size = 0;
423      else
424         size = wclique(g, w, list);
425      return size;
426}
427
428/**********************************************************************/
429
430void scg_delete_graph(SCG *g)
431{     /* delete cliqued graph */
432      dmp_delete_pool(g->pool);
433      xfree(g->i_ptr);
434      xfree(g->j_ptr);
435      xfree(g->c_ptr);
436      xfree(g->v_ptr);
437      xfree(g->flag);
438      xfree(g);
439      return;
440}
441
442/* eof */