OSDipBlockBearcatSolver.cpp
Go to the documentation of this file.
1 /* $Id: OSBlockBearcatSolver.cpp 3038 2009-11-07 11:43:44Z kmartin $ */
16 #include "OSErrorClass.h"
17 #include "OSDataStructures.h"
18 
19 #include <sstream>
20 using std::ostringstream;
21 
22 
23 
24 
26  m_osinstance(NULL) {
27  std::cout << "INSIDE OSDipBlockBearcatSolver CONSTRUCTOR" << std::endl;
28 }//end OSDipBlockBearcatSolver
29 
30 
32  std::cout << "INSIDE OSDipBlockBearcatSolver CONSTRUCTOR" << std::endl;
33  std::cout << "whichBlock = " << m_whichBlock<< std::endl;
34  try{
37 
39 
40  m_demand = NULL;
41  m_u = NULL;
42  m_v = NULL;
43  m_g = NULL;
44  m_px = NULL;
45  m_tx =NULL;
46  m_varIdx = NULL;
47 
48  //now get data
49  std::vector<SolverOption*> solverOptions;
50  std::vector<SolverOption*>::iterator vit;
51  std::vector<int >demand;
52 
53 
54  solverOptions = m_osoption->getSolverOptions("Dip");
55  //iterate over the vector
56 
57  int tmpVal;
58 
59  for (vit = solverOptions.begin(); vit != solverOptions.end(); vit++) {
60 
61 
62  //std::cout << (*vit)->name << std::endl;
63 
64  //(*vit)->name.compare("initialCol") == 0
65  //if(rowNames[ i3].find("routeCapacity(1)") == string::npos )
66 
67  if( (*vit)->name.find("numHubs") != std::string::npos){
68 
69  std::istringstream buffer( (*vit)->value);
70  buffer >> m_numHubs;
71  std::cout << "numHubs = " << m_numHubs << std::endl;
72 
73  }else{
74 
75  if((*vit)->name.find("numNodes") != std::string::npos){
76 
77  std::istringstream buffer( (*vit)->value);
78  buffer >> m_numNodes;
79  std::cout << "numNodes = " << m_numNodes << std::endl;
80 
81  }else{
82  if((*vit)->name.find("totalDemand") != std::string::npos){
83 
84  std::istringstream buffer( (*vit)->value);
85  buffer >> m_totalDemand;
86  std::cout << "m_totalDemand = " << m_totalDemand << std::endl;
87 
88  }else{
89  if((*vit)->name.find("minDemand") != std::string::npos){
90 
91  std::istringstream buffer( (*vit)->value);
92  buffer >> m_minDemand;
93  std::cout << "m_minDemand = " << m_minDemand << std::endl;
94 
95  }else{
96  if( (*vit)->name.find("demand") != std::string::npos ){
97 
98  std::istringstream buffer( (*vit)->value);
99  buffer >> tmpVal;
100  demand.push_back( tmpVal);
101  //std::cout << "demand = " << tmpVal << std::endl;
102 
103  }
104  }
105  }
106  }
107  }//end if on solver options
108 
109  }//end for loop on options
110 
111  //now fill in demand
112  m_demand = new int[ m_numNodes];
113 
114  m_varIdx = new int[ m_numNodes];
115 
116  m_u = new double*[ m_numNodes];
117  m_v = new double*[ m_numNodes];
118  m_g = new double*[ m_numNodes];
119 
120  m_px = new int*[ m_numNodes];
121  m_tx = new int*[ m_numNodes];
122 
123  if(m_numNodes != (int)demand.size( ) ) throw ErrorClass("inconsistent number of demand nodes");
124  int i;
125  for (i = 0; i < m_numNodes; i++) {
126 
127  m_demand[ i] = demand[i];
128  m_u[ i] = new double[ m_totalDemand + 1];
129  m_v[ i] = new double[ m_totalDemand + 1];
130  m_g[ i] = new double[ m_numNodes];
131 
132  m_px[ i] = new int[ m_totalDemand + 1];
133  m_tx[ i] = new int[ m_totalDemand + 1];
134 
135 
136  }
137 
138 
139  } catch (const ErrorClass& eclass) {
140 
141  throw ErrorClass(eclass.errormsg);
142 
143  }
144 
145 }//end OSDipBlockBearcatSolver
146 
148 
149  std::cout << "INSIDE ~OSDipBlockBearcatSolver()" << std::endl;
150 
151 
152  std::vector<IndexValuePair*>::iterator vit;
153 
154  for (vit = m_primalVals.begin(); vit != m_primalVals.end(); vit++) {
155 
156  delete *vit;
157  }
158 
159  //delete data structures for arrays used in calculating minimum reduced cost
160  int i;
161 
162  for(i = 0; i < m_numNodes; i++){
163 
164  delete[] m_u[i];
165  m_u[i] = NULL;
166 
167  delete[] m_v[i];
168  m_v[i] = NULL;
169 
170  delete[] m_g[i];
171  m_g[i] = NULL;
172 
173  delete[] m_px[i];
174  m_px[i] = NULL;
175 
176  delete[] m_tx[i];
177  m_tx[i] = NULL;
178  }
179 
180  delete[] m_u;
181  m_u= NULL;
182 
183  delete[] m_v;
184  m_v= NULL;
185 
186  delete[] m_g;
187  m_g= NULL;
188 
189  delete[] m_px;
190  m_px= NULL;
191 
192  delete[] m_tx;
193  m_tx= NULL;
194 
195 
196 
197  if(m_demand != NULL) delete[] m_demand;
198 
199  if(m_varIdx != NULL) delete[] m_varIdx;
200  //if(m_osrlreader != NULL) delete m_osrlreader;
201 }//end ~OSDipBlockBearcatSolver
202 
203 void OSDipBlockBearcatSolver::solve(double *cost, std::vector<IndexValuePair*> *solIndexValPair, double *optVal){
204 
205  try{
206 
207  int i;
208 
209 
210 
211  struct IndexValuePair *primalValPair;
212 
213  std::vector<IndexValuePair*>::iterator vit;
214 
215  for (vit = m_primalVals.begin(); vit != m_primalVals.end(); vit++) {
216 
217  delete *vit;
218  }
219  m_primalVals.clear();
220 
221 
222  //now, lets find the minimum reduced cost
223 
224  *optVal = OSDBL_MAX;
225 
226  //std::cout << "Number of Var " << m_numberOfVar << std::endl;
227 
228  int l ;
229 
230  int bestl;
231  int kountVar;
232  double testVal;
233 
234  for(l = m_minDemand; l <= m_totalDemand; l++){
235 
236  testVal = qrouteCost(m_whichBlock, l, cost, &kountVar);
237  if( testVal < *optVal){
238  *optVal = testVal;
239  bestl = l;
240  }
241 
242 
243 
244  }
245 
246  *optVal = qrouteCost(m_whichBlock, bestl, cost, &kountVar);
247  std::cout << "best reduced cost = " << *optVal << std::endl;
248 
249  std::map<int, int> indexCount;
250 
251  std::map<int, int>::iterator mit;
252 
253 
254  for(i = 0; i < kountVar; i++){
255 
256  if( indexCount.find( m_varIdx[ i]) == indexCount.end() ){
257 
258  indexCount[ m_varIdx[ i]] = 1;
259 
260  }else{
261 
262  indexCount[ m_varIdx[ i]] += 1;
263 
264  }
265 
266 
267  }
268 
269  for (mit = indexCount.begin(); mit != indexCount.end(); mit++){
270 
271  //std::cout << "Variable Index " << mit->first << " Count = " << mit->second << std::endl;
272 
273  //std::cout << "Variable Name = " << m_osinstance->getVariableNames()[ mit->first] << std::endl;
274 
275  // get unique indexes
276  primalValPair = new IndexValuePair();
277 
278  primalValPair->value = mit->second;
279  primalValPair->idx = mit->first;
280 
281 
282  m_primalVals.push_back( primalValPair);
283  }
284 
285  /*** old way could lead to duplicate indexes
286  for (i = 0; i < kountVar; i++){
287 
288  //std::cout << "Variable Number = " << m_varIdx[ i] << std::endl;
289  std::cout << "Variable Name = " << m_osinstance->getVariableNames()[ m_varIdx[ i]] << std::endl;
290 
291  // get unique indexes
292  primalValPair = new IndexValuePair();
293 
294  primalValPair->value = 1.0;
295  primalValPair->idx = m_varIdx[ i];
296 
297 
298  m_primalVals.push_back( primalValPair);
299 
300 
301  }
302  */
303 
304  *solIndexValPair = m_primalVals;
305 
306 
307  } catch (const ErrorClass& eclass) {
308 
309  throw ErrorClass(eclass.errormsg);
310 
311  }
312 
313 }//end solve
314 
315 
316 void OSDipBlockBearcatSolver::solve(double *cost, std::string *osrl){
317 
318 
319  try{
320  //not implemented
321 
322  } catch (const ErrorClass& eclass) {
323 
324  throw ErrorClass(eclass.errormsg);
325 
326  }
327 
328 }//end solve
329 
330 
331 
332 double OSDipBlockBearcatSolver::qrouteCost(const int& k, const int& l, double* c, int* kountVar){
333 
334  //critical -- nodes 0, ..., m_numNodes - 1 are the hub nodes
335  // we are doing the calculation for hub k, k <= m_numNodes - 1
336  //l should be the total demand on the network -- we are NOT using 0 based counting
337  double rcost;
338  rcost = OSDBL_MAX;
339 
340  //start of the cost vector for hub k plus move to start of l demands
341  // now move the pointer up
342  //int startPnt = k*totalDemand*(m_numNodes*m_numNodes - m_numNodes) + (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
343 
344  int startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
345  c+= startPnt ;
346 
347  *kountVar = 0;
348  int bestLastNode = 0;
349  int i;
350  int j;
351  int l1;
352  int l2;
353  //for(i = 0; i < 20; i++){
354  // std::cout << "i = " << i << " c[i] = " << *(c + i) << std::endl ;
355  //}
356 
357 
358 
359  // l is the total demand on this network
360  //address of node (j, i) is j*(m_numNodes-1) + i when i < j
361  //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
362 
363  //
364  // initialize
365 
366  for(i = m_numHubs; i < m_numNodes; i++){
367 
368 
369  for(l1 = m_minDemand; l1 <= l; l1++){ //l-1 is total demand on network
370 
371  m_u[i][l1] = OSDBL_MAX;
372  m_v[i][l1] = OSDBL_MAX;
373 
374  if(l1 == *(m_demand + i) ){
375 
376  m_px[i][l1] = k;
377  // want the cost for arc (k, i)
378  // note: k < i so use that formula
379  m_u[i][l1] = *(c + k*(m_numNodes - 1) + i - 1); // put in correct cost
380 
381 
382  }
383  }
384  }
385  //end initialize
386  //
387 
388  //if l = minDemand we visit only one node, let's get the reduced cost in this case
389  if(l == m_minDemand){
390 
391  for(i = m_numHubs; i < m_numNodes; i++){
392 
393 
394  if( m_u[i][l] + *(c + i*(m_numNodes-1) + k ) < rcost){
395 
396  rcost = m_u[i][l] + *(c + i*(m_numNodes-1) + k );
397 
398  //std::cout << " m_u[i][l2] = " << m_u[i][l2] << std::endl;
399 
400  //std::cout << " *(c + i*(m_numNodes-1) + k ) = " << *(c + i*(m_numNodes-1) + k ) << std::endl;
401  //push node back
402  bestLastNode = i;
403  }
404 
405  }
406 
407  //go from node bestLastNode to node k
408  //node bestLastNode is a higher number than k
409  *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
410  *(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1) + bestLastNode - 1;
411 
412 
413  return rcost;
414  }//end if on l == minDemand
415 
416 // now calculate values for demand 2 or greater
417  //address of node (j, i) is j*(m_numNodes-1) + i when i < j
418  //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
419  // we start l2 at 2 since demand must be at least 1
420  // change to min demand + 1
421  int lowerVal = m_minDemand + 1;
422  for(l2 = lowerVal; l2 <= l; l2++){// loop over possible demand values assuming we have already gone to at least one node
423 
424  for(i = m_numHubs; i < m_numNodes; i++) { //we are finding least cost to node i
425 
426 
427  if( *(m_demand + i) < l2 ){ // kipp < OR <= ????
428 
429  for(j = m_numHubs; j < i; j++){ //we are coming from node j
430 
431 
432 
433  //calculate g -- l2 - d[ i] is the demand trough and including node j
434  l1 = l2 - *(m_demand + i);
435 
436  if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
437 
438 
439  m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i - 1) ;
440 
441 
442 
443 
444  }else{
445 
446  m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i - 1) ;
447 
448 
449 
450  }
451 
452  // update u and the pointer for p
453 
454  if(m_g[j][i] < m_u[i][l2] ){
455 
456  m_u[i][l2] = m_g[j][i];
457  m_px[i][l2] = j;
458 
459  }
460 
461 
462 
463  }//end of first for loop on j, j < i
464 
465 
466  for(j = i + 1; j < m_numNodes; j++){ //we are coming from node j
467 
468 
469  //calculate g -- l2 - d[ i] is the demand trough and including node j
470  l1 = l2 - *(m_demand + i);
471 
472  if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
473 
474 
475  m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i ) ;
476 
477 
478  }else{
479 
480  m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i ) ;
481 
482  }
483 
484  // update u and the pointer for p
485 
486  if(m_g[j][i] < m_u[i][l2] ){
487 
488  m_u[i][l2] = m_g[j][i];
489  m_px[i][l2] = j;
490 
491  }
492 
493 
494  }//end of second for loop on j, j > i
495 
496 
497  //now calculate the second best solution and point
498  //right now m_px[i][l2] points to the best j node
499 
500  for(j =m_numHubs; j < m_numNodes; j++){ //we are coming from node j
501 
502  if(j != i){
503 
504  if( (m_g[j][i] < m_v[i][l2] ) && (m_px[i][l2] != j) ){ // kipp the && gives the second best
505 
506  //if(g(j, i) < v(i, l2),
507 
508  m_v[i][l2] = m_g[j][i];
509  m_tx[i][l2] = j;
510 
511 
512  }
513 
514  }
515 
516 
517  }//end second best calculation, anothe for loop on j
518 
519  //now if l2 = l we are done
520  if(l2 == l ){
521 
522  if( m_u[i][l2] + *(c + i*(m_numNodes-1) + k ) < rcost){
523 
524  rcost = m_u[i][l2] + *(c + i*(m_numNodes-1) + k );
525 
526  bestLastNode = i;
527  }
528 
529  }
530 
531 
532  }//end if on demand less than l2
533 
534 
535  }//i loop
536 
537 
538  }//l2 loop
539 
540 
541  //std::cout << "best Last Node = " << bestLastNode << std::endl;
542 
543  // now get the path that gives the reduced cost
544 
545  int currentNode;
546  int successorNode;
547  int lvalue;
548 
549  //initialize
550  // we work backwords from bestLastNode
551  //in our recursion we recurse on the currentNode and assume
552  //it is in the optimal path
553 
554  //node bestLastNode is a higher number than k
555  *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
556 
557 
558 
559 
560  //by successor, I mean node after current node in the path
561  successorNode = k;
562  currentNode = bestLastNode;
563  //the lvalue is the demand through the currentNode
564  lvalue = l ;
565 
566 
567  while(currentNode != k){
568 
569 
570  if( m_px[ currentNode][ lvalue ] != successorNode){
571 
572 
573 
574  //update nodes
575  successorNode = currentNode;
576  currentNode = m_px[ currentNode][ lvalue ];
577 
578 
579  if(currentNode - successorNode > 0){
580  //below diagonal
581  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
582 
583  }else{
584  //above diagonal
585  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
586 
587  }
588 
589 
590  }else{ //take second best
591 
592 
593  //update nodes
594  successorNode = currentNode;
595  currentNode = m_tx[ currentNode][ lvalue ];
596 
597  if(currentNode - successorNode > 0){
598  //below diagonal
599  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
600 
601  }else{
602  //above diagonal
603  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
604 
605  }
606 
607  }
608 
609  //update lvalue
610  lvalue = lvalue - *(m_demand + successorNode);
611 
612 
613 
614  }
615 
616 
617  //go from node k to bestLastNode -- already done in loop above
618  //*(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1) + currentNode - 1;
619 
620 
621 
622 
623  return rcost;
624 
625 }//end qroute
626 
virtual void solve(double *cost, std::vector< IndexValuePair * > *solIndexValPair, double *optVal)
int m_whichBlock
m_whichBlock is the index of the subproblem we are working with
double value
value is a double that holds the value of the entity
Definition: OSGeneral.h:638
int getVariableNumber()
Get number of variables.
A commonly used structure holding an index-value pair.
Definition: OSGeneral.h:630
std::string errormsg
errormsg is the error that is causing the exception to be thrown
Definition: OSErrorClass.h:42
std::vector< SolverOption * > getSolverOptions(std::string solver_name)
Get the options associated with a given solver.
Definition: OSOption.cpp:4508
OSOption * osoption
std::vector< IndexValuePair * > m_primalVals
static char * j
Definition: OSdtoa.cpp:3622
The Option Class.
Definition: OSOption.h:3564
OSDipBlockBearcatSolver()
Default Constructor.
void fint fint * k
const double OSDBL_MAX
Definition: OSParameters.h:93
int idx
idx holds the index of an entity (such as a variable, constraint, objective) that is part of a sparse...
Definition: OSGeneral.h:635
The in-memory representation of an OSiL instance..
Definition: OSInstance.h:2262
OSInstance * osinstance
double qrouteCost(const int &k, const int &l, double *c, int *kountVar)
~OSDipBlockBearcatSolver()
Default destructor.
used for throwing exceptions.
Definition: OSErrorClass.h:31
real c