graph_gcg.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program */
4 /* GCG --- Generic Column Generation */
5 /* a Dantzig-Wolfe decomposition based extension */
6 /* of the branch-cut-and-price framework */
7 /* SCIP --- Solving Constraint Integer Programs */
8 /* */
9 /* Copyright (C) 2010-2018 Operations Research, RWTH Aachen University */
10 /* Zuse Institute Berlin (ZIB) */
11 /* */
12 /* This program is free software; you can redistribute it and/or */
13 /* modify it under the terms of the GNU Lesser General Public License */
14 /* as published by the Free Software Foundation; either version 3 */
15 /* of the License, or (at your option) any later version. */
16 /* */
17 /* This program is distributed in the hope that it will be useful, */
18 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
19 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
20 /* GNU Lesser General Public License for more details. */
21 /* */
22 /* You should have received a copy of the GNU Lesser General Public License */
23 /* along with this program; if not, write to the Free Software */
24 /* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
25 /* */
26 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
27 
33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34 #include <algorithm>
35 #include <iostream>
36 #include <cmath>
37 #include <set>
38 #include "graph_gcg.h"
39 
40 #ifdef WITH_GSL
41 #include <gsl/gsl_errno.h>
42 #endif
43 
44 using namespace std;
45 
46 namespace gcg {
47 
48 
49 GraphGCG::GraphGCG()
50 {
51  undirected = true;
52  locked = false;
53  initialized = false;
54  nodes = vector<int>(0);
55 #ifdef WITH_GSL
56  adj_matrix_sparse = gsl_spmatrix_alloc(1,1);
57  working_adj_matrix = NULL;
58 #endif
59  //adj_matrix = new vector<vector<double>>(0, vector<double>(0));
60 }
61 
62 GraphGCG::GraphGCG(int _n_nodes, bool _undirected)
63 {
64  undirected = _undirected;
65  locked = false;
66  initialized = true;
67  assert(_n_nodes >= 0);
68  nodes = vector<int>(_n_nodes, 0);
69 #ifdef WITH_GSL
70  adj_matrix_sparse = gsl_spmatrix_alloc(_n_nodes, _n_nodes);
71  working_adj_matrix = NULL;
72 #else
73  adj_matrix = vector<vector<double>>(_n_nodes, vector<double>(_n_nodes, 0.0));
74 #endif
75 }
76 
77 GraphGCG::~GraphGCG()
78 {
79  vector<int>().swap(this->nodes);
80 #ifndef WITH_GSL
81  vector<vector<double>>().swap(adj_matrix);
82  assert((int)adj_matrix.size() == 0);
83 #else
84  gsl_spmatrix_free(adj_matrix_sparse);
85  /*if(working_adj_matrix != NULL)
86  gsl_spmatrix_free(working_adj_matrix);*/
87 #endif
88  for(auto edge: edges)
89  {
90  if (edge != NULL)
91  delete edge;
92  }
93  initialized = false;
94 }
95 
96 SCIP_RETCODE GraphGCG::addNNodes(int _n_nodes)
97 {
98  if(locked || initialized)
99  {
100  return SCIP_ERROR;
101  }
102 #ifdef WITH_GSL
103  if(adj_matrix_sparse == NULL)
104  adj_matrix_sparse = gsl_spmatrix_alloc(_n_nodes, _n_nodes);
105  // add diagonal because of MCL algorithm
106  for(auto i = 0; i < _n_nodes; i++)
107  {
108  gsl_spmatrix_set(adj_matrix_sparse, i, i, 1.0);
109  }
110 #else
111  assert(adj_matrix.size() == 0);
112  adj_matrix = vector<vector<double>>(_n_nodes, vector<double>(_n_nodes, 0.0));
113 
114 #endif
115  nodes = vector<int>(_n_nodes, 0);
116  initialized = true;
117  return SCIP_OKAY;
118 }
119 
120 SCIP_RETCODE GraphGCG::addNNodes(int _n_nodes, vector<int> weights)
121 {
122  auto res = addNNodes(_n_nodes);
123  nodes = vector<int>(weights);
124  return res;
125 }
126 
127 
128 int GraphGCG::getNNodes()
129 {
130  if(!initialized) return 0;
131 #ifdef WITH_GSL
132  assert(nodes.size() == adj_matrix_sparse->size1);
133  assert(adj_matrix_sparse->size1 == adj_matrix_sparse->size2);
134  return adj_matrix_sparse->size1;
135 #else
136  assert(nodes.size() == adj_matrix.size());
137  return (int)adj_matrix.size();
138 #endif
139 }
140 
141 
142 #ifdef WITH_GSL
143 gsl_spmatrix* GraphGCG::getAdjMatrix()
144 {
145  cout << "Return adj_matrix_sparse..." << endl;
146  return adj_matrix_sparse;
147 }
148 
149 
150 void GraphGCG::expand(int factor)
151 {
152  assert(working_adj_matrix->sptype == (size_t)1);
153  const double alpha = 1.0;
154  gsl_spmatrix *tmp;
155  for(int i = 1; i < factor; i++)
156  {
157  tmp = gsl_spmatrix_alloc_nzmax(working_adj_matrix->size1, working_adj_matrix->size2, working_adj_matrix->nz, GSL_SPMATRIX_CCS);
158  gsl_spblas_dgemm(alpha, working_adj_matrix, working_adj_matrix, tmp);
159  gsl_spmatrix_free(working_adj_matrix);
160  working_adj_matrix = gsl_spmatrix_alloc_nzmax(tmp->size1, tmp->size2, tmp->nz, GSL_SPMATRIX_CCS);
161  gsl_spmatrix_memcpy(working_adj_matrix, tmp);
162  gsl_spmatrix_free(tmp);
163  }
164 }
165 
166 
167 // inflate columns and normalize them
168 void GraphGCG::inflate(double factor)
169 {
170  assert(working_adj_matrix->sptype == (size_t)1);
171  size_t i = 0;
172  while(i < working_adj_matrix->nz)
173  {
174  working_adj_matrix->data[i] = pow(working_adj_matrix->data[i], factor);
175  i++;
176  }
177  colL1Norm();
178 }
179 
180 
181 // normalize columns, but remove values below 1e-6
182 void GraphGCG::colL1Norm()
183 {
184  assert(working_adj_matrix->sptype == (size_t)1);
185  // normalize the columns to sum up to 1
186  for(int col = 0; col < (int)working_adj_matrix->size2; col++)
187  {
188  double col_sum = 0.0;
189  size_t begin_ = working_adj_matrix->p[col];
190  const size_t end_ = working_adj_matrix->p[col+1];
191  while(begin_ < end_)
192  {
193  col_sum += working_adj_matrix->data[begin_++];
194  }
195  begin_ = working_adj_matrix->p[col];
196  while(begin_ < end_)
197  {
198  working_adj_matrix->data[begin_] /= col_sum;
199  begin_++;
200  }
201  }
202 
203  // We need this part to remove very very small values, otherwise GSL behaves strangely
204  double threshold = 1e-8;
205  gsl_spmatrix *tmp = gsl_spmatrix_alloc(working_adj_matrix->size1, working_adj_matrix->size2);
206 
207  for (int col = 0; col < (int)working_adj_matrix->size1; col++)
208  {
209  size_t begin_ = working_adj_matrix->p[col];
210  const size_t end_ = working_adj_matrix->p[col+1];
211  while(begin_ < end_)
212  {
213  if(working_adj_matrix->data[begin_] > threshold)
214  {
215  gsl_spmatrix_set(tmp, working_adj_matrix->i[begin_], col, working_adj_matrix->data[begin_]);
216  }
217  begin_++;
218  }
219  }
220 
221  gsl_spmatrix *tmp_comp = gsl_spmatrix_compcol(tmp);
222  gsl_spmatrix_free(working_adj_matrix);
223  working_adj_matrix = gsl_spmatrix_alloc_nzmax(tmp_comp->size1, tmp_comp->size2, tmp_comp->nz, GSL_SPMATRIX_CCS);
224  gsl_spmatrix_memcpy(working_adj_matrix, tmp_comp);
225  gsl_spmatrix_free(tmp);
226  gsl_spmatrix_free(tmp_comp);
227 }
228 
229 
230 // remove values below 1e-3 or 5*1e-4 or 1e-4 and then normalize
231 void GraphGCG::prune()
232 {
233  double threshold_start = 1e-3;
234  gsl_spmatrix *tmp = gsl_spmatrix_alloc(working_adj_matrix->size1, working_adj_matrix->size2);
235 
236  double ave_col_sum = 0.0;
237  for (int col = 0; col < (int)working_adj_matrix->size1; col++)
238  {
239  size_t begin_ = working_adj_matrix->p[col];
240  const size_t end_ = working_adj_matrix->p[col+1];
241  while(begin_ < end_)
242  {
243  if(working_adj_matrix->data[begin_] > threshold_start)
244  {
245  gsl_spmatrix_set(tmp, working_adj_matrix->i[begin_], col, working_adj_matrix->data[begin_]);
246  ave_col_sum += working_adj_matrix->data[begin_];
247  }
248  begin_++;
249  }
250  }
251  ave_col_sum /= (double)(tmp->size1);
252 
253  double thresholds[2] = {5*1e-4, 1e-4};
254  for(int i = 0; i < 2; i++)
255  {
256  if(ave_col_sum > 0.85)
257  {
258  break;
259  }
260 
261  gsl_spmatrix_set_zero(tmp);
262  ave_col_sum = 0.0;
263  for (int col = 0; col < (int)working_adj_matrix->size1; col++)
264  {
265  size_t begin_ = working_adj_matrix->p[col];
266  const size_t end_ = working_adj_matrix->p[col+1];
267  while(begin_ < end_)
268  {
269  if(working_adj_matrix->data[begin_] > thresholds[i])
270  {
271  gsl_spmatrix_set(tmp, working_adj_matrix->i[begin_], col, working_adj_matrix->data[begin_]);
272  ave_col_sum += working_adj_matrix->data[begin_];
273  }
274  begin_++;
275  }
276  }
277  ave_col_sum /= (double)(tmp->size1);
278  }
279 
280  gsl_spmatrix *tmp_comp = gsl_spmatrix_compcol(tmp);
281  gsl_spmatrix_free(working_adj_matrix);
282  working_adj_matrix = gsl_spmatrix_alloc_nzmax(tmp_comp->size1, tmp_comp->size2, tmp_comp->nz, GSL_SPMATRIX_CCS);
283  gsl_spmatrix_memcpy(working_adj_matrix, tmp_comp);
284 
285  gsl_spmatrix_free(tmp);
286  gsl_spmatrix_free(tmp_comp);
287 
288  colL1Norm();
289 }
290 
291 
292 // checks if A*A - A stays unchanged
293 bool GraphGCG::stopMCL(int iter)
294 {
295  if(iter > 6)
296  {
297  gsl_spmatrix *tmpSquare = gsl_spmatrix_alloc_nzmax(working_adj_matrix->size1, working_adj_matrix->size2, working_adj_matrix->nz, GSL_SPMATRIX_CCS);
298  gsl_spmatrix *tmpMinus = gsl_spmatrix_alloc_nzmax(working_adj_matrix->size1, working_adj_matrix->size2, working_adj_matrix->nz, GSL_SPMATRIX_CCS);
299  gsl_spmatrix *res = gsl_spmatrix_alloc_nzmax(working_adj_matrix->size1, working_adj_matrix->size2, working_adj_matrix->nz, GSL_SPMATRIX_CCS);
300  gsl_spmatrix_memcpy(tmpMinus, working_adj_matrix);
301  gsl_spmatrix_scale(tmpMinus, -1.0);
302  double alpha = 1.0;
303  gsl_spblas_dgemm(alpha, working_adj_matrix, working_adj_matrix, tmpSquare);
304  gsl_spmatrix_add(res, tmpSquare, tmpMinus);
305  double min, max;
306  gsl_set_error_handler_off();
307  int status = gsl_spmatrix_minmax(res, &min, &max);
308  if(status == GSL_EINVAL)
309  {
310  min = 0.0;
311  max = 0.0;
312  }
313 
314  gsl_spmatrix_free(res);
315  gsl_spmatrix_free(tmpSquare);
316  gsl_spmatrix_free(tmpMinus);
317  if(abs(max - min) < 1e-8)
318  {
319  return true;
320  }
321  }
322  return false;
323 }
324 
325 vector<int> GraphGCG::getClustersMCL()
326 {
327  assert(working_adj_matrix->size1 == working_adj_matrix->size2);
328  vector<int> res(working_adj_matrix->size1, 0);
329  map<int, int> row_to_label;
330  set<int> nzrows;
331  vector<int> labels(working_adj_matrix->size1);
332 
333  for(int i = 0; i < (int)working_adj_matrix->nz; i++)
334  {
335  nzrows.insert(working_adj_matrix->i[i]);
336  }
337  int c = 0;
338  for(auto nzrow: nzrows)
339  {
340  row_to_label[nzrow] = c;
341  c++;
342  }
343 
344  map<int, vector<int> > clust_map;
345  // all data that belongs to that row belongs to its cluster
346  for (int col = 0; col < (int)working_adj_matrix->size2; col++)
347  {
348  size_t begin_ = working_adj_matrix->p[col];
349  const size_t end_ = working_adj_matrix->p[col+1];
350  while(begin_ < end_)
351  {
352  labels[col] = row_to_label[working_adj_matrix->i[begin_]];
353  clust_map[row_to_label[working_adj_matrix->i[begin_]]].push_back(col); // python version
354  begin_++;
355  }
356  }
357 
358  for(auto item : clust_map)
359  {
360  for(auto val : item.second)
361  {
362  //cout << val << ", ";
363  labels[val] = item.first;
364  }
365  }
366  // enumerations of labels may be broken (e.g. we can have labels 0,1,4,5), so we have to fix it to 0,1,2,3
367  set<int> existinglabels;
368  map<int, int> fixlabels;
369  for(int i = 0; i < (int)labels.size(); i++)
370  {
371  existinglabels.insert(labels[i]);
372  }
373 
374  c = 0;
375  for(auto existinglabel: existinglabels)
376  {
377  fixlabels[existinglabel] = c;
378  c++;
379  }
380  for(int i = 0; i < (int)labels.size(); i++)
381  {
382  labels[i] = fixlabels[labels[i]];
383  }
384  return labels;
385 }
386 
387 
388 void GraphGCG::initMCL()
389 {
390  working_adj_matrix = gsl_spmatrix_alloc_nzmax(adj_matrix_sparse->size1, adj_matrix_sparse->size2, adj_matrix_sparse->nz, GSL_SPMATRIX_CCS);
391  gsl_spmatrix_memcpy(working_adj_matrix, adj_matrix_sparse);
392 }
393 
394 
395 void GraphGCG::clearMCL()
396 {
397  gsl_spmatrix_free(working_adj_matrix);
398 }
399 
400 #else
401 vector<vector<double>> GraphGCG::getAdjMatrix()
402 {
403  return adj_matrix;
404 }
405 #endif
406 
407 
408 // TODO: we can use edges.size() here
409 int GraphGCG::getNEdges()
410 {
411  int n_edges = 0;
412 #ifdef WITH_GSL
413  if(initialized)
414  n_edges = adj_matrix_sparse->nz - adj_matrix_sparse->size1;
415  else
416  n_edges = 0;
417 #else
418  for(int i = 0; i < (int)adj_matrix.size(); i++)
419  {
420  n_edges += getNNeighbors(i);
421  }
422 #endif
423  if(undirected)
424  {
425  assert(n_edges % 2 == 0);
426  n_edges = (int) ( n_edges / 2.0);
427  }
428  assert(n_edges == (int)edges.size());
429  return n_edges;
430 }
431 
432 SCIP_Bool GraphGCG::isEdge(int node_i, int node_j)
433 {
434  assert(node_i >= 0);
435  assert(node_j >= 0);
436 #ifdef WITH_GSL
437  if(gsl_spmatrix_get(adj_matrix_sparse, node_i, node_j) != 0.0)
438  {
439  return 1;
440  }
441 #else
442  if(adj_matrix[node_i][node_j] != 0.0)
443  {
444  return 1;
445  }
446 #endif
447  return 0;
448 }
449 
450 int GraphGCG::getNNeighbors(int node)
451 {
452  assert(node >= 0);
453  int n_neighbors;
454 #ifdef WITH_GSL
455  if(!initialized) return 0;
456  assert(adj_matrix_sparse->sptype == (size_t)1);
457  assert(node < (int)adj_matrix_sparse->size2);
458  const size_t begin_ = adj_matrix_sparse->p[node];
459  const size_t end_ = adj_matrix_sparse->p[node+1];
460  n_neighbors = (int)(end_ - begin_);
461  if(gsl_spmatrix_get(adj_matrix_sparse, node, node) != 0.0)
462  {
463  n_neighbors -=1 ;
464  }
465 #else
466  assert(node < (int)adj_matrix.size());
467  n_neighbors = count_if(adj_matrix[node].begin(), adj_matrix[node].end(), [](double i) {return i != 0.0;});
468  // remove itself as a neighbor
469  if(adj_matrix[node][node] != 0.0)
470  {
471  n_neighbors -= 1;
472  }
473 #endif
474  return n_neighbors;
475 }
476 
477 vector<int> GraphGCG::getNeighbors(int node)
478 {
479  vector<int> res;
480 #ifdef WITH_GSL
481  if(!initialized || !locked) return res;
482  assert(adj_matrix_sparse->sptype == (size_t)1);
483  assert(node < (int)adj_matrix_sparse->size2);
484  size_t begin_ = adj_matrix_sparse->p[node];
485  size_t end_ = adj_matrix_sparse->p[node+1];
486  res.resize(end_ - begin_);
487  size_t curr = 0;
488  bool self_found = false;
489  while(begin_ < end_){
490  if((int)adj_matrix_sparse->i[begin_] == node) self_found = true;
491  res[curr++] = adj_matrix_sparse->i[begin_++];
492  }
493  assert(curr == res.size());
494  //res = row_inds;
495  if(self_found)
496  res.erase(remove(res.begin(), res.end(), node), res.end());
497 #else
498  for(int i = 0; i < (int)adj_matrix[node].size(); i++)
499  {
500  if(adj_matrix[node][i] != 0.0 && node != i)
501  {
502  res.push_back(i);
503  }
504  }
505 #endif
506  return res;
507 }
508 
509 vector<pair<int, double> > GraphGCG::getNeighborWeights(int node)
510 {
511  vector<pair<int, double> > res;
512 #ifdef WITH_GSL
513  if(!initialized || !locked) return res;
514  assert(adj_matrix_sparse->sptype == (size_t)1);
515  assert(node < (int)adj_matrix_sparse->size2);
516  size_t begin_ = adj_matrix_sparse->p[node];
517  size_t end_ = adj_matrix_sparse->p[node+1];
518  res.resize(end_ - begin_);
519  size_t curr = 0;
520  bool self_found = false;
521  int found_pos = 0;
522  while(begin_ < end_){
523  if((int)adj_matrix_sparse->i[begin_] == node)
524  {
525  self_found = true;
526  found_pos = curr;
527  }
528  double value = adj_matrix_sparse->data[begin_];
529  int row = (int)adj_matrix_sparse->i[begin_];
530  res[curr++] = pair<int, double>(row, value);
531  begin_++;
532  }
533  assert(curr == res.size());
534  if(self_found)
535  {
536  res.erase(res.begin() + found_pos);
537  }
538 #else
539  for(int i = 0; i < (int)adj_matrix[node].size(); i++)
540  {
541  if(adj_matrix[node][i] != 0.0 && node != i)
542  {
543  res.push_back(make_pair(i, adj_matrix[node][i]));
544  }
545  }
546 #endif
547 
548 
549  return res;
550 }
551 
553 SCIP_RETCODE GraphGCG::addNode(int node, int weight)
554 {
555 
556  if(locked)
557  {
558  return SCIP_ERROR;
559  }
560  int next_id = (int)nodes.size();
561  assert(node == next_id);
562 
563 #ifdef WITH_GSL
564  if(adj_matrix_sparse == NULL)
565  adj_matrix_sparse = gsl_spmatrix_alloc(1, 1);
566  // add diagonal because of MCL algorithm
567  gsl_spmatrix_set(adj_matrix_sparse, next_id, next_id, 1.0);
568 
569  //return SCIP_NOTIMPL;
570 #else
571 
572  // add new column
573  for(size_t i = 0; i < adj_matrix.size(); i++)
574  {
575  adj_matrix[i].push_back(0.0);
576  }
577  // add new row
578  adj_matrix.push_back(vector<double>(next_id+1));
579 
580  // add self loop
581  adj_matrix[next_id][next_id] = 1.0;
582  assert(adj_matrix.size() == adj_matrix[0].size());
583 #endif
584  initialized = true;
585  nodes.push_back(weight);
586  return SCIP_OKAY;
587 }
588 
589 SCIP_RETCODE GraphGCG::addNode()
590 {
591 #ifdef WITH_GSL
592  return addNode((int)adj_matrix_sparse->size2, 0);
593 #else
594  return addNode((int)adj_matrix.size(), 0);
595 #endif
596 }
597 
598 SCIP_RETCODE GraphGCG::deleteNode(int node)
599 {
600  if(locked)
601  {
602  return SCIP_ERROR;
603  }
604  // Would be very inefficient and is not necessary
605  return SCIP_INVALIDCALL;
606 }
607 
608 SCIP_RETCODE GraphGCG::addEdge(int node_i, int node_j)
609 {
610  return addEdge(node_i, node_j, 1);
611 }
612 
613 SCIP_RETCODE GraphGCG::addEdge(int node_i, int node_j, double weight)
614 {
615  if(locked || !initialized || node_i == node_j)
616  {
617  return SCIP_ERROR;
618  }
619  bool new_edge = true;
620  assert(weight >= 0.0);
621  assert(node_i >= 0);
622  assert(node_j >= 0);
623 #ifdef WITH_GSL
624  assert(node_i < (int)adj_matrix_sparse->size2);
625  assert(node_j < (int)adj_matrix_sparse->size2);
626  if (gsl_spmatrix_get(adj_matrix_sparse, node_i, node_j) != 0.0)
627  {
628  new_edge = false;
629  }
630  if(new_edge)
631  {
632  gsl_spmatrix_set(adj_matrix_sparse, node_i, node_j, weight);
633  if(undirected)
634  {
635  gsl_spmatrix_set(adj_matrix_sparse, node_j, node_i, weight);
636  }
637  }
638 #else
639  assert(node_i < (int)adj_matrix.size());
640  assert(node_j < (int)adj_matrix.size());
641  (adj_matrix[node_i][node_j] == 0.0) ? new_edge = true : new_edge = false;
642  adj_matrix[node_i][node_j] = weight;
643  if(undirected)
644  {
645  adj_matrix[node_j][node_i] = weight;
646  }
647 #endif
648  if(new_edge && node_i != node_j)
649  {
650  int first, second;
651  first = node_i < node_j ? node_i : node_j;
652  second = node_i == first ? node_j : node_i;
653  EdgeGCG* edge1 = new EdgeGCG(first, second, weight);
654  edges.push_back(edge1);
655  }
656 
657  return SCIP_OKAY;
658 }
659 
660 SCIP_RETCODE GraphGCG::setEdge(int node_i, int node_j, double weight)
661 {
662  if(locked || !initialized)
663  {
664  return SCIP_ERROR;
665  }
666  assert(weight >= 0.0);
667  assert(node_i >= 0);
668  assert(node_j >= 0);
669 #ifdef WITH_GSL
670  assert(node_i < (int)adj_matrix_sparse->size2);
671  assert(node_j < (int)adj_matrix_sparse->size2);
672  gsl_spmatrix_set(adj_matrix_sparse, node_i, node_j, weight);
673  if(undirected)
674  {
675  gsl_spmatrix_set(adj_matrix_sparse, node_j, node_i, weight);
676  }
677 #else
678  assert(node_i < (int)adj_matrix.size());
679  assert(node_j < (int)adj_matrix.size());
680  adj_matrix[node_i][node_j] = weight;
681  if(undirected)
682  {
683  adj_matrix[node_j][node_i] = weight;
684  }
685 #endif
686  if(node_i != node_j)
687  for(auto edge: edges)
688  {
689  if((edge->src == node_i && edge->dest == node_j) || (edge->src == node_j && edge->dest == node_i))
690  {
691  edge->weight = weight;
692  }
693  }
694  return SCIP_OKAY;
695 }
696 
697 SCIP_RETCODE GraphGCG::deleteEdge(int node_i, int node_j)
698 {
699  return setEdge(node_i, node_j, 0);
700 }
701 
702 int GraphGCG::graphGetWeights(int node)
703 {
704  return nodes[node];
705 }
706 
707 double GraphGCG::getEdgeWeight(int node_i, int node_j)
708 {
709  assert(node_i >= 0);
710  assert(node_j >= 0);
711  double weight;
712 #ifdef WITH_GSL
713  assert(node_i < (int)adj_matrix_sparse->size1);
714  assert(node_j < (int)adj_matrix_sparse->size1);
715  weight = gsl_spmatrix_get(adj_matrix_sparse, node_i, node_j);
716 #else
717  assert(node_i < (int)adj_matrix.size());
718  assert(node_j < (int)adj_matrix.size());
719  //assert(adj_matrix[node_i][node_j] == 0);
720  weight = adj_matrix[node_i][node_j];
721 #endif
722  return weight;
723 }
724 
725 SCIP_RETCODE GraphGCG::flush()
726 {
727  locked = true;
728 #ifdef WITH_GSL
729  gsl_spmatrix *adj_matrix_sparse_tmp = gsl_spmatrix_compcol(adj_matrix_sparse);
730  gsl_spmatrix_free(adj_matrix_sparse);
731  adj_matrix_sparse = adj_matrix_sparse_tmp;
732 #endif
733  return SCIP_OKAY;
734 }
735 
736 SCIP_RETCODE GraphGCG::normalize()
737 {
738  double scaler = 0.0;
739 #ifdef WITH_GSL
740  double min, max;
741  gsl_spmatrix_minmax(adj_matrix_sparse, &min, &max);
742  scaler = max;
743  gsl_spmatrix_scale(adj_matrix_sparse, double(1.0/scaler));
744 
745 #else
746  for(int i = 0; i < (int)adj_matrix.size(); i++)
747  {
748  vector<double>::iterator tmp = max_element(adj_matrix[i].begin(), adj_matrix[i].end());
749  double curr_max = *tmp;
750  if( curr_max > scaler )
751  {
752  scaler = curr_max;
753  }
754  }
755 
756  for(int i = 0; i < (int)adj_matrix.size(); i++)
757  {
758  for(int j = 0; j < (int)adj_matrix[i].size(); j++)
759  {
760  adj_matrix[i][j] /= scaler;
761  assert(adj_matrix[i][j] <= 1);
762  assert(adj_matrix[i][j] >= 0);
763  }
764  }
765 #endif
766  return SCIP_OKAY;
767 }
768 
769 double GraphGCG::getEdgeWeightPercentile(double q)
770 {
771  double res = -1;
772  vector<double> all_weights;
773 #ifdef WITH_GSL
774  all_weights = vector<double>(adj_matrix_sparse->data, adj_matrix_sparse->data + adj_matrix_sparse->nz);
775 #else
776  for(int i = 0; i < (int)adj_matrix.size(); i++)
777  {
778  for(int j = 0; j < (int)adj_matrix[i].size(); j++)
779  {
780  if(adj_matrix[i][j] != 0.0)
781  {
782  all_weights.push_back(adj_matrix[i][j]);
783  }
784  }
785  }
786 #endif
787  sort(all_weights.begin(), all_weights.end());
788 
789  int upper_pos = (int) floor((q/100.0)*all_weights.size());
790  int lower_pos = (int) ceil((q/100.0)*all_weights.size());
791  if(upper_pos != lower_pos)
792  {
793  res = (all_weights[upper_pos] + all_weights[lower_pos]) / 2.0;
794  }
795  else
796  {
797  res = all_weights[lower_pos];
798  }
799  return res;
800 }
801 
802 
803 SCIP_RETCODE GraphGCG::getEdges(vector<void *>& _edges)
804 {
805  _edges.resize(edges.size());
806  //_edges = static_cast< vector<void *> >(edges);
807  for(int i = 0; i < (int)edges.size(); i++)
808  {
809  assert(edges[i]->src < (int)nodes.size());
810  assert(edges[i]->dest < (int)nodes.size());
811  assert(edges[i]->src != edges[i]->dest);
812 
813  (_edges[i]) = (void *)(edges[i]);
814  }
815  return SCIP_OKAY;
816 }
817 
818 
819 // Compare two edges according to their weights.
820 // Used in sort() for sorting an array of edges
821 int GraphGCG::edgeComp(const EdgeGCG* a, const EdgeGCG* b)
822 {
823  return (a->src < b->src) || (a->src == b->src && a->weight < b->weight);
824 }
825 
826 
827 
828 
829 }
STL namespace.
Implementation of the graph which supports both node and edge weights.
double weight
Definition: graph_gcg.h:52