Scippy

SCIP

Solving Constraint Integer Programs

HeurFarthestInsert.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2017 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License. */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file HeurFarthestInsert.cpp
17  * @brief farthest insert - combinatorial heuristic for TSP
18  * @author Timo Berthold
19  */
20 
21 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
22 
23 
24 #include <iostream>
25 #include <cassert>
26 
27 #include "objscip/objscip.h"
28 #include "GomoryHuTree.h"
29 #include "HeurFarthestInsert.h"
30 #include "ProbDataTSP.h"
31 
32 using namespace tsp;
33 using namespace std;
34 
35 /** method finding the edge going from the node with id index1 to the node with id index2 */
37  GRAPHNODE* nodes, /**< all nodes of the graph */
38  int index1, /**< id of the node where the searched edge starts */
39  int index2 /**< id of the node where the searched edge ends */
40  )
41 {
42  GRAPHEDGE* edge = nodes[index1].first_edge;
43 
44  // regard every outgoing edge of node index1 and stop if adjacent to node index2
45  while( edge != NULL )
46  {
47  if( edge->adjac->id == index2 )
48  break;
49  edge = edge->next;
50  }
51 
52  return edge;
53 }
54 
55 /** method updating the distances of the nodes to the tour after having inserted one node with id index */
57  GRAPHNODE* nodes, /**< all nodes of the graph */
58  double* dist, /**< array with current distances of all nodes to the subtour */
59  int idx /**< id of the inserted node */
60  )
61 {
62  GRAPHEDGE* edge = nodes[idx].first_edge;
63 
64  // regard all outgoing edges of the node and update,
65  // if the length and therefore the distance of the adjacent is smaller
66  // and the edge is not fixed to 0.
67  while( edge != NULL )
68  {
69  if( dist[edge->adjac->id] > edge->length && SCIPvarGetUbGlobal(edge->var) != 0.0 )
70  dist[edge->adjac->id] = edge->length;
71  edge = edge->next;
72  }
73 }
74 
75 
76 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
77 SCIP_DECL_HEURFREE(HeurFarthestInsert::scip_free)
78 {
79  return SCIP_OKAY;
80 } /*lint !e715*/
81 
82 /** initialization method of primal heuristic (called after problem was transformed) */
83 SCIP_DECL_HEURINIT(HeurFarthestInsert::scip_init)
84 {
85  ProbDataTSP* probdata = dynamic_cast<ProbDataTSP*>(SCIPgetObjProbData(scip));
86  graph_ = probdata->getGraph(); /*lint !e613*/
87  capture_graph(graph_);
88  return SCIP_OKAY;
89 } /*lint !e715*/
90 
91 /** deinitialization method of primal heuristic (called before transformed problem is freed) */
92 SCIP_DECL_HEUREXIT(HeurFarthestInsert::scip_exit)
93 {
94  release_graph(&graph_);
95 
96  return SCIP_OKAY;
97 } /*lint !e715*/
98 
99 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin)
100  *
101  * This method is called when the presolving was finished and the branch and bound process is about to begin.
102  * The primal heuristic may use this call to initialize its branch and bound specific data.
103  *
104  */
105 SCIP_DECL_HEURINITSOL(HeurFarthestInsert::scip_initsol)
106 {
107  return SCIP_OKAY;
108 } /*lint !e715*/
109 
110 /** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed)
111  *
112  * This method is called before the branch and bound process is freed.
113  * The primal heuristic should use this call to clean up its branch and bound data.
114  */
115 SCIP_DECL_HEUREXITSOL(HeurFarthestInsert::scip_exitsol)
116 {
117  return SCIP_OKAY;
118 } /*lint !e715*/
119 
120 /** execution method of primal heuristic
121  *
122  * Searches for feasible primal solutions. The method is called in the node processing loop.
123  *
124  * possible return values for *result:
125  * - SCIP_FOUNDSOL : at least one feasible primal solution was found
126  * - SCIP_DIDNOTFIND : the heuristic searched, but did not find a feasible solution
127  * - SCIP_DIDNOTRUN : the heuristic was skipped
128  * - SCIP_DELAYED : the heuristic was skipped, but should be called again as soon as possible, disregarding
129  * its frequency
130  */
131 SCIP_DECL_HEUREXEC(HeurFarthestInsert::scip_exec)
132 {
133  int nnodes = graph_->nnodes; /*lint !e613*/
134  int nedges = graph_->nedges; /*lint !e613*/
135 
136  SCIP_Bool hasFixedEdges = FALSE;
137  for(int e = 0; e < nedges; ++e)
138  {
139  GRAPHEDGE* edge = &(graph_->edges[e]); /*lint !e613*/
140  if( SCIPvarGetLbGlobal(edge->var) == 1.0 )
141  {
142  hasFixedEdges = true;
143  break;
144  }
145  }
146 
147  // no longer need "SCIPgetNRuns(scip) > 1" since we now respect fixed variables after restart
148  if( nnodes < 3 || hasFixedEdges )
149  *result = SCIP_DIDNOTRUN;
150  else
151  {
152  bool* subtour;
153  int i;
154  double* dist;
155 
156  GRAPHNODE* startnode;
157  GRAPHNODE* node;
158  GRAPHEDGE* edge;
159 
160  GRAPHEDGE** bestedges; // will contain the best insertion of a given node into a subtour
161  GRAPHEDGE** edges; // will contain some insertion of a given node into a subtour
162  GRAPHEDGE** successor; // stores the successor of a node in the current subtour
163  GRAPHNODE* nodes = graph_->nodes; /*lint !e613*/
164 
165  for( i = 0; i < nnodes; i++ )
166  assert( i == nodes[i].id );
167 
168  // memory allocation
169  SCIP_CALL( SCIPallocBufferArray(scip, &subtour, nnodes) ); /*lint !e530*/
170  SCIP_CALL( SCIPallocBufferArray(scip, &dist, nnodes) ); /*lint !e530*/
171  SCIP_CALL( SCIPallocBufferArray(scip, &successor, nnodes) ); /*lint !e530*/
172  SCIP_CALL( SCIPallocBufferArray(scip, &edges, 3) ); /*lint !e530*/
173  SCIP_CALL( SCIPallocBufferArray(scip, &bestedges, 3) ); /*lint !e530*/
174 
175  BMSclearMemoryArray(subtour, nnodes);
176  for( i = 0; i < nnodes; i++ )
177  dist[i] = DBL_MAX;
178 
179  // building up a 3-circle, only using edges not fixed to 0
180  SCIP_Bool foundThreeCircle = FALSE;
181  for(int u = 0; u < nnodes - 2 && !foundThreeCircle; ++u)
182  {
183  for(int v = u + 1; v < nnodes - 1 && !foundThreeCircle; ++v)
184  {
185  GRAPHEDGE * uv = findEdge(nodes, u, v);
186  assert(uv != NULL);
187  if( SCIPvarGetUbGlobal(uv->var) == 0.0 )
188  continue;
189  for(int w = v + 1; (w < nnodes) && !foundThreeCircle; ++w) /*lint !e845*/
190  {
191  GRAPHEDGE * vw = findEdge(nodes, v, w);
192  assert(vw != NULL);
193  GRAPHEDGE * wu = findEdge(nodes, w, u);
194  assert(wu != NULL);
195  if( SCIPvarGetUbGlobal(vw->var) == 0.0 || SCIPvarGetUbGlobal(wu->var) == 0.0 )
196  continue;
197  else {
198  foundThreeCircle = true;
199 
200  subtour[u] = true;
201  dist[u] = 0.0;
202  updateDistances(nodes, dist, u);
203  subtour[v] = true;
204  dist[v] = 0.0;
205  updateDistances(nodes, dist, v);
206  subtour[w] = true;
207  dist[w] = 0.0;
208  updateDistances(nodes, dist, w);
209  successor[u] = uv;
210  successor[v] = vw;
211  successor[w] = wu;
212  } // foundThreeCircle with no fixed variables
213  } // for w
214  } // for v
215  } // for u
216 
217  if( !foundThreeCircle )
218  {
219  *result = SCIP_DIDNOTFIND;
220  }
221  else
222  {
223  double maxmin;
224  double minval;
225  int newnodeindex;
226 
227  SCIP_Bool couldNotInsert = FALSE;
228 
229  // widen the subtour by one node each step until you have a complete tour, actually the farthest insert heuritic
230  int subtourlength = 3;
231  for(; subtourlength < nnodes; subtourlength++ )
232  {
233  // find the node with the maximal distance to the tour
234  maxmin = 0.0;
235  newnodeindex = -1;
236  for( i = 0; i < nnodes; i++)
237  {
238  if( (maxmin < dist[i] && dist[i] != DBL_MAX) || (maxmin == dist[i] && !subtour[i]) ) /*lint !e777*/
239  {
240  maxmin = dist[i];
241  newnodeindex = i;
242  }
243  }
244  if(newnodeindex == -1)
245  {
246  couldNotInsert = TRUE;
247  break;
248  }
249 
250  // find connection to one node in the tour
251  BMSclearMemoryArray(bestedges, 3);
252  edge = nodes[newnodeindex].first_edge;
253  startnode = NULL;
254 
255  while( edge != NULL )
256  {
257  if( subtour[edge->adjac->id] && SCIPvarGetUbGlobal(edge->var) != 0.0 )
258  break;
259  edge = edge->next;
260  }
261 
262  assert(edge != NULL);
263  assert(subtour[edge->adjac->id]);
264 
265  // find best insertion of the new node by trying to replace any edge connecting by the two edges connecting
266  // its end node with the new node
267  minval = DBL_MAX;
268  edges[0] = edge;
269  startnode = edge->adjac;
270  node = startnode;
271 
272  // succeed to the next edge in the subtour
273  do
274  {
275  edges[1] = successor[node->id];
276  edges[2] = findEdge(nodes, edges[1]->adjac->id, newnodeindex);
277  assert( edges[2] != NULL );
278 
279  // check, whether you have found a better (feasible) insertion
280  if( edges[0]->back->length - edges[1]->length + edges[2]->back->length < minval
281  && SCIPvarGetUbGlobal(edges[0]->var) != 0.0
282  && SCIPvarGetUbGlobal(edges[2]->var) != 0.0 )
283  {
284  minval = edges[0]->back->length - edges[1]->length + edges[2]->back->length;
285  for( i = 0; i < 3; i++ )
286  bestedges[i] = edges[i];
287  }
288  node = edges[1]->adjac;
289  edges[0] = edges[2]->back;
290  }
291  while( node != startnode);
292 
293  if( minval == DBL_MAX ) /*lint !e777*/
294  {
295  couldNotInsert = TRUE;
296  break;
297  }
298  else
299  {
300  // bestedges should contain a 3-cycle (modulo orientation) connecting new node with two incident ones of the tour
301  assert(bestedges[0]->adjac->id == bestedges[1]->back->adjac->id);
302  assert(bestedges[1]->adjac->id == bestedges[2]->back->adjac->id);
303  assert(bestedges[2]->adjac->id == bestedges[0]->back->adjac->id);
304  assert(subtour[bestedges[0]->adjac->id]);
305  assert(subtour[bestedges[1]->adjac->id]);
306  assert(bestedges[2]->adjac->id == newnodeindex);
307  assert(!subtour[newnodeindex]);
308 
309  // now officially insert the new node into the tour
310  successor[bestedges[0]->adjac->id] = bestedges[0]->back;
311  successor[bestedges[2]->adjac->id] = bestedges[2]->back;
312  dist[newnodeindex] = 0.0;
313  subtour[newnodeindex] = true;
314  updateDistances(nodes, dist, newnodeindex);
315  } // minval < DBL_MAX
316  } // for subtourlength
317 
318  if(couldNotInsert)
319  {
320  *result = SCIP_DIDNOTFIND;
321  }
322  else
323  {
324  SCIP_SOL* sol;
325  SCIP_Bool success;
326 
327  // now create a solution out of the edges stored in successor and try to add it to SCIP
328  SCIP_CALL( SCIPcreateSol (scip, &sol, heur) );
329  for( i = 0; i < nnodes; i++ )
330  {
331  SCIP_CALL( SCIPsetSolVal(scip, sol, successor[i]->var, 1.0) );
332  }
333  SCIP_CALL( SCIPtrySol(scip, sol, FALSE, FALSE, FALSE, FALSE, FALSE, &success) );
334  if( success )
335  *result = SCIP_FOUNDSOL;
336  else
337  *result = SCIP_DIDNOTFIND;
338  SCIP_CALL( SCIPfreeSol(scip, &sol) );
339  } // couldNotInsert == FALSE
340  } // foundThreeCircle == TRUE
341 
342  // free all local memory
343  SCIPfreeBufferArray(scip, &bestedges);
344  SCIPfreeBufferArray(scip, &edges);
345  SCIPfreeBufferArray(scip, &successor);
346  SCIPfreeBufferArray(scip, &dist);
347  SCIPfreeBufferArray(scip, &subtour);
348  }
349 
350  return SCIP_OKAY;
351 } /*lint !e715*/
352 
353 /** clone method which will be used to copy a objective plugin */
354 SCIP_DECL_HEURCLONE(scip::ObjCloneable* HeurFarthestInsert::clone) /*lint !e665*/
355 {
356  return new HeurFarthestInsert(scip);
357 }
struct GraphEdge * next
Definition: GomoryHuTree.h:59
double length
Definition: GomoryHuTree.h:57
SCIP_DECL_HEURFREE(HeurFarthestInsert::scip_free)
SCIP_Real SCIPvarGetLbGlobal(SCIP_VAR *var)
Definition: var.c:17166
void updateDistances(GRAPHNODE *nodes, double *dist, int idx)
#define FALSE
Definition: def.h:64
#define TRUE
Definition: def.h:63
SCIP_DECL_HEUREXITSOL(HeurFarthestInsert::scip_exitsol)
generator for global cuts in undirected graphs
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip.h:21937
farthest insert - combinatorial heuristic for TSP
SCIP_VAR * var
Definition: GomoryHuTree.h:64
SCIP_Real SCIPvarGetUbGlobal(SCIP_VAR *var)
Definition: var.c:17176
GRAPHNODE * adjac
Definition: GomoryHuTree.h:62
struct GraphEdge * back
Definition: GomoryHuTree.h:60
struct GraphEdge * first_edge
Definition: GomoryHuTree.h:43
#define NULL
Definition: lpi_spx1.cpp:137
C++ wrapper classes for SCIP.
#define SCIP_CALL(x)
Definition: def.h:306
GRAPH * getGraph()
Definition: ProbDataTSP.h:111
C++ problem data for TSP.
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip.h:21925
SCIP_RETCODE SCIPsetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var, SCIP_Real val)
Definition: scip.c:37867
#define SCIP_Bool
Definition: def.h:61
SCIP_DECL_HEURCLONE(scip::ObjCloneable *HeurFarthestInsert::clone)
SCIP_RETCODE SCIPfreeSol(SCIP *scip, SCIP_SOL **sol)
Definition: scip.c:37631
SCIP_RETCODE SCIPtrySol(SCIP *scip, SCIP_SOL *sol, SCIP_Bool printreason, SCIP_Bool completely, SCIP_Bool checkbounds, SCIP_Bool checkintegrality, SCIP_Bool checklprows, SCIP_Bool *stored)
Definition: scip.c:39749
scip::ObjProbData * SCIPgetObjProbData(SCIP *scip)
SCIP_DECL_HEURINIT(HeurFarthestInsert::scip_init)
SCIP_DECL_HEUREXEC(HeurFarthestInsert::scip_exec)
SCIP_DECL_HEURINITSOL(HeurFarthestInsert::scip_initsol)
Definition of base class for all clonable classes.
Definition: objcloneable.h:38
SCIP_DECL_HEUREXIT(HeurFarthestInsert::scip_exit)
GRAPHEDGE * findEdge(GRAPHNODE *nodes, int index1, int index2)
void capture_graph(GRAPH *gr)
void release_graph(GRAPH **gr)
#define nnodes
Definition: gastrans.c:65
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:85
SCIP_RETCODE SCIPcreateSol(SCIP *scip, SCIP_SOL **sol, SCIP_HEUR *heur)
Definition: scip.c:37005