Scippy

SCIP

Solving Constraint Integer Programs

ReaderTSP.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-2019 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 scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file ReaderTSP.cpp
17  * @brief C++ file reader for TSP data files
18  * @author Timo Berthold
19  */
20 
21 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
22 
23 #include <iostream>
24 #include <fstream>
25 #include <string>
26 #include <sstream>
27 
28 #include "objscip/objscip.h"
29 
30 #include "scip/cons_linear.h"
31 #include <math.h>
32 
33 #include "ReaderTSP.h"
34 #include "ProbDataTSP.h"
35 #include "ConshdlrSubtour.h"
36 #include "GomoryHuTree.h"
37 
38 using namespace tsp;
39 using namespace scip;
40 using namespace std;
41 
42 #define NINT(x) (floor(x+0.5))
43 
44 /** parses the node list */
45 void ReaderTSP::getNodesFromFile(
46  std::ifstream& filedata, /**< filestream containing the data to extract */
47  double* x_coords, /**< double array to be filled with the x-coordinates of the nodes */
48  double* y_coords, /**< same for y-coordinates */
49  GRAPH* graph /**< the graph which is to be generated by the nodes */
50  )
51 {
52  int i = 0;
53  int nodenumber;
54  GRAPHNODE* node = &(graph->nodes[0]);
55 
56  // extract every node out of the filestream
57  while ( i < graph->nnodes && !filedata.eof() )
58  {
59  filedata >> nodenumber >> x_coords[i] >> y_coords[i];
60 
61  // assign everything
62  node->id = i;
63  if( nodenumber-1 != i)
64  cout<<"warning: nodenumber <" <<nodenumber<<"> does not match its index in node list <"<<i+1
65  <<">. Node will get number "<<i+1<<" when naming variables and constraints!"<<endl;
66  node->x = x_coords[i];
67  node->y = y_coords[i];
68  node->first_edge = NULL;
69  node++;
70  i++;
71  }
72  assert( i == graph->nnodes );
73 } /*lint !e1762*/
74 
75 /** adds a variable to both halfedges and captures it for usage in the graph */
76 SCIP_RETCODE ReaderTSP::addVarToEdges(
77  SCIP* scip, /**< SCIP data structure */
78  GRAPHEDGE* edge, /**< an edge of the graph */
79  SCIP_VAR* var /**< variable corresponding to that edge */
80  )
81 {
82  assert(scip != NULL);
83  assert(edge != NULL);
84  assert(var != NULL);
85 
86  /* add variable to forward edge and capture it for usage in graph */
87  edge->var = var;
88  SCIP_CALL( SCIPcaptureVar(scip, edge->var) );
89 
90  /* two parallel halfedges have the same variable,
91  * add variable to backward edge and capture it for usage in graph */
92  edge->back->var = edge->var;
93  SCIP_CALL( SCIPcaptureVar(scip, edge->back->var) );
94 
95  return SCIP_OKAY;
96 }
97 
98 /** method asserting, that the file has had the correct format and everything was set correctly */
99 bool ReaderTSP::checkValid(
100  GRAPH* graph, /**< the constructed graph, schould not be NULL */
101  std::string name, /**< the name of the file */
102  std::string type, /**< the type of the problem, should be "TSP" */
103  std::string edgeweighttype, /**< type of the edgeweights, should be "EUC_2D", "MAX_2D", "MAN_2D",
104  * "ATT", or "GEO" */
105  int nnodes /**< dimension of the problem, should at least be one */
106  )
107 {
108  // if something seems to be strange, it will be reported, that file was not valid
109  if( nnodes < 1 )
110  {
111  cout << "parse error in file <" << name << "> dimension should be greater than 0"<< endl ;
112  return false;
113  }
114  if (type != "TSP" )
115  {
116  cout << "parse error in file <" << name << "> type should be TSP" << endl;
117  return false;
118  }
119  if ( !( edgeweighttype == "EUC_2D" || edgeweighttype == "MAX_2D" || edgeweighttype == "MAN_2D"
120  || edgeweighttype == "GEO" || edgeweighttype == "ATT") )
121  {
122  cout << "parse error in file <" << name
123  << "> unknown weight type, should be EUC_2D, MAX_2D, MAN_2D, ATT, or GEO" << endl;
124  return false;
125  }
126  if( graph == NULL)
127  {
128  cout << "error while reading file <" << name << ">, graph is uninitialized. ";
129  cout << "Probably NODE_COORD_SECTION is missing" << endl;
130  return false;
131  }
132  return true;
133 } /*lint !e1762*/
134 
135 
136 /** destructor of file reader to free user data (called when SCIP is exiting) */
137 SCIP_DECL_READERFREE(ReaderTSP::scip_free)
138 {
139  return SCIP_OKAY;
140 } /*lint !e715*/
141 
142 /** problem reading method of reader
143  *
144  * possible return values for *result:
145  * - SCIP_SUCCESS : the reader read the file correctly and created an appropritate problem
146  * - SCIP_DIDNOTRUN : the reader is not responsible for given input file
147  *
148  * If the reader detected an error in the input file, it should return with RETCODE SCIP_READERR or SCIP_NOFILE.
149  */
150 SCIP_DECL_READERREAD(ReaderTSP::scip_read)
151 {
152  SCIP_RETCODE retcode;
153 
154  GRAPH* graph = NULL;
155  GRAPHNODE* node;
156  GRAPHNODE* nodestart; // the two incident nodes of an edge
157  GRAPHNODE* nodeend;
158  GRAPHEDGE* edgeforw; // two converse halfedges
159  GRAPHEDGE* edgebackw;
160  GRAPHEDGE* edge;
161 
162  double* x_coords = NULL; // arrays of coordinates of the nodes
163  double* y_coords = NULL;
164 
165 #ifdef SCIP_DEBUG
166  double** weights = NULL;
167 #endif
168 
169  double x; // concrete coordinates
170  double y;
171 
172  int nnodes = 0;
173  int nedges = 0;
174  int i;
175  int j;
176 
177  string name = "MY_OWN_LITTLE_TSP";
178  string token;
179  string type = "TSP";
180  string edgeweighttype = "EUC_2D";
181 
182  retcode = SCIP_OKAY;
183  *result = SCIP_DIDNOTRUN;
184 
185  // open the file
186  ifstream filedata(filename);
187  if( !filedata )
188  return SCIP_READERROR;
189  filedata.clear();
190 
191  // read the first lines of information
192  filedata >> token;
193  while( !filedata.eof() )
194  {
195  if( token == "NAME:" )
196  filedata >> name;
197  else if( token == "NAME" )
198  filedata >> token >> name;
199  else if( token == "TYPE:" )
200  filedata >> type;
201  else if( token == "TYPE" )
202  filedata >> token >> type;
203  else if( token == "DIMENSION:" )
204  {
205  filedata >> nnodes;
206  nedges = nnodes * (nnodes-1);
207  }
208  else if( token == "DIMENSION" )
209  {
210  filedata >> token >> nnodes;
211  nedges = nnodes * (nnodes-1);
212  }
213  else if( token == "EDGE_WEIGHT_TYPE:" )
214  filedata >> edgeweighttype;
215  else if( token == "EDGE_WEIGHT_TYPE" )
216  filedata >> token >> edgeweighttype;
217  else if( token == "NODE_COORD_SECTION" || token == "DISPLAY_DATA_SECTION" )
218  {
219  // there should be some nodes to construct a graph
220  if( nnodes < 1 )
221  {
222  retcode = SCIP_READERROR;
223  break;
224  }
225  // the graph is created and filled with nodes
226  else if( create_graph(nnodes, nedges, &graph) )
227  {
228  assert(x_coords == NULL);
229  assert(y_coords == NULL);
230 
231  x_coords = new double[nnodes];
232  y_coords = new double[nnodes];
233  getNodesFromFile(filedata, x_coords, y_coords, graph);
234  }
235  else
236  {
237  retcode = SCIP_NOMEMORY;
238  break;
239  }
240  }
241  else if( token == "COMMENT:" || token == "COMMENT" ||
242  token == "DISPLAY_DATA_TYPE" || token == "DISPLAY_DATA_TYPE:" )
243  (void) getline( filedata, token );
244  else if( token == "EOF" )
245  break;
246  else if( token == "" )
247  ;
248  else
249  {
250  cout << "parse error in file <" << name << "> unknown keyword <" << token << ">" << endl;
251  return SCIP_READERROR;
252  }
253  filedata >> token;
254  }// finished parsing the input
255 
256  // check whether the input data was valid
257  if( !checkValid(graph, name, type, edgeweighttype, nnodes) )
258  retcode = SCIP_READERROR;
259 
260  assert(graph != NULL);
261 
262  if( retcode == SCIP_OKAY )
263  {
264  edgeforw = &( graph->edges[0] ); /*lint !e613*/
265  edgebackw= &( graph->edges[nedges/2] ); /*lint !e613*/
266 
267 #ifdef SCIP_DEBUG
268  weights = new double* [nnodes];
269  for( i = 0; i < nnodes; ++i )
270  weights[i] = new double[nnodes];
271 #endif
272 
273  // construct all edges in a complete digraph
274  for( i = 0; i < nnodes; i++ )
275  {
276  nodestart = &graph->nodes[i]; /*lint !e613*/
277  for( j = i+1; j < nnodes; j++ )
278  {
279  nodeend = &graph->nodes[j]; /*lint !e613*/
280 
281  // construct two 'parallel' halfedges
282  edgeforw->adjac = nodeend;
283  edgebackw->adjac = nodestart;
284  edgeforw->back = edgebackw;
285  edgebackw->back = edgeforw;
286 
287  // calculate the Euclidean / Manhattan / Maximum distance of the two nodes
288  x = x_coords[(*nodestart).id] - x_coords[(*nodeend).id]; /*lint !e613*/
289  y = y_coords[(*nodestart).id] - y_coords[(*nodeend).id]; /*lint !e613*/
290  if( edgeweighttype == "EUC_2D")
291  edgeforw->length = sqrt( x*x + y*y );
292  else if( edgeweighttype == "MAX_2D")
293  edgeforw->length = MAX( ABS(x), ABS(y) );
294  else if( edgeweighttype == "MAN_2D")
295  edgeforw->length = ABS(x) + ABS(y);
296  else if( edgeweighttype == "ATT")
297  edgeforw->length = ceil( sqrt( (x*x+y*y)/10.0 ) );
298  else if( edgeweighttype == "GEO")
299  {
300  const double pi = 3.141592653589793;
301  double rads[4];
302  double coords[4];
303  double degs[4];
304  double mins[4];
305  double euler[3];
306  int k;
307 
308  coords[0] = x_coords[(*nodestart).id]; /*lint !e613*/
309  coords[1] = y_coords[(*nodestart).id]; /*lint !e613*/
310  coords[2] = x_coords[(*nodeend).id]; /*lint !e613*/
311  coords[3] = y_coords[(*nodeend).id]; /*lint !e613*/
312 
313  for( k = 0; k < 4; k++ )
314  {
315  degs[k] = coords[k] >= 0 ? floor(coords[k]) : ceil(coords[k]);
316  mins[k] = coords[k] - degs[k];
317  rads[k] = pi*(degs[k]+5.0*mins[k]/3.0)/180.0;
318  }
319 
320  euler[0] = cos(rads[1]-rads[3]);
321  euler[1] = cos(rads[0]-rads[2]);
322  euler[2] = cos(rads[0]+rads[2]);
323  edgeforw->length = floor(6378.388 * acos(0.5*((1.0+euler[0])*euler[1]-(1.0-euler[0])*euler[2]))+1.0);
324  }
325 
326  // in TSP community, it is common practice to round lengths to next integer
327  if( round_lengths_ )
328  edgeforw->length = NINT(edgeforw->length);
329 
330  edgebackw->length = edgeforw->length;
331 #ifdef SCIP_DEBUG
332  weights[i][j] = edgeforw->length;
333  weights[j][i] = edgebackw->length;
334 #endif
335 
336  // insert one of the halfedges into the edge list of the node
337  if (nodestart->first_edge == NULL)
338  {
339  nodestart->first_edge = edgeforw;
340  nodestart->first_edge->next = NULL;
341  }
342  else
343  {
344  edgeforw->next = nodestart->first_edge;
345  nodestart->first_edge = edgeforw;
346  }
347 
348  // dito
349  if (nodeend->first_edge == NULL)
350  {
351  nodeend->first_edge = edgebackw;
352  nodeend->first_edge->next = NULL;
353  }
354  else
355  {
356  edgebackw->next = nodeend->first_edge;
357  nodeend->first_edge = edgebackw;
358  }
359 
360  edgeforw++;
361  edgebackw++;
362  }
363  }
364  }
365 
366  delete[] y_coords;
367  delete[] x_coords;
368 
369  if( retcode != SCIP_OKAY )
370  {
371 #ifdef SCIP_DEBUG
372  if( weights != NULL )
373  {
374  for( i = 0; i < nnodes; i++ )
375  {
376  delete[] weights[i];
377  }
378  delete[] weights;
379  }
380 #endif
381  return retcode;
382  }
383 
384 #ifdef SCIP_DEBUG
385  printf("Matrix:\n");
386  for( i = 0; i < nnodes; i++ )
387  {
388  for( j = 0; j < nnodes; j++ )
389  printf(" %4.0f ",weights[i][j]);
390  printf("\n");
391  delete[] weights[i];
392  }
393  delete[] weights;
394 #endif
395 
396  // create the problem's data structure
397  SCIP_CALL( SCIPcreateObjProb(scip, name.c_str(), new ProbDataTSP(graph), TRUE) );
398 
399  // add variables to problem and link them for parallel halfedges
400  for( i = 0; i < nedges/2; i++ )
401  {
402  SCIP_VAR* var;
403 
404  stringstream varname;
405  edge = &graph->edges[i]; /*lint !e613*/
406 
407 /**! [SnippetTSPVariableCreation] */
408 
409  // the variable is named after the two nodes connected by the edge it represents
410  varname << "x_e_" << edge->back->adjac->id+1 << "-" << edge->adjac->id+1;
411  SCIP_CALL( SCIPcreateVar(scip, &var, varname.str().c_str(), 0.0, 1.0, edge->length,
413 
414  /* add variable to SCIP and to the graph */
415  SCIP_CALL( SCIPaddVar(scip, var) );
416  SCIP_CALL( addVarToEdges(scip, edge, var) );
417 
418  /* release variable for the reader. */
419  SCIP_CALL( SCIPreleaseVar(scip, &var) );
420 
421 /**! [SnippetTSPVariableCreation] */
422 
423  }
424 
425  /* add all n node degree constraints */
426  if( nnodes >= 2 )
427  {
428 
429  for( i = 0, node = &(graph->nodes[0]); i < nnodes; i++, node++ ) /*lint !e613*/
430  {
431 /**! [SnippetTSPDegreeConstraintCreation] */
432  SCIP_CONS* cons;
433  stringstream consname;
434  consname << "deg_con_v" << node->id+1;
435 
436  // a new degree constraint is created, named after a node
437  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, consname.str().c_str(), 0, NULL, NULL, 2.0, 2.0,
439 
440  edge = node->first_edge;
441  // sum up the values of all adjacent edges
442  while( edge != NULL )
443  {
444  SCIP_CALL( SCIPaddCoefLinear(scip, cons, edge->var, 1.0) );
445  edge = edge->next;
446  }
447 
448  // add the constraint to SCIP
449  SCIP_CALL( SCIPaddCons(scip, cons) );
450  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
451 /**! [SnippetTSPDegreeConstraintCreation] */
452  }
453  }
454 
455 /**! [SnippetTSPNosubtourConstraintCreation] */
456 
457  /* last, we need a constraint forbidding subtours */
458  SCIP_CONS* cons;
459  SCIP_CALL( SCIPcreateConsSubtour(scip, &cons, "subtour", graph,
460  FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE ) );
461  SCIP_CALL( SCIPaddCons(scip, cons) );
462  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
463 
464 /**! [SnippetTSPNosubtourConstraintCreation] */
465 
466  release_graph(&graph);
467  *result = SCIP_SUCCESS;
468 
469  return SCIP_OKAY;
470 } /*lint !e715*/
471 
472 /** problem writing method of reader; NOTE: if the parameter "genericnames" is TRUE, then
473  * SCIP already set all variable and constraint names to generic names; therefore, this
474  * method should always use SCIPvarGetName() and SCIPconsGetName();
475  *
476  * possible return values for *result:
477  * - SCIP_SUCCESS : the reader read the file correctly and created an appropritate problem
478  * - SCIP_DIDNOTRUN : the reader is not responsible for given input file
479  *
480  * If the reader detected an error in the writing to the file stream, it should return
481  * with RETCODE SCIP_WRITEERROR.
482  */
483 SCIP_DECL_READERWRITE(ReaderTSP::scip_write)
484 {
485  *result = SCIP_DIDNOTRUN;
486 
487  return SCIP_OKAY;
488 } /*lint !e715*/
struct GraphEdge * next
Definition: GomoryHuTree.h:58
SCIP_RETCODE SCIPcreateConsLinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, SCIP_VAR **vars, SCIP_Real *vals, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool initial, SCIP_Bool separate, SCIP_Bool enforce, SCIP_Bool check, SCIP_Bool propagate, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool dynamic, SCIP_Bool removable, SCIP_Bool stickingatnode)
#define NULL
Definition: def.h:253
double length
Definition: GomoryHuTree.h:56
Definition: grph.h:57
SCIP_Bool create_graph(int n, int m, GRAPH **gr)
SCIP_RETCODE SCIPcaptureVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_var.c:1217
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
SCIP_DECL_READERFREE(ReaderTSP::scip_free)
Definition: ReaderTSP.cpp:137
SCIPInterval cos(const SCIPInterval &x)
#define FALSE
Definition: def.h:73
#define TRUE
Definition: def.h:72
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
SCIP_DECL_READERWRITE(ReaderTSP::scip_write)
Definition: ReaderTSP.cpp:483
generator for global cuts in undirected graphs
SCIP_VAR ** x
Definition: circlepacking.c:54
SCIP_VAR * var
Definition: GomoryHuTree.h:63
SCIP_RETCODE SCIPcreateConsSubtour(SCIP *scip, SCIP_CONS **cons, const char *name, GRAPH *graph, SCIP_Bool initial, SCIP_Bool separate, SCIP_Bool enforce, SCIP_Bool check, SCIP_Bool propagate, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool dynamic, SCIP_Bool removable)
SCIPInterval sqrt(const SCIPInterval &x)
Definition: pqueue.h:28
GRAPHNODE * adjac
Definition: GomoryHuTree.h:61
struct GraphEdge * back
Definition: GomoryHuTree.h:59
struct GraphEdge * first_edge
Definition: GomoryHuTree.h:43
C++ wrapper classes for SCIP.
#define SCIP_CALL(x)
Definition: def.h:365
C++ problem data for TSP.
C++ constraint handler for TSP subtour elimination constraints.
SCIP_RETCODE SCIPcreateVar(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype, SCIP_Bool initial, SCIP_Bool removable, SCIP_DECL_VARDELORIG((*vardelorig)), SCIP_DECL_VARTRANS((*vartrans)), SCIP_DECL_VARDELTRANS((*vardeltrans)), SCIP_DECL_VARCOPY((*varcopy)), SCIP_VARDATA *vardata)
Definition: scip_var.c:104
#define NINT(x)
Definition: ReaderTSP.cpp:42
Constraint handler for linear constraints in their most general form, .
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1667
C++ file reader for TSP data files.
#define MAX(x, y)
Definition: def.h:222
double x
Definition: GomoryHuTree.h:35
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1251
SCIP_RETCODE SCIPcreateObjProb(SCIP *scip, const char *name, scip::ObjProbData *objprobdata, SCIP_Bool deleteobject)
SCIP_VAR ** y
Definition: circlepacking.c:55
SCIP_DECL_READERREAD(ReaderTSP::scip_read)
Definition: ReaderTSP.cpp:150
double y
Definition: GomoryHuTree.h:36
int edges
Definition: grph.h:92
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2765
void release_graph(GRAPH **gr)
#define nnodes
Definition: gastrans.c:65
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1109
#define ABS(x)
Definition: def.h:218