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