My Project
OSBearcatSolverXij.cpp
Go to the documentation of this file.
1 /* $Id: OSBearcatSolverXij.cpp 3038 2009-11-07 11:43:44Z kmartin $ */
14 #include "OSBearcatSolverXij.h"
15 
16 #include "OSErrorClass.h"
17 #include "OSDataStructures.h"
18 
19 #include "OSInstance.h"
20 #include "OSCoinSolver.h"
21 #include "OSConfig.h"
22 #include "OSResult.h"
23 
24 #include "OSiLReader.h"
25 #include "OSiLWriter.h"
26 #include "OSoLReader.h"
27 #include "OSoLWriter.h"
28 #include "OSrLReader.h"
29 #include "OSrLWriter.h"
30 #include "OSInstance.h"
31 #include "OSFileUtil.h"
32 #include "OSStringUtil.h"
33 
34 #include "CoinTime.hpp"
35 
36 #include "ClpFactorization.hpp"
37 #include "ClpNetworkMatrix.hpp"
38 #include "OsiClpSolverInterface.hpp"
39 
40 #include <iostream>
41 
42 #include <algorithm>
43 
44 
45 #ifdef HAVE_CMATH
46 # include <cmath>
47 #else
48 # ifdef HAVE_MATH_H
49 # include <math.h>
50 # else
51 # error "don't have header file for math"
52 # endif
53 #endif
54 
55 #ifdef HAVE_CTIME
56 # include <ctime>
57 #else
58 # ifdef HAVE_TIME_H
59 # include <time.h>
60 # else
61 # error "don't have header file for time"
62 # endif
63 #endif
64 
65 using std::ostringstream;
66 
67 
69  std::cout << "INSIDE OSBearcatSolverXij CONSTRUCTOR with OSOption argument" << std::endl;
70 }//end default OSBearcatSolverXij constructor
71 
73  std::cout << "INSIDE OSBearcatSolverXij CONSTRUCTOR with OSOption argument" << std::endl;
74 
75  m_hubPoint = NULL;
76 
79 
80  m_upperBoundL = NULL;
83 
84  m_u = NULL;
85  m_v = NULL;
86  m_g = NULL;
87  m_px = NULL;
88  m_tx =NULL;
89  m_varIdx = NULL;
90 
91  m_optL = NULL;
92  m_optD = NULL;
93  m_vv = NULL;
94  m_vvpnt = NULL;
95 
96  m_demand = NULL;
97  m_cost = NULL;
98 
99  m_rc = NULL;
100 
101  m_routeCapacity = NULL;
102  m_routeMinPickup = NULL;
103 
105 
106  m_multiCommodCutLimit = 250;
107  m_numMultCuts = 0;
108 
109 
110 }//end OSBearcatSolverXijDestructor
111 
113 
114  int k;
115  int i;
116  int l;
117  int tmpVal;
118 
119  try{
120 
121 
122  //get all the parameter values
124  // we need to get cost data from option file
125  if(m_cost == NULL) throw ErrorClass("Option file did not contain cost data");
126  //now get the variable index map
128 
129 
131 
132  m_upperBoundL = new int[ m_numHubs];
133  m_lowerBoundL = new int[ m_numHubs];
134 
135  m_hubPoint = new int[ m_numHubs];
136  //get the permutation of the hubs
137  permuteHubs();
138 
139 
140  for(k = 0; k < m_numHubs; k++){
141 
142  //adjust routeMinPickup
143  tmpVal = 0;
144  for(i = 0; i < m_numHubs; i++)
145  if( i != k) tmpVal += m_routeCapacity[ i];
146 
147  m_lowerBoundL[ k] = std::max( m_routeMinPickup[ k], (m_totalDemand - tmpVal) ) ;
148 
149  m_upperBoundL[ k] = m_routeCapacity[ k];
150 
151  //set m_upperBoundL cannot exceed total demand
154 
155  }
156 
157  //m_varIdx = new int[ m_numNodes];
158  // I think we can make this the max of m_upperBoundL
159  // over the k
160  m_varIdx = new int[ m_upperBoundLMax + 1];
161 
162 
163  m_u = new double*[ m_numNodes];
164  m_v = new double*[ m_numNodes];
165  m_g = new double*[ m_numNodes];
166 
167  m_px = new int*[ m_numNodes];
168  m_tx = new int*[ m_numNodes];
169 
170 
171 
181  for (i = 0; i < m_numNodes; i++) {
182 
183  //kipp we can use the biggest possible m_upperBoundL
184  m_u[ i] = new double[ m_upperBoundLMax + 1];
185  m_v[ i] = new double[ m_upperBoundLMax + 1];
186 
187 
188 
189  for(l = 0; l <= m_upperBoundLMax; l++){
190 
191  m_u[ i][l] = OSDBL_MAX;
192  m_v[ i][l] = OSDBL_MAX;
193  }
194 
195  m_g[ i] = new double[ m_numNodes];
196  m_px[ i] = new int[ m_upperBoundLMax + 1];
197  m_tx[ i] = new int[m_upperBoundLMax + 1];
198 
199 
200  }
201 
202 
203  //outer dynamic programming arrays
204  m_optL = new int[ m_numHubs];
205  m_optD = new int[ m_numHubs];
206 
207  m_vv = new double*[ m_numHubs];
208  m_vvpnt = new int*[ m_numHubs];
209  m_rc = new double*[ m_numHubs];
210 
211  for (k = 0; k < m_numHubs; k++) {
212 
213 
214  m_vv[ k] = new double[ m_totalDemand + 1];
215  m_vvpnt[ k] = new int[ m_totalDemand + 1];
216 
217  //allocate memory for the reduced cost vector.
218  //assume order is k, l, i, j
219  //assume order is l, i, j
220  m_rc[ k] = new double[ m_upperBoundL[ k]*(m_numNodes*m_numNodes - m_numNodes)];
221 
222 
223  }
224 
225  m_optValHub = new double[ m_numHubs];
226 
227  m_variableNames = new std::string[ m_numNodes*(m_numNodes - 1)];
228 
230 
231  //these are constraints that say we must be incident to each
232  //non-hub node -- there are m_numNodes - m_numHubs of these
233  m_pntAmatrix = new int[ m_numNodes - m_numHubs + 1];
234  //the variables -- the variable space we are in is x_{ij} NOT
235  // x_{ijk} -- a bit tricky
236  //m_Amatrix[ j] is a variable index -- this logic works
237  //since the Amatrix coefficient is 1 -- we don't need a value
238  //it indexes variable that points into the node
239  m_Amatrix = new int[ (m_numNodes - m_numHubs)*(m_numNodes - 1) ];
240  createAmatrix();
241 
242  //this has size of the number of x variables
243  //int numVar = m_numNodes*m_numNodes - m_numHubs ;
244  int numVar = m_numNodes*m_numNodes - m_numNodes ;
245  m_tmpScatterArray = new int[ numVar ];
246  for(i = 0; i < numVar; i++){
247 
248  m_tmpScatterArray[ i] = 0;
249 
250  }
251 
252  //New column arrays
253  m_newColumnNonz = new int[ m_numHubs] ; //at most one column per Hub
254  m_costVec = new double[ m_numHubs];
255  m_newColumnRowIdx = new int*[ m_numHubs];
256  m_newColumnRowValue = new double*[ m_numHubs];
257 
258  for (k = 0; k < m_numHubs; k++) {
259  //size of arrray is maximum number of constraints
260  // 1) coupling + 2) tour breaking + 3) branching
261  m_newColumnRowValue[ k] = new double[ m_maxBmatrixCon + m_numNodes];
262  m_newColumnRowIdx[ k] = new int[ m_maxBmatrixCon + m_numNodes];
263 
264  }
265 
266 
267 
268  //New row arrays
269  //m_newRowNonz = new int[ m_numHubs] ; //at most one cut per Hub
270  //m_newRowColumnIdx = new int*[ m_numHubs] ; //at most one cut per Hub
271  //m_newRowColumnValue = new double*[ m_numHubs] ; //at most one cut per Hub
272  //m_newRowUB = new double[ m_numHubs]; //at most one cut per Hub
273  //m_newRowLB = new double[ m_numHubs]; //at most one cut per Hub
274 
275 
276  m_newRowNonz = new int[ 100] ; //at most one cut per Hub
277  m_newRowColumnIdx = new int*[ 100] ; //at most one cut per Hub
278  m_newRowColumnValue = new double*[ 100] ; //at most one cut per Hub
279  m_newRowUB = new double[ 100]; //at most one cut per Hub
280  m_newRowLB = new double[ 100]; //at most one cut per Hub
281 
282 
283  //for now, the number of columns will be m_maxMasterColumns
284  //for (k = 0; k < m_numHubs; k++) {
285  for (k = 0; k < 100; k++) {
286 
287 
288  m_newRowColumnValue[ k] = new double[ m_maxMasterColumns];
289  m_newRowColumnIdx[ k] = new int[ m_maxMasterColumns];
290 
291  }
292 
293  //new array for keeping track of convexity rows
295  //new arrays for branches
296 
298  branchCutValues = new double[ m_maxMasterColumns];
299 
300  m_thetaPnt = new int[ m_maxMasterColumns + 1];
301  for(i = 0; i <= m_maxMasterColumns; i++){
302  m_thetaPnt[ i] = 0;
303  }
304  m_thetaCost = new double[ m_maxMasterColumns];
305  m_thetaIndex = new int[ m_maxThetaNonz];
306  m_numThetaVar = 0;
307  m_numThetaNonz = 0;
309 
310 
311  //the number of cuts will be at most nubmer of tour
312  // breaking constraints + braching variable cuts
313  // branching variable constraints = m_numNodes*(m_numNodes - 1)
314  m_pntBmatrix = new int[ m_maxBmatrixCon];
315  // number of nonzeros in the Bmatrix
316  m_BmatrixIdx = new int[ m_maxBmatrixNonz];
318  for(i = 0; i < m_maxBmatrixCon; i++) m_BmatrixRowIndex[ i] = -1;
319 
320  // number of nonzeros in the Bmatrix
321  m_BmatrixVal = new double[ m_maxBmatrixNonz];
322 
323 
324  m_numBmatrixCon = 0;
325  m_numBmatrixNonz = 0;
327 
328 ;
329 
330 
331  m_separationIndexMap = new int[m_numNodes*(m_numNodes - 1)];
332 
333  for(i = 0; i < m_numNodes*(m_numNodes - 1); i++){
334 
336 
337  }
338 
339 
340 
341 
342  //kipp -- move this later
344  for(k = 0; k < m_numHubs; k++){
345 
347 
348  }
349 
350  } catch (const ErrorClass& eclass) {
351 
352  throw ErrorClass(eclass.errormsg);
353 
354  }
355 
356 
357 }//end initializeDataStructures
358 
359 
361 
362  std::cout << "INSIDE ~OSBearcatSolverXij DESTRUCTOR" << std::endl;
363 
364 
365 
366  //delete data structures for arrays used in calculating minimum reduced cost
367  int i;
368 
369  delete[] m_routeCapacity;
370  m_routeCapacity = NULL;
371 
372 
373  delete[] m_routeMinPickup;
374  m_routeMinPickup = NULL;
375 
376  for(i = 0; i < m_numNodes; i++){
377 
378 
379 
380  delete[] m_v[i];
381  delete[] m_g[i];
382  delete[] m_px[i];
383  delete[] m_tx[i];
384  delete[] m_u[i];
385 
386 
387  }
388 
389  delete[] m_u;
390  m_u= NULL;
391 
392  delete[] m_v;
393  m_v= NULL;
394 
395  delete[] m_g;
396  m_g= NULL;
397 
398  delete[] m_px;
399  m_px= NULL;
400 
401  delete[] m_tx;
402  m_tx= NULL;
403 
404 
405 
406  if(m_demand != NULL){
407  //std::cout << "I AM DELETING m_demand" << std::endl;
408  delete[] m_demand;
409  }
410 
411 
412  if(m_varIdx != NULL){
413  delete[] m_varIdx;
414  m_varIdx = NULL;
415 
416  }
417 
418  if(m_cost != NULL ){
419  delete[] m_cost;
420  m_cost = NULL;
421  }
422 
423  for(i = 0; i < m_numHubs; i++){
424 
425  delete[] m_vv[i];
426  delete[] m_vvpnt[i];
427  delete[] m_rc[ i];
428 
429 
430  }
431  delete[] m_optL;
432  m_optL = NULL;
433  delete[] m_optD;
434  m_optD = NULL;
435  delete[] m_vv;
436  m_vv = NULL;
437  delete[] m_vvpnt;
438  m_vvpnt = NULL;
439 
440  delete[] m_rc;
441  m_rc = NULL;
442 
443  delete[] m_upperBoundL;
444  m_upperBoundL = NULL;
445 
446  delete[] m_lowerBoundL;
447  m_lowerBoundL = NULL;
448 
449  delete[] m_hubPoint;
450  m_hubPoint = NULL;
451 
452  delete[] m_optValHub;
453  m_optValHub = NULL;
454 
455  delete[] m_variableNames;
456  m_variableNames = NULL;
457 
458  delete[] m_pntAmatrix;
459  m_pntAmatrix = NULL;
460 
461  delete[] m_Amatrix;
462  m_Amatrix = NULL;
463 
464  delete[] m_tmpScatterArray;
465  m_tmpScatterArray = NULL;
466 
467  delete[] m_newColumnNonz ;
468  m_newColumnNonz = NULL;
469  delete[] m_costVec ;
470  m_costVec = NULL;
471 
472  for(i = 0; i < m_numHubs; i++){
473 
474  delete[] m_newColumnRowIdx[i];
475  delete[] m_newColumnRowValue[i];
476  }
477 
478  delete[] m_newColumnRowIdx;
479  m_newColumnRowIdx = NULL;
480 
481  delete[] m_newColumnRowValue;
482  m_newColumnRowValue = NULL;
483 
484 
485  delete[] m_convexityRowIndex;
486  m_convexityRowIndex = NULL;
487 
488 
489  //getCut arrays
490  delete[] m_newRowNonz;
491  m_newRowNonz = NULL;
492 
493  delete[] m_newRowUB ;
494  m_newRowUB = NULL;
495 
496  delete[] m_newRowLB ;
497  m_newRowLB = NULL;
498 
499  //garbage collection on row arrays
500  for (i = 0; i < m_numHubs; i++) {
501 
502  delete[] m_newRowColumnValue[ i];
503  delete[] m_newRowColumnIdx[ i];
504 
505  }
506 
507  delete[] m_newRowColumnIdx;
508  m_newRowColumnIdx = NULL;
509 
510  delete[] m_newRowColumnValue;
511  m_newRowColumnValue = NULL;
512 
513 
514  delete[] branchCutIndexes ;
515  branchCutIndexes = NULL;
516 
517  delete[] branchCutValues ;
518  branchCutValues = NULL;
519 
520 
521  delete[] m_thetaPnt;
522  m_thetaPnt = NULL;
523 
524  delete[] m_thetaIndex;
525  m_thetaIndex = NULL;
526 
527 
528  delete[] m_thetaCost;
529  m_thetaCost = NULL;
530 
531 
532  delete[] m_pntBmatrix ;
533  m_pntBmatrix = NULL;
534 
535  delete[] m_BmatrixIdx ;
536  m_BmatrixIdx = NULL;
537 
538  delete[] m_BmatrixVal ;
539  m_BmatrixVal = NULL;
540 
541  delete[] m_BmatrixRowIndex ;
542  m_BmatrixRowIndex = NULL;
543 
544 
545 
546 
547  delete[] m_separationIndexMap;
548  m_separationIndexMap = NULL;
549 
550  delete m_separationClpModel;
551  m_separationClpModel = NULL;
552 
553  delete m_osinstanceSeparation;
554  m_osinstanceSeparation = NULL;
555 
556  std::vector<CoinSolver*>::iterator vit;
557  for(vit = m_multCommodCutSolvers.begin(); vit < m_multCommodCutSolvers.end(); vit++){
558 
559  delete *vit;
560  }
561 
562 }//end ~OSBearcatSolverXij
563 
564 
565 
566 
567 
568 
569 
570 void OSBearcatSolverXij::getOptL( double** c) {
571 
572  //initialize the first HUB
573 
574  int k;
575  int d;
576  int d1;
577  int kountVar;
578  double testVal;
579  int l;
580  //int startPntInc;
581  double trueMin;
582 
583  bool isFeasible;
584  isFeasible = false;
585 
586  kountVar = 0;
587  //startPntInc = m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes);
588  //m_hubPoint[k] is the pointer
589 
590  m_vv[ m_hubPoint[0] ][ 0] = 0;
591  for(d = 1; d <= m_totalDemand; d++){
592 
593  m_vv[ m_hubPoint[0] ][ d] = OSDBL_MAX;
594 
595  }
596  //initialize up to last hub
597  for(k = 1; k < m_numHubs - 1; k++){
598  for(d = 0; d <= m_totalDemand; d++){
599 
600  m_vv[ m_hubPoint[ k] ][ d] = OSDBL_MAX;
601 
602  }
603  }
604 
605 
606  //now loop over the other HUBS
607 
608  int dlower;
609  dlower = 0;
610 
611  for(k = 1; k < m_numHubs; k++){
612 
613  dlower += m_lowerBoundL[ m_hubPoint[k - 1] ];
614 
615  //kipp make d the min demand for the previous routes
616  for(d = dlower; d <= m_totalDemand; d++){
617 
618  m_vv[ m_hubPoint[ k] ][ d] = OSDBL_MAX;
619 
620  //d1 is the state variable at stage k -1
621  //for(d1 = 0; d1 <= m_totalDemand; d1++){
622  for(d1 = 0; d1 <= d; d1++){
623 
624  l = d - d1;
625 
626  //std::cout << "L = " << l << " m_upperBoundL[ k - 1] " << m_upperBoundL[ k - 1] << std::endl;
627  //kipp make m_upperBoundL the route capapcity
628  if( (m_vv[ m_hubPoint[ k - 1] ][ d1] < OSDBL_MAX) && (l <= m_upperBoundL[ m_hubPoint[ k - 1] ]) && (l >= m_lowerBoundL[ m_hubPoint[k - 1] ]) ){
629 
630  //std::cout << "k - 1 = " << k - 1 << " L = " << l << " m_upperBoundL[ k - 1] " << m_upperBoundL[ k - 1] << std::endl;
631  // l was the decision at state d1 in stage k-1
632  // l + d1 brings us to state d at stage k
633  // d is the total carried on routes 0 -- k-1
634 
635  testVal = qrouteCost( m_hubPoint[ k - 1], l, c[ m_hubPoint[ k - 1] ], &kountVar);
636 
637  //std::cout << "L = " << l << std::endl;
638  //std::cout << "testVal " << testVal << std::endl;
639 
640  if( m_vv[ m_hubPoint[ k-1] ][ d1] + testVal < m_vv[ m_hubPoint[ k] ][ d] ){
641 
642  m_vv[ m_hubPoint[ k] ][ d] = m_vv[ m_hubPoint[ k-1] ][ d1] + testVal;
643  //now point to the best way to get to d
644  m_vvpnt[ m_hubPoint[ k] ][ d] = d1;
645 
646  }
647 
648 
649  }
650 
651  }
652 
653  }
654 
655  //c+= startPntInc ;
656 
657  }// end for on k
658 
659  trueMin = OSDBL_MAX;
660  //we now enter the last stage through the other hubs
661  // have satisfied total demand d
662 
663 
664  //if (m_numHubs > 1) dlower = 1;
665 
666  //std::cout << " dlower = " << dlower << " m_totalDemand = " << m_totalDemand << std::endl;
667 
668  for(d = dlower; d < m_totalDemand; d++){
669 
670  //std::cout << "m_vv[ m_numHubs - 1 ][ d] " << m_vv[ m_numHubs - 1 ][ d] << std::endl;
671  l = m_totalDemand - d;
672 
673  if(m_vv[ m_hubPoint[ m_numHubs - 1] ][ d] < OSDBL_MAX && l <= m_upperBoundL[ m_hubPoint[ m_numHubs - 1] ] && l >= m_lowerBoundL[ m_hubPoint[ m_numHubs - 1] ]){
674 
675  //must execute this loop at least once
676 
677  //std::cout << "LL = " << l << std::endl;
678 
679  isFeasible = true;
680 
681 
682  testVal = qrouteCost(m_hubPoint[m_numHubs -1] , l, c[ m_hubPoint[ m_numHubs -1] ], &kountVar);
683 
684  //std::cout << "l = " << l << std::endl;
685  //std::cout << "testVal = " << testVal << std::endl;
686 
687  if(m_vv[ m_hubPoint[ m_numHubs - 1] ][ d] + testVal < trueMin){
688 
689  trueMin = m_vv[ m_hubPoint[ m_numHubs -1] ][ d] + testVal;
690  m_optD[ m_hubPoint[ m_numHubs -1] ] = d;
691  m_optL[ m_hubPoint[ m_numHubs -1] ] = l;
692 
693  }
694 
695 
696  }
697  }
698 
699  //std::cout << "TRUE MIN = " << trueMin << std::endl;
700 
701  if( isFeasible == false){
702 
703  std::cout << "NOT ENOUGH CAPACITY " << std::endl;
704  for(k = 0; k < m_numHubs; k++) std::cout << " k perm = " << m_hubPoint[ k ]<< std::endl;
705  throw ErrorClass( "NOT ENOUGH CAPACITY ");
706  }
707 
708  k = m_numHubs - 1;
709 
710  while( k - 1 >= 0) {
711 
712  m_optD[ m_hubPoint[ k - 1] ] = m_vvpnt[ m_hubPoint[ k] ][ m_optD[ m_hubPoint[ k] ] ];
713 
714  m_optL[ m_hubPoint[ k - 1] ] = m_optD[ m_hubPoint[ k ] ] - m_optD[ m_hubPoint[ k - 1] ] ;
715 
716  k--;
717 
718  }
719 
720 }//end getOptL
721 
722 
723 
724 
725 
726 
727 double OSBearcatSolverXij::qrouteCost(const int& k, const int& l, const double* c, int* kountVar){
728 
729  //critical -- nodes 0, ..., m_numNodes - 1 are the hub nodes
730  // we are doing the calculation for hub k, k <= m_numNodes - 1
731  //l should be the total demand on the network -- we are NOT using 0 based counting
732  double rcost;
733  rcost = OSDBL_MAX;
734 
735 
736 
737  if(l <= 0){
738 
739  std::cout << "LVALUE NEGATIVE OR ZERO " << l << std::endl;
740  exit( 1);
741  }
742 
743 
744 
745  if(l > m_upperBoundL[ k]){
746 
747  std::cout << "LVALUE BIGGER THAN UPPER BOUND " << l << std::endl;
748  exit( 1);
749  }
750 
751  //start of the cost vector for hub k plus move to start of l demands
752  // now move the pointer up
753  //int startPnt = m_upperBoundL[ k]*(m_numNodes*m_numNodes - m_numNodes) + (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
754 
755  int startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
756  c+= startPnt ;
757 
758 
759 
760  *kountVar = 0;
761  int bestLastNode;
762  //initialize
763  bestLastNode = OSINT_MAX;
764  int i;
765  int j;
766  int l1;
767  int l2;
768  //for(i = 0; i < 20; i++){
769  // std::cout << "i = " << i << " c[i] = " << *(c + i) << std::endl ;
770  //}
771 
772 
773 
774  // l is the total demand on this network
775  //address of node (j, i) is j*(m_numNodes-1) + i when i < j
776  //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
777 
778  //
779  // initialize
780 
781 
782  for(i = m_numHubs; i < m_numNodes; i++){
783 
784 
785  for(l1 = m_minDemand; l1 <= l; l1++){ //l-1 is total demand on network
786  //std::cout << "HUB : " << k << " i = " << i << " l1 " << l1 << std::endl;
787  //std::cout << "m_upperBoundLMax: " << m_upperBoundLMax << std::endl;
788  m_u[i][l1] = OSDBL_MAX;
789  //std::cout << "DONE: " << " i = " << i << " l1 " << l1 << std::endl;
790  m_v[i][l1] = OSDBL_MAX;
791  m_px[i][l1] = -1; //a node we don't have
792 
793 
794  if(l1 == *(m_demand + i) ){//this should be valid even if demand is zero
795 
796  m_px[i][l1] = k;
797  // want the cost for arc (k, i)
798  // note: k < i so use that formula
799  m_u[i][l1] = *(c + k*(m_numNodes - 1) + i - 1); // put in correct cost
800 
801 
802 
803  }
804 
805  }
806  }
807  //end initialize
808 
809 
810 
811  //
812 
813  //if l = minDemand we visit only one node, let's get the reduced cost in this case
814  if(l == m_minDemand){
815 
816  for(i = m_numHubs; i < m_numNodes; i++){
817 
818 
819  if( m_u[i][l] + *(c + i*(m_numNodes-1) + k ) < rcost){
820 
821  rcost = m_u[i][l] + *(c + i*(m_numNodes-1) + k );
822 
823  //std::cout << " m_u[i][l2] = " << m_u[i][l2] << std::endl;
824 
825  //std::cout << " *(c + i*(m_numNodes-1) + k ) = " << *(c + i*(m_numNodes-1) + k ) << std::endl;
826  //push node back
827  bestLastNode = i;
828  }
829 
830  }
831 
832  //go from node bestLastNode to node k
833  //node bestLastNode is a higher number than k
834  *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
835  *(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1) + bestLastNode - 1;
836 
837 
838  return rcost;
839  }//end if on l == minDemand
840 
841 
842 // now calculate values for demand 2 or greater
843  //address of node (j, i) is j*(m_numNodes-1) + i when i < j
844  //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
845  // we start l2 at 2 since demand must be at least 1
846  // change to min demand + 1
847  //int lowerVal = m_minDemand + 1;
848  //lowerVal = m_minDemand;
849  //for(l2 = lowerVal; l2 <= l; l2++){// loop over possible demand values assuming we have already gone to at least one node
850  for(l2 = m_minDemand + 1; l2 <= l; l2++){// loop over possible demand values assuming we have already gone to at least one node
851 
852  for(i = m_numHubs; i < m_numNodes; i++) { //we are finding least cost to node i
853 
854 
855  if( *(m_demand + i) < l2 ){ // kipp < OR <= ????
856 
857  for(j = m_numHubs; j < i; j++){ //we are coming from node j
858 
859 
860 
861  //calculate g -- l2 - d[ i] is the demand trough and including node j
862  l1 = l2 - *(m_demand + i);
863 
864  if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
865 
866 
867  m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i - 1) ;
868 
869 
870 
871 
872  }else{
873 
874  m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i - 1) ;
875 
876 
877 
878  }
879 
880  // update u and the pointer for p
881 
882  if(m_g[j][i] < m_u[i][l2] ){
883 
884  m_u[i][l2] = m_g[j][i];
885  m_px[i][l2] = j;
886 
887  }
888 
889 
890 
891  }//end of first for loop on j, j < i
892 
893 
894  for(j = i + 1; j < m_numNodes; j++){ //we are coming from node j
895 
896 
897  //calculate g -- l2 - d[ i] is the demand trough and including node j
898  l1 = l2 - *(m_demand + i);
899 
900  if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
901 
902 
903  m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i ) ;
904 
905 
906  }else{
907 
908  m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i ) ;
909 
910  }
911 
912  // update u and the pointer for p
913 
914  if(m_g[j][i] < m_u[i][l2] ){
915 
916  m_u[i][l2] = m_g[j][i];
917  m_px[i][l2] = j;
918 
919  }
920 
921 
922  }//end of second for loop on j, j > i
923 
924 
925  //now calculate the second best solution and point
926  //right now m_px[i][l2] points to the best j node
927 
928  for(j =m_numHubs; j < m_numNodes; j++){ //we are coming from node j
929 
930  if(j != i){
931 
932  if( (m_g[j][i] < m_v[i][l2] ) && (m_px[i][l2] != j) ){ // kipp the && gives the second best
933 
934  //if( m_g[j][i] < m_v[i][l2] ) {
935 
936  m_v[i][l2] = m_g[j][i];
937  m_tx[i][l2] = j;
938 
939 
940  }
941 
942  }
943 
944 
945  }//end second best calculation, another for loop on j
946 
947  //now if l2 = l we are done
948  if(l2 == l ){
949 
950  if( m_u[i][l2] + *(c + i*(m_numNodes-1) + k ) < rcost){
951 
952  rcost = m_u[i][l2] + *(c + i*(m_numNodes-1) + k );
953 
954  bestLastNode = i;
955  }
956 
957  }
958 
959 
960  }//end if on demand less than l2
961 
962 
963  }//i loop
964 
965 
966  }//l2 loop
967 
968  //
969 
970  //std::cout << "best Last Node = " << bestLastNode << std::endl;
971 
972  // now get the path that gives the reduced cost
973 
974 
975  int currentNode;
976  int successorNode;
977  int lvalue;
978 
979  //initialize
980  // we work backwords from bestLastNode
981  //in our recursion we recurse on the currentNode and assume
982  //it is in the optimal path
983 
984  //node bestLastNode is a higher number than k
985  *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
986 
987 
988  //if we are here we should have a value for bestLastNode
989  //if not return infinity
990  if( bestLastNode == OSINT_MAX) return OSDBL_MAX;
991 
992  //by successor, I mean node after current node in the path
993  successorNode = k;
994  currentNode = bestLastNode;
995  //the lvalue is the demand through the currentNode
996  lvalue = l ;
997 
998 
999 
1000 
1001  while(currentNode != k){
1002  //std::cout << "currentNode = " << currentNode << " " << "lvalue " << lvalue << std::endl;
1003  if(*kountVar == m_upperBoundLMax + 1) return OSDBL_MAX;
1004  if( m_px[ currentNode][ lvalue ] != successorNode){
1005 
1006 
1007 
1008  //update nodes
1009  successorNode = currentNode;
1010  currentNode = m_px[ currentNode][ lvalue ];
1011 
1012 
1013  if(currentNode - successorNode > 0){
1014  //below diagonal
1015 
1016  //std::cout << "startPnt " << startPnt << " currentNode " << currentNode << " successorNode " << successorNode << " Kount " << *kountVar << std::endl;
1017 
1018  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
1019 
1020 
1021  }else{
1022  //above diagonal
1023 
1024  //std::cout << "startPnt " << startPnt << " currentNode " << currentNode << " successorNode " << successorNode << " Kount " << *kountVar << std::endl;
1025 
1026  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
1027 
1028  }
1029 
1030 
1031  }else{ //take second best
1032 
1033 
1034  //update nodes
1035  successorNode = currentNode;
1036  currentNode = m_tx[ currentNode][ lvalue ];
1037 
1038  if(currentNode - successorNode > 0){
1039  //below diagonal
1040 
1041  //std::cout << "startPnt " << startPnt << " currentNode " << currentNode << " successorNode " << successorNode << " Kount " << *kountVar << std::endl;
1042  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
1043 
1044  }else{
1045  //above diagonal
1046  //std::cout << "startPnt " << startPnt << " currentNode " << currentNode << " successorNode " << successorNode << " Kount " << *kountVar << std::endl;
1047  *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
1048 
1049  }
1050 
1051  }
1052 
1053  //update lvalue
1054  //if ( *(m_demand + successorNode) == 0) lvalue = lvalue - 1;
1055  //else lvalue = lvalue - *(m_demand + successorNode);
1056  lvalue = lvalue - *(m_demand + successorNode);
1057 
1058 
1059 
1060  }//end while
1061 
1062 
1063  //go from node k to bestLastNode -- already done in loop above
1064  //*(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1) + currentNode - 1;
1065 
1066 
1067  return rcost;
1068 
1069 }//end qroute
1070 
1071 
1072 
1073 
1074 void OSBearcatSolverXij::getColumns(const double* yA, const int numARows,
1075  const double* yB, const int numBRows,
1076  int &numNewColumns, int* &numNonzVec, double* &costVec,
1077  int** &rowIdxVec, double** &valuesVec, double &lowerBound)
1078 {
1079 
1080 //first strip of the phi dual values and then the convexity row costs
1081 
1082  int i;
1083  int j;
1084  int numCoulpingConstraints;
1085  numCoulpingConstraints = m_numNodes - m_numHubs;
1086 
1087  int numVar;
1088  numVar = m_numNodes*m_numNodes - m_numHubs;
1089  int numNonz;
1090 
1091  try{
1092 
1093 
1094 
1095  if(numARows != m_numNodes) throw ErrorClass("inconsistent row count in getColumns");
1096 
1097 
1098 
1099  //get the reduced costs
1100  calcReducedCost( yA, yB );
1101 
1102  int kountVar;
1103  double testVal;
1104  testVal = 0;
1105  int k;
1106  int startPntInc;
1107  int rowCount;
1108  double rowValue;
1109 
1111 
1112  double cpuTime;
1113  double start = CoinCpuTime();
1114  getOptL( m_rc);
1115  cpuTime = CoinCpuTime() - start;
1116  std::cout << "DYNAMIC PROGRAMMING CPU TIME " << cpuTime << std::endl;
1117  m_lowerBnd = 0.0;
1118  for(k = 0; k < m_numHubs; k++){
1119 
1120 
1121  //startPntInc = k*m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes) + (m_optL[ k] - 1)*(m_numNodes*m_numNodes - m_numNodes);
1122  startPntInc = (m_optL[ k] - 1)*(m_numNodes*m_numNodes - m_numNodes);
1123 
1124  std::cout << " whichBlock = " << k << " L = " << m_optL[ k] << std::endl;
1125 
1126  testVal += m_optL[ k];
1127 
1128  kountVar = 0;
1129 
1131  m_optValHub[ k] = qrouteCost(k, m_optL[ k], m_rc[ k], &kountVar);
1132 
1133  m_optValHub[ k] -= yA[ k + numCoulpingConstraints];
1134 
1135  std::cout << "Best Reduced Cost Hub " << k << " = " << m_optValHub[ k] << std::endl;
1136  m_lowerBnd += m_optValHub[ k];
1137 
1138  //loop over the rows, scatter each row and figure
1139  //out the column coefficients in this row
1140  //first scatter the sparse array m_varIdx[ j]
1141 
1142  m_costVec[ k] = 0.0;
1143 
1144  for(j = 0; j < kountVar; j++){
1145 
1146 
1147  //we are counting the NUMBER of times the variable used
1148  //the same variable can appear more than once in m_varIdx
1149  m_tmpScatterArray[ m_varIdx[ j] - startPntInc ] += 1;
1150 
1151 
1152  // is variable m_varIdx[ j] - startPntInc in this row
1153 
1154  m_costVec[ k] += m_cost[ m_varIdx[ j] - startPntInc ];
1155 
1156  }
1157 
1158 
1159 
1160  numNonz = 0;
1161  //multiply the sparse array by each A matrix constraint
1162  for(i = 0; i < numCoulpingConstraints; i++){
1163 
1164  rowCount = 0;
1165 
1166  for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
1167 
1168  //m_Amatrix[ j] is a variable index -- this logic works
1169  //since the Amatrix coefficient is 1 -- we don't need a value
1170  //it indexes variable that points into the node
1171  rowCount += m_tmpScatterArray[ m_Amatrix[ j] ];
1172 
1173 
1174  }
1175 
1176  if(rowCount > 0){
1177 
1178  //std::cout << "GAIL NODE " << i + m_numHubs << " COUNT = " << rowCount << std::endl;
1179  m_newColumnRowIdx[ k][ numNonz] = i;
1180  m_newColumnRowValue[ k][ numNonz] = rowCount;
1181  numNonz++;
1182  }
1183 
1184 
1185  }//end loop on coupling constraints
1186 
1187  //add a 1 in the convexity row
1188  m_newColumnRowIdx[ k][ numNonz] = m_numNodes - m_numHubs + k;
1189  m_newColumnRowValue[ k][ numNonz++] = 1.0;
1190 
1191 
1192 
1193  //now multiply the sparse array by each B-matrix constraint
1194 
1195  for(i = 0; i < m_numBmatrixCon; i++){
1196  //if the row corresponds to a multi-commodity row then m_BmatrixRowIndex[ i] = k
1197  if(m_BmatrixRowIndex[ i] == -1 || m_BmatrixRowIndex[ i] == k ){
1198 
1199  //rowCount = 0;
1200  rowValue = 0;
1201 
1202  for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
1203 
1204  //m_BmatrixIdx[ j] is a variable index -- this logic works
1205  //since the Bmatrix coefficient is 1 -- we don't need a value
1206  //it indexes variable that points into the node
1207  //rowCount += m_tmpScatterArray[ m_BmatrixIdx[ j] ];
1208  //now assume coefficients not necessarily 1
1209 
1210 
1211  rowValue += m_tmpScatterArray[ m_BmatrixIdx[ j] ]*m_BmatrixVal[ j];
1212 
1213 
1214  }
1215 
1216  if( rowValue > m_osDecompParam.zeroTol || rowValue < -m_osDecompParam.zeroTol){
1217 
1218 
1219  m_newColumnRowIdx[ k][ numNonz] = i + m_numNodes;
1220  m_newColumnRowValue[ k][ numNonz++] = rowValue;
1221 
1222  }
1223 
1224  }
1225 
1226 
1227  }
1228 
1229  m_newColumnNonz[ k] = numNonz;
1230 
1231 
1232  // std::cout << std::endl << std::endl;
1233  // std::cout << "HERE IS COLUMN " << m_numThetaVar << std::endl;
1234 
1235  //zero out the scatter vector and store the generated column
1236  for(j = 0; j < kountVar; j++){
1237 
1238  //std::cout << m_variableNames[ m_varIdx[ j] - startPntInc] << std::endl;
1239 
1240  m_thetaIndex[ m_numThetaNonz++ ] = m_varIdx[ j] - startPntInc ;
1241  m_tmpScatterArray[ m_varIdx[ j] - startPntInc ] = 0;
1242 
1243  // is variable m_varIdx[ j] - startPntInc in this row
1244 
1245  }
1246 
1247  //std::cout << std::endl << std::endl;
1248 
1249 
1250  intVarSet.insert ( std::pair<int,double>( m_numThetaVar, 1.0) );
1252  m_costVec[ k] = m_optL[ k]*m_costVec[ k];
1253  m_thetaCost[ m_numThetaVar++ ] = m_costVec[ k];
1255 
1256 
1257 
1258 
1259 
1260  //stuff for debugging
1261  //*****************************//
1294  //**************************//
1295  //end debugging stuff
1296 
1297 
1298  }//end of loop on hubs
1299 
1300 
1301 
1302  numNonzVec = m_newColumnNonz;
1303  costVec = m_costVec;
1304  rowIdxVec = m_newColumnRowIdx;
1305  valuesVec = m_newColumnRowValue;
1306  std::cout << "Lower Bound = " << m_lowerBnd << std::endl;
1307 
1308 
1309  if(testVal != m_totalDemand) {
1310 
1311  std::cout << "TOTAL DEMAND = " << m_totalDemand << std::endl;
1312  std::cout << "Test Value = " << testVal << std::endl;
1313  throw ErrorClass( "inconsistent demand calculation" );
1314  }
1315 
1316 
1317 
1318 
1319 
1320 
1321  } catch (const ErrorClass& eclass) {
1322 
1323  throw ErrorClass(eclass.errormsg);
1324 
1325  }
1326 
1327 
1328  //set method parameters
1329  numNewColumns = m_numHubs;
1330  lowerBound = m_lowerBnd;
1331 
1332  std::cout << "LEAVING GET COLUMNS" << std::endl;
1333  return;
1334 
1335 }//end getColumns
1336 
1337 
1338 
1633 
1635 
1636  std::cout << "Executing OSBearcatSolverXij::getInitialRestrictedMaster( )" << std::endl;
1637 
1638  //this master will have m_numNodes artificial variables
1639  int numVarArt;
1640  //numVarArt = 2*m_numNodes;
1641  numVarArt = m_numNodes;
1642 
1643  int numXVar;
1644  numXVar = m_numNodes*m_numNodes - m_numNodes;
1645 
1646  double* xVar = NULL;
1647  xVar = new double[ numXVar];
1648 
1649  int i;
1650  int k;
1651  std::string testFileName;
1652  std::string osil;
1653 
1654  std::map<int, std::map<int, std::vector<int> > >::iterator mit;
1655  std::map<int, std::vector<int> >::iterator mit2;
1656  std::vector<int>::iterator vit;
1657 
1658  m_osinstanceMaster = NULL;
1659  //add linear constraint coefficients
1660  //number of values will nodes.size() the coefficients in the node constraints
1661  //plus coefficients in convexity constraints which is the number of varaibles
1662  int kountNonz;
1663  int kount;
1664  m_numberOfSolutions = 1;
1665  int numThetaVar = m_numberOfSolutions*m_numHubs;
1666  //the extra is for the artificial variables
1667  double *values = new double[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar + numVarArt];
1668  int *indexes = new int[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar + numVarArt] ;
1669  int *starts = new int[ numThetaVar + 1 + numVarArt];
1670  kount = 0;
1671  starts[ 0] = 0;
1672  int startsIdx;
1673  startsIdx = 0;
1674  startsIdx++;
1675 
1676  double tspObjValue;
1677  double demandSum;
1678 
1679  SparseVector *objcoeff = NULL;
1680 
1681  for(i = 0; i < numVarArt; i++){
1683  m_thetaPnt[ m_numThetaVar++] = 0;
1684 
1685  }
1686  try {
1687 
1688  //start building the restricted master here
1690  m_osinstanceMaster->setInstanceDescription("The Initial Restricted Master");
1691 
1692  // first the variables
1693  // we have numVarArt] artificial variables
1695 
1696  // now add the objective function
1698  //add m_numNodes artificial variables
1699  objcoeff = new SparseVector( m_numberOfSolutions*m_numHubs + numVarArt);
1700 
1701  // now the constraints
1703 
1704  //add the artificial variables -- they must be first in the model
1705 
1706  int varNumber;
1707  varNumber = 0;
1708  std::string masterVarName;
1709  kountNonz = 0;
1710 
1711  for(i = 0; i < m_numNodes; i++){
1712 
1713  objcoeff->indexes[ varNumber ] = varNumber ;
1714  //if obj too large we get the following error
1715  //Assertion failed: (fabs(obj[i]) < 1.0e25), function createRim,
1716  //file ../../../Clp/src/ClpSimplex.cpp, l
1717  //objcoeff->values[ varNumber ] = OSDBL_MAX;
1718  //objcoeff->values[ varNumber ] = 1.0e24;
1719  objcoeff->values[ varNumber ] = m_osDecompParam.artVarCoeff;
1720 
1721  m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt("AP", i ) ,
1722  0, 1.0, 'C');
1723 
1724  values[ kountNonz] = 1;
1725  indexes[ kountNonz++] = i ;
1726  starts[ startsIdx++] = kountNonz;
1727 
1728  }
1729 
1730  m_bestIPValue = 0;
1731  mit = m_initSolMap.find( 0);
1732  if(mit == m_initSolMap.end() ) throw ErrorClass("Problem finding first solution in m_initSolMap");
1733 
1734  CoinSolver *solver = NULL;
1735  solver = getTSP(m_numNodes, m_cost);
1736 
1737 
1738  for(k = 0; k < m_numHubs; k++){
1739  //locate hub k
1740  mit2 = mit->second.find( k);
1741  if(mit2 == mit->second.end() ) throw ErrorClass("Problem finding hub k in the solution map");
1742 
1743  tspObjValue = getRouteDistance(m_numNodes, mit2->first, solver,
1744  mit2->second, xVar);
1745  demandSum = 0;
1746  //get demands
1747  for ( vit = mit2->second.begin() ; vit != mit2->second.end(); vit++ ) {
1748 
1749  demandSum += m_demand[ *vit];
1750  values[ kountNonz] = 1;
1751  indexes[ kountNonz++] = *vit - m_numHubs ;
1752 
1753  }
1754 
1755  //update theta indexes
1756 
1757  for(i = 0; i < numXVar; i++){
1758 
1759  //std::cout << "xVar = " << xVar[ i] << std::endl;
1760  if(xVar[ i] > .1) {
1761  m_thetaIndex[ m_numThetaNonz++ ] = i;
1762  //std::cout << m_variableNames[ i] << std::endl;
1763  }
1764  }
1765  //exit( 1);
1766  //now for convexity row k
1767 
1768  values[ kountNonz] = 1;
1769  indexes[ kountNonz++] = m_numNodes - m_numHubs + k ;
1770 
1771 
1772 
1773 
1775  m_thetaCost[ m_numThetaVar++ ] = demandSum*tspObjValue;
1777 
1778  masterVarName = makeStringFromInt("theta(", k);
1779  masterVarName += makeStringFromInt(",", 0);
1780  masterVarName += ")";
1781  intVarSet.insert ( std::pair<int,double>(varNumber, 1.0) );
1782  m_osinstanceMaster->addVariable(varNumber++, masterVarName, 0, 1, 'C');
1783 
1784  std::cout << "Optimal Objective Value = " << demandSum*tspObjValue << std::endl;
1785 
1786 
1787  objcoeff->indexes[ k + numVarArt ] = k + numVarArt ;
1788  objcoeff->values[ k + numVarArt ] = demandSum*tspObjValue;
1789 
1790  m_bestIPValue += demandSum*tspObjValue;
1791 
1792  std::cout << "m_bestIPValue = " << m_bestIPValue << std::endl;
1793  starts[ startsIdx++] = kountNonz;
1794 
1795  }//end of index on k
1796  std::cout << " m_numThetaVar " << m_numThetaVar << std::endl;
1797  //add the constraints
1798  //add the row saying we must visit each node
1799 
1800  for( i = 0; i < m_numNodes - m_numHubs ; i++){
1801 
1802  m_osinstanceMaster->addConstraint(i, makeStringFromInt("visitNode_", i + m_numHubs) , 1.0, 1.0, 0);
1803  }
1804 
1805  kount = 0;
1806 
1807  //add the convexity row
1808  for( i = m_numNodes - m_numHubs; i < m_numNodes ; i++){
1809 
1810  m_osinstanceMaster->addConstraint(i, makeStringFromInt("convexityRowRoute_", kount++ ) , 1.0, 1.0, 0);
1811  }
1812 
1813  m_osinstanceMaster->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
1814 
1815 
1816  //add the linear constraints coefficients
1818  values, 0, kountNonz - 1, indexes, 0, kountNonz - 1, starts, 0, startsIdx);
1819 
1820 
1821  std::cout << m_osinstanceMaster->printModel( ) << std::endl;
1822  std::cout << "NONZ = " << kountNonz << std::endl;
1823 
1824 
1825  delete objcoeff;
1826  objcoeff = NULL;
1827  delete[] xVar;
1828  xVar = NULL;
1829  delete solver->osinstance;
1830  delete solver;
1831 
1832 
1833 
1834 
1835  } catch (const ErrorClass& eclass) {
1836  std::cout << std::endl << std::endl << std::endl;
1837 
1838  if(objcoeff != NULL) delete objcoeff;
1839  if(xVar != NULL) delete xVar;
1840  // Problem with the parser
1841  throw ErrorClass(eclass.errormsg);
1842  }
1843 
1844  return m_osinstanceMaster;
1845 }//end generateInitialRestrictedMaster2
1846 
1847 
1848 
1850 
1851 
1852  std::cout << "Executing getOptions(OSOption *osoption)" << std::endl;
1853  //get any options relevant to OSColGenApp
1854  try{
1855 
1856  int i;
1857 
1858 
1859  std::vector<SolverOption*> solverOptions;
1860  std::vector<SolverOption*>::iterator vit;
1861  std::vector<int>::iterator vit2;
1862  std::vector<int >demand;
1863  std::vector<std::string >nodeName;
1864  std::vector<int >routeCapacity;
1865  std::vector<int >routeMinPickup;
1866  std::vector<double >arcCost;
1867  std::vector<double >::iterator vit3;
1868  std::vector<std::string>::iterator vit4;
1869 
1870 
1871  m_numberOfSolutions = 0;
1872  solverOptions = osoption->getSolverOptions("routeSolver");
1873  if (solverOptions.size() == 0) throw ErrorClass( "options for routeSolver not available");
1874  //iterate over the vector
1875 
1876  int tmpVal;
1877  double tmpCost;
1878 
1879 
1880  std::string routeString; //variable for parsing a category option
1881  std::string solutionString; //variable for parsing a category option
1882  std::string::size_type pos1; //variable for parsing a category option
1883  std::string::size_type pos2; //variable for parsing a category option
1884  std::string::size_type pos3; //variable for parsing a category option
1885 
1886 
1887  std::map<int, std::map<int, std::vector<int> > >::iterator mit;
1888  std::map<int, std::vector<int> >::iterator mit2;
1889  int solutionNumber;
1890  int routeNumber;
1891 
1892 
1893  for (vit = solverOptions.begin(); vit != solverOptions.end(); vit++) {
1894 
1895 
1896  //std::cout << (*vit)->name << std::endl;
1897 
1898  //(*vit)->name.compare("initialCol") == 0
1899  //if(rowNames[ i3].find("routeCapacity(1)") == string::npos )
1900 
1901  if( (*vit)->name.find("numHubs") != std::string::npos){
1902 
1903 
1904  std::istringstream hubBuffer( (*vit)->value);
1905  hubBuffer >> m_numHubs;
1906  std::cout << "numHubs = " << m_numHubs << std::endl;
1907 
1908  }else{
1909 
1910  if((*vit)->name.find("numNodes") != std::string::npos){
1911 
1912 
1913  std::istringstream numNodesBuffer( (*vit)->value);
1914  numNodesBuffer >> m_numNodes;
1915  std::cout << "numNodes = " << m_numNodes << std::endl;
1916 
1917  }else{
1918  if((*vit)->name.find("totalDemand") != std::string::npos){
1919 
1920 
1921  std::istringstream totalDemandBuffer( (*vit)->value);
1922  totalDemandBuffer >> m_totalDemand;
1923  std::cout << "m_totalDemand = " << m_totalDemand << std::endl;
1924 
1925  }else{
1926  if((*vit)->name.find("routeMinPickup") != std::string::npos){
1927 
1928 
1929  std::istringstream routeMinPickupBuffer( (*vit)->value);
1930  routeMinPickupBuffer >> tmpVal;
1931  routeMinPickup.push_back( tmpVal);
1932  //std::cout << "m_minDemand = " << tmpVal << std::endl;
1933 
1934  }else{
1935  if( (*vit)->name.find("demand") != std::string::npos ){
1936 
1937 
1938  std::istringstream demandBuffer( (*vit)->value);
1939  demandBuffer >> tmpVal;
1940  if(tmpVal <= 0 && demand.size() > (unsigned int) m_numHubs) throw ErrorClass("must have strictly positive demand");
1941  if(tmpVal < m_minDemand && demand.size() > (unsigned int) m_numHubs ) m_minDemand = tmpVal;
1942  demand.push_back( tmpVal);
1943  nodeName.push_back( (*vit)->description);
1944  //std::cout << "demand = " << tmpVal << std::endl;
1945 
1946  }else{
1947  if((*vit)->name.find("routeCapacity") != std::string::npos ){
1948  std::istringstream routeCapacityBuffer( (*vit)->value);
1949  routeCapacityBuffer >> tmpVal;
1950  routeCapacity.push_back( tmpVal);
1951  std::cout << "m_routeCapacity = " << tmpVal << std::endl;
1952 
1953  }else{
1954 
1955  if((*vit)->name.find("osilFile") != std::string::npos ){
1956  m_initOSiLFile = (*vit)->value;
1957  std::cout << "m_initOSiLFile = " << m_initOSiLFile << std::endl;
1958 
1959  }else{
1960 
1961  if( (*vit)->name.find("restrictedMasterSolution") != std::string::npos ){
1962  //std::istringstream buffer( (*vit)->value);
1963  //buffer >> m_numberOfSolutions;
1964 
1965  //get the block number and solution number
1966  //first get routeString and soluionString
1967  //parse the category string base on :
1968  pos1 = (*vit)->category.find( ":");
1969  if(pos1 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
1970 
1971  //solutionString = (*vit)->category.substr( pos1 + 1, pos2 - pos1 - 1);
1972  solutionString = (*vit)->category.substr( 0, pos1);
1973  routeString = (*vit)->category.substr( pos1 + 1);
1974 
1975  pos2 = solutionString.find( "n");
1976  if(pos2 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
1977 
1978  std::istringstream solutionBuffer( solutionString.substr( pos2 + 1) );
1979  solutionBuffer >> solutionNumber;
1980  //std::cout << "solution number = " << solutionNumber << std::endl;
1981 
1982 
1983  pos3 = routeString.find( "e");
1984  if(pos3 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
1985  std::istringstream routeBuffer( routeString.substr( pos3 + 1) );
1986  routeBuffer >> routeNumber;
1987  //std::cout << "route number = " << routeNumber << std::endl;
1988  std::istringstream nodeBuffer( (*vit)->value);
1989  nodeBuffer >> tmpVal;
1990 
1991  mit = m_initSolMap.find( solutionNumber );
1992 
1993  if( mit != m_initSolMap.end() ){
1994  // we have solution from before
1995  // do we have a new route?
1996 
1997  mit2 = mit->second.find( routeNumber);
1998 
1999  if(mit2 != mit->second.end() ){
2000  // we have a route from before and solution from before
2001 
2002 
2003  mit2->second.push_back( tmpVal);
2004 
2005 
2006  }else{
2007 
2008  //we have a new route, but old solution
2009  std::vector<int> tmpVec;
2010  tmpVec.push_back( tmpVal) ;
2011  mit->second.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
2012 
2013 
2014  }
2015 
2016  }else{
2017  // we have a new solution
2018  std::vector<int> tmpVec;
2019  tmpVec.push_back( tmpVal) ;
2020 
2021  std::map<int, std::vector<int> > tmpMap;
2022  tmpMap.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
2023  m_initSolMap.insert( std::pair<int, std::map<int, std::vector<int> > >(solutionNumber, tmpMap) ) ;
2024 
2025  }
2026  }//if on restricted master solution
2027  else{
2028  if( (*vit)->name.find("maxMasterColumns") != std::string::npos){
2029 
2030 
2031  std::istringstream maxMasterColumns( (*vit)->value);
2032  maxMasterColumns >> m_maxMasterColumns;
2033  std::cout << "m_maxMasterColumn = " << m_maxMasterColumns << std::endl;
2034 
2035  }else{
2036  if( (*vit)->name.find("maxThetaNonz") != std::string::npos){
2037 
2038  std::istringstream maxThetaNonz( (*vit)->value);
2039  maxThetaNonz >> m_maxThetaNonz;
2040  std::cout << "m_maxThetaNonz = " << m_maxThetaNonz << std::endl;
2041 
2042  }else{
2043  if( (*vit)->name.find("use1OPTstart") != std::string::npos){
2044  m_use1OPTstart = false;
2045  if ( (*vit)->value.find("true") != std::string::npos ) m_use1OPTstart = true;
2046  std::cout << "m_use1OPTstart = " << m_use1OPTstart << std::endl;
2047 
2048  }else{
2049  if( (*vit)->name.find("maxBmatrixCon") != std::string::npos ){
2050 
2051  std::istringstream maxBmatrixCon( (*vit)->value);
2052  maxBmatrixCon >> m_maxBmatrixCon;
2053  std::cout << "m_maxBmatrixCon = " << m_maxBmatrixCon << std::endl;
2054 
2055  }else{
2056  if( (*vit)->name.find("maxBmatrixNonz") != std::string::npos ){
2057 
2058  std::istringstream maxBmatrixNonz( (*vit)->value);
2059  maxBmatrixNonz >> m_maxBmatrixNonz;
2060  std::cout << "m_maxBmatrixNonz = " << m_maxBmatrixNonz << std::endl;
2061 
2062 
2063  }else{
2064 
2065  if( (*vit)->name.find("cost") != std::string::npos ){
2066 
2067 
2068  std::istringstream costBuffer( (*vit)->value);
2069  costBuffer >> tmpCost;
2070  arcCost.push_back( tmpCost);
2071 
2072  }
2073  }
2074  }
2075  }
2076  }
2077  }
2078  }
2079  }
2080  }
2081  }
2082  }
2083  }
2084  }
2085  }//end if on solver options
2086 
2087  }//end for loop on options
2088 
2089 
2090  //now fill in route capacities
2091  i = 0;
2092  m_routeCapacity = new int[ m_numHubs];
2093  if( (unsigned int)m_numHubs != routeCapacity.size( ) ) throw ErrorClass("inconsistent number of HUBS");
2094  for (vit2 = routeCapacity.begin(); vit2 != routeCapacity.end(); vit2++) {
2095 
2096  *(m_routeCapacity + i++) = *vit2;
2097 
2098  std::cout << "ROUTE CAP = " << *vit2 << std::endl;
2099 
2100  }
2101  routeCapacity.clear();
2102 
2103 
2104  //now fill in route min pickups
2105  i = 0;
2106  m_routeMinPickup = new int[ m_numHubs];
2107  if( (unsigned int)m_numHubs != routeMinPickup.size( ) ) throw ErrorClass("inconsistent number of HUBS");
2108  //initialize
2109  for(int k = 0; k < m_numHubs; k++){
2110  m_routeMinPickup[ k] = 1;
2111  }
2112  for (vit2 = routeMinPickup.begin(); vit2 != routeMinPickup.end(); vit2++) {
2113 
2114  *(m_routeMinPickup + i++) = *vit2;
2115 
2116  }
2117  routeMinPickup.clear();
2118 
2119 
2120  //now fill in demand
2121  i = 0;
2122  m_demand = new int[ m_numNodes];
2123  if( (unsigned int)m_numNodes != demand.size( ) )
2124  throw ErrorClass("inconsistent number of demand nodes");
2125  for (vit2 = demand.begin(); vit2 != demand.end(); vit2++) {
2126 
2127  *(m_demand + i++) = *vit2;
2128 
2129  }
2130  demand.clear();
2131 
2132  //now fill in node names
2133  i = 0;
2134  m_nodeName = new std::string[ m_numNodes];
2135 
2136  for (vit4 = nodeName.begin(); vit4 != nodeName.end(); vit4++) {
2137 
2138  *(m_nodeName + i++) = *vit4;
2139 
2140  }
2141  nodeName.clear();
2142 
2143 
2144 
2145  //now fill in costs
2146  m_cost = NULL;
2147  m_costSetInOption = false;
2148  if(arcCost.size() != (unsigned int)(m_numNodes*m_numNodes - m_numNodes) )
2149  throw ErrorClass("input cost vector not consistent with number of nodes");
2150  if(arcCost.size() >= 1){
2151  m_costSetInOption = true;
2152  i = 0;
2153  m_cost = new double[ m_numNodes*m_numNodes - m_numNodes ];
2154  for (vit3 = arcCost.begin(); vit3 != arcCost.end(); vit3++) {
2155 
2156  *(m_cost + i++) = *vit3;
2157 
2158  }
2159  arcCost.clear();
2160  }
2161 
2162  //kipp -- fill in numberOfRestricedMasterSolutions from map size
2164 
2165 
2166  } catch (const ErrorClass& eclass) {
2167 
2168  throw ErrorClass(eclass.errormsg);
2169 
2170  }
2171 
2172 }//end getOptions
2173 
2174 
2175 
2176 void OSBearcatSolverXij::getCutsTheta(const double* theta, const int numTheta,
2177  int &numNewRows, int* &numNonz, int** &colIdx,
2178  double** &values, double* &rowLB, double* &rowUB) {
2179  //critical -- the variables that come in the theta variables
2180  //not the x variables, we must convert to x, find a cut in x-space
2181  //and then convert back to theta
2182 
2183  int i;
2184  int j;
2185  int k;
2186  int index;
2187  int rowKount;
2188  int tmpKount;
2189  int indexAdjust = m_numNodes - m_numHubs;
2190  double* tmpRhs;
2191  int numSepRows = m_osinstanceSeparation->getConstraintNumber() ;
2192 
2193  tmpRhs = new double[ numSepRows ];
2194 
2195  for(i = 0; i < numSepRows; i++){
2196 
2197  tmpRhs[ i] = 0;
2198  }
2199 
2200  try{
2202  //m_numNodes is the number of artificial variables
2203  if(numTheta != m_numThetaVar ) throw
2204  ErrorClass("number of master varibles in OSBearcatSolverXij::getCuts inconsistent");
2205 
2206  //for(i = 0; i < numTheta; i++){
2207 
2208  //std::cout << "numTheta = " << numTheta << std::endl;
2209  //std::cout << "m_numThetaVar = " << m_numThetaVar - 1 << std::endl;
2210 
2211  //exit( 1);
2212 
2213  for(i = 0; i < numTheta; i++){
2214 
2215  //get a postive theta
2216  if(theta[ i] > m_osDecompParam.zeroTol){
2217 
2218  //get the xij indexes associated with this variable
2219  for(j = m_thetaPnt[ i ]; j < m_thetaPnt[ i + 1 ]; j++ ){
2220 
2221  //get the xij index
2222 
2223 
2224 
2225  rowKount = m_separationIndexMap[ m_thetaIndex[ j] ];
2226 
2227  //std::cout << "rowKount = " << rowKount <<std::endl;
2228 
2229  if(rowKount < OSINT_MAX ){
2230 
2231  tmpRhs[ rowKount] -= theta[ i];
2232 
2233  }
2234 
2235  }
2236  }
2237  }
2238 
2239 
2240  // don't adjust the kludge row
2241 
2242  for(i = indexAdjust; i < numSepRows - 1; i++){
2243 
2244  if(-tmpRhs[ i] > 1 + m_osDecompParam.zeroTol ){
2245  // quick and dirty does x_{ij} + x_{ji} <= 1 cut off anything
2246  //std::cout << " tmpRhs[ i] = " << tmpRhs[ i] << std::endl;
2247  //which variable is this
2248  //kipp this an inefficient way of finding i and j -- improve this
2249  int tmpKount = indexAdjust;
2250  for(int i1 = m_numHubs; i1 < m_numNodes; i1++){
2251 
2252 
2253 
2254  for(int j1 = i1+1; j1 < m_numNodes; j1++){
2255 
2256  if(tmpKount == i){
2257 
2258  //std::cout << "i = " << i1 << std::endl;
2259  //std::cout << "j = " << j1 << std::endl;
2260  //okay generate a cut that says
2261  // x(i1,j1) + x(j1, i1) << 1
2263  //get index for i1,j1
2264  m_BmatrixIdx[ m_numBmatrixNonz++ ] = i1*(m_numNodes - 1) + j1 - 1 ;
2265 
2267  //get index for j1,i1
2268  m_BmatrixIdx[ m_numBmatrixNonz++ ] = j1*(m_numNodes - 1) + i1 ;
2269  m_numBmatrixCon++;
2271 
2272  numNewRows = 1;
2273 
2274  m_newRowNonz[ 0] = 0;
2275  m_newRowUB[ 0] = 1;
2276  m_newRowLB[ 0] = 0;
2277 
2278  //now we have to get the theta column indexes
2279  //scatter the constraint in the x - variables
2280 
2281  for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2282  j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2283 
2284 
2285  std::cout << " m_BmatrixIdx[ j] " << m_BmatrixIdx[ j] << std::endl ;
2286 
2287  m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 1;
2288 
2289  }
2290 
2291 
2292 
2293 
2294  for(k = 0; k < m_numThetaVar ; k++){
2295 
2296  //get the xij indexes in this column
2297  tmpKount = 0;
2298 
2299 
2300  for(j = m_thetaPnt[k]; j < m_thetaPnt[k + 1] ; j++){
2301 
2302  if(m_tmpScatterArray[ m_thetaIndex[ j] ] > 0 ){ //count number of xij for theta_i
2303 
2304  std::cout << " Variable " << m_variableNames[ m_thetaIndex[ j] ] << std::endl;
2305 
2306  tmpKount++;
2307 
2308  }
2309 
2310  }//end loop on j
2311 
2312  if(tmpKount > 0){
2313  //theta_i has a nonzero coefficient in this row
2314 
2315  m_newRowColumnIdx[0][ m_newRowNonz[ 0] ] = k ;
2316  //m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
2317  m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
2318 
2319 
2320  }
2321 
2322  }//end loop on k
2323 
2324 
2325  //zero out the scatter array again
2326 
2327  //::cout << " m_numBmatrixCon - 1 " << m_numBmatrixCon - 1 << std::endl;
2328  //std::cout << " m_pntBmatrix[ m_numBmatrixCon - 1] " << m_pntBmatrix[ m_numBmatrixCon - 1] << std::endl ;
2329  //std::cout << " m_pntBmatrix[ m_numBmatrixCon ] " << m_pntBmatrix[ m_numBmatrixCon ] << std::endl ;
2330  for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2331  j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2332 
2333  m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 0;
2334 
2335  }
2336 
2337  numNonz = m_newRowNonz;
2338  colIdx = m_newRowColumnIdx;
2339  values = m_newRowColumnValue;
2340  rowUB = m_newRowUB;
2341  rowLB = m_newRowLB;
2342 
2343  delete[] tmpRhs;
2344  tmpRhs = NULL;
2345 
2346  //we found a row, add the corresponding artificial variables
2347  //to the transformation matrix
2348  m_numThetaVar++;
2351  //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz;
2352  //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz;
2353 
2354  return;
2355 
2356  } //end loop on if tmpKount
2357 
2358  tmpKount++;
2359 
2360  }//loop on j1
2361 
2362  }//loop on i1
2363 
2364 
2365  }//end if on tmpRHS
2366 
2367  m_separationClpModel->setRowUpper(i, tmpRhs[ i] );
2368  m_separationClpModel->setRowLower(i, tmpRhs[ i] );
2369 
2370  }//end loop on i
2371 
2372 
2373  //std::cout << m_osinstanceSeparation->printModel() << std::endl;
2374 
2375 
2376  std::vector<int> dualIdx;
2377  std::vector<int>::iterator vit1;
2378  std::vector<int>::iterator vit2;
2379 
2380  //if the objective function value is greater than zero we have a cut
2381  //the cut is based on the nodes with dual value - 1
2382 
2383  for(k = 0; k < indexAdjust; k++){
2384  //std::cout << std::endl << std::endl;
2385 
2386 
2387  m_separationClpModel->setRowUpper(k, 0.0);
2388  //we don't need output
2389 
2390  m_separationClpModel->setLogLevel( 0);
2391 
2392  m_separationClpModel->primal();
2393 
2394  if(m_separationClpModel->getObjValue() > m_osDecompParam.zeroTol){
2395  std::cout << "DOING SEPARATION FOR NODE " << k + m_numHubs << std::endl;
2396  std::cout << "SEPERATION OBJ VALUE = " << m_separationClpModel->getObjValue() << std::endl;
2397  numNewRows = 1;
2398 
2399  for(i = 0; i < m_numNodes - m_numHubs ; i++){
2400  //std::cout << m_osinstanceSeparation->getConstraintNames()[ i] << " = " << m_separationClpModel->getRowPrice()[ i] << std::endl;
2401  if( m_separationClpModel->getRowPrice()[ i] - m_osDecompParam.zeroTol <= -1) dualIdx.push_back( i) ;
2402  }
2403 
2404  for (vit1 = dualIdx.begin(); vit1 != dualIdx.end(); vit1++) {
2405 
2406  i = *vit1 + m_numHubs;
2407 
2408  for (vit2 = dualIdx.begin(); vit2 != dualIdx.end(); vit2++) {
2409 
2410  j = *vit2 + m_numHubs;
2411 
2412  if( i > j ){
2413 
2414  index = i*(m_numNodes -1) + j;
2415  std::cout << "CUT VARIABLE = " << m_variableNames[ index ] <<std::endl;
2417  m_BmatrixIdx[ m_numBmatrixNonz++ ] = index ;
2418 
2419  }else{
2420 
2421  if( i < j ){
2422 
2423  index = i*(m_numNodes -1) + j - 1;
2424  std::cout << "CUT VARIABLE = " << m_variableNames[ index ] <<std::endl;
2426  m_BmatrixIdx[ m_numBmatrixNonz++ ] = index ;
2427 
2428  }
2429  }
2430 
2431  }//end for on vit2
2432  }//end for on vit1
2433 
2434  //add the tour-breaking cut
2435  m_numBmatrixCon++;
2437 
2438  // multiply the transformation matrix times this cut to get the cut in theta space
2439  // do the usual trick and scatter m_BmatrixIdx into a dense vector
2440 
2441  //reset
2442  // don't adjust the kludge row
2443  for(i = indexAdjust; i < numSepRows - 1; i++){
2444 
2445  m_separationClpModel->setRowUpper(i, 0.0 );
2446  m_separationClpModel->setRowLower(i, 0.0 );
2447 
2448 
2449  }
2450  m_separationClpModel->setRowUpper(k, 1.0);
2451  delete[] tmpRhs;
2452  tmpRhs = NULL;
2453 
2454 
2455  m_newRowNonz[ 0] = 0;
2456  m_newRowUB[ 0] = dualIdx.size() - 1;
2457  m_newRowLB[ 0] = 0;
2458 
2459  dualIdx.clear();
2460 
2461  //now we have to get the theata column indexes
2462  //scatter the constraint in the x - variables
2463 
2464  for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2465  j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2466 
2467  m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 1;
2468 
2469  }
2470 
2471 
2472 
2473 
2474  for(i = 0; i < m_numThetaVar ; i++){
2475 
2476  //get the xij indexes in this column
2477  tmpKount = 0;
2478  for(j = m_thetaPnt[i]; j < m_thetaPnt[i + 1] ; j++){
2479 
2480  if(m_tmpScatterArray[ m_thetaIndex[ j] ] > 0 ){ //count number of xij for theta_i
2481 
2482  tmpKount++;
2483 
2484  }
2485 
2486  }//end loop on j
2487 
2488  if(tmpKount > 0){
2489  //theta_i has a nonzero coefficient in this row
2490 
2491  m_newRowColumnIdx[0][ m_newRowNonz[ 0] ] = i ;
2492 
2493  m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
2494 
2495 
2496  }
2497 
2498  }//end loop on i
2499 
2500 
2501  //zero out the scatter array again
2502 
2503  for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2504  j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2505 
2506  m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 0;
2507 
2508  }
2509 
2510 
2511 
2512  numNonz = m_newRowNonz;
2513  colIdx = m_newRowColumnIdx;
2514  values = m_newRowColumnValue;
2515  rowUB = m_newRowUB;
2516  rowLB = m_newRowLB;
2517  m_numThetaVar++;
2520  //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz; //first artificial variable
2521  //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz; //second artificial varaible
2522 
2523  return;
2524 
2525 
2526 
2527  }//end if on obj value
2528  m_separationClpModel->setRowUpper(k, 1.0);
2529  dualIdx.clear();
2530 
2531  }//end loop on k
2532 
2533  //if we are here there was no cut
2534  //reset
2535  // don't adjust the kludge row
2536  for(i = indexAdjust; i < numSepRows - 1; i++){
2537 
2538  m_separationClpModel->setRowUpper(i, 0.0 );
2539  m_separationClpModel->setRowLower(i, 0.0 );
2540 
2541 
2542  }
2543 
2544  delete[] tmpRhs;
2545  tmpRhs = NULL;
2546 
2547  } catch (const ErrorClass& eclass) {
2548 
2549  throw ErrorClass(eclass.errormsg);
2550 
2551  }
2552 
2553 
2554 
2555 }//end getCutsTheta
2556 
2557 
2558 void OSBearcatSolverXij::getCutsMultiCommod(const double* theta, const int numTheta,
2559  int &numNewRows, int* &numNonz, int** &colIdx,
2560  double** &values, double* &rowLB, double* &rowUB) {
2561  //critical -- the variables that come in the theta variables
2562  //not the x variables, we must convert to x, find a cut in x-space
2563  //and then convert back to theta
2564 
2565  numNewRows = 0;
2566 
2567 
2568  //m_convexityRowIndex
2569  int i;
2570  int j;
2571  int k;
2572  int ivalue;
2573  int jvalue;
2574  int thetaIdx;
2575  int inodenum;
2576  int jnodenum;
2577  int j1;
2578  int j2;
2579  int kount;
2580  double wVal;
2581  double uVal;
2582  bool foundCut;
2583  double objVal;
2584 
2585  int ckijIndex;
2586  int numNonHubs;
2587  numNonHubs = m_numNodes - m_numHubs;
2588 
2589  int numVar;
2590  double rowValue;
2591 
2592  double* scatterValues;
2593  int numXijVar = m_numNodes*m_numNodes - m_numNodes;
2594  scatterValues = new double[ numXijVar ];
2595  for(i = 0; i < numXijVar; i++ )scatterValues[ i] = 0;
2596 
2597  double tmpCoeff;
2598 
2599  double *wcoeff = NULL;
2600  wcoeff = new double[ numNonHubs];
2601  CoinSolver *solver;
2602 
2603  std::cout << std::endl << std::endl;
2604  std::cout << "INSIDE getCutsMultiCommod " << std::endl;
2605  std::cout << std::endl << std::endl;
2606 
2607  std::map<int, std::vector<std::pair<int,double> > > xVarMap;
2608  std::map<int, std::vector<std::pair<int,double> > >::iterator mit;
2609  std::vector<std::pair<int,double> >::iterator vit;
2610 
2611  std::map<std::pair<int, int>,int>xVarIndexMap;
2612  std::pair<int,int> indexPair;
2613  ostringstream cutString;
2614 
2615 
2616 
2617 
2618 
2619 
2620  //construct index map
2648  //intVarSet.insert ( std::pair<int,double>(mit->second, 1.0) );
2649 
2650  try{
2651 
2652  for(i = 0; i < numTheta; i++){
2653  xVarMap[ m_convexityRowIndex[ i] ] ;
2654  //get a postive theta
2655  if(theta[ i] > m_osDecompParam.zeroTol){
2656 
2657  //get the xij indexes associated with this variable
2658  for(j = m_thetaPnt[ i ]; j < m_thetaPnt[ i + 1 ]; j++ ){
2659 
2660  mit = xVarMap.find( m_convexityRowIndex[ i]) ;
2661 
2662  if(mit != xVarMap.end() ){
2663 
2664  mit->second.push_back( std::pair<int, double>( m_thetaIndex[ j], theta[ i]) );
2665 
2666  }
2667  }
2668  }
2669  }
2670 
2671 
2672 
2673  //get a cut for each hub
2674  for(k = 0; k < m_numHubs; k++){
2675 
2676 
2677  mit = xVarMap.find( k) ;
2678 
2679  solver = m_multCommodCutSolvers[ k];
2680 
2681  numVar = solver->osiSolver->getNumCols();
2682  for(i = 0; i < numVar; i++ ) solver->osiSolver->setObjCoeff( i, 0);
2683 
2684  for(i = 0; i < numNonHubs; i++) wcoeff[ i ] = 0;
2685 
2686  if(mit != xVarMap.end() ){
2687 
2688  std::cout << "CONVEXITY ROW " << " = " << mit->first << std::endl;
2689 
2690  for(vit = mit->second.begin(); vit < mit->second.end(); vit++){
2691 
2692  //std::cout << m_variableNames[ vit->first ] << " = " << vit->second << std::endl;
2693 
2694  ivalue = vit->first /(m_numNodes - 1);
2695 
2696  jvalue = vit->first - ivalue*(m_numNodes - 1);
2697 
2698  if( jvalue >= ivalue ){
2699  //std::cout << " i NODE NUMBER = " << ivalue << std::endl;
2700  //std::cout << " j NODE NUMBER = " << jvalue + 1 << std::endl;
2701  inodenum = ivalue;
2702  jnodenum = jvalue + 1;
2703 
2704 
2705  }else{
2706  //std::cout << " i NODE NUMBER = " << ivalue << std::endl;
2707  //std::cout << " j NODE NUMBER = " << jvalue << std::endl;
2708 
2709  inodenum = ivalue;
2710  jnodenum = jvalue;
2711  }
2712 
2713  if(jnodenum != k){
2714 
2715  //std::cout << "GAIL " << m_nodeName[ inodenum] << " " << m_nodeName[ jnodenum] << std::endl;
2716  wcoeff[ jnodenum - m_numHubs ] += vit->second;
2717 
2718  if( inodenum == k ){
2719  ckijIndex = jnodenum - m_numHubs;
2720  } else{
2721  //inodenum and jnodenum NOT hub nodes
2722  inodenum -= m_numHubs;
2723  jnodenum -= m_numHubs;
2724 
2725 
2726  if( jnodenum > inodenum) ckijIndex = inodenum*(numNonHubs - 1) + jnodenum - 1 ;
2727  else ckijIndex = inodenum*(numNonHubs - 1) + jnodenum ;
2728 
2729  //increment by the number (k, i) variables --- there are numNonHUbs
2730  ckijIndex += numNonHubs ;
2731 
2732  }
2733 
2734 
2735  ckijIndex += numNonHubs;
2736 
2737  tmpCoeff = solver->osiSolver->getObjCoefficients()[ ckijIndex] ;
2738  //std::cout << " HONDA " << "cijkIndex = " << ckijIndex << std::endl;
2739 
2740  tmpCoeff = tmpCoeff - vit->second;
2741  if( ckijIndex > numVar - 1) throw ErrorClass( "problem with ckijIndex calculation");
2742  //std::cout << "cijkIndex = " << ckijIndex << " tmp Coeff After = " << tmpCoeff << " vit->second = "<< vit->second << std::endl;
2743  solver->osiSolver->setObjCoeff( ckijIndex, tmpCoeff );
2744  }
2745 
2746  }//end looping over xkij
2747 
2748  }//end if on mit
2749  foundCut = false;
2750  //loop over s to get a cut
2751  for(i = 0; i < numNonHubs; i++){
2752 
2753  //set s^{\hat} coefficient
2754  solver->osiSolver->setObjCoeff( i, wcoeff[ i ] );
2755 
2756  //solve the LP
2757  solver->solve();
2758  //solver->osiSolver->initialSolve();
2759  //kipp -- change the hard coding
2760  //we have a cut
2761  if(solver->osiSolver->getObjValue() > .1){
2762  objVal = 0;
2763  std::cout << "Separation Obj Value = " <<
2764  solver->osiSolver->getObjValue() << " Node " << m_nodeName[i + m_numHubs] << std::endl;
2765  //if(k == 1)solver->osiSolver->writeLp( "tmpLP--April17");
2766  foundCut = true;
2767  m_newRowNonz[ numNewRows ] = 0;
2768  //get the solution, let's see what the cut is in x(i, j) variables
2769  //empty the buffer
2770  cutString.str("");
2771  cutString << "";
2772  //first get the coefficients on x(i, shat)
2773  // this is the sum of w(shat) and the u(i, shat)
2774  kount = numNonHubs;
2775 
2776  wVal = solver->osiSolver->getColSolution()[ i];
2777  objVal += wVal*solver->osiSolver->getObjCoefficients()[ i];
2778 
2779  if(wVal < m_osDecompParam.zeroTol )throw ErrorClass("problem with wVal in generating a multicommodity cut");
2780  //get the u(k,j2) variables
2781  //j1 = k
2782  for(j2 = m_numHubs; j2 < m_numNodes; j2++){ //j1 = k
2783 
2784  indexPair.first = k;
2785  indexPair.second = j2;
2786 
2787  uVal = solver->osiSolver->getColSolution()[ kount];
2788  objVal += uVal*solver->osiSolver->getObjCoefficients()[ kount];
2789 
2790  //if(solver->osiSolver->getObjCoefficients()[ kount] < -0.001 ) std::cout << m_variableNames[ m_xVarIndexMap[ indexPair] ]<< " valueeeee = " << uVal << " coeff " << solver->osiSolver->getObjCoefficients()[ kount] << " kount " << kount << std::endl;
2791 
2792 
2793  if (j2 == (i + m_numHubs) ){
2794 
2795  if( (wVal - uVal) > m_osDecompParam.zeroTol || (uVal - wVal) > m_osDecompParam.zeroTol ){
2796  //variable (wVal - uVal)*x(k, shat)
2797  cutString << " =";
2798  cutString << (wVal - uVal);
2799  cutString << "*" ;
2800  //cutString << m_nodeName[ k] ;
2801  //cutString << "," ;
2802  //cutString << m_nodeName[ i + m_numHubs] ; //i is indexing a non-hub
2803  //cutString << ") " ;
2804 
2805  m_BmatrixVal[ m_numBmatrixNonz ] = (wVal - uVal)/wVal ;
2806  //get index for k, j2
2807  //m_BmatrixIdx[ m_numBmatrixNonz++ ] = k*(m_numNodes - 1) + j2 - 1 ;
2808  if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
2809  std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
2810  throw ErrorClass( "index mapping problem in generating multicommodity cut");
2811  }else{
2812  m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
2813  cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
2814  }
2815 
2816  }
2817 
2818  }else{
2819 
2820  if( (-uVal > m_osDecompParam.zeroTol) || (uVal > m_osDecompParam.zeroTol) ){
2821  //variable -uVal*x(j1, j2)
2822  cutString << " ";
2823  cutString << - uVal;
2824  cutString << "*" ;
2825  //cutString << m_nodeName[ k ];
2826  //cutString << "," ;
2827  //cutString << m_nodeName[ j2 ];
2828  //cutString << ") " ;
2829 
2830 
2831  m_BmatrixVal[ m_numBmatrixNonz ] = -uVal/wVal ;
2832  //get index for k, j2
2833  //m_BmatrixIdx[ m_numBmatrixNonz++ ] = k*(m_numNodes - 1) + j2 - 1 ;
2834 
2835  if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
2836  std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
2837  throw ErrorClass( "index mapping problem in generating multicommodity cut");
2838  }else{
2839  m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
2840  cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
2841  }
2842  }
2843 
2844  }
2845 
2846  kount++;
2847  }
2848 
2849  //get the u(j1, j2) variables for j1 not a hub
2850  for(j1 = m_numHubs; j1 < m_numNodes; j1++){ //<= since we include hub
2851 
2852  //j1 = 0 corresponds to the hub
2853 
2854  //for(j2 = j1 + 1; j2 < m_numNodes; j2++){
2855  for(j2 = m_numHubs; j2 < j1; j2++){
2856 
2857  indexPair.first = j1;
2858  indexPair.second = j2;
2859 
2860  uVal = solver->osiSolver->getColSolution()[ kount];
2861  objVal += uVal*solver->osiSolver->getObjCoefficients()[ kount];
2862 
2863  //if(solver->osiSolver->getObjCoefficients()[ kount] < -0.001 ) std::cout << m_variableNames[ m_xVarIndexMap[ indexPair] ] << " valueeeee = " << uVal << " coeff " << solver->osiSolver->getObjCoefficients()[ kount] << " kount " << kount << std::endl;
2864 
2865  if (j2 == (i + m_numHubs) ){
2866 
2867 
2868  if( (wVal - uVal) > m_osDecompParam.zeroTol || (uVal - wVal) > m_osDecompParam.zeroTol ){
2869  //variable (wVal - uVal)*x(j1, shat)
2870  cutString << " +";
2871  cutString << (wVal - uVal) ;
2872  cutString << "*" ;
2873  //cutString << m_nodeName[ j1] ;
2874  //cutString << "," ;
2875  //cutString << m_nodeName[ i + m_numHubs];
2876  //cutString << ") " ;
2877  m_BmatrixVal[ m_numBmatrixNonz ] = (wVal - uVal)/wVal ;
2878  //get index for j1, j2 with j1 < j2
2879  //m_BmatrixIdx[ m_numBmatrixNonz++ ] = j1*(m_numNodes - 1) + j2 - 1 ;
2880 
2881  if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
2882  std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
2883  throw ErrorClass( "index mapping problem in generating multicommodity cut");
2884  }else{
2885  m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
2886  cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
2887  }
2888 
2889  }
2890 
2891  }else{
2892 
2893  if( (-uVal > m_osDecompParam.zeroTol) || (uVal > m_osDecompParam.zeroTol) ){
2894  //variable -uVal*x(j1, j2)
2895  cutString << " ";
2896  cutString << - uVal;
2897  cutString << "*" ;
2898  //cutString << m_nodeName[ j1 ];
2899  //cutString << "," ;
2900  //cutString << m_nodeName[ j2 ];
2901  //cutString << ") " ;
2902 
2903  m_BmatrixVal[ m_numBmatrixNonz ] = -uVal/wVal ;
2904  //get index for j1, j2 with j1 < j2
2905  //m_BmatrixIdx[ m_numBmatrixNonz++ ] = j1*(m_numNodes - 1) + j2 - 1 ;
2906 
2907  if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
2908  std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
2909  throw ErrorClass( "index mapping problem in generating multicommodity cut");
2910  }else{
2911  m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
2912  cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
2913  }
2914  }
2915  }
2916  kount++;
2917  }
2918 
2919  //for(j2 = m_numHubs; j2 < j1; j2++){
2920  for(j2 = j1 + 1; j2 < m_numNodes; j2++){
2921 
2922  indexPair.first = j1;
2923  indexPair.second = j2;
2924 
2925  uVal = solver->osiSolver->getColSolution()[ kount];
2926  objVal += uVal*solver->osiSolver->getObjCoefficients()[ kount];
2927 
2928  //if(solver->osiSolver->getObjCoefficients()[ kount] < -0.001) std::cout << m_variableNames[ m_xVarIndexMap[ indexPair] ] << " valueeeee = " << uVal << " coeff " << solver->osiSolver->getObjCoefficients()[ kount] << " kount " << kount << std::endl;
2929 
2930  if (j2 == (i + m_numHubs) ){
2931 
2932  if( (wVal - uVal) > m_osDecompParam.zeroTol || (uVal - wVal) > m_osDecompParam.zeroTol ){
2933  //variable (wVal - uVal)*x(j1, shat)
2934  cutString << " +";
2935  cutString << (wVal - uVal);
2936  cutString << "*" ;
2937  //cutString << m_nodeName[ j1 ];
2938  //cutString << "," ;
2939  //cutString << m_nodeName[ i+ m_numHubs] ;
2940  //cutString << ") " ;
2941  m_BmatrixVal[ m_numBmatrixNonz ] = (wVal - uVal)/wVal ;
2942  //get index for j1, j2 with j1 > j2
2943  //m_BmatrixIdx[ m_numBmatrixNonz++ ] = j1*(m_numNodes - 1) + j2 ;
2944 
2945  if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
2946  std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
2947  throw ErrorClass( "index mapping problem in generating multicommodity cut");
2948  }else{
2949  m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
2950  cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
2951  }
2952  }
2953 
2954  }else{
2955 
2956  if( (-uVal > m_osDecompParam.zeroTol) || (uVal > m_osDecompParam.zeroTol) ){
2957  //variable -uVal*x(j1, j2)
2958  cutString << " ";
2959  cutString << - uVal;
2960  cutString << "*" ;
2961  //cutString << m_nodeName[ j1 ];
2962  //cutString << "," ;
2963  //cutString << m_nodeName[ j2] ;
2964  //cutString << ") " ;
2965  m_BmatrixVal[ m_numBmatrixNonz ] = -uVal/wVal ;
2966  //get index for j1, j2 with j1 > j2
2967  //m_BmatrixIdx[ m_numBmatrixNonz++ ] = j1*(m_numNodes - 1) + j2 ;
2968  if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
2969  std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
2970  throw ErrorClass( "index mapping problem in generating multicommodity cut");
2971  }else{
2972  m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
2973  cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
2974  }
2975  }
2976  }
2977  kount++;
2978  }
2979  }//end loop on j1
2980  cutString << std::endl;
2981  //std::cout << cutString.str() << " kount = " << kount << std::endl;
2982  //std::cout << "OPTIMAL OBJECTIVE FUNCTION VALUE = " << objVal << std::endl;
2983  //std::cout << "OPTIMAL W VALUE = " << wVal << std::endl;
2984 
2985  //now add the cut
2986  //
2987  m_numBmatrixCon++;
2989  //we have taken care of the B-matrix -- the Xij space
2990  //now take care the theta indexes and values
2991 
2992  //scatter the constraint in the x - variables, we
2993  //are scattering the B matrix constraint just found
2994 
2995  for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2996  j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2997 
2998  m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 1;
2999  scatterValues[ m_BmatrixIdx[ j] ] = m_BmatrixVal[ j];
3000  //std::cout << m_variableNames[ m_BmatrixIdx[ j] ] << " xkij cut coeff = " << m_BmatrixVal[ j] << std::endl;
3001 
3002  }
3003 
3004  //done with scatter operation
3005 
3006 
3007 
3008  if(numTheta != m_numThetaVar)throw
3009  ErrorClass( "Inconsistent Number of theta variables in multicommondity cut separation problem" );
3010 
3011 
3012  for(thetaIdx = 0; thetaIdx < m_numThetaVar ; thetaIdx++){
3013  //make sure we only consider thetaIdx in convexity row k
3014 
3015  if(m_convexityRowIndex[ thetaIdx] == k){
3016 
3017 
3018  rowValue = 0;
3019  for(j = m_thetaPnt[ thetaIdx]; j < m_thetaPnt[ thetaIdx + 1] ; j++){
3020 
3021  //std::cout << "row value = " << rowValue << " m_tmpScatterArray[ j ] " << m_tmpScatterArray[ j ] << std::endl;
3022 
3023  rowValue += m_tmpScatterArray[ m_thetaIndex[ j] ]*scatterValues[ m_thetaIndex[ j] ];
3024  //std::cout << "theta index " << thetaIdx << " " << m_variableNames[ m_thetaIndex[ j] ] << " = " << m_BmatrixVal[ j] << " row value = " << rowValue << std::endl;
3025  }
3026 
3027  if(rowValue > m_osDecompParam.zeroTol || -rowValue > m_osDecompParam.zeroTol){
3028  //std::cout << "numNewRows = " << numNewRows << " m_newRowNonz[ numNewRows ] " << m_newRowNonz[ numNewRows ] << std::endl;
3029  m_newRowColumnIdx[ numNewRows][ m_newRowNonz[ numNewRows ] ] = thetaIdx ;
3030  m_newRowColumnValue[ numNewRows][ m_newRowNonz[ numNewRows]++ ] = rowValue;
3031  }
3032 
3033 
3034  }//end of if on convexity row
3035 
3036 
3037  }//end loop on thetaIdx
3038 
3039  m_newRowLB[ numNewRows] = -OSDBL_MAX;
3040  m_newRowUB[ numNewRows] = 0;
3041  numNewRows++;
3042 
3043 
3044  //zero out the scatter array again
3045  for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
3046  j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
3047 
3048  m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 0;
3049  scatterValues[ m_BmatrixIdx[ j] ] = 0;
3050 
3051  }
3052  //
3053  //done adding cut
3054 
3055  //kipp -- don't forget to increment the artificial variable since a cut is added
3056  //
3057 
3058 
3059 
3062 
3063  }//end iff on positive obj value
3064  //set objcoefficient back to 0
3065  solver->osiSolver->setObjCoeff( i, 0 );
3066 
3067  //we have a cut so break from the loop
3068  //if(foundCut == true) break;
3069  }//end loop over i
3070  std::cout << std::endl << std::endl;
3071 
3072  //reset the coefficients
3073 
3074  for(i = 0; i < numVar; i++ ) solver->osiSolver->setObjCoeff( i, 0 );
3075 
3076 
3077  }//end loop over k
3078 
3079  delete[] wcoeff;
3080  wcoeff = NULL;
3081 
3082  delete[] scatterValues;
3083  scatterValues = NULL;
3084 
3085  numNonz = m_newRowNonz;
3086  colIdx = m_newRowColumnIdx;
3087  values = m_newRowColumnValue;
3088  rowUB = m_newRowUB;
3089  rowLB = m_newRowLB;
3090 
3091  for(i = 0; i < numNewRows; i++){
3092 
3093  //we found a row, add the corresponding artificial variables
3094  //to the transformation matrix
3095  m_numThetaVar++;
3098  }
3099 
3100  return;
3101 
3102 
3103  } catch (const ErrorClass& eclass) {
3104  if(wcoeff != NULL){
3105  delete[] wcoeff;
3106  wcoeff = NULL;
3107  }
3108 
3109 
3110 
3111  if(scatterValues != NULL) {
3112  delete[] scatterValues;
3113  scatterValues = NULL;
3114  }
3115 
3116 
3117  throw ErrorClass(eclass.errormsg);
3118 
3119  }
3120 
3121 
3122 }//end getCutsMultiCommod
3123 
3124 void OSBearcatSolverXij::getCutsX(const double* x, const int numX,
3125  int &numNewRows, int* &numNonz, int** &colIdx,
3126  double** &values, double* &rowLB, double* &rowUB) {
3127  //critical -- we are assuming that the size of x is going to be
3128  // m_numNodes*(m_numNodes - 1)
3129 
3130  int i;
3131  int j;
3132  int k;
3133  int index;
3134  int rowKount;
3135 
3136 
3137  int indexAdjust = m_numNodes - m_numHubs;
3138  double* tmpRhs;
3139  int numSepRows = m_osinstanceSeparation->getConstraintNumber() ;
3140 
3141  tmpRhs = new double[ numSepRows ];
3142 
3143  for(i = 0; i < numSepRows; i++){
3144 
3145  tmpRhs[ i] = 0;
3146  }
3147 
3148  try{
3150 
3151  for(i = 0; i < numX; i++){
3152 
3153  //get a postive theta
3154  if(x[ i] > m_osDecompParam.zeroTol){
3155 
3156  //the row index for x_{ij}
3157  rowKount = m_separationIndexMap[ i ];
3158 
3159  if(rowKount < OSINT_MAX ){
3160 
3161  tmpRhs[ rowKount] -= x[ i];
3162 
3163  }
3164 
3165  }
3166  }// end i loop
3167 
3168  for(i = indexAdjust; i < numSepRows - 1; i++){
3169 
3170  if(-tmpRhs[ i] > 1 + m_osDecompParam.zeroTol ){
3171 
3172  // quick and dirty does x_{ij} + x_{ji} <= 1 cut off anything
3173  //std::cout << " tmpRhs[ i] = " << tmpRhs[ i] << std::endl;
3174  //which variable is this
3175  //kipp this an inefficient way of finding i and j -- improve this
3176  int tmpKount = indexAdjust;
3177  for(int i1 = m_numHubs; i1 < m_numNodes; i1++){
3178 
3179  for(int j1 = i1+1; j1 < m_numNodes; j1++){
3180 
3181  if(tmpKount == i){
3182 
3183  numNewRows = 1;
3184 
3185  m_newRowNonz[ 0] = 2;
3186  m_newRowUB[ 0] = 1;
3187  m_newRowLB[ 0] = 0;
3188 
3189  m_newRowColumnIdx[ 0][ 0 ] = i1*(m_numNodes - 1) + j1 - 1;
3190  m_newRowColumnIdx[ 0][ 1 ] = j1*(m_numNodes - 1) + i1;
3191  m_newRowColumnValue[ 0][ 0] = 1;
3192  m_newRowColumnValue[ 0][ 1] = 1;
3193 
3194  numNonz = m_newRowNonz;
3195  colIdx = m_newRowColumnIdx;
3196  values = m_newRowColumnValue;
3197  rowUB = m_newRowUB;
3198  rowLB = m_newRowLB;
3199 
3200  delete[] tmpRhs;
3201  tmpRhs = NULL;
3202  return;
3203 
3204 
3205 
3206  }
3207 
3208  tmpKount++;
3209 
3210  }// end loop on j1
3211 
3212  }//end loop on i1
3213 
3214 
3215  }//end if on tmpRHS
3216 
3217  m_separationClpModel->setRowUpper(i, tmpRhs[ i] );
3218  m_separationClpModel->setRowLower(i, tmpRhs[ i] );
3219 
3220  }//end loop on i
3221 
3222 
3223  //std::cout << m_osinstanceSeparation->printModel() << std::endl;
3224 
3225 
3226  std::vector<int> dualIdx;
3227  std::vector<int>::iterator vit1;
3228  std::vector<int>::iterator vit2;
3229 
3230  //if the objective function value is greater than zero we have a cut
3231  //the cut is based on the nodes with dual value - 1
3232 
3233  for(k = 0; k < indexAdjust; k++){
3234  std::cout << std::endl << std::endl;
3235 
3236 
3237  m_separationClpModel->setRowUpper(k, 0.0);
3238 
3239 
3240  m_separationClpModel->primal();
3241 
3242  if(m_separationClpModel->getObjValue() > m_osDecompParam.zeroTol){
3243  std::cout << "DOING SEPARATION FOR NODE " << k + m_numHubs << std::endl;
3244  std::cout << "SEPERATION OBJ = " << m_separationClpModel->getObjValue() << std::endl;
3245  numNewRows = 1;
3246  m_newRowNonz[ 0] = 0;
3247  m_newRowLB[ 0] = 0;
3248 
3249  for(i = 0; i < m_numNodes - m_numHubs ; i++){
3250  //std::cout << m_osinstanceSeparation->getConstraintNames()[ i] << " = " << m_separationClpModel->getRowPrice()[ i] << std::endl;
3251  if( m_separationClpModel->getRowPrice()[ i] - m_osDecompParam.zeroTol <= -1) dualIdx.push_back( i) ;
3252  }
3253 
3254  for (vit1 = dualIdx.begin(); vit1 != dualIdx.end(); vit1++) {
3255 
3256  i = *vit1 + m_numHubs;
3257 
3258  for (vit2 = dualIdx.begin(); vit2 != dualIdx.end(); vit2++) {
3259 
3260  j = *vit2 + m_numHubs;
3261 
3262  if( i > j ){
3263 
3264  index = i*(m_numNodes -1) + j;
3265  std::cout << "CUT VARIABLE = " << m_variableNames[ index] <<std::endl;
3266  m_newRowColumnValue[ 0][ m_newRowNonz[ 0] ] = 1.0;
3267  m_newRowColumnIdx[ 0][ m_newRowNonz[ 0]++ ] = index;
3268 
3269  }else{
3270 
3271  if( i < j ){
3272 
3273  index = i*(m_numNodes -1) + j - 1;
3274  std::cout << "CUT VARIABLE = " << m_variableNames[ index] <<std::endl;
3275  m_newRowColumnValue[ 0][ m_newRowNonz[ 0] ] = 1.0;
3276  m_newRowColumnIdx[ 0][ m_newRowNonz[ 0]++ ] = index;
3277 
3278  }
3279  }
3280 
3281  }//end for on vit2
3282  }//end for on vit1
3283 
3284 
3285  m_newRowUB[ 0] = dualIdx.size() - 1;
3286 
3287  dualIdx.clear();
3288  //reset
3289  // don't adjust the kludge row
3290  for(i = indexAdjust; i < numSepRows - 1; i++){
3291 
3292  m_separationClpModel->setRowUpper(i, 0.0 );
3293  m_separationClpModel->setRowLower(i, 0.0 );
3294 
3295 
3296  }
3297  m_separationClpModel->setRowUpper(k, 1.0);
3298  delete[] tmpRhs;
3299  tmpRhs = NULL;
3300 
3301 
3302  numNonz = m_newRowNonz;
3303  colIdx = m_newRowColumnIdx;
3304  values = m_newRowColumnValue;
3305  rowUB = m_newRowUB;
3306  rowLB = m_newRowLB;
3307 
3308  return;
3309 
3310 
3311 
3312  }//end if on obj value
3313  m_separationClpModel->setRowUpper(k, 1.0);
3314  dualIdx.clear();
3315 
3316  }//end loop on k
3317 
3318  //if we are here there was no cut
3319  //reset
3320  // don't adjust the kludge row
3321  for(i = indexAdjust; i < numSepRows - 1; i++){
3322 
3323  m_separationClpModel->setRowUpper(i, 0.0 );
3324  m_separationClpModel->setRowLower(i, 0.0 );
3325 
3326 
3327  }
3328 
3329  delete[] tmpRhs;
3330  tmpRhs = NULL;
3331 
3332  } catch (const ErrorClass& eclass) {
3333 
3334  throw ErrorClass(eclass.errormsg);
3335 
3336  }
3337 
3338 
3339 }//end getCutsX
3340 
3341 
3342 void OSBearcatSolverXij::calcReducedCost( const double* yA, const double* yB){
3343 
3344  int k;
3345  int i;
3346  int j;
3347  int l;
3348  int kount;
3349 
3350  int tmpVal;
3351  tmpVal = m_numNodes - 1;
3352 
3353  for(k = 0; k < m_numHubs; k++){
3354  kount = 0;
3355 
3356  for(l = 1; l <= m_upperBoundL[ k]; l++){
3357 
3358 
3359  for(i = 0; i< m_numNodes; i++){
3360 
3361 
3362 
3363  for(j = 0; j < i; j++){
3364 
3365  if(j < m_numHubs){
3366 
3367  m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j ] ;
3368 
3369  }else{
3370 
3371  m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j ] - yA[ j - m_numHubs] ;
3372  }
3373 
3374 
3375  }
3376 
3377 
3378 
3379  for(j = i + 1; j < m_numNodes; j++){
3380 
3381 
3382  if(j < m_numHubs){
3383 
3384  m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j - 1 ];
3385 
3386  } else {
3387 
3388 
3389  m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j - 1 ] - yA[ j - m_numHubs ];
3390 
3391  }
3392 
3393  }
3394 
3395 
3396  }
3397 
3398 
3399  }//end l loop
3400 
3401 
3402  }//end hub loop
3403 
3404  //now adjust for tour breaking constraints
3405 
3406 
3407  int startPnt ;
3408 
3409  for(i = 0; i < m_numBmatrixCon; i++){
3410 
3411  //get the xij
3412 
3413  for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
3414 
3415 
3416 
3417  for(k = 0; k < m_numHubs; k++){
3418 
3419 
3420  for(l = 1; l <= m_upperBoundL[ k]; l++){
3421 
3422  //startPnt = k*m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes) + (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
3423  startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
3424 
3425  if(m_BmatrixRowIndex[ i] == -1 || m_BmatrixRowIndex[ i] == k ) m_rc[ k][ startPnt + m_BmatrixIdx[ j] ] -= yB[ i]*m_BmatrixVal[ j];
3426 
3427  }
3428 
3429  }
3430 
3431 
3432  }
3433 
3434  }
3435 
3436 }//end calcReducedCost
3437 
3438 
3440 
3441  int i;
3442  int j;
3443  int kount;
3444 
3445  kount = 0;
3446 
3447  for(i = 0; i< m_numNodes; i++){
3448 
3449  //if we have (i, j) where j is hub then do not subtract off phi[ j]
3450  for(j = 0; j < i; j++){
3451 
3452  if(m_nodeName[ i] == "") m_variableNames[ kount] = makeStringFromInt("x(" , i);
3453  else m_variableNames[ kount] = "x(" + m_nodeName[ i];
3454  if(m_nodeName[ i] == "") m_variableNames[ kount] += makeStringFromInt( "," , j);
3455  else m_variableNames[ kount] += "," + m_nodeName[ j];
3456  m_variableNames[ kount] += ")";
3457  //std::cout << "GAIL VARIABLE NAME " << m_variableNames[ kount] << std::endl;
3458 
3459  kount++;
3460 
3461  }
3462 
3463  for(j = i + 1; j < m_numNodes; j++){
3464 
3465  if(m_nodeName[ i] == "") m_variableNames[ kount] = makeStringFromInt("x(" , i);
3466  else m_variableNames[ kount] = "x(" + m_nodeName[ i];
3467  if(m_nodeName[ i] == "") m_variableNames[ kount] += makeStringFromInt( "," , j);
3468  else m_variableNames[ kount] += "," + m_nodeName[ j];
3469  m_variableNames[ kount] += ")";
3470 
3471  //std::cout << "GAIL VARIABLE NAME " << m_variableNames[ kount] << std::endl;
3472  kount++;
3473 
3474  }
3475 
3476 
3477  }
3478 }//end createVariableNames
3479 
3481 
3482  //arrays for the coupling constraint matrix
3483  //this is in the x variable space, not theta
3484  //int* m_pntAmatrix;
3485  //int* m_Amatrix;
3486 
3487 
3488  int i;
3489  int j;
3490  int numNonz;
3491 
3492  //loop over nodes
3493  m_pntAmatrix[ 0] = 0;
3494  numNonz = 0;
3495 
3496  for(j = m_numHubs; j < m_numNodes; j++){
3497 
3498 
3499  for(i = 0; i < j; i++){
3500 
3501  m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j - 1 ;
3502 
3503  }
3504 
3505  for(i = j + 1; i < m_numNodes; i++){
3506 
3507  m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j ;
3508 
3509  }
3510 
3511  m_pntAmatrix[ j - m_numHubs + 1] = numNonz;
3512 
3513  }
3514 
3515  /*
3516  for(i = 0; i < m_numNodes - m_numHubs; i++){
3517 
3518  for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
3519 
3520  std::cout << m_variableNames[ m_Amatrix[ j ] ] << std::endl;
3521 
3522  }
3523 
3524  }
3525  */
3526 
3527 }//end createAmatrix
3528 
3529 void OSBearcatSolverXij::pauHana( std::vector<int> &m_zOptIndexes,
3530  std::vector<double> &m_zRootLPx_vals,
3531  int numNodes, int numColsGen, std::string message){
3532 
3533  std::cout << std::endl;
3534  std::cout << " PAU HANA TIME! " << std::endl;
3535 
3536  double cost;
3537  cost = 0;
3538  std::vector<int>::iterator vit;
3539  try{
3540  int i;
3541  int j;
3542 
3543 
3544  //we better NOT have any artifical variables positive
3545  //for(i = 0; i < numVarArt ; i++){
3546  //
3547  // if(theta[ i] > m_osDecompParam.zeroTol) throw ErrorClass("we have a positive artificial variable");
3548  //}
3549 
3550  //for(i = 0; i < m_numThetaVar ; i++){
3551 
3552  //cost += theta[ i]*m_thetaCost[ i ];
3553  //std::cout << "COLUMN VALUE = " << theta[ i] << std::endl;
3554 
3555  //}
3556 
3557  double routeDemand;
3558 
3559 
3560  int ivalue;
3561  int jvalue;
3562 
3563  for(vit = m_zOptIndexes.begin() ; vit != m_zOptIndexes.end(); vit++){
3564 
3565  i = *vit;
3566  std::cout << "x variables for column " << i << " CONVEXITY ROW = "<< m_convexityRowIndex[ i] << std::endl;
3567 
3568 
3569  //cost += m_thetaCost[ i ];
3570  routeDemand = 0;
3571 
3572  for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
3573 
3574 
3575  //std::cout << "INDEX = " << m_thetaIndex[ j] << std::endl;
3576  std::cout << m_variableNames[ m_thetaIndex[ j] ] << " = " << 1 << " DISTANCE = " << m_cost[ m_thetaIndex[ j] ] << std::endl;
3577 
3578  ivalue = m_thetaIndex[ j] /(m_numNodes - 1);
3579 
3580  jvalue = m_thetaIndex[ j] - ivalue*(m_numNodes - 1);
3581 
3582  if( jvalue >= ivalue ){
3583  //std::cout << " i NODE NUMBER = " << ivalue << std::endl;
3584  //std::cout << " j NODE NUMBER = " << jvalue + 1 << std::endl;
3585  //std::cout << " DEMAND = = " << m_demand[ jvalue + 1] << std::endl;
3586  routeDemand += m_demand[ jvalue + 1];
3587 
3588 
3589  }else{
3590  //std::cout << " i NODE NUMBER = " << ivalue << std::endl;
3591  //std::cout << " j NODE NUMBER = " << jvalue << std::endl;
3592  //std::cout << " DEMAND = = " << m_demand[ jvalue ] << std::endl;
3593  routeDemand += m_demand[ jvalue ];
3594  }
3595  }
3596 
3597  std::cout << "route demand = " << routeDemand << std::endl << std::endl;
3598  std::cout << "distance for this route " << m_thetaCost[ i ] / routeDemand << std::endl;
3599 
3600  }
3601 
3602 
3603  //std::cout << "cost = " << cost << std::endl << std::endl;
3604 
3605  std::cout << std::endl << std::endl;
3606  std::cout << message << std::endl;
3607  std::cout << "LINEAR PROGRAMMING RELAXATION VALUE = " << m_rootLPValue << std::endl;
3608  std::cout << "NONLINEAR RELAXATION VALUE = " << calcNonlinearRelax( m_zRootLPx_vals) << std::endl;
3609  std::cout << "LOWER BOUND VALUE = " << m_bestLPValue << std::endl;
3610  std::cout << "FINAL BEST IP SOLUTION VALUE = " << m_bestIPValue << std::endl;
3611  std::cout << "NUMBER OF COLUMNS IN FINAL MASTER = " << m_numThetaVar << std::endl;
3612  //std::cout << "NUMBER OF GENERATED COLUMNS = " << m_numThetaVar - 2*m_numNodes - 2*m_numBmatrixCon << std::endl;
3613  //the original master has m_numHubs + m_numNodes columns
3614  std::cout << "NUMBER OF GENERATED COLUMNS = " << numColsGen << std::endl;
3615  std::cout << "NUMBER OF GENERATED CUTS = " << m_numBmatrixCon << std::endl;
3616  std::cout << "NUMBER OF NODES = " << numNodes << std::endl;
3617  std::cout << " PAU!!!" << std::endl;
3618 
3619  std::cout << std::endl << std::endl;
3620 
3621  std::cout << std::endl << std::endl;
3622  }catch (const ErrorClass& eclass) {
3623 
3624  throw ErrorClass(eclass.errormsg);
3625 
3626  }
3627 
3628 }//end pauHana -- no pun intended
3629 
3630 double OSBearcatSolverXij::calcNonlinearRelax( std::vector<double> &m_zRootLPx_vals){
3631 
3632 
3633  std::vector<double>::iterator dvit;
3634  std::vector<double>::iterator dvit2;
3635  int index = 0;
3636  int i;
3637  int j;
3638  int ivalue;
3639  int jvalue;
3640  double *hubDemand = NULL;
3641  double *hubCost = NULL;
3642  hubDemand = new double[m_numHubs ];
3643  hubCost = new double[ m_numHubs];
3644 
3645  double objVal = 0.0;
3646  double objValAux = 0.0;
3647 
3648  std::map<int, std::vector<double> > extPointDemands;
3649  std::map<int, std::vector<double> > extPointCosts;
3650  std::map<int, std::vector<double> > extPointValues;
3651 
3652  std::map<int, std::vector<double> >::iterator mit;
3653 
3654  double tmpDemand;
3655  double tmpCost;
3656 
3657  for(i = 0; i < m_numHubs; i++){
3658  hubDemand[ i] = 0;
3659  hubCost[ i] = 0;
3660 
3661  }
3662 
3663  try{
3664  for (dvit = m_zRootLPx_vals.begin(); dvit < m_zRootLPx_vals.end(); dvit++ ){
3665 
3666  if(*dvit > m_osDecompParam.zeroTol){
3667  //std::cout << std::endl;
3668  //std::cout << "LP VALUE = " << *dvit << std::endl;
3669  //std::cout << "x variables for column " << index << " CONVEXITY ROW = "<< m_convexityRowIndex[ index] << std::endl;
3670 
3671  tmpDemand = 0;
3672  tmpCost = 0;
3673 
3674  for(j = m_thetaPnt[ index]; j < m_thetaPnt[ index + 1] ; j++){
3675 
3676 
3677  //std::cout << "INDEX = " << m_thetaIndex[ j] << std::endl;
3678  //std::cout << m_variableNames[ m_thetaIndex[ j] ] << " = " << 1 << " DISTANCE = " << m_cost[ m_thetaIndex[ j] ] << std::endl;
3679  hubCost[ m_hubPoint[ m_convexityRowIndex[ index] ] ] += m_cost[ m_thetaIndex[ j] ]*( *dvit);
3680  tmpCost += m_cost[ m_thetaIndex[ j] ];
3681 
3682  ivalue = m_thetaIndex[ j] /(m_numNodes - 1);
3683  jvalue = m_thetaIndex[ j] - ivalue*(m_numNodes - 1);
3684 
3685  if( jvalue >= ivalue ){
3686  //std::cout << " i NODE NUMBER = " << ivalue << std::endl;
3687  //std::cout << " j NODE NUMBER = " << jvalue + 1 << std::endl;
3688  //std::cout << " DEMAND = = " << m_demand[ jvalue + 1] << std::endl;
3689  hubDemand[ m_hubPoint[ m_convexityRowIndex[ index] ] ] += m_demand[ jvalue + 1]*( *dvit);
3690  tmpDemand += m_demand[ jvalue + 1];
3691 
3692 
3693 
3694  }else{
3695  //std::cout << " i NODE NUMBER = " << ivalue << std::endl;
3696  //std::cout << " j NODE NUMBER = " << jvalue << std::endl;
3697  //std::cout << " DEMAND = = " << m_demand[ jvalue ] << std::endl;
3698  hubDemand[ m_hubPoint[ m_convexityRowIndex[ index] ] ] += m_demand[ jvalue ]*( *dvit);
3699  tmpDemand += m_demand[ jvalue];
3700 
3701  }
3702  }
3703 
3704  //std::cout << "TOTAL DEMAND = = " << tmpDemand << " TOTAL COST = " << tmpCost << std::endl;
3705  extPointDemands[ m_convexityRowIndex[ index] ].push_back( tmpDemand);
3706  extPointCosts[ m_convexityRowIndex[ index] ].push_back( tmpCost);
3707  extPointValues[ m_convexityRowIndex[ index] ].push_back(*dvit);
3708 
3709  }//end if
3710 
3711  //get x values
3712 
3713  index++;
3714 
3715  }
3716 
3717  int mapSize;
3718  objValAux = 0;
3719  for (i = 0; i < m_numHubs; i++){
3720  //std::cout << "HUB DEMAND " << hubDemand[ m_hubPoint[ i]] << std::endl;
3721  //std::cout << "HUB COST " << hubCost[ m_hubPoint[ i]] << std::endl;
3722  objVal += hubDemand[ m_hubPoint[ i] ]*hubCost[ m_hubPoint[ i]];
3723  tmpDemand = 0;
3724  tmpCost = 0;
3725 
3726  mapSize = extPointDemands[ m_hubPoint[ i] ].size();
3727  //std::cout << std::endl;
3728  //std::cout << " HUB NUBMER = " << m_hubPoint[ i] << std::endl;
3729 
3730  for (j = 0; j < mapSize; j++){
3731 
3732  tmpDemand += extPointDemands[ m_hubPoint[ i] ][j]*extPointValues[ m_hubPoint[ i] ][j];
3733  tmpCost += extPointCosts[ m_hubPoint[ i] ][j]*extPointValues[ m_hubPoint[ i] ][j];
3734 
3735  //std::cout << " DEMAND = " << extPointDemands[ m_hubPoint[ i] ][j] <<
3736  // " COST " << extPointCosts[ m_hubPoint[ i] ][j] << std::endl;
3737  objValAux +=
3738  extPointCosts[ m_hubPoint[ i] ][j]*extPointDemands[ m_hubPoint[ i] ][j]*extPointValues[ m_hubPoint[ i] ][j];
3739 
3740  }
3741 
3742  //std::cout << std::endl;
3743  //std::cout << "HUB DEMAND 2 " << tmpDemand << std::endl;
3744  //std::cout << "HUB COST 2 " << tmpCost << std::endl;
3745 
3746  }
3747 
3748  std::cout << "AUXILIARY VALUE COST = " << objValAux << std::endl;
3749 
3750  //garbage collection
3751  delete[] hubDemand;
3752  hubDemand = NULL;
3753  delete[] hubCost;
3754  hubCost = NULL;
3755 
3756 
3757  return objVal;
3758 
3759  }catch (const ErrorClass& eclass) {
3760 
3761  //garbage collection
3762  delete[] hubDemand;
3763  hubDemand = NULL;
3764  delete[] hubCost;
3765  hubCost = NULL;
3766 
3767 
3768  throw ErrorClass(eclass.errormsg);
3769 
3770  }
3771 
3772 
3773 }// end calcNonlinearRelax
3774 
3775 
3777 
3778 
3779 
3780 
3781  m_osinstanceSeparation = NULL;
3782 
3783  //add linear constraint coefficients
3784  //number of values will nodes.size() the coefficients in the node constraints
3785  //plus coefficients in convexity constraints which is the number of varaibles
3786  int kountNonz;
3787  int kount;
3788  int startsIdx;
3789  //we build these on nodes that do not include the hubs
3790  int numYvar = (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1);
3791  int numVvar = m_numNodes - m_numHubs;
3792  //the plus 1 is for the kludge row
3793  int numCon = (m_numNodes - m_numHubs) + (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1)/2 + 1;
3794  double *values = new double[ 2*numYvar + 2*numVvar];
3795  int *indexes = new int[ 2*numYvar + 2*numVvar];
3796  int *starts = new int[ numYvar + numVvar + 1];
3797  starts[ 0] = 0;
3798  startsIdx = 0;
3799  startsIdx++;
3800  kountNonz = 0;
3801  int i;
3802  int j;
3803 
3804 
3805  std::string separationVarName;
3806  std::string separationConName;
3807 
3808  try {
3809 
3811 
3812  //start building the separation instance
3813 
3814  m_osinstanceSeparation->setInstanceDescription("The Tour Breaking Separation Problem");
3815 
3816 
3817  // now the constraints
3819 
3820 
3821  //add the node rows
3822  for( i = 0; i < m_numNodes - m_numHubs ; i++){
3823 
3824  m_osinstanceSeparation->addConstraint(i, makeStringFromInt("nodeRow_", i+ m_numHubs ) , 0.0, 1.0, 0);
3825 
3826  }
3827 
3828  //add the variable rows rows
3829 
3830  int rowKounter;
3831  rowKounter = m_numNodes - m_numHubs;
3832 
3833  for(i = m_numHubs; i < m_numNodes; i++){
3834 
3835 
3836 
3837  for(j = i+1; j < m_numNodes; j++){
3838  separationConName = makeStringFromInt("Row_(", i);
3839  separationConName += makeStringFromInt(",", j);
3840  separationConName += ")";
3841 
3842  m_osinstanceSeparation->addConstraint(rowKounter++, separationConName , 0, 0, 0);
3843  }
3844 
3845  }
3846 
3847  // the klude row so we have +/-1 in every column
3848  m_osinstanceSeparation->addConstraint(rowKounter++, "kludgeRow" , 0, m_numNodes, 0);
3849 
3850  // the variables
3851  m_osinstanceSeparation->setVariableNumber( numYvar + numVvar);
3852 
3853 
3854 
3855  std::cout << "NUMBER OF VARIABLES SET = " << numYvar + numVvar << std::endl;
3856  //add the v variables
3857  for(i = 0; i < numVvar; i++){
3858 
3859  separationVarName = makeStringFromInt("v", i + m_numHubs);
3860 
3861  m_osinstanceSeparation->addVariable(i, separationVarName, 0, 1, 'C');
3862 
3863  values[ kountNonz ] = -1.0;
3864  indexes[ kountNonz ] = i;
3865  kountNonz++;
3866 
3867  values[ kountNonz ] = 1.0;
3868  indexes[ kountNonz ] = rowKounter - 1;
3869  kountNonz++;
3870 
3871 
3872 
3873  starts[ startsIdx++ ] = kountNonz;
3874 
3875 
3876  }
3877  //add the y variables
3878  kount = numVvar;
3879 
3880  int i1;
3881  int j1;
3882  int kountCon;
3883  kountCon = m_numNodes - m_numHubs;
3884 
3885  for(i1 = 0; i1 < m_numNodes - m_numHubs; i1++){
3886 
3887 
3888  i = i1 + m_numHubs;
3889 
3890  for(j1 = i1 + 1; j1 < m_numNodes - m_numHubs; j1++){
3891 
3892 
3893  j = j1 + m_numHubs;
3894 
3895  separationVarName = makeStringFromInt("y(", i);
3896  separationVarName += makeStringFromInt(",", j);
3897  separationVarName += ")";
3898  m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
3899 
3900  //map the variable to row kountCon
3901 
3902  // i < j case
3903  m_separationIndexMap[ i*(m_numNodes - 1) + (j - 1) ] = kountCon;
3904 
3905  values[ kountNonz ] = 1.0;
3906  indexes[ kountNonz ] = i1;
3907  kountNonz++;
3908 
3909  values[ kountNonz ] = -1.0;
3910  indexes[ kountNonz ] = kountCon ;
3911  kountNonz++;
3912 
3913  starts[ startsIdx++ ] = kountNonz;
3914 
3915 
3916 
3917 
3918  separationVarName = makeStringFromInt("y(", j );
3919  separationVarName += makeStringFromInt(",", i);
3920  separationVarName += ")";
3921  m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
3922 
3923  values[ kountNonz ] = 1.0;
3924  indexes[ kountNonz ] = j1;
3925  kountNonz++;
3926 
3927  // i < j case
3928  m_separationIndexMap[ j*(m_numNodes - 1) + i ] = kountCon;
3929 
3930  values[ kountNonz ] = -1.0;
3931  indexes[ kountNonz ] = kountCon ;
3932  kountNonz++;
3933 
3934  starts[ startsIdx++ ] = kountNonz;
3935 
3936 
3937  kountCon++;
3938 
3939 
3940  }
3941 
3942  }
3943 
3944  std::cout << "NUMBER OF VARIABLES ADDED = " << kount << std::endl;
3945 
3946  // now add the objective function
3948  SparseVector *objcoeff = new SparseVector( numVvar);
3949 
3950 
3951  for(i = 0; i < numVvar; i++){
3952 
3953  objcoeff->indexes[ i] = i;
3954  objcoeff->values[ i] = 1.0;
3955 
3956  }
3957 
3958 
3959 
3960 
3961 
3962  m_osinstanceSeparation->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
3963  //now for the nonzeros
3964  //add the linear constraints coefficients
3966  values, 0, kountNonz - 1, indexes, 0, kountNonz - 1, starts, 0, startsIdx);
3967 
3968 
3969 
3970  //std::cout << m_osinstanceSeparation->printModel( ) << std::endl;
3971  //
3972  delete objcoeff;
3973 
3974 
3975  // now create the Clp model
3976 
3977 
3978  //below is temporty see if we can setup as a Clp network problem
3979  CoinPackedMatrix* matrix;
3980  bool columnMajor = true;
3981  double maxGap = 0;
3982  matrix = new CoinPackedMatrix(
3983  columnMajor, //Column or Row Major
3988  columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->indexes : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->indexes, //Pointer to start of minor dimension indexes -- change to allow for row storage
3990  0, 0, maxGap );
3991 
3992  ClpNetworkMatrix network( *matrix);
3993 
3994  m_separationClpModel = new ClpSimplex();
3995 
3996  //if( m_osinstanceSeparation->getObjectiveMaxOrMins()[0] == "min") osiSolver->setObjSense(1.0);
3997  m_separationClpModel->setOptimizationDirection( 1);
4002  );
4003 
4004  m_separationClpModel->factorization()->maximumPivots(200 + m_separationClpModel->numberRows() / 100);
4005 
4006 
4007  delete matrix;
4008 
4009  }catch (const ErrorClass& eclass) {
4010 
4011  throw ErrorClass(eclass.errormsg);
4012 
4013  }
4014 
4015  return NULL;
4016 }//end getSeparationInstance
4017 
4018 
4019 
4020 int OSBearcatSolverXij::getBranchingVar(const double* theta, const int numThetaVar ) {
4021 
4022  int varIdx;
4023  varIdx = -1;
4024  int i;
4025  int j;
4026  int numVar = m_numNodes*m_numNodes - m_numHubs ;
4027 
4028  double from1Distance;
4029  double from0Distance;
4030  double fraction;
4031  double minFraction;
4032 
4033  double *xvalues;
4034 
4035 
4036  xvalues = new double[ numVar];
4037  for(i = 0; i < numVar; i++){
4038  xvalues[ i] = 0;
4039  }
4040 
4041  try{
4042  if(numThetaVar != m_numThetaVar) throw ErrorClass("inconsistent number of variables in getBranchingVar");
4043  //loop over the fractional thetas
4044  for(i = 0; i < m_numThetaVar; i++){
4045 
4046  if( ( theta[ i ] > m_osDecompParam.zeroTol ) && ( theta[ i ] < 1 - m_osDecompParam.zeroTol ) ){
4047 
4048  for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
4049 
4050  xvalues[ m_thetaIndex[ j] ] += theta[ i ] ;
4051 
4052  }
4053 
4054  }
4055 
4056 
4057  }
4058 
4059  //let's branch on a variable in and out of hub first
4060  minFraction = 1.0;
4061  //ideally we find minFraction very close to .5
4062 
4063  for(i = 0; i < m_numHubs; i++){
4064 
4065  for( j = 0; j < i; j++){
4066 
4067  //j < i so the index is i*(m_numNodes - 1) + j
4068  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
4069  from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
4070  fraction = std::max(from1Distance, from0Distance);
4071  //try to find fractional variable that is the closest to .5
4072  if(fraction < minFraction){
4073 
4074  minFraction = fraction;
4075  varIdx = i*(m_numNodes - 1) + j;
4076  }
4077 
4078  }
4079 
4080  for(j = i + 1; j < m_numNodes; j++){
4081 
4082  //j < i so the index is i*(m_numNodes - 1) + j - 1
4083  //j < i so the index is i*(m_numNodes - 1) + j
4084  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
4085  from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
4086  fraction = std::max(from1Distance, from0Distance);
4087  //try to find fractional variable that is the closest to .5
4088  if(fraction < minFraction) {
4089 
4090  minFraction = fraction;
4091  varIdx = i*(m_numNodes - 1) + j - 1;
4092  }
4093 
4094 
4095  }
4096 
4097  }
4098 
4099  //if we have a candidate among arcs in/out of hubs, take it
4100 
4101  if(minFraction > 1 - m_osDecompParam.zeroTol){
4102 
4103  for(i = m_numHubs; i < m_numNodes; i++){
4104 
4105  for( j = 0; j < i; j++){
4106 
4107  //j < i so the index is i*(m_numNodes - 1) + j
4108  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
4109  from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
4110  fraction = std::max(from1Distance, from0Distance);
4111  //try to find fractional variable that is the closest to .5
4112  if(fraction < minFraction) {
4113 
4114  minFraction = fraction;
4115  varIdx = i*(m_numNodes - 1) + j ;
4116  }
4117 
4118  }
4119 
4120  for(j = i + 1; j < m_numNodes; j++){
4121 
4122  //j < i so the index is i*(m_numNodes - 1) + j - 1
4123  //j < i so the index is i*(m_numNodes - 1) + j
4124  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
4125  from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
4126  fraction = std::max(from1Distance, from0Distance);
4127  //try to find fractional variable that is the closest to .5
4128  if(fraction < minFraction) {
4129 
4130  minFraction = fraction;
4131  varIdx = i*(m_numNodes - 1) + j - 1;
4132  }
4133 
4134  }
4135 
4136  }
4137 
4138  }//end of if on minFraction
4139 
4140 
4141  delete[] xvalues;
4142 
4143  xvalues = NULL;
4144 
4145  return varIdx;
4146 
4147  }catch (const ErrorClass& eclass) {
4148 
4149  delete[] xvalues;
4150  xvalues = NULL;
4151 
4152  throw ErrorClass(eclass.errormsg);
4153 
4154  }
4155 
4156 
4157 }//end getBranchingVar Dense
4158 
4159 
4160 
4161 int OSBearcatSolverXij::getBranchingVar(const int* thetaIdx, const double* theta,
4162  const int numThetaVar) {
4163 
4164  int varIdx;
4165  varIdx = -1;
4166  int i;
4167  int j;
4168  int numVar = m_numNodes*m_numNodes - m_numHubs ;
4169  double from1Distance;
4170  double from0Distance;
4171  double fraction;
4172  double minFraction;
4173 
4174  double *xvalues;
4175 
4176 
4177  xvalues = new double[ numVar];
4178  for(i = 0; i < numVar; i++){
4179  xvalues[ i] = 0;
4180  }
4181 
4182  try{
4183  //loop over the fractional thetas
4184  //working with a sparse matrix
4185  for(i = 0; i < numThetaVar; i++){
4186 
4187  if( ( theta[ i ] > m_osDecompParam.zeroTol ) && ( theta[ i ] < 1 - m_osDecompParam.zeroTol ) ){
4188 
4189  for(j = m_thetaPnt[ thetaIdx[ i] ]; j < m_thetaPnt[ thetaIdx[ i] + 1] ; j++){
4190 
4191  xvalues[ m_thetaIndex[ j] ] += theta[ i ] ;
4192 
4193  }
4194 
4195  }
4196 
4197 
4198  }
4199 
4200 
4201  //let's branch on a variable in and out of hub first
4202  minFraction = 1.0;
4203  //ideally we find minFraction very close to .5
4204 
4205  for(i = 0; i < m_numHubs; i++){
4206 
4207  for( j = 0; j < i; j++){
4208 
4209  //j < i so the index is i*(m_numNodes - 1) + j
4210  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
4211  from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
4212  fraction = std::max(from1Distance, from0Distance);
4213  //try to find fractional variable that is the closest to .5
4214  if(fraction < minFraction){
4215 
4216  minFraction = fraction;
4217  varIdx = i*(m_numNodes - 1) + j;
4218  }
4219 
4220  }
4221 
4222  for(j = i + 1; j < m_numNodes; j++){
4223 
4224  //j < i so the index is i*(m_numNodes - 1) + j - 1
4225  //j < i so the index is i*(m_numNodes - 1) + j
4226  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
4227  from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
4228  fraction = std::max(from1Distance, from0Distance);
4229  //try to find fractional variable that is the closest to .5
4230  if(fraction < minFraction) {
4231 
4232  minFraction = fraction;
4233  varIdx = i*(m_numNodes - 1) + j - 1;
4234  }
4235 
4236 
4237  }
4238 
4239  }
4240 
4241  //if we have a candidate among arcs in/out of hubs, take it
4242 
4243  std::cout << "MIN FRACTION = " << minFraction << std::endl;
4244 
4245  if(minFraction > 1 - m_osDecompParam.zeroTol){
4246 
4247  for(i = m_numHubs; i < m_numNodes; i++){
4248 
4249 
4250 
4251  for( j = 0; j < i; j++){
4252 
4253  //j < i so the index is i*(m_numNodes - 1) + j
4254  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
4255  from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
4256  fraction = std::max(from1Distance, from0Distance);
4257  //try to find fractional variable that is the closest to .5
4258  if(fraction < minFraction) {
4259 
4260  minFraction = fraction;
4261  varIdx = i*(m_numNodes - 1) + j ;
4262  }
4263 
4264  }
4265 
4266  for(j = i + 1; j < m_numNodes; j++){
4267 
4268  //j < i so the index is i*(m_numNodes - 1) + j - 1
4269  //j < i so the index is i*(m_numNodes - 1) + j
4270  from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
4271  from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
4272  fraction = std::max(from1Distance, from0Distance);
4273  //try to find fractional variable that is the closest to .5
4274  if(fraction < minFraction) {
4275 
4276  minFraction = fraction;
4277  varIdx = i*(m_numNodes - 1) + j - 1;
4278  }
4279 
4280  }
4281 
4282  }
4283 
4284  }//end of if on minFraction
4285 
4286  //zero out the scatter array
4287 
4288  delete[] xvalues;
4289  xvalues = NULL;
4290 
4291  return varIdx;
4292 
4293  }catch (const ErrorClass& eclass) {
4294 
4295  delete[] xvalues;
4296  xvalues = NULL;
4297 
4298  throw ErrorClass(eclass.errormsg);
4299 
4300  }
4301 
4302 
4303 }//end getBranchingVar Sparse
4304 
4305 
4306 void OSBearcatSolverXij::getBranchingCut(const double* thetaVar, const int numThetaVar,
4307  const std::map<int, int> &varConMap, int &varIdx, int &numNonz,
4308  int* &indexes, double* &values) {
4309 
4310  //get a branching variable
4311  int i;
4312  int j;
4313  int kount;
4314  numNonz = 0;
4315  //keep numNonz at zero if there is no cut
4316  //there will be no cut if the xij is in conVarMap
4317 
4318  try{
4319 
4320  if(numThetaVar != m_numThetaVar) throw ErrorClass("inconsistent number of variables in getBranchingCut");
4321 
4322 
4323  varIdx = getBranchingVar(thetaVar, numThetaVar );
4324  std::cout << "VARIABLE INDEX = " << varIdx << std::endl;
4325  std::cout << "Branching on Variable: " << m_variableNames[ varIdx] << std::endl;
4326 
4327  //if this variable is in the map, then we just return with the index,
4328  //if this variable is NOT in the map then we add a cut
4329 
4330  if( varConMap.find( varIdx) == varConMap.end() ){
4331 
4332  for(i = 0; i < m_numThetaVar; i++){
4333 
4334  kount = 0;
4335 
4336  for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
4337 
4338  if ( m_thetaIndex[ j] == varIdx) kount++ ;
4339 
4340  }
4341 
4342  //count is the number times variable i appears in the constraint
4343 
4344  if(kount > 0){
4345 
4346  branchCutIndexes[ numNonz] = i;
4347  branchCutValues[ numNonz++] = kount ;
4348 
4349  }
4350 
4351  }
4352 
4353 
4354  //add varIdx cut to B matrix
4356  m_BmatrixIdx[ m_numBmatrixNonz++] = varIdx;
4357  m_numBmatrixCon++;
4359 
4360  //make sure to add artificial variables
4361  //of course they have no nonzero elements in
4362  //the transformation matrix
4363  m_numThetaVar++;
4366  //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz; //first artificial variable
4367  //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz; //second artificial variable
4368 
4369 
4370  }//end of if on checking for map membership
4371 
4372  //set return arguments
4373 
4374  indexes = branchCutIndexes;
4375  values = branchCutValues;
4376 
4377  return;
4378 
4379  }catch (const ErrorClass& eclass) {
4380 
4381  throw ErrorClass(eclass.errormsg);
4382 
4383  }
4384 
4385 }//end getBranchingCut dense
4386 
4387 
4388 void OSBearcatSolverXij::getBranchingCut(const int* thetaIdx, const double* thetaVar,
4389  const int numThetaVar, const std::map<int, int> &varConMap,
4390  int &varIdx, int &numNonz, int* &indexes, double* &values) {
4391 
4392  //get a branching variable
4393  int i;
4394  int j;
4395  int kount;
4396  numNonz = 0;
4397  //keep numNonz at zero if there is no cut
4398  //there will be no cut if the xij is in conVarMap
4399 
4400  try{
4401 
4402 
4403 
4404  varIdx = getBranchingVar(thetaIdx, thetaVar, numThetaVar );
4405 
4406  std::cout << "Branching on Variable: " << m_variableNames[ varIdx] << std::endl;
4407 
4408  //if this variable is in the map, then we just return with the index,
4409  //if this variable is NOT in the map then we add a cut
4410 
4411  if( varConMap.find( varIdx) == varConMap.end() ){
4412 
4413 
4414 
4415 
4416 
4417 
4418  for(i = 0; i < m_numThetaVar; i++){
4419 
4420  kount = 0;
4421 
4422  for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
4423 
4424  if ( m_thetaIndex[ j] == varIdx) kount++ ;
4425 
4426  }
4427 
4428  //count is the number times variable i appears in the constraint
4429 
4430  if(kount > 0){
4431 
4432  branchCutIndexes[ numNonz] = i;
4433  branchCutValues[ numNonz++] = kount ;
4434 
4435  }
4436 
4437  }
4438 
4439 
4440  //add varIdx cut to B matrix
4442  m_BmatrixIdx[ m_numBmatrixNonz++] = varIdx;
4443  m_numBmatrixCon++;
4445 
4446  //make sure to add artificial variables
4447  //of course they have no nonzero elements in
4448  //the transformation matrix
4449  m_numThetaVar++;
4452  //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz; //first artificial variable
4453  //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz; // second artificial variable
4454 
4455 
4456  }//end of if on checking for map membership
4457 
4458  //set return arguments
4459 
4460  indexes = branchCutIndexes;
4461  values = branchCutValues;
4462 
4463  return;
4464 
4465  }catch (const ErrorClass& eclass) {
4466 
4467  throw ErrorClass(eclass.errormsg);
4468 
4469  }
4470 
4471 }//end getBranchingCut sparse
4472 
4473 
4475 
4476  int k;
4477  double *routeDemand;
4478  routeDemand = NULL;
4479  routeDemand = new double[ m_numHubs];
4480  std::map<int, std::vector<int> > routeMap;
4481  std::vector<int>::iterator vit;
4482  bool isFeasible;
4483 
4484 
4485  isFeasible = true;
4486  try{
4487  //see if we have a feasible solution
4488  routeMap = m_initSolMap[ 0];
4489 
4490  for(k = 0; k < m_numHubs; k++){
4491 
4492  routeDemand[ k] = 0;
4493  for(vit = routeMap[k].begin(); vit!=routeMap[k].end(); ++vit){
4494 
4495 
4496  routeDemand[ k] += m_demand[ *vit];
4497  }
4498 
4499  if(routeDemand[k] > m_routeCapacity[ k] ){
4500 
4501  isFeasible = false;
4502  break;
4503 
4504  }
4505  }
4506 
4507  delete[] routeDemand;
4508  routeDemand = NULL;
4509 
4510  if(isFeasible == false) getFeasibleSolution();
4511  //call the OneOPT heuristic;
4512  OneOPT();
4513 
4514 
4515  }catch (const ErrorClass& eclass) {
4516 
4517  if(routeDemand == NULL){
4518 
4519  delete[] routeDemand;
4520  routeDemand = NULL;
4521  }
4522 
4523  throw ErrorClass(eclass.errormsg);
4524 
4525  }
4526 
4527 
4528 }//end getInitialSolution
4529 
4530 
4531 void OSBearcatSolverXij::resetMaster( std::map<int, int> &inVars, OsiSolverInterface *si){
4532 
4533  int i;
4534  int j;
4535 
4536  int kount;
4537  int numNonz;
4538 
4539  std::map<int, int>::iterator mit;
4540  //temporarily create memory for the columns we keep
4541  int numVars = inVars.size();
4542  int numVarArt;
4543  //there 2*m_numNodes in the A matrix
4544  //there are m_numBmatrixCon B matrix constraints
4545  //numVarArt = 2*m_numNodes + 2*m_numBmatrixCon;
4546  numVarArt = m_numNodes + m_numBmatrixCon;
4547 
4548  //arrays for the new osinstance
4549  std::vector<double> valuesVec;
4550  double *values = NULL;
4551 
4552  std::vector<int> indexesVec;
4553  int *indexes = NULL ;
4554 
4555  int *starts = new int[ numVars + 1 + numVarArt];
4556 
4557  int startsIdx;
4558 
4559 
4560 
4561  //temporay holders
4562  int* thetaPntTmp;
4563  int* thetaIndexTmp;
4564  int* tmpConvexity = new int[ m_numThetaVar];
4565 
4566  //get the number of nonzeros that we need
4567  numNonz = 0;
4568 
4569  for(mit = inVars.begin(); mit != inVars.end(); mit++){
4570 
4571  numNonz += m_thetaPnt[mit->first + 1 ] - m_thetaPnt[ mit->first ];
4572  }
4573 
4574  //allocate the memory
4575  thetaPntTmp = new int[ numVars + 1];
4576  thetaIndexTmp = new int[ numNonz];
4577 
4578 
4579  //error check
4580  for(mit = inVars.begin(); mit != inVars.end(); mit++){
4581 
4582 
4583  //std::cout << "VARIABLE INDEX = " << mit->first << " OBJ COEF = " << si->getObjCoefficients()[ mit->first ] << std::endl;
4584  if( m_convexityRowIndex[ mit->first] == -1) throw ErrorClass( "we have an artificial variable in reset master");
4585 
4586 
4587  }
4588 
4589  //fill in the temporary arrays
4590  kount = 0;
4591  numNonz = 0;
4592  thetaPntTmp[ kount] = 0;
4593 
4594  for(mit = inVars.begin(); mit != inVars.end(); mit++){
4595 
4596  //std::cout << "mit->first " << mit->first << " mit->second " << mit->second << std::endl;
4597 
4598  kount++;
4599 
4600  for(i = m_thetaPnt[ mit->first ]; i < m_thetaPnt[mit->first + 1 ]; i++){
4601 
4602  thetaIndexTmp[ numNonz++] = m_thetaIndex[ i];
4603 
4604  //std::cout << "Column = " << mit->first << " Variable " << m_variableNames[ m_thetaIndex[ i] ] << std::endl;
4605 
4606  }
4607 
4608  thetaPntTmp[ kount] = numNonz;
4609 
4610  //std::cout << "kount = " << kount << " thetaPntTmp[ kount] = " << thetaPntTmp[ kount] << std::endl;
4611  //readjust numbering to take into account artificial variables
4612  //mit->second += numVarArt;
4613  //kipp -- double check calculation below
4614  inVars[ mit->first] = numVarArt + kount - 1 ;
4615 
4616  }
4617 
4618  //std::cout << "kount = " << kount << std::endl;
4619  //std::cout << "numVars = " << numVars << std::endl;
4620 
4621 
4622 
4623  //copy old values of m_convexityRowIndex
4624  for(i = 0; i < m_numThetaVar; i++) tmpConvexity[ i] = m_convexityRowIndex[ i];
4625 
4626  //reset the theta pointers
4627  //first the artificial variables
4628  m_numThetaVar = 0;
4629  m_numThetaNonz = 0;
4630  for(i = 0; i < numVarArt; i++){
4631 
4633  m_thetaPnt[ m_numThetaVar++] = 0;
4634 
4635 
4636  }
4637  //now fill in the other pointers from the temp arrarys
4638  //std::cout << "Number of artificial variables = " << numVarArt << std::endl;
4639  intVarSet.clear();
4640  for(mit = inVars.begin(); mit != inVars.end(); mit++){
4641 
4642 
4643  intVarSet.insert ( std::pair<int,double>(mit->second, 1.0) );
4644 
4645  //std::cout << " m_numThetaVar = " << m_numThetaVar << " m_numThetaNonz = " << m_numThetaNonz << std::endl;
4646  //std::cout << "Variable number " << mit->first << " OBJ coefficient = " << si->getObjCoefficients()[ mit->first] << std::endl;
4647 
4648  m_convexityRowIndex[ m_numThetaVar] = tmpConvexity[ mit->first];
4649 
4651 
4652  for(j = thetaPntTmp[ mit->second - numVarArt]; j < thetaPntTmp[ mit->second - numVarArt + 1 ]; j++){
4653 
4654 
4655  m_thetaIndex[ m_numThetaNonz ] = thetaIndexTmp[ m_numThetaNonz] ;
4656  m_numThetaNonz++;
4657 
4658  }
4659 
4660  }
4661 
4663  //std::cout << " number of art vars = " << numVarArt << std::endl;
4664  //std::cout << " m_numThetaVar = " << m_numThetaVar << std::endl;
4665  //std::cout << " m_numThetaNonz = " << m_numThetaNonz << std::endl;
4666  //done with the transformation matrix
4667 
4668 
4669 
4670  //
4671  //write old master --- just for testing
4672  si->writeLp( "gailTest" );
4673 
4674  //now create the formulation
4675 
4676  //first get each column of the new master
4677  //first take care of the artificial variables
4678  numNonz = 0;
4679  startsIdx = 0;
4680  starts[ startsIdx++] = numNonz;
4681 
4682  for(i = 0; i < numVarArt; i++){
4683  numNonz++;
4684  starts[ startsIdx++] = numNonz;
4685  valuesVec.push_back( 1.0);
4686  indexesVec.push_back( i);
4687 
4688  }
4689 
4690 
4691  int rowCount;
4692 
4693  int numAmatrixRows;
4694  numAmatrixRows = m_numNodes - m_numHubs;
4695 
4696  for(mit = inVars.begin(); mit != inVars.end(); mit++){
4697 
4698  //std::cout << "CONVEXITY ROW = " << m_convexityRowIndex[ mit->second] << std::endl;
4699  valuesVec.push_back( 1.0);
4700  indexesVec.push_back( numAmatrixRows + m_convexityRowIndex[ mit->second] );
4701  //increment numNonz by 1 for the convexity row
4702  numNonz++;
4703 
4704  for(j = m_thetaPnt[ mit->second ]; j < m_thetaPnt[ mit->second + 1 ]; j++){
4705 
4707 
4708  //std::cout << "Column = " << mit->second << " Variable " << m_variableNames[ m_thetaIndex[ j] ] << std::endl;
4709 
4710  }
4711 
4712 
4713 
4714  //multiply the sparse array by each A matrix constraint
4715  for(i = 0; i < numAmatrixRows; i++){
4716 
4717  rowCount = 0;
4718 
4719  for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
4720 
4721  //m_Amatrix[ j] is a variable index -- this logic works
4722  //since the Amatrix coefficient is 1 -- we don't need a value
4723  //it indexes variable that points into the node
4724  rowCount += m_tmpScatterArray[ m_Amatrix[ j] ];
4725 
4726 
4727  }
4728 
4729  if(rowCount > 0){
4730 
4731  numNonz++;
4732 
4733  //std::cout << "Column = " << mit->second << " Nonzero in A marix row " << i << " Value = " << rowCount << std::endl;
4734  valuesVec.push_back( rowCount);
4735  indexesVec.push_back( i);
4736 
4737 
4738  }
4739 
4740 
4741  }//end loop on coupling constraints
4742 
4743 
4744  //multiply the sparse array by each B matrix constraint
4745  for(i = 0; i < m_numBmatrixCon; i++){
4746 
4747  rowCount = 0;
4748 
4749  for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
4750 
4751  //m_Amatrix[ j] is a variable index -- this logic works
4752  //since the Amatrix coefficient is 1 -- we don't need a value
4753  //it indexes variable that points into the node
4754  rowCount += m_tmpScatterArray[ m_BmatrixIdx[ j] ];
4755 
4756 
4757  }
4758 
4759  if(rowCount > 0){
4760  numNonz++;
4761 
4762  //std::cout << "Column = " << mit->first << " Nonzero in B matrix row " << i + m_numNodes<< " Value = " << rowCount << std::endl;
4763 
4764  valuesVec.push_back( rowCount);
4765  indexesVec.push_back( i + m_numNodes);
4766  }
4767 
4768 
4769  }//end loop on B matrix constraints
4770 
4771 
4772  //zero out the scatter array
4773 
4774  for(j = m_thetaPnt[ mit->second ]; j < m_thetaPnt[ mit->second + 1 ]; j++){
4775 
4776  m_tmpScatterArray[ m_thetaIndex[ j] ] = 0;
4777 
4778  }
4779 
4780  starts[ startsIdx++] = numNonz;
4781 
4782  }
4783 
4784 
4785  //for(i = 0; i < startsIdx; i++) std::cout << "starts[ i] = " << starts[ i] << std::endl;
4786  values = new double[ numNonz];
4787  indexes = new int[ numNonz];
4788 
4789  if( (unsigned int)numNonz != valuesVec.size() ) throw ErrorClass("dimension problem in reset");
4790  if( (unsigned int)numNonz != indexesVec.size() ) throw ErrorClass("dimension problem in reset");
4791 
4792  for(i = 0; i < numNonz; i++){
4793 
4794  values[ i] = valuesVec[i];
4795  indexes[ i] = indexesVec[i];
4796 
4797  }
4798 
4799 
4800 
4801  //construct the new master
4802  //create an OSInstance from the tmp arrays
4803  // delete the old m_osinstanceMaster
4804 
4805  delete m_osinstanceMaster;
4806  m_osinstanceMaster = NULL;
4807 
4808  //start building the restricted master here
4810  m_osinstanceMaster->setInstanceDescription("The Restricted Master");
4811 
4812  // first the variables
4813  // we have numVarArt] artificial variables
4814  m_osinstanceMaster->setVariableNumber( numVars + numVarArt );
4815 
4816  // now add the objective function
4817  //m_osinstanceMaster->setObjectiveNumber( 1);
4818  //add m_numNodes artificial variables
4819  SparseVector *objcoeff = new SparseVector( numVars + numVarArt);
4820 
4821 
4822 
4823  // now the constraints
4825 
4826 
4827  //add the artificial variables first
4828 
4829  int varNumber;
4830  varNumber = 0;
4831 
4832 
4833  //define the artificial variables
4834  for(i = 0; i < numVarArt; i++){
4835 
4836  objcoeff->indexes[ varNumber ] = varNumber ;
4837 
4838  objcoeff->values[ varNumber ] = m_osDecompParam.artVarCoeff;
4839 
4840  m_thetaCost[ varNumber] = m_osDecompParam.artVarCoeff;
4841 
4842  m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt("x", i ) ,
4843  0, 1.0, 'C');
4844 
4845 
4846  }
4847 
4848  // now the theta variables
4849  kount = 0;
4850  for(mit = inVars.begin(); mit != inVars.end(); mit++){
4851 
4852  objcoeff->indexes[ varNumber ] = varNumber ;
4853 
4854  objcoeff->values[ varNumber ] = si->getObjCoefficients()[ mit->first] ;
4855 
4856  m_thetaCost[ varNumber] = si->getObjCoefficients()[ mit->first];
4857 
4858  m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt("x", kount + numVarArt ) ,
4859  0, 1.0, 'C');
4860 
4861  kount++;
4862 
4863 
4864 
4865  }
4866 
4867 
4868 
4869  for(i = 0; i < m_numNodes; i++){
4870 
4872  1.0, 1.0, 0);
4873 
4874  }
4875 
4876 
4877  for(i = m_numNodes; i < m_numBmatrixCon + m_numNodes; i++){
4878 
4880  si->getRowLower()[ i], si->getRowUpper()[ i], 0);
4881 
4882 
4883  }
4884 
4885 
4886  // now add the objective function
4888  m_osinstanceMaster->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
4889 
4890  //add the linear constraints coefficients
4892  values, 0, numNonz - 1, indexes, 0, numNonz - 1, starts, 0, startsIdx);
4893 
4894 
4895  //std::cout << m_osinstanceMaster->printModel( ) << std::endl;
4896 
4897  //garbage collection
4898  delete[] tmpConvexity;
4899  tmpConvexity = NULL;
4900  delete[] thetaPntTmp;
4901  thetaPntTmp = NULL;
4902  delete[] thetaIndexTmp;
4903  thetaIndexTmp = NULL;
4904  delete objcoeff;
4905  objcoeff = NULL;
4906 }//end resetMaster
4907 
4908 
4909 
4910 double OSBearcatSolverXij::getRouteDistance(int numNodes, int hubIndex,
4911  CoinSolver* solver, std::vector<int> zk, double* xVar){
4912 
4913  int i;
4914  int j;
4915  int numVar;
4916  numVar = numNodes*numNodes - numNodes;
4917 
4918  int kount;
4919  double objValue;
4920  std::vector<int>::iterator vit;
4921  std::map<std::pair<int, int>, int > cutMap;
4922  std::map<std::pair<int, int>, int >::iterator mit;
4923  std::vector<IndexValuePair*> primalValPair;
4924  //
4925  //stuff for cut generation
4926  //
4927 
4928  //zero this array out
4929  for(i = 0; i < numVar; i++) xVar[ i] = 0;
4930  //getRows function call return parameters
4931  int numNewRows;
4932  int* numRowNonz = NULL;
4933  int** colIdx = NULL;
4934  double** rowValues = NULL ;
4935  double* rowLB;
4936  double* rowUB;
4937  //end of getRows function call return parameters
4938  //
4939  //end of cut generation stuff
4940  OsiSolverInterface *si = solver->osiSolver;
4941  try{
4942 
4943  //adjust the RHS of the hubIndex constraints
4944  si->setRowUpper(hubIndex, 1.0);
4945  si->setRowUpper(hubIndex + numNodes, 1.0);
4946  si->setRowLower(hubIndex, 1.0);
4947  si->setRowLower(hubIndex + numNodes, 1.0);
4948  //
4949  //adjust the lower bounds on the variables
4950 
4951  kount = zk.size();
4952 
4953  int idx1;
4954  int idx2;
4955 
4956  for(i = 0; i < kount ; i++){
4957 
4958  //adjust the RHS of the hubIndex constraints
4959  si->setRowUpper(zk[ i], 1.0);
4960  si->setRowUpper(zk[ i] + numNodes, 1.0);
4961  si->setRowLower(zk[ i], 1.0);
4962  si->setRowLower(zk[ i] + numNodes, 1.0);
4963 
4964  idx1 = hubIndex;
4965  idx2 = zk[i ];
4966 
4967  //get variable number for (idx1, idx2)
4968  if( idx1 > idx2 ){
4969  si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
4970  }else{
4971 
4972  if( idx1 < idx2 ){
4973 
4974  si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
4975 
4976  }
4977  }
4978 
4979  idx1 = zk[ i];
4980  idx2 = hubIndex;
4981 
4982  //get variable number for (idx1, idx2)
4983  if( idx1 > idx2 ){
4984 
4985  si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
4986 
4987  }else{
4988 
4989  if( idx1 < idx2 ){
4990 
4991  si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
4992 
4993  }
4994  }
4995 
4996  for(j = i + 1; j < kount; j++){
4997 
4998  idx1 = zk[i];
4999  idx2 = zk[j];
5000 
5001  //get variable number for (idx1, idx2)
5002  if( idx1 > idx2 ){
5003 
5004  si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
5005 
5006  }else{
5007 
5008  if( idx1 < idx2 ){
5009 
5010  si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
5011  }
5012  }
5013 
5014 
5015  idx1 = zk[j];
5016  idx2 = zk[i];
5017 
5018  //get variable number for (idx1, idx2)
5019  if( idx1 > idx2 ){
5020  si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
5021  }else{
5022 
5023  if( idx1 < idx2 ){
5024 
5025  si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
5026  }
5027  }
5028  }
5029  }
5030 
5031  solver->solve();
5032 
5033  //throw an exception if we don't get an optimal solution
5034 
5035  if(solver->osresult->getSolutionStatusType( 0 ) != "optimal" ) throw ErrorClass("Trouble with integer program used to construct initial master");
5036  objValue = solver->osresult->getObjValue(0, 0) ;
5037  // now get the X values -- use these to get a cut
5038  primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
5039 
5040  //loop over routes again to create master objective coefficients
5041  bool isCutAdded;
5042  isCutAdded = true;
5043  while(isCutAdded == true){
5044 
5045  isCutAdded = false;
5046 
5047  for(i = 0; i < numVar; i++) xVar[ primalValPair[ i]->idx ] = primalValPair[ i]->value;
5048 
5049  numNewRows = 0;
5050  getCutsX(xVar, numVar, numNewRows, numRowNonz,
5051  colIdx,rowValues, rowLB, rowUB);
5052 
5053  if(numNewRows >= 1){
5054  isCutAdded = true;
5055  std::cout << "WE HAVE A CUT " << std::endl;
5056  std::cout << "EXPRESS CUT IN X(I, J) SPACE" << std::endl;
5057  //for(i = 0; i < numRowNonz[ 0]; i++) std::cout << m_variableNames[ colIdx[0][ i] ] << std::endl;
5058 
5059  for(i = 0; i < numNewRows; i++){
5060 
5061  std::cout << "CUT UPPER BOUND = " << rowUB[ i] << std::endl;
5062  si->addRow(numRowNonz[ i], colIdx[ i], rowValues[ i], rowLB[ i], rowUB[ i] ) ;
5063  }
5064 
5065  std::cout << "A CUT WAS ADDED, CALL SOLVE AGAIN" << std::endl;
5066  solver->solve();
5067  if(solver->osresult->getSolutionStatusType( 0 ) != "optimal" ) throw ErrorClass("Trouble with integer program used to construct initial master");
5068  primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
5069  std::cout << "New Solution Status = " << solver->osresult->getSolutionStatusType( 0 ) << std::endl;
5070  std::cout << "Optimal Objective Value = " << solver->osresult->getObjValue(0, 0) << std::endl;
5071 
5072  // zero out xVar
5073  //for(i = 0; i < numVar; i++) xVar[ i ] = 0;
5074 
5075  }//end if on numNewRows >= 1
5076 
5077  }//end while on isCutAdded
5078 
5079  objValue = solver->osresult->getObjValue(0, 0) ;
5080  // now get the X values -- use these to get a cut
5081  primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
5082 
5083  //for(i = 0; i < numVar; i++)
5084  // if( primalValPair[ i]->value > .5) std::cout << m_variableNames[ primalValPair[ i]->idx ] << std::endl;
5085  std::cout << "Objective Function Value = " << objValue << std::endl;
5086 
5087  //reset the upper bounds
5088  for(i = 0; i < numVar; i++) si->setColUpper(i, 0);
5089  for(i = 0; i < 2*numNodes; i++) si->setRowUpper(i, 0);
5090  for(i = 0; i < 2*numNodes; i++) si->setRowLower(i, 0);
5091 
5092  return objValue;
5093 
5094 
5095  } catch (const ErrorClass& eclass) {
5096 
5097  std::cout << std::endl << std::endl << std::endl;
5098 
5099  if( xVar != NULL)
5100  delete xVar;
5101  // Problem with the parser
5102  throw ErrorClass(eclass.errormsg);
5103  }
5104 
5105 }//end getRouteDistance
5106 
5107 
5108 
5109 CoinSolver* OSBearcatSolverXij::getTSP(int numNodes, double* cost){
5110 
5111 
5112  int i;
5113  int j;
5114  int numVar;
5115  numVar = numNodes*numNodes - numNodes;
5116  int numNonz;
5117  int kount;
5118 
5119  std::vector<int>::iterator vit;
5120  double* rhsVec;
5121  CoinSolver *solver = NULL;
5122  std::map<std::pair<int, int>, int > cutMap;
5123  std::map<std::pair<int, int>, int >::iterator mit;
5124 
5125  //numCuts is the number of cuts of the form Xij + Xji <= 1
5126  int numCuts;
5127 
5128  rhsVec = new double[ 2*numNodes];
5129 
5130  for(i = 0; i < 2*numNodes; i++){
5131 
5132  //this will then get changed to 1
5133  //when we know assignments
5134  rhsVec[ i] = 0;
5135  }
5136 
5137 
5138 
5139  //now count and get variables in Xij + Xji <=1 cut
5140  numCuts = 0;
5141  std::pair <int,int> pairIJ;
5142  std::pair <int,int> pairJI;
5143 
5144  for(i = 0; i < numNodes - 1; i++){
5145 
5146 
5147  for(j = i + 1; j < numNodes; j++){
5148 
5149 
5150  pairIJ = std::make_pair(i, j);
5151  pairJI = std::make_pair(j, i);
5152 
5153  cutMap[pairIJ ] = 2*numNodes + numCuts;
5154  cutMap[pairJI ] = 2*numNodes + numCuts;
5155  numCuts++;
5156 
5157  }
5158  }
5159 
5160 
5161 
5162  OSInstance *osinstance = new OSInstance();
5163  try{
5164 
5165  osinstance->setInstanceSource("generated from OSBearcatSoleverXij");
5166  osinstance->setInstanceDescription("Find the minimum tour for a route");
5167  osinstance->setVariableNumber( numVar);
5168 
5169  for(i = 0; i < numVar; i++){
5170 
5171  osinstance->addVariable(i, m_variableNames[ i], 0, 0, 'B');
5172 
5173  }
5174  // now add the objective function
5175  osinstance->setObjectiveNumber( 1);
5176 
5177  // now the coefficient
5178  SparseVector *objcoeff = new SparseVector( numVar);
5179 
5180  for(i = 0; i < numVar; i++){
5181 
5182  objcoeff->indexes[ i] = i;
5183  objcoeff->values[ i] = cost[ i];
5184 
5185  }
5186 
5187  osinstance->addObjective(-1, "minimizeDistance", "min", 0.0, 1.0, objcoeff);
5188  objcoeff->bDeleteArrays = true;
5189  delete objcoeff;
5190 
5191  osinstance->setConstraintNumber( 2*numNodes + numCuts);
5192  //bool addConstraint(int index, string name, double lowerBound, double upperBound, double constant);
5193 
5194  for(i = 0; i < numNodes; i++){
5195 
5196  osinstance->addConstraint(i, makeStringFromInt("enter_node_", i) , rhsVec[i], rhsVec[i], 0);
5197 
5198  }
5199 
5200  for(i = numNodes; i < 2*numNodes; i++){
5201 
5202  osinstance->addConstraint(i, makeStringFromInt("leave_node_", i - numNodes) , rhsVec[i], rhsVec[i], 0);
5203 
5204  }
5205 
5206  //now the Xij cuts
5207 
5208  for(i = 0; i < numCuts; i++){
5209 
5210  osinstance->addConstraint(2*numNodes + i, makeStringFromInt("XijCut_", i) , 0, 1, 0);
5211 
5212  }
5213 
5214  //now the A matrix
5215  //note: each Xij + Xji <= 1 has two nonzero elements
5216  numNonz = numVar*numVar - numVar + 2*numCuts;
5217 
5218 
5219  double *values = new double[ numNonz];
5220  int *rowIndexes = new int[ numNonz];
5221  int *starts = new int[ numVar + 1];
5222 
5223 
5224  kount = 0;
5225  numNonz = 0;
5226  starts[ kount++] = 0;
5227  for(i = 0; i < numNodes; i++){
5228 
5229  for(j = 0; j < i; j++){
5230  //we have xij, j < i
5231 
5232  rowIndexes[ numNonz] = j; //enter node j
5233  values[ numNonz++] = 1.0;
5234 
5235  rowIndexes[ numNonz] = numNodes + i; //leave node i
5236  values[ numNonz++] = 1.0;
5237 
5238 
5239  pairIJ = std::make_pair(i, j);
5240  mit = cutMap.find( pairIJ);
5241  if(mit != cutMap.end() ){
5242 
5243  rowIndexes[ numNonz] = mit->second;
5244  values[ numNonz++] = 1.0;
5245  }
5246 
5247 
5248  starts[ kount++] = numNonz;
5249 
5250 
5251  }
5252 
5253  for(j = i+1; j < numNodes; j++){
5254  //we have xij, j > i
5255 
5256  rowIndexes[ numNonz] = j; //enter node j
5257  values[ numNonz++] = 1.0;
5258 
5259  rowIndexes[ numNonz] = numNodes + i; //leave node i
5260  values[ numNonz++] = 1.0;
5261 
5262 
5263  pairIJ = std::make_pair(i, j);
5264  mit = cutMap.find( pairIJ);
5265  if(mit != cutMap.end() ){
5266 
5267  rowIndexes[ numNonz] = mit->second;
5268  values[ numNonz++] = 1.0;
5269  }
5270 
5271 
5272  starts[ kount++] = numNonz;
5273 
5274  }
5275 
5276 
5277  }
5278 
5279  osinstance->setLinearConstraintCoefficients(numNonz, true, values, 0, numNonz - 1,
5280  rowIndexes, 0, numNonz - 1, starts, 0, numVar);
5281 
5282  //std::cout << osinstance->printModel() << std::endl;
5283 
5284 
5285  solver = new CoinSolver();
5286  solver->sSolverName ="cbc";
5287  solver->osinstance = osinstance;
5288  solver->osoption = m_osoption;
5289  solver->buildSolverInstance();
5290  solver->osoption = m_osoption;
5291 
5292  delete[] rhsVec;
5293  rhsVec = NULL;
5294 
5295  return solver;
5296 
5297 
5298 
5299  } catch (const ErrorClass& eclass) {
5300 
5301  std::cout << std::endl << std::endl << std::endl;
5302 
5303  if(rhsVec != NULL){
5304  delete[] rhsVec;
5305  rhsVec = NULL;
5306  }
5307 
5308  // Problem with the parser
5309  throw ErrorClass(eclass.errormsg);
5310  }
5311 
5312 
5313 }//end getTSP
5314 
5315 
5317 
5318  int k;
5319  int kprime;
5320  double *routeCost = NULL;
5321  int *routeDemand = NULL;
5322  double *xVar = NULL;
5323  int numXVar = m_numNodes*m_numNodes - m_numNodes;
5324  double routeCostInf;
5325 
5326  routeCostInf = OSINT_MAX;
5327  double feasMult = 10000; //factor by which we multiply feasibility improvement
5328 
5329  routeCost = new double[m_numHubs];
5330  routeDemand = new int[m_numHubs];
5331  xVar = new double[ numXVar];
5332 
5333  //get current capacities
5334 
5335  std::map<int, std::vector<int> > routeMap;
5336  std::vector<int> tmpVec;
5337  std::vector<int>::iterator vit;
5338  std::vector<int>::iterator vit2;
5339 
5340  routeMap = m_initSolMap[ 0];
5341  CoinSolver *solver = NULL;
5342 
5343  double totalCost;
5344  bool foundMoBetta; //set to true if improved
5345  bool foundLocalImprovement;
5346 
5347  try{
5348 
5349  solver = getTSP(m_numNodes, m_cost);
5350 
5351  totalCost = 0;
5352  for(k = 0; k < m_numHubs; k++){
5353 
5354  routeDemand[ k] = 0;
5355  for(vit = routeMap[k].begin(); vit!=routeMap[k].end(); ++vit){
5356 
5357  std::cout << "node i = " << *vit << " demand " << m_demand[ *vit] << std::endl;
5358  routeDemand[ k] += m_demand[ *vit];
5359  }
5360 
5361  if(routeDemand[k] <= m_routeCapacity[ k] ){
5362 
5363  routeCost[ k] = getRouteDistance(m_numNodes, k, solver,
5364  routeMap[k], xVar)*routeDemand[ k];
5365 
5366 
5367  std::cout << "rout = " << k << " cost " << routeCost[ k] << std::endl;
5368  totalCost += routeCost[ k];
5369  }else{
5370  std::cout << "rout demand " << routeDemand[k] << " route capacity = " << m_routeCapacity[ k] << std::endl;
5371  routeCost[ k] = routeCostInf;
5372  totalCost += routeCost[ k];
5373 
5374  }
5375 
5376 
5377  }
5378 
5379  //now loop as long as there is improvement
5380  foundMoBetta = true;
5381  double improvement;
5382 
5383  double tmpCostk;
5384  double tmpCostkPrime;
5385  double optCostk;
5386  double optCostkPrime;
5387  double tmpVal;
5388  int tmpIdx;
5389  std::vector<int>::iterator vitIdx;
5390 
5391  while(foundMoBetta == true){
5392  foundMoBetta = false;
5393 
5394  for(k = 0; k < m_numHubs; k++){
5395 
5396  foundLocalImprovement = false;
5397 
5398  for(kprime = 0; kprime < m_numHubs; kprime++){
5399 
5400  if(kprime != k && routeCost[ kprime] < routeCostInf){
5401 
5402  //try to move a node from route k to route kprime
5403  improvement = 0;
5404  optCostk = routeCostInf;
5405  optCostkPrime = routeCostInf;
5406 
5407  for(vit = routeMap[k].begin(); vit!=routeMap[k].end(); ++vit){//looping over nodes
5408 
5409  foundLocalImprovement = false;
5410 
5411  //consider making a switch if feasible
5412  if( m_demand[ *vit] + routeDemand[ kprime] <= m_routeCapacity[ kprime] ){
5413 
5414  tmpCostk = routeCostInf;
5415  tmpCostkPrime = routeCostInf;
5416 
5417 
5418  //make a switch
5419  //add *vit to route kprime and take away from k
5420  routeMap[ kprime].push_back( *vit);
5421  //calculate cost for route kprime
5422  tmpCostkPrime = getRouteDistance(m_numNodes, kprime, solver,
5423  routeMap[kprime], xVar)*(routeDemand[ kprime] + m_demand[ *vit] );
5424 
5425  //now switch back
5426  routeMap[ kprime].pop_back( );
5427 
5428 
5429  //kipp -- handle both infinite and finite
5430  if(routeCost[k] == routeCostInf ){
5431 
5432  //std::cout << "WE HAVE INFINITY CASE" << std::endl;
5433 
5434 
5435  //don't bother to solve TSP for route k if we are still infinite
5436  if( routeDemand[ k] - m_demand[ *vit] <= m_routeCapacity[ k]) {
5437 
5438  for(vit2 = routeMap[k].begin(); vit2 != routeMap[k].end(); vit2++){
5439 
5440  if(vit != vit2) tmpVec.push_back( *vit2);
5441 
5442  }
5443 
5444  tmpCostk = getRouteDistance(m_numNodes, k, solver,
5445  tmpVec, xVar)*(routeDemand[ k] - m_demand[ *vit] );
5446 
5447  tmpVec.clear();
5448 
5449  }
5450 
5451  tmpVal = std::max(*vit, routeDemand[ k] - m_demand[ *vit])*feasMult
5452  + ( routeCost[kprime] - tmpCostkPrime);
5453 
5454  if( tmpVal > improvement) {
5455  foundLocalImprovement = true;
5456  improvement = tmpVal;
5457  tmpIdx = *vit;
5458  vitIdx = vit;
5459  optCostk = tmpCostk;
5460  optCostkPrime = tmpCostkPrime;
5461 
5462  }
5463 
5464 
5465  }else{// not infinite cost
5466 
5467  //std::cout << "WE HAVE FINITE CASE" << std::endl;
5468 
5469  //evaluate cost on route k
5470  for(vit2 = routeMap[k].begin(); vit2 != routeMap[k].end(); vit2++){
5471 
5472  if(vit != vit2) tmpVec.push_back( *vit2);
5473 
5474  }
5475 
5476  tmpCostk = getRouteDistance(m_numNodes, k, solver,
5477  tmpVec, xVar)*(routeDemand[ k] - m_demand[ *vit] );
5478 
5479  tmpVec.clear();
5480 
5481  if( ( (routeCost[k] + routeCost[kprime]) -
5482  ( tmpCostkPrime + tmpCostk ) ) > improvement ) {
5483 
5484  improvement = (routeCost[k] + routeCost[kprime]) -
5485  ( tmpCostkPrime + tmpCostk );
5486 
5487  foundLocalImprovement = true;
5488  tmpIdx = *vit;
5489  vitIdx = vit;
5490  optCostk = tmpCostk;
5491  optCostkPrime = tmpCostkPrime;
5492 
5493 
5494 
5495  }
5496  }
5497  }//end if on capacity test
5498 
5499  if(foundLocalImprovement == true)break;
5500  }// looping over nodes in route k
5501 
5502  //do updates here if we found an improvement
5503 
5504  //make switch on best found
5505  if( foundLocalImprovement == true) {
5506 
5507 
5508  std::cout << "index = " << tmpIdx << std::endl;
5509  //add tmpIdx to route kPrime
5510  routeMap[ kprime].push_back( tmpIdx);
5511  //adjust route demand
5512  routeDemand[ kprime] += m_demand[ tmpIdx];
5513  //adjust route cost
5514  routeCost[ kprime] = optCostkPrime;
5515 
5516  //std::cout << "kprime route cost = " << routeCost[ kprime] << std::endl;
5517  //std::cout << "kprime route demand = " << routeDemand[ kprime] << std::endl;
5518 
5519 
5520  //delete tmpIdx to route kPrime
5521  std::cout << "tmpIdx = " << tmpIdx << std::endl;
5522  std::cout << "vitIdx = " << *vitIdx << std::endl;
5523  routeMap[ k].erase( vitIdx);
5524  //adjust rouet demand
5525  routeDemand[ k] -= m_demand[ tmpIdx];
5526  //adjust route cost
5527  routeCost[ k] = optCostk;
5528 
5529  //std::cout << "k route cost = " << routeCost[ k] << std::endl;
5530  //std::cout << "k route demand = " << routeDemand[ k] << std::endl;
5531 
5532 
5533 
5534  foundMoBetta = true;
5535 
5536  }//if on OSDBL_MAX
5537  }//close if on k != kprime
5538  }//loop on kprime
5539  }//loop on k
5540 
5541 
5542  }//loop on while
5543 
5544 
5545  //summarize
5546  totalCost = 0;
5547  for(k = 0; k < m_numHubs; k++){
5548 
5549  std::cout << std::endl << std::endl;
5550 
5551  std::cout << "ROUTE " << k << " Total Demand = " << routeDemand[ k] << std::endl;
5552  std::cout << "ROUTE " << k << " Total Cost = " << routeCost[ k] << std::endl;
5553  std::cout << "ROUTE " << k << " Nodes " << std::endl;
5554 
5555  //for(vit = routeMap[ k].begin(); vit != routeMap[k].end(); vit++){
5556  // std::cout << *vit << std::endl;
5557  //}
5558  totalCost += routeCost[ k];
5559  }
5560 
5561  std::cout << "Total Cost = " << totalCost << std::endl;
5562  //get the solution
5563  m_initSolMap[ 0] = routeMap;
5564 
5565 
5566  //garbage collection
5567  delete[] routeCost;
5568  routeCost = NULL;
5569  delete[] routeDemand;
5570  routeDemand = NULL;
5571  delete[] xVar;
5572  xVar = NULL;
5573  delete solver->osinstance;
5574  delete solver;
5575  //exit( 1);
5576  if( totalCost >= routeCostInf )return false;
5577  else return true;
5578 
5579 
5580  } catch (const ErrorClass& eclass) {
5581 
5582  std::cout << std::endl << std::endl << std::endl;
5583 
5584  if(routeCost != NULL){
5585  delete[] routeCost;
5586  routeCost = NULL;
5587  }
5588 
5589  if(routeDemand != NULL){
5590  delete[] routeDemand;
5591  routeDemand = NULL;
5592  }
5593 
5594  if(xVar != NULL){
5595  delete[] xVar;
5596  xVar = NULL;
5597  }
5598 
5599 
5600 
5601  // Problem with the parser
5602  throw ErrorClass(eclass.errormsg);
5603 }
5604 
5605 }//1OPT
5606 
5607 
5608 
5609 
5610 
5612 
5613 
5614  int i;
5615  int j;
5616  int numVar;
5617  int numNonHubs;
5618  numNonHubs = m_numNodes - m_numHubs;
5619  //first the w varibles
5620  numVar = numNonHubs;
5621  //u(i, j) varibles with i = hubIndex
5622  numVar += numNonHubs;
5623  //u(i, j) variable where i, j are not hubs
5624  numVar += numNonHubs*numNonHubs - numNonHubs;
5625  int numNonz;
5626  int kount;
5627  int numCon;
5628  CoinSolver *solver = NULL;
5629  SparseVector *objcoeff = NULL;
5630 
5631  numCon = numNonHubs + (numNonHubs*numNonHubs - numNonHubs) + 1;
5632 
5633 
5634 
5635  OSInstance *osinstance = new OSInstance();
5636  try{
5637 
5638  osinstance->setInstanceSource("generated from OSBearcatSoleverXij");
5639  osinstance->setInstanceDescription("Instance for finding a multicommodity cut");
5640  osinstance->setVariableNumber( numVar);
5641 
5642  numVar = 0;
5643 
5644  for(i = m_numHubs; i < m_numNodes; i++){
5645 
5646  if(m_nodeName[ i] != "")
5647  osinstance->addVariable(numVar++, "w[" + m_nodeName[ i] +"]", 0, OSDBL_MAX, 'C');
5648  else
5649  osinstance->addVariable(numVar++, makeStringFromInt("w[", i) +"]", 0, OSDBL_MAX, 'C');
5650 
5651  m_tmpVarMap.insert( std::pair<int,std::string>(numVar, "w[" + m_nodeName[ i] +"]" ) );
5652 
5653  }
5654 
5655 
5656 
5657  for(j = m_numHubs; j < m_numNodes; j++){
5658 
5659  if(m_nodeName[ hubIndex ] != "" && m_nodeName[ j] != "")
5660  osinstance->addVariable(numVar++, "u[" + m_nodeName[ hubIndex] + "," + m_nodeName[ j] +"]", 0, OSDBL_MAX, 'C');
5661  else
5662  osinstance->addVariable(numVar++, makeStringFromInt("u[", hubIndex) + makeStringFromInt(",", j) +"]", 0, OSDBL_MAX, 'C');
5663 
5664  m_tmpVarMap.insert( std::pair<int,std::string>(numVar, m_nodeName[ hubIndex] + "," + m_nodeName[ j] +"]") );
5665 
5666 
5667  }
5668 
5669  for(i = m_numHubs; i < m_numNodes; i++){
5670 
5671 
5672  for(j = m_numHubs; j < i; j++){
5673 
5674 
5675 
5676  if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
5677  osinstance->addVariable(numVar++, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", 0, OSDBL_MAX, 'C');
5678  else
5679  osinstance->addVariable(numVar++, makeStringFromInt("u[", i) + makeStringFromInt(",", j) +"]", 0, OSDBL_MAX, 'C');
5680 
5681 
5682  m_tmpVarMap.insert( std::pair<int,std::string>(numVar, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]") );
5683 
5684  }
5685 
5686  for(j = i + 1; j < m_numNodes; j++){
5687 
5688  if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
5689  osinstance->addVariable(numVar++, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", 0, OSDBL_MAX, 'C');
5690  else
5691  osinstance->addVariable(numVar++, makeStringFromInt("u[", i) + makeStringFromInt(",", j) +"]", 0, OSDBL_MAX, 'C');
5692 
5693  m_tmpVarMap.insert( std::pair<int,std::string>(numVar, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]") );
5694 
5695 
5696  }
5697 
5698  }
5699 
5700 
5701  // now add the objective function
5702  osinstance->setObjectiveNumber( 1);
5703 
5704  // now the coefficient
5705 
5706  objcoeff = new SparseVector( numVar);
5707 
5708  for(i = 0; i < numVar; i++){
5709 
5710  objcoeff->indexes[ i] = i;
5711  objcoeff->values[ i] = 0;
5712 
5713  }
5714 
5715  osinstance->addObjective(-1, "cutViolation", "max", 0.0, 1.0, objcoeff);
5716  objcoeff->bDeleteArrays = true;
5717  delete objcoeff;
5718 
5719  osinstance->setConstraintNumber( numCon );
5720  //bool addConstraint(int index, string name, double lowerBound, double upperBound, double constant);
5721  numCon = 0;
5722  for(i = m_numHubs; i < m_numNodes; i++){
5723 
5724  if(m_nodeName[ hubIndex] != "" && m_nodeName[ i] != "")
5725  osinstance->addConstraint(numCon++, "dualCon[" + m_nodeName[ hubIndex] + "," + m_nodeName[ i] + "]", -OSDBL_MAX, 0, 0);
5726  else
5727  osinstance->addConstraint(numCon++, makeStringFromInt("dualCon[", hubIndex) + makeStringFromInt(",", i) +"]", -OSDBL_MAX, 0, 0);
5728  }
5729 
5730 
5731  for(i = m_numHubs; i < m_numNodes; i++){
5732 
5733 
5734  for(j = m_numHubs; j < i; j++){
5735 
5736  if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
5737  osinstance->addConstraint(numCon++, "dualCon[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", -OSDBL_MAX, 0, 0);
5738  else
5739  osinstance->addConstraint(numCon++, makeStringFromInt("dualCon[", i) + makeStringFromInt(",", j) +"]", -OSDBL_MAX, 0, 0);
5740 
5741 
5742  }
5743 
5744  for(j = i + 1; j < m_numNodes; j++){
5745 
5746  if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
5747  osinstance->addConstraint(numCon++, "dualCon[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", -OSDBL_MAX, 0, 0);
5748  else
5749  osinstance->addConstraint(numCon++, makeStringFromInt("dualCon[", i) + makeStringFromInt(",", j) +"]", -OSDBL_MAX, 0, 0);
5750 
5751 
5752  }
5753 
5754  }
5755 
5756  double upperBound;
5757  upperBound = numVar - numNonHubs ;
5758  upperBound = 100;
5759  osinstance->addConstraint(numCon++, "boundConstraint", -OSDBL_MAX, upperBound, 0);
5760 
5761  //now the A matrix
5762  numCon = numNonHubs + (numNonHubs*numNonHubs - numNonHubs) + 1;
5763  numNonz = 2*numNonHubs;
5764  numNonz += 3*(numNonHubs*numNonHubs - numNonHubs);
5765  numNonz += (numNonHubs*numNonHubs - numNonHubs) + numNonHubs;
5766 
5767 
5768  double *values = new double[ numNonz];
5769  int *columnIndexes = new int[ numNonz];
5770  //store by row major
5771  int *starts = new int[ numCon + 1];
5772 
5773 
5774  kount = 0;
5775  numNonz = 0;
5776  starts[ kount++] = 0;
5777 
5778 
5780 
5781 
5782  int uijKount;
5783  uijKount = numNonHubs;
5784 
5785  for(j = m_numHubs; j < m_numNodes; j++){
5786 
5787  //-u(k, j) + w(j) =l= 0;
5788  //variable w(j)
5789  columnIndexes[ numNonz] = j - m_numHubs ;
5790  values[ numNonz++] = 1.0;
5791 
5792  //variable -u(k, j)
5793  columnIndexes[ numNonz] = uijKount ;
5794  values[ numNonz++] = -1.0;
5795 
5796  starts[ kount++] = numNonz;
5797  uijKount++;
5798  }
5799 
5800 
5801  for(i = m_numHubs; i < m_numNodes; i++){
5802 
5803 
5804  for(j = m_numHubs; j < i; j++){
5805 
5806  //-u(i, j) - w(i) + w(j) =l= 0;
5807  //variable w(i)
5808  columnIndexes[ numNonz] = i - m_numHubs ;
5809  values[ numNonz++] = -1.0;
5810 
5811  //variable w(j)
5812  columnIndexes[ numNonz] = j - m_numHubs ;
5813  values[ numNonz++] = 1.0;
5814 
5815  //variable -u(i, j)
5816  columnIndexes[ numNonz] = uijKount ;
5817  values[ numNonz++] = -1.0;
5818 
5819 
5820  starts[ kount++] = numNonz;
5821  uijKount++;
5822 
5823 
5824  }
5825 
5826  for(j = i + 1; j < m_numNodes; j++){
5827 
5828  //-u(i, j) - w(i) + w(j) =l= 0;
5829  //variable w(i)
5830  columnIndexes[ numNonz] = i - m_numHubs ;
5831  values[ numNonz++] = -1.0;
5832 
5833  //variable w(j)
5834  columnIndexes[ numNonz] = j - m_numHubs ;
5835  values[ numNonz++] = 1.0;
5836 
5837  //variable -u(i, j)
5838  columnIndexes[ numNonz] = uijKount ;
5839  values[ numNonz++] = -1.0;
5840 
5841 
5842  starts[ kount++] = numNonz;
5843  uijKount++;
5844 
5845 
5846  }
5847 
5848  }
5849 
5850  //the last row
5851  for(i = numNonHubs; i < numVar; i++ ){
5852 
5853  //variable u(i, j)
5854  columnIndexes[ numNonz] = i ;
5855  values[ numNonz++] = 1.0;
5856 
5857  }
5858 
5859  starts[ kount++] = numNonz;
5860  osinstance->setLinearConstraintCoefficients(numNonz, false, values, 0, numNonz - 1,
5861  columnIndexes, 0, numNonz - 1, starts, 0, numVar);
5862 
5863  //std::cout << osinstance->printModel() << std::endl;
5864 
5865 
5866  solver = new CoinSolver();
5867  solver->sSolverName ="clp";
5868  solver->osinstance = osinstance;
5869  solver->buildSolverInstance();
5870  solver->osoption = m_osoption;
5871 
5872 
5873  return solver;
5874 
5875 
5876 
5877  } catch (const ErrorClass& eclass) {
5878 
5879  if( objcoeff != NULL ){
5880  delete objcoeff;
5881  objcoeff = NULL;
5882  }
5883  // Problem with the parser
5884  throw ErrorClass(eclass.errormsg);
5885  }
5886 
5887 
5888 }//end getMultiCommodInstance
5889 
5890 
5891 
5893 
5894  int i;
5895  int k;
5896  int numVar;
5897  int numNonHubs;
5898  numNonHubs = m_numNodes - m_numHubs;
5899  //the xki varibles
5900  numVar = m_numHubs*numNonHubs;
5901 
5902  int numNonz;
5903  int kount;
5904 
5905  CoinSolver *solver = NULL;
5906  SparseVector *objcoeff = NULL;
5907 
5908  std::pair<int,int> indexPair1;
5909  std::pair<int,int> indexPair2;
5910 
5911  std::map<int, std::vector<int> > routeMap;
5912 
5913  OSInstance *osinstance = new OSInstance();
5914  try{
5915 
5916  osinstance->setInstanceSource("generated from OSBearcatSoleverXij");
5917  osinstance->setInstanceDescription("Generalized Assignment Instance for finding a feasible solution");
5918  osinstance->setVariableNumber( numVar);
5919 
5920  numVar = 0;
5921  // xki = 1 if hub k serves node i
5922  for(k = 0; k < m_numHubs; k++){
5923 
5924  for(i = m_numHubs; i < m_numNodes; i++){
5925 
5926  if( m_nodeName[ k] != "" && m_nodeName[ i] != "")
5927  osinstance->addVariable(numVar++, "x(" + m_nodeName[ k] + "," + m_nodeName[ i] +")", 0, 1, 'B');
5928  else
5929  osinstance->addVariable(numVar++, makeStringFromInt("x(" ,k) + makeStringFromInt(",", i) +")", 0, 1, 'B');
5930 
5931  }
5932 
5933  }
5934 
5935  // now add the objective function
5936  osinstance->setObjectiveNumber( 1);
5937 
5938  // now the coefficient
5939 
5940  objcoeff = new SparseVector( numVar);
5941 
5942  kount = 0;
5943  for(k = 0; k < m_numHubs; k++){
5944 
5945  indexPair1.first = k;
5946  indexPair2.second = k;
5947 
5948  for(i = m_numHubs; i < m_numNodes; i++){
5949 
5950  indexPair1.second = i;
5951  indexPair2.first = i;
5952  objcoeff->indexes[ kount] = kount;
5953 
5954  if( (m_xVarIndexMap.find( indexPair1) == m_xVarIndexMap.end() ) ||
5955  (m_xVarIndexMap.find( indexPair2) == m_xVarIndexMap.end() ) ){
5956  throw ErrorClass("Index problem in generalized assignment problem to find feasible solution");
5957  }
5958 
5959  objcoeff->values[ kount++] = (m_cost[ m_xVarIndexMap[ indexPair1] ] +
5960  m_cost[ m_xVarIndexMap[ indexPair2] ])*m_demand[i] ;
5961 
5962  }
5963 
5964  }
5965 
5966  osinstance->addObjective(-1, "feasibililtyObj", "min", 0.0, 1.0, objcoeff);
5967  objcoeff->bDeleteArrays = true;
5968  delete objcoeff;
5969  objcoeff = NULL;
5970 
5971  osinstance->setConstraintNumber( m_numNodes );
5972  //bool addConstraint(int index, string name, double lowerBound, double upperBound, double constant);
5973  //
5974  for(k = 0; k < m_numHubs; k++){
5975 
5976  if(m_nodeName[ k] != "" )
5977  osinstance->addConstraint(k, "capacityCon[" + m_nodeName[ k] + "]", 0, m_routeCapacity[ k], 0);
5978  else
5979  osinstance->addConstraint(k, makeStringFromInt("dualCon[", k) +"]", 0, m_routeCapacity[ k], 0);
5980  }
5981 
5982 
5983  for(i = m_numHubs; i < m_numNodes; i++){
5984 
5985  if(m_nodeName[ i] != "" )
5986  osinstance->addConstraint(i, "assingCon[" + m_nodeName[ i] +"]", 1, 1, 0);
5987  else
5988  osinstance->addConstraint(i, makeStringFromInt("assignCon[", i) +"]", 1, 1, 0);
5989 
5990  }
5991 
5992 
5993  //now the A matrix
5994 
5995  numNonz = 2*numVar;
5996 
5997  //store by column major
5998  double *values = new double[ numNonz];
5999  int *rowIndexes = new int[ numNonz];
6000  int *starts = new int[ numVar + 1];
6001 
6002 
6003  kount = 0;
6004  numNonz = 0;
6005  starts[ kount++] = 0;
6006 
6007 
6009 
6010 
6011  for(k = 0; k < m_numHubs; k++){
6012 
6013 
6014  for(i = m_numHubs; i < m_numNodes; i++){
6015 
6016 
6017  rowIndexes[ numNonz] = k ;
6018  values[ numNonz++] = m_demand[ i];
6019 
6020  rowIndexes[ numNonz] = i ;
6021  values[ numNonz++] = 1.0;
6022 
6023  starts[ kount++] = numNonz;
6024 
6025  }
6026  }
6027 
6028 
6029  osinstance->setLinearConstraintCoefficients(numNonz, true, values, 0, numNonz - 1,
6030  rowIndexes, 0, numNonz - 1, starts, 0, numVar);
6031 
6032  //std::cout << osinstance->printModel() << std::endl;
6033 
6034  //
6035 
6036  solver = new CoinSolver();
6037  solver->sSolverName ="cbc";
6038  solver->osinstance = osinstance;
6039  solver->buildSolverInstance();
6040  solver->osoption = m_osoption;
6041  solver->solve();
6042  //now lets get the solution
6043  //first make sure we are optimal
6044  OSResult *osresult;
6045  osresult = solver->osresult;
6046  std::string solStatus;
6047  std::vector<IndexValuePair*> primalValPair;
6048  int vecSize;
6049  double optSolValue;
6050  // the argument is the solution index
6051  solStatus = osresult->getSolutionStatusType( 0 );
6052  // if solStatus is optimal get the optimal solution value
6053  if( solStatus.find("ptimal") != std::string::npos ){
6054  //first index is objIdx, second is solution index
6055  optSolValue = osresult->getOptimalObjValue( -1, 0);
6056  std::cout << "OPTIMAL SOLUTION VALUE " << optSolValue << std::endl;
6057  }else{
6058  throw ErrorClass("There is no feasible solution to this problem!");
6059  }
6060 
6061  primalValPair = osresult->getOptimalPrimalVariableValues( 0);
6062  vecSize = primalValPair.size();
6063 
6064 
6065  kount = 0;
6066  for(k = 0; k < m_numHubs; k++){
6067 
6068  indexPair1.first = k;
6069  std::cout << std::endl << std::endl;
6070  for(i = m_numHubs; i < m_numNodes; i++){
6071 
6072  indexPair1.second = i;
6073  if(primalValPair[ kount ]->value > m_osDecompParam.zeroTol) {
6074 
6075  std::cout << "Variable = " << m_variableNames[ m_xVarIndexMap[ indexPair1] ]
6076  << " value = " << primalValPair[ kount ]->value << std::endl;
6077 
6078  routeMap[k].push_back( i);
6079  }
6080 
6081  kount++;
6082  }
6083 
6084 
6085  }
6086  //exit( 1);
6087  m_initSolMap[ 0] = routeMap;
6088  delete solver;
6089  solver = NULL;
6090 
6091 
6092  } catch (const ErrorClass& eclass) {
6093 
6094  if( objcoeff != NULL ){
6095  delete objcoeff;
6096  objcoeff = NULL;
6097  }
6098  // Problem with the parser
6099  throw ErrorClass(eclass.errormsg);
6100  }
6101 
6102 
6103 }//end getFeasibleSolution
6104 
6106 
6107  int i;
6108  int j;
6109  int kount;
6110  std::pair<int,int> indexPair;
6111  //construct index map
6112  kount = 0;
6113  for(i = 0; i < m_numNodes; i++){
6114 
6115  for(j = 0; j < i; j++){ //j < i
6116 
6117  indexPair.first = i;
6118  indexPair.second = j;
6119  m_xVarIndexMap[ indexPair] = kount;
6120  kount++;
6121  }
6122 
6123  for(j = i + 1; j < m_numNodes; j++){ // i < j
6124 
6125  indexPair.first = i;
6126  indexPair.second = j;
6127  m_xVarIndexMap[ indexPair] = kount;
6128  kount++;
6129  }
6130  }
6131  //end construct map
6132 
6133 }//end getVariableIndexMap
6134 
6135 
6137 
6138  int k1;
6139  int k2;
6140 
6141  double tmpVal;
6142  double *tmpCap;
6143  int tmpIdx;
6144  tmpCap = new double[ m_numHubs];
6145 
6146  for(k1 = 0; k1 < m_numHubs; k1++) tmpCap[ k1] = m_routeCapacity[ k1]; //initialize capacities
6147  for(k1 = 0; k1 < m_numHubs; k1++) m_hubPoint[ k1] = k1; //initialize capacities
6148 
6149  for(k1 = 0; k1 < m_numHubs - 1; k1++){
6150 
6151 
6152  for(k2 = k1 + 1; k2 < m_numHubs; k2++){
6153 
6154  if( tmpCap[ k2 ] < tmpCap[ k1 ] ){ //make switch
6155 
6156 
6157  tmpVal = tmpCap[ k1 ];
6158  tmpCap[ k1 ] = tmpCap[ k2 ];
6159  tmpCap[ k2 ] = tmpVal;
6160 
6161  tmpIdx = m_hubPoint[ k1];
6162  m_hubPoint[ k1] = m_hubPoint[ k2];
6163  m_hubPoint[ k2] = tmpIdx;
6164 
6165  }
6166 
6167  }// end k2 loop
6168  }// end k1 loop
6169 
6170  for(k1 = 0; k1 < m_numHubs; k1++) std::cout << "m_hubPoint = " << m_hubPoint[ k1] << std::endl;
6171  for(k1 = 0; k1 < m_numHubs; k1++) std::cout << "tmp Cap = " << tmpCap[ k1] << std::endl;
6172  for(k1 = 0; k1 < m_numHubs; k1++) std::cout << "hub capacity = " << m_routeCapacity[ m_hubPoint[ k1] ]<< std::endl;
6173  //exit( 1);
6174  delete[] tmpCap;
6175  tmpCap = NULL;
6176 
6177 }
6178 
6179 
6180 
6181 
6182 
6183 
6184 
OSOption * osoption
std::string makeStringFromInt(std::string theString, int theInt)
#define d1
Implements a solve method for the Coin solvers.
Definition: OSCoinSolver.h:38
virtual void buildSolverInstance()
The implementation of the corresponding virtual function.
OsiSolverInterface * osiSolver
osiSolver is the osi solver object – in this case clp, glpk, cbc, cplex, symphony or dylp
Definition: OSCoinSolver.h:93
virtual void solve()
The implementation of the corresponding virtual function.
std::string sSolverName
sSolverName is the name of the Coin solver used, e.g.
OSInstance * osinstance
osinstance holds the problem instance in-memory as an OSInstance object
OSOption * osoption
osoption holds the solver options in-memory as an OSOption object
OSResult * osresult
osresult holds the solution or result of the model in-memory as an OSResult object
used for throwing exceptions.
Definition: OSErrorClass.h:32
std::string errormsg
errormsg is the error that is causing the exception to be thrown
Definition: OSErrorClass.h:42
int * m_demand
m_demand is the vector of node demands
std::map< int, std::string > m_tmpVarMap
double getRouteDistance(int numNodes, int hubIndex, CoinSolver *solver, std::vector< int > zk, double *xVar)
call this method to get a minimum TSP tour for a given assignment of nodes to routes INPUT: int numNo...
double * m_cost
the distance/cost vectors
double calcNonlinearRelax(std::vector< double > &m_zRootLPx_vals)
calculate the nonlinear relaxation value for an LP solution
std::string * m_nodeName
m_nodeName is the vector of node names
int * m_upperBoundL
largest possible L-value on a route – this will be the minimum of m_routeCapacity and total demand
void getCutsMultiCommod(const double *thetaVar, const int numThetaVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)
This is the routine that generates the multi-item cuts.
std::string m_initOSiLFile
int * m_BmatrixRowIndex
m_BmatrixRowIndex holds the index of the convexity row that the constraint corresponds to,...
CoinSolver * getTSP(int numNodes, double *cost)
call this method to get a TSP instance
int * m_lowerBoundL
smallest possible L-value on a route for now this will equal
int * m_hubPoint
m_hubPoint[ k] points to the the k'th hub that we use in getOptL
ClpSimplex * m_separationClpModel
int * m_separationIndexMap
m_separationIndexMap maps the variable index into the appropriate row in the separation problem for t...
double qrouteCost(const int &k, const int &l, const double *c, int *kountVar)
kipp – document
void getVariableIndexMap()
this method will populate: std::map<std::pair<int, int>,int>m_xVarIndexMap this gives us
int * m_routeMinPickup
the minimum number of students that we pickup on a route this can vary with the route/hub
double ** m_rc
the reduced cost vector for each k, we asssume order is (l, i, j)
int m_maxThetaNonz
m_maxMasterNonz is the maximumn number of nonzero elements we allow in the transformation matrix betw...
virtual OSInstance * getInitialRestrictedMaster()
OSInstance* OSBearcatSolverXij::getInitialRestrictedMaster( ){.
~OSBearcatSolverXij()
Default Destructor.
void getInitialSolution()
generate an intitial feasible solution in theta space for the initial master
std::map< std::pair< int, int >, int > m_xVarIndexMap
m_xVarIndexMap takes a variable indexed by (i,j) and returns the index of the variable in one dimensi...
void getOptions(OSOption *osoption)
virtual void resetMaster(std::map< int, int > &inVars, OsiSolverInterface *si)
INPUT:
int m_upperBoundLMax
largest possible L-value over all routes
virtual void initializeDataStructures()
allocate memory and initialize arrays
virtual void pauHana(std::vector< int > &m_zOptIndexes, std::vector< double > &m_zRootLPx_vals, int numNodes, int numColsGen, std::string message)
OSInstance * m_osinstanceSeparation
std::map< int, std::map< int, std::vector< int > > > m_initSolMap
the index on the outer map is on the solution number, the index on the inner map indexes the route nu...
void permuteHubs()
this method will calculate a permuation of the hubs so that they are in ascending order,...
virtual void getBranchingCut(const double *thetaVar, const int numThetaVar, const std::map< int, int > &varConMap, int &varIdx, int &numNonz, int *&indexes, double *&values)
RETURN VALUES: varIdx – the variable number x_{ij} for branching numNonz – number of theta indexes in...
int getBranchingVar(const double *theta, const int numThetaVar)
RETURN VALUES: return the integer index of a fractional x_{ij} variable.
int m_minDemand
m_minDemand is the value of the minimum demand node – it is not the minimum demand that must be carri...
std::vector< CoinSolver * > m_multCommodCutSolvers
m_multCommodCutSolvers is a vector of solvers, one solver for each hub, used to find multicommodity f...
OSInstance * getSeparationInstance()
void getOptL(double **c)
bool m_costSetInOption
m_costSetInOption is true if the costs are set using the OSOption file
void getFeasibleSolution()
call this method to get generate an instance of the generalized assignment problem and find a feasibl...
int * m_routeCapacity
the route capacity – bus seating limit this can vary with the route/hub
OSBearcatSolverXij()
Default Constructor.
bool OneOPT()
try and find a feasible solution, return false if solution not feasible
bool m_use1OPTstart
if m_use1OPTstart is true we use the option file to fix the nodes to hubs found by SK's 1OPT heuristi...
int * m_convexityRowIndex
m_convexityRowIndex holds the index of the convexity row that the theta columns are in.
virtual void getColumns(const double *yA, const int numARows, const double *yB, const int numBRows, int &numNewColumns, int *&numNonz, double *&cost, int **&rowIdx, double **&values, double &lowerBound)
RETURN VALUES: int numNewColumns – number of new columns generated int* numNonz – number of nonzeros ...
void getCutsX(const double *xVar, const int numXVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)
RETURN VALUES: int numNewRows – number of new rows generated int* numNonz – number of nonzeros in eac...
void getCutsTheta(const double *thetaVar, const int numThetaVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)
RETURN VALUES: int numNewRows – number of new rows generated int* numNonz – number of nonzeros in eac...
CoinSolver * getMultiCommodInstance(int hubIndex)
call this method to get an instance that is used to generate a multicommodity cut
void calcReducedCost(const double *yA, const double *yB)
calculate the reduced costs c – input of the objective function costs yA – dual values on node assign...
double artVarCoeff
artVarCoeff is the coefficient of the artificial variable in the objective function
Definition: OSDecompParam.h:55
double zeroTol
we terminate column generation when the reduced costs are not smaller than zeroTol
Definition: OSDecompParam.h:50
OSDecompParam m_osDecompParam
share the parameters with the decomposition solver
int m_maxBmatrixCon
m_maxBmatrixCon is the maximum number of B matrix constraints it is the number of tour breaking const...
double m_rootLPValue
double m_bestIPValue
int m_numNodes
m_numNodes is the number of nodes (both pickup and hub) in the model
int m_maxMasterRows
m_maxMasterColumns is the maximumn number of rows we allow in the master, in this application it is e...
int m_numBmatrixCon
m_numBmatrixCon is the number of constraints in B - 1, we have the -1 because: m_pntBmatrix[ k] point...
int m_numHubs
m_numHubs is the number of hubs/routes
std::string * m_variableNames
std::set< std::pair< int, double > > intVarSet
intVarSet holds and std::pair where the first element is the index of an integer variable and the sec...
OSInstance * m_osinstanceMaster
OSOption * m_osoption
int m_maxMasterColumns
m_maxMasterColumns is the maximumn number of columns we allow in the master
int m_maxBmatrixNonz
m_maxBmatrixNonz is the maximum number of nonzero elements in the B matrix constraints
double * m_BmatrixVal
double m_bestLPValue
The in-memory representation of an OSiL instance..
Definition: OSInstance.h:2263
double * getConstraintLowerBounds()
Get constraint lower bounds.
double * getVariableUpperBounds()
Get variable upper bounds.
bool setConstraintNumber(int number)
set the number of constraints.
bool addVariable(int index, std::string name, double lowerBound, double upperBound, char type)
add a variable.
std::string printModel()
Print the infix representation of the problem.
bool addConstraint(int index, std::string name, double lowerBound, double upperBound, double constant)
add a constraint.
int getConstraintNumber()
Get number of constraints.
bool bConstraintsModified
bConstraintsModified is true if the constraints data has been modified.
Definition: OSInstance.h:2298
int getLinearConstraintCoefficientNumber()
Get number of specified (usually nonzero) linear constraint coefficient values.
SparseMatrix * getLinearConstraintCoefficientsInRowMajor()
Get linear constraint coefficients in row major.
SparseMatrix * getLinearConstraintCoefficientsInColumnMajor()
Get linear constraint coefficients in column major.
bool setInstanceSource(std::string source)
set the instance source.
double ** getDenseObjectiveCoefficients()
getDenseObjectiveCoefficients.
bool setLinearConstraintCoefficients(int numberOfValues, bool isColumnMajor, double *values, int valuesBegin, int valuesEnd, int *indexes, int indexesBegin, int indexesEnd, int *starts, int startsBegin, int startsEnd)
set linear constraint coefficients
double * getVariableLowerBounds()
Get variable lower bounds.
int getVariableNumber()
Get number of variables.
bool setInstanceDescription(std::string description)
set the instance description.
bool addObjective(int index, std::string name, std::string maxOrMin, double constant, double weight, SparseVector *objectiveCoefficients)
add an objective.
bool setObjectiveNumber(int number)
set the number of objectives.
bool setVariableNumber(int number)
set the number of variables.
double * getConstraintUpperBounds()
Get constraint upper bounds.
The Option Class.
Definition: OSOption.h:3565
std::vector< SolverOption * > getSolverOptions(std::string solver_name)
Get the options associated with a given solver.
Definition: OSOption.cpp:4508
The Result Class.
Definition: OSResult.h:2549
std::string getSolutionStatusType(int solIdx)
Get the [i]th optimization solution status type, where i equals the given solution index.
Definition: OSResult.cpp:2051
double getOptimalObjValue(int objIdx, int solIdx)
Get one solution's optimal objective value.
Definition: OSResult.cpp:3065
std::vector< IndexValuePair * > getOptimalPrimalVariableValues(int solIdx)
Get one solution of optimal primal variable values.
Definition: OSResult.cpp:2215
double getObjValue(int solIdx, int objIdx)
Definition: OSResult.cpp:3050
int * indexes
indexes holds an integer array of rowIdx (or colIdx) elements in coefMatrix (AMatrix).
Definition: OSGeneral.h:258
int * starts
starts holds an integer array of start elements in coefMatrix (AMatrix), which points to the start of...
Definition: OSGeneral.h:252
double * values
values holds a double array of value elements in coefMatrix (AMatrix), which contains nonzero element...
Definition: OSGeneral.h:264
a sparse vector data structure
Definition: OSGeneral.h:123
double * values
values holds a double array of nonzero values.
Definition: OSGeneral.h:164
int * indexes
indexes holds an integer array of indexes whose corresponding values are nonzero.
Definition: OSGeneral.h:159
bool bDeleteArrays
bDeleteArrays is true if we delete the arrays in garbage collection set to true by default
Definition: OSGeneral.h:149
const int OSINT_MAX
Definition: OSParameters.h:94
const double OSDBL_MAX
Definition: OSParameters.h:93
OSResult * osresult