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