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