My Project
OSBearcatSolverXkij.cpp
Go to the documentation of this file.
1/* $Id: OSBearcatSolverXkij.cpp 3038 2009-11-07 11:43:44Z kmartin $ */
14#include "OSBearcatSolverXkij.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
33#include "CoinTime.hpp"
34
35#include "ClpFactorization.hpp"
36#include "ClpNetworkMatrix.hpp"
37#include "OsiClpSolverInterface.hpp"
38
39
40#ifdef HAVE_CMATH
41# include <cmath>
42#else
43# ifdef HAVE_MATH_H
44# include <math.h>
45# else
46# error "don't have header file for math"
47# endif
48#endif
49
50#ifdef HAVE_CTIME
51# include <ctime>
52#else
53# ifdef HAVE_TIME_H
54# include <time.h>
55# else
56# error "don't have header file for time"
57# endif
58#endif
59
60
61
62std::string makeStringFromInt2(std::string theString, int theInt);
63
64
66 std::cout << "INSIDE OSBearcatSolverXkij CONSTRUCTOR with OSOption argument" << std::endl;
67}//end default OSBearcatSolverXkij constructor
68
70 std::cout << "INSIDE OSBearcatSolverXkij CONSTRUCTOR with OSOption argument" << std::endl;
71
72
75
76 m_upperBoundL = NULL;
79
80 m_u = NULL;
81 m_v = NULL;
82 m_g = NULL;
83 m_px = NULL;
84 m_tx =NULL;
85 m_varIdx = NULL;
86
87 m_optL = NULL;
88 m_optD = NULL;
89 m_vv = NULL;
90 m_vvpnt = NULL;
91
92 m_demand = NULL;
93 m_cost = NULL;
94
95 m_rc = NULL;
96
97 m_routeCapacity = NULL;
98 m_routeMinPickup = NULL;
99
101
102}//end OSBearcatSolverXkijDestructor
103
105
106 int k;
107 int i;
108 int l;
109
110 try{
111
112
113 //get all the parameter values
115
117
118 m_upperBoundL = new int[ m_numHubs];
119 m_lowerBoundL = new int[ m_numHubs];
120
121 for(k = 0; k < m_numHubs; k++){
122
125 //set m_upperBoundL cannot exceed total demand
128
129 }
130
131 //m_varIdx = new int[ m_numNodes];
132 // I think we can make this the max of m_upperBoundL
133 // over the k
134 m_varIdx = new int[ m_upperBoundLMax + 1];
135
136
137 m_u = new double*[ m_numNodes];
138 m_v = new double*[ m_numNodes];
139 m_g = new double*[ m_numNodes];
140
141 m_px = new int*[ m_numNodes];
142 m_tx = new int*[ m_numNodes];
143
144
145
155 for (i = 0; i < m_numNodes; i++) {
156
157 //kipp we can use the biggest possible m_upperBoundL
158 m_u[ i] = new double[ m_upperBoundLMax + 1];
159 m_v[ i] = new double[ m_upperBoundLMax + 1];
160
161
162
163 for(l = 0; l <= m_upperBoundLMax; l++){
164
165 m_u[ i][l] = OSDBL_MAX;
166 m_v[ i][l] = OSDBL_MAX;
167 }
168
169 m_g[ i] = new double[ m_numNodes];
170 m_px[ i] = new int[ m_upperBoundLMax + 1];
171 m_tx[ i] = new int[m_upperBoundLMax + 1];
172
173
174 }
175
176
177 //outer dynamic programming arrays
178 m_optL = new int[ m_numHubs];
179 m_optD = new int[ m_numHubs];
180
181 m_vv = new double*[ m_numHubs];
182 m_vvpnt = new int*[ m_numHubs];
183 m_cost = new double*[ m_numHubs];
184 m_rc = new double*[ m_numHubs];
185
186 for (k = 0; k < m_numHubs; k++) {
187
188
189 m_vv[ k] = new double[ m_totalDemand + 1];
190 m_vvpnt[ k] = new int[ m_totalDemand + 1];
191 m_cost[ k] = new double[ m_numNodes*m_numNodes - m_numNodes];
192
193 //allocate memory for the reduced cost vector.
194 //assume order is k, l, i, j
195 //assume order is l, i, j
196 m_rc[ k] = new double[ m_upperBoundL[ k]*(m_numNodes*m_numNodes - m_numNodes)];
197
198
199 }
200
201 m_optValHub = new double[ m_numHubs];
202
203 m_variableNames = new string[ m_numNodes*(m_numNodes - 1)];
204
206
207 //these are constraints that say we must be incident to each
208 //non-hub node -- there are m_numNodes - m_numHubs of these
209 m_pntAmatrix = new int[ m_numNodes - m_numHubs + 1];
210 //the variables -- the variable space we are in is x_{ij} NOT
211 // x_{ijk} -- a bit tricky
212 //m_Amatrix[ j] is a variable index -- this logic works
213 //since the Amatrix coefficient is 1 -- we don't need a value
214 //it indexes variable that points into the node
215 m_Amatrix = new int[ (m_numNodes - m_numHubs)*(m_numNodes - 1) ];
217
218 //this has size of the number of x variables
219 int numVar = m_numNodes*m_numNodes - m_numHubs ;
220 m_tmpScatterArray = new int[ numVar ];
221 for(i = 0; i < numVar; i++){
222
223 m_tmpScatterArray[ i] = 0;
224
225 }
226
227 //New column arrays
228 m_newColumnNonz = new int[ m_numHubs] ; //at most one column per Hub
229 m_costVec = new double[ m_numHubs];
230 m_newColumnRowIdx = new int*[ m_numHubs];
231 m_newColumnRowValue = new double*[ m_numHubs];
232
233 for (k = 0; k < m_numHubs; k++) {
234 //size of arrray is maximum number of constraints
235 // 1) coupling + 2) tour breaking + 3) branching
238
239 }
240
241
242
243 //New row arrays
244 m_newRowNonz = new int[ m_numHubs] ; //at most one cut per Hub
245 m_newRowColumnIdx = new int*[ m_numHubs] ; //at most one cut per Hub
246 m_newRowColumnValue = new double*[ m_numHubs] ; //at most one cut per Hub
247 m_newRowUB = new double[ m_numHubs]; //at most one cut per Hub
248 m_newRowLB = new double[ m_numHubs]; //at most one cut per Hub
249
250
251 //for now, the number of columns will be m_maxMasterColumns
252 for (k = 0; k < m_numHubs; k++) {
253
254 m_newRowColumnValue[ k] = new double[ m_maxMasterColumns];
256
257 }
258
259 //new array for keeping track of convexity rows
261 //new arrays for branches
262
264 branchCutValues = new double[ m_maxMasterColumns];
265
266 m_thetaPnt = new int[ m_maxMasterColumns + 1];
267 for(i = 0; i <= m_maxMasterColumns; i++){
268 m_thetaPnt[ i] = 0;
269 }
270 m_thetaCost = new double[ m_maxMasterColumns];
271 m_thetaIndex = new int[ m_maxThetaNonz];
272 m_numThetaVar = 0;
273 m_numThetaNonz = 0;
275
276
277 //the number of cuts will be at most nubmer of tour
278 // breaking constraints + braching variable cuts
279 // branching variable constraints = m_numNodes*(m_numNodes - 1)
280 m_pntBmatrix = new int[ m_maxBmatrixCon];
281 // number of nonzeros in the Bmatrix
282 m_Bmatrix = new int[ m_maxBmatrixNonz];
283 m_numBmatrixCon = 0;
286
287;
288
289
291
292 for(i = 0; i < m_numNodes*(m_numNodes - 1); i++){
293
295
296 }
297
298
299
300
301 //kipp -- move this later
303 } catch (const ErrorClass& eclass) {
304
305 throw ErrorClass(eclass.errormsg);
306
307 }
308
309
310}//end initializeDataStructures
311
312
314
315 std::cout << "INSIDE ~OSBearcatSolverXkij DESTRUCTOR" << std::endl;
316
317
318
319 //delete data structures for arrays used in calculating minimum reduced cost
320 int i;
321
322 delete[] m_routeCapacity;
323 m_routeCapacity = NULL;
324
325
326 delete[] m_routeMinPickup;
327 m_routeMinPickup = NULL;
328
329 for(i = 0; i < m_numNodes; i++){
330
331
332
333 delete[] m_v[i];
334 delete[] m_g[i];
335 delete[] m_px[i];
336 delete[] m_tx[i];
337 delete[] m_u[i];
338
339
340 }
341
342 delete[] m_u;
343 m_u= NULL;
344
345 delete[] m_v;
346 m_v= NULL;
347
348 delete[] m_g;
349 m_g= NULL;
350
351 delete[] m_px;
352 m_px= NULL;
353
354 delete[] m_tx;
355 m_tx= NULL;
356
357
358
359 if(m_demand != NULL){
360 //std::cout << "I AM DELETING m_demand" << std::endl;
361 delete[] m_demand;
362 }
363
364
365 if(m_varIdx != NULL) delete[] m_varIdx;
366
367 for(i = 0; i < m_numHubs; i++){
368
369 delete[] m_vv[i];
370 delete[] m_vvpnt[i];
371 delete[] m_cost[ i];
372 delete[] m_rc[ i];
373
374
375 }
376 delete[] m_optL;
377 m_optL = NULL;
378 delete[] m_optD;
379 m_optD = NULL;
380 delete[] m_vv;
381 m_vv = NULL;
382 delete[] m_vvpnt;
383 m_vvpnt = NULL;
384
385 delete[] m_cost;
386 m_cost = NULL;
387
388 delete[] m_rc;
389 m_rc = NULL;
390
391 delete[] m_upperBoundL;
392 m_upperBoundL = NULL;
393
394 delete[] m_lowerBoundL;
395 m_lowerBoundL = NULL;
396
397
398 delete[] m_optValHub;
399 m_optValHub = NULL;
400
401 delete[] m_variableNames;
402 m_variableNames = NULL;
403
404 delete[] m_pntAmatrix;
405 m_pntAmatrix = NULL;
406
407 delete[] m_Amatrix;
408 m_Amatrix = NULL;
409
410 delete[] m_tmpScatterArray;
411 m_tmpScatterArray = NULL;
412
413 delete[] m_newColumnNonz ;
414 m_newColumnNonz = NULL;
415 delete[] m_costVec ;
416 m_costVec = NULL;
417
418 for(i = 0; i < m_numHubs; i++){
419
420 delete[] m_newColumnRowIdx[i];
421 delete[] m_newColumnRowValue[i];
422 }
423
424 delete[] m_newColumnRowIdx;
425 m_newColumnRowIdx = NULL;
426
427 delete[] m_newColumnRowValue;
428 m_newColumnRowValue = NULL;
429
430
431 delete[] convexityRowIndex;
432 convexityRowIndex = NULL;
433
434
435 //getCut arrays
436 delete[] m_newRowNonz;
437 m_newRowNonz = NULL;
438
439 delete[] m_newRowUB ;
440 m_newRowUB = NULL;
441
442 delete[] m_newRowLB ;
443 m_newRowLB = NULL;
444
445 //garbage collection on row arrays
446 for (i = 0; i < m_numHubs; i++) {
447
448 delete[] m_newRowColumnValue[ i];
449 delete[] m_newRowColumnIdx[ i];
450
451 }
452
453 delete[] m_newRowColumnIdx;
454 m_newRowColumnIdx = NULL;
455
456 delete[] m_newRowColumnValue;
457 m_newRowColumnValue = NULL;
458
459
460 delete[] branchCutIndexes ;
461 branchCutIndexes = NULL;
462
463 delete[] branchCutValues ;
464 branchCutValues = NULL;
465
466
467 delete[] m_thetaPnt;
468 m_thetaPnt = NULL;
469
470 delete[] m_thetaIndex;
471 m_thetaIndex = NULL;
472
473
474 delete[] m_thetaCost;
475 m_thetaCost = NULL;
476
477
478 delete[] m_pntBmatrix ;
479 m_pntBmatrix = NULL;
480
481 delete[] m_Bmatrix ;
482 m_Bmatrix = NULL;
483
484
485
486
487 delete[] m_separationIndexMap;
489
492
495
496}//end ~OSBearcatSolverXkij
497
498
499
500
501
502
503
505
506 //initialize the first HUB
507
508 int k;
509 int d;
510 int d1;
511 int kountVar;
512 double testVal;
513 int l;
514 //int startPntInc;
515 double trueMin;
516
517 bool isFeasible;
518 isFeasible = false;
519
520 kountVar = 0;
521 //startPntInc = m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes);
522
523 m_vv[ 0][ 0] = 0;
524 for(d = 1; d <= m_totalDemand; d++){
525
526 m_vv[ 0][ d] = OSDBL_MAX;
527
528 }
529 //initialize up to last hub
530 for(k = 1; k < m_numHubs - 1; k++){
531 for(d = 0; d <= m_totalDemand; d++){
532
533 m_vv[ k][ d] = OSDBL_MAX;
534
535 }
536 }
537
538
539 //now loop over the other HUBS
540
541 int dlower;
542 dlower = 0;
543
544 for(k = 1; k < m_numHubs; k++){
545
546 dlower += m_lowerBoundL[ k - 1];
547
548 //kipp make d the min demand for the previous routes
549 for(d = dlower; d <= m_totalDemand; d++){
550
551 m_vv[ k][ d] = OSDBL_MAX;
552
553 //d1 is the state variable at stage k -1
554 for(d1 = 0; d1 <= m_totalDemand; d1++){
555
556 l = d - d1;
557 //kipp make m_upperBoundL the route capapcity
558 if( (m_vv[ k - 1][ d1] < OSDBL_MAX) && (l <= m_upperBoundL[ k - 1]) && (l >= m_lowerBoundL[ k - 1]) ){
559
560
561 // l was the decision at state d1 in stage k-1
562 // l + d1 brings us to state d at stage k
563 // d is the total carried on routes 0 -- k-1
564
565 testVal = qrouteCost(k - 1, l, c[ k - 1], &kountVar);
566
567 //std::cout << "L = " << l << std::endl;
568 //std::cout << "testVal " << testVal << std::endl;
569
570 if( m_vv[ k-1][ d1] + testVal < m_vv[ k][ d] ){
571
572 m_vv[ k][ d] = m_vv[ k-1][ d1] + testVal;
573 //now point to the best way to get to d
574 m_vvpnt[ k][ d] = d1;
575
576 }
577
578
579 }
580
581 }
582
583 }
584
585 //c+= startPntInc ;
586
587 }// end for on k
588
589 trueMin = OSDBL_MAX;
590 //we now enter the last stage through the other hubs
591 // have satisfied total demand d
592
593
594 //if (m_numHubs > 1) dlower = 1;
595
596 for(d = dlower; d < m_totalDemand; d++){
597
598 //std::cout << "m_vv[ m_numHubs - 1 ][ d] " << m_vv[ m_numHubs - 1 ][ d] << std::endl;
599 l = m_totalDemand - d;
600
601 if(m_vv[ m_numHubs - 1 ][ d] < OSDBL_MAX && l <= m_upperBoundL[ m_numHubs - 1] && l >= m_lowerBoundL[ m_numHubs - 1]){
602
603 //must execute this loop at least once
604
605 //std::cout << "LL = " << l << std::endl;
606
607 isFeasible = true;
608
609
610 testVal = qrouteCost(m_numHubs -1 , l, c[ m_numHubs -1], &kountVar);
611
612 //std::cout << "l = " << l << std::endl;
613 //std::cout << "testVal = " << testVal << std::endl;
614
615 if(m_vv[ m_numHubs - 1][ d] + testVal < trueMin){
616
617 trueMin = m_vv[ m_numHubs -1][ d] + testVal;
618 m_optD[ m_numHubs -1 ] = d;
619 m_optL[ m_numHubs -1 ] = l;
620
621 }
622
623
624 }
625 }
626
627 //std::cout << "TRUE MIN = " << trueMin << std::endl;
628
629 if( isFeasible == false){
630
631 std::cout << "NOT ENOUGH CAPACITY " << std::endl;
632 throw ErrorClass( "NOT ENOUGH CAPACITY ");
633 }
634
635 k = m_numHubs -1;
636
637 while( k - 1 >= 0) {
638
639 m_optD[ k - 1 ] = m_vvpnt[ k][ m_optD[ k ] ];
640
641 m_optL[ k - 1 ] = m_optD[ k ] - m_optD[ k - 1 ] ;
642
643 //std::cout << "k = " << k << std::endl;
644 //std::cout << "m_optD[ k ] = " << m_optD[ k ] << std::endl;
645 //std::cout << "m_optD[ k -1 ] " << m_optD[ k - 1 ] << std::endl;
646
647 k--;
648
649
650 }
651
652}//end getOptL
653
654
655
656
657
658
659double OSBearcatSolverXkij::qrouteCost(const int& k, const int& l, const double* c, int* kountVar){
660
661 //critical -- nodes 0, ..., m_numNodes - 1 are the hub nodes
662 // we are doing the calculation for hub k, k <= m_numNodes - 1
663 //l should be the total demand on the network -- we are NOT using 0 based counting
664 double rcost;
665 rcost = OSDBL_MAX;
666
667
668
669 if(l < 0){
670
671 std::cout << "LVALUE NEGATIVE " << l << std::endl;
672 exit( 1);
673 }
674
675
676
677 if(l > m_upperBoundL[ k]){
678
679 std::cout << "LVALUE BIGGER THAN UPPER BOUND " << l << std::endl;
680 exit( 1);
681 }
682
683 //start of the cost vector for hub k plus move to start of l demands
684 // now move the pointer up
685 //int startPnt = m_upperBoundL[ k]*(m_numNodes*m_numNodes - m_numNodes) + (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
686
687 int startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
688 c+= startPnt ;
689
690
691
692 *kountVar = 0;
693 int bestLastNode;
694 //initialize
695 bestLastNode = OSINT_MAX;
696 int i;
697 int j;
698 int l1;
699 int l2;
700 //for(i = 0; i < 20; i++){
701 // std::cout << "i = " << i << " c[i] = " << *(c + i) << std::endl ;
702 //}
703
704
705
706 // l is the total demand on this network
707 //address of node (j, i) is j*(m_numNodes-1) + i when i < j
708 //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
709
710 //
711 // initialize
712
713
714
715 for(i = m_numHubs; i < m_numNodes; i++){
716
717
718 for(l1 = m_minDemand; l1 <= l; l1++){ //l-1 is total demand on network
719
720 m_u[i][l1] = OSDBL_MAX;
721 m_v[i][l1] = OSDBL_MAX;
722 m_px[i][l1] = -1; //a node we don't have
723 if(l1 == *(m_demand + i) ){
724
725 m_px[i][l1] = k;
726 // want the cost for arc (k, i)
727 // note: k < i so use that formula
728 m_u[i][l1] = *(c + k*(m_numNodes - 1) + i - 1); // put in correct cost
729
730
731
732 }
733 }
734 }
735 //end initialize
736
737 //
738
739 //if l = minDemand we visit only one node, let's get the reduced cost in this case
740 if(l == m_minDemand){
741
742 for(i = m_numHubs; i < m_numNodes; i++){
743
744
745 if( m_u[i][l] + *(c + i*(m_numNodes-1) + k ) < rcost){
746
747 rcost = m_u[i][l] + *(c + i*(m_numNodes-1) + k );
748
749 //std::cout << " m_u[i][l2] = " << m_u[i][l2] << std::endl;
750
751 //std::cout << " *(c + i*(m_numNodes-1) + k ) = " << *(c + i*(m_numNodes-1) + k ) << std::endl;
752 //push node back
753 bestLastNode = i;
754 }
755
756 }
757
758 //go from node bestLastNode to node k
759 //node bestLastNode is a higher number than k
760 *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
761 *(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1) + bestLastNode - 1;
762
763
764
765 return rcost;
766 }//end if on l == minDemand
767
768
769// now calculate values for demand 2 or greater
770 //address of node (j, i) is j*(m_numNodes-1) + i when i < j
771 //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
772 // we start l2 at 2 since demand must be at least 1
773 // change to min demand + 1
774 int lowerVal = m_minDemand + 1;
775 for(l2 = lowerVal; l2 <= l; l2++){// loop over possible demand values assuming we have already gone to at least one node
776
777 for(i = m_numHubs; i < m_numNodes; i++) { //we are finding least cost to node i
778
779
780 if( *(m_demand + i) < l2 ){ // kipp < OR <= ????
781
782 for(j = m_numHubs; j < i; j++){ //we are coming from node j
783
784
785
786 //calculate g -- l2 - d[ i] is the demand trough and including node j
787 l1 = l2 - *(m_demand + i);
788
789 if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
790
791
792 m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i - 1) ;
793
794
795
796
797 }else{
798
799 m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i - 1) ;
800
801
802
803 }
804
805 // update u and the pointer for p
806
807 if(m_g[j][i] < m_u[i][l2] ){
808
809 m_u[i][l2] = m_g[j][i];
810 m_px[i][l2] = j;
811
812 }
813
814
815
816 }//end of first for loop on j, j < i
817
818
819 for(j = i + 1; j < m_numNodes; j++){ //we are coming from node j
820
821
822 //calculate g -- l2 - d[ i] is the demand trough and including node j
823 l1 = l2 - *(m_demand + i);
824
825 if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
826
827
828 m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i ) ;
829
830
831 }else{
832
833 m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i ) ;
834
835 }
836
837 // update u and the pointer for p
838
839 if(m_g[j][i] < m_u[i][l2] ){
840
841 m_u[i][l2] = m_g[j][i];
842 m_px[i][l2] = j;
843
844 }
845
846
847 }//end of second for loop on j, j > i
848
849
850 //now calculate the second best solution and point
851 //right now m_px[i][l2] points to the best j node
852
853 for(j =m_numHubs; j < m_numNodes; j++){ //we are coming from node j
854
855 if(j != i){
856
857 if( (m_g[j][i] < m_v[i][l2] ) && (m_px[i][l2] != j) ){ // kipp the && gives the second best
858
859 //if( m_g[j][i] < m_v[i][l2] ) {
860
861 m_v[i][l2] = m_g[j][i];
862 m_tx[i][l2] = j;
863
864
865 }
866
867 }
868
869
870 }//end second best calculation, another for loop on j
871
872 //now if l2 = l we are done
873 if(l2 == l ){
874
875 if( m_u[i][l2] + *(c + i*(m_numNodes-1) + k ) < rcost){
876
877 rcost = m_u[i][l2] + *(c + i*(m_numNodes-1) + k );
878
879 bestLastNode = i;
880 }
881
882 }
883
884
885 }//end if on demand less than l2
886
887
888 }//i loop
889
890
891 }//l2 loop
892
893
894 //std::cout << "best Last Node = " << bestLastNode << std::endl;
895
896 // now get the path that gives the reduced cost
897
898
899 int currentNode;
900 int successorNode;
901 int lvalue;
902
903 //initialize
904 // we work backwords from bestLastNode
905 //in our recursion we recurse on the currentNode and assume
906 //it is in the optimal path
907
908 //node bestLastNode is a higher number than k
909 *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
910
911
912 //if we are here we should have a value for bestLastNode
913 //if not return infinity
914 if( bestLastNode == OSINT_MAX) return OSDBL_MAX;
915
916 //by successor, I mean node after current node in the path
917 successorNode = k;
918 currentNode = bestLastNode;
919 //the lvalue is the demand through the currentNode
920 lvalue = l ;
921
922
923 while(currentNode != k){
924 //std::cout << "currentNode = " << currentNode << " " << "lvalue " << lvalue << std::endl;
925 if( m_px[ currentNode][ lvalue ] != successorNode){
926
927
928
929 //update nodes
930 successorNode = currentNode;
931 currentNode = m_px[ currentNode][ lvalue ];
932
933
934 if(currentNode - successorNode > 0){
935 //below diagonal
936
937 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
938
939
940 }else{
941 //above diagonal
942
943 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
944
945 }
946
947
948 }else{ //take second best
949
950
951 //update nodes
952 successorNode = currentNode;
953 currentNode = m_tx[ currentNode][ lvalue ];
954
955 if(currentNode - successorNode > 0){
956 //below diagonal
957 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
958
959 }else{
960 //above diagonal
961 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
962
963 }
964
965 }
966
967 //update lvalue
968 lvalue = lvalue - *(m_demand + successorNode);
969
970
971
972 }//end while
973
974
975 //go from node k to bestLastNode -- already done in loop above
976 //*(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1) + currentNode - 1;
977
978
979 return rcost;
980
981}//end qroute
982
983
984
985
986void OSBearcatSolverXkij::getColumns(const double* yA, const int numARows,
987 const double* yB, const int numBRows,
988 int &numNewColumns, int* &numNonzVec, double* &costVec,
989 int** &rowIdxVec, double** &valuesVec, double &lowerBound)
990{
991
992//first strip of the phi dual values and then the convexity row costs
993
994 int i;
995 int j;
996 int numCoulpingConstraints;
997 numCoulpingConstraints = m_numNodes - m_numHubs;
998
999 int numVar;
1000 numVar = m_numNodes*m_numNodes - m_numHubs;
1001 int numNonz;
1002
1003 try{
1004
1005
1006
1007 if(numARows != m_numNodes) throw ErrorClass("inconsistent row count in getColumns");
1008
1009
1010
1011 //get the reduced costs
1012 calcReducedCost( yA, yB );
1013
1014 int kountVar;
1015 double testVal;
1016 testVal = 0;
1017 int k;
1018 int startPntInc;
1019 int rowCount;
1020
1021
1023
1024 double cpuTime;
1025 double start = CoinCpuTime();
1026 getOptL( m_rc);
1027 cpuTime = CoinCpuTime() - start;
1028 std::cout << "DYNAMIC PROGRSMMING CPU TIME " << cpuTime << std::endl;
1029 m_lowerBnd = 0.0;
1030 for(k = 0; k < m_numHubs; k++){
1031
1032
1033 //startPntInc = k*m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes) + (m_optL[ k] - 1)*(m_numNodes*m_numNodes - m_numNodes);
1034 startPntInc = (m_optL[ k] - 1)*(m_numNodes*m_numNodes - m_numNodes);
1035
1036 std::cout << " whichBlock = " << k << " L = " << m_optL[ k] << std::endl;
1037
1038 testVal += m_optL[ k];
1039
1040 kountVar = 0;
1041
1043 m_optValHub[ k] = qrouteCost(k, m_optL[ k], m_rc[ k], &kountVar);
1044
1045 m_optValHub[ k] -= yA[ k + numCoulpingConstraints];
1046
1047 std::cout << "Best Reduced Cost Hub " << k << " = " << m_optValHub[ k] << std::endl;
1048 m_lowerBnd += m_optValHub[ k];
1049
1050 //loop over the rows, scatter each row and figure
1051 //out the column coefficients in this row
1052 //first scatter the sparse array m_varIdx[ j]
1053
1054 m_costVec[ k] = 0.0;
1055
1056 for(j = 0; j < kountVar; j++){
1057
1058
1059 //we are counting the NUMBER of times the variable used
1060 //the same variable can appear more than once in m_varIdx
1061 m_tmpScatterArray[ m_varIdx[ j] - startPntInc ] += 1;
1062
1063 // is variable m_varIdx[ j] - startPntInc in this row
1064
1065 m_costVec[ k] += m_cost[k][ m_varIdx[ j] - startPntInc ];
1066
1067 }
1068
1069
1070
1071 numNonz = 0;
1072 //multiply the sparse array by each A matrix constraint
1073 for(i = 0; i < numCoulpingConstraints; i++){
1074
1075 rowCount = 0;
1076
1077 for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
1078
1079 //m_Amatrix[ j] is a variable index -- this logic works
1080 //since the Amatrix coefficient is 1 -- we don't need a value
1081 //it indexes variable that points into the node
1082 rowCount += m_tmpScatterArray[ m_Amatrix[ j] ];
1083
1084
1085 }
1086
1087 if(rowCount > 0){
1088
1089 //std::cout << "GAIL NODE " << i + m_numHubs << " COUNT = " << rowCount << std::endl;
1090 m_newColumnRowIdx[ k][ numNonz] = i;
1091 m_newColumnRowValue[ k][ numNonz] = rowCount;
1092 numNonz++;
1093 }
1094
1095
1096 }//end loop on coupling constraints
1097
1098 //add a 1 in the convexity row
1099 m_newColumnRowIdx[ k][ numNonz] = m_numNodes - m_numHubs + k;
1100 m_newColumnRowValue[ k][ numNonz++] = 1.0;
1101
1102
1103
1104 //now multiply the sparse array by each B-matrix constraint
1105
1106 for(i = 0; i < m_numBmatrixCon; i++){
1107
1108 rowCount = 0;
1109
1110 for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
1111
1112 //m_Bmatrix[ j] is a variable index -- this logic works
1113 //since the Bmatrix coefficient is 1 -- we don't need a value
1114 //it indexes variable that points into the node
1115 rowCount += m_tmpScatterArray[ m_Bmatrix[ j] ];
1116
1117
1118 }
1119
1120 if(rowCount > 0){
1121
1122
1123 m_newColumnRowIdx[ k][ numNonz] = i + m_numNodes;
1124 m_newColumnRowValue[ k][ numNonz++] = rowCount;
1125
1126 }
1127
1128
1129 }
1130
1131 m_newColumnNonz[ k] = numNonz;
1132
1133
1134
1135 //zero out the scatter vector and store the generated column
1136 for(j = 0; j < kountVar; j++){
1137
1138
1139 m_thetaIndex[ m_numThetaNonz++ ] = m_varIdx[ j] - startPntInc ;
1140 m_tmpScatterArray[ m_varIdx[ j] - startPntInc ] = 0;
1141
1142 // is variable m_varIdx[ j] - startPntInc in this row
1143
1144 }
1145
1146
1147 intVarSet.insert ( std::pair<int,double>( m_numThetaVar, 1.0) );
1149 m_costVec[ k] = m_optL[ k]*m_costVec[ k];
1152
1153
1154
1155
1156
1157 //stuff for debugging
1158 //*****************************//
1191 //**************************//
1192 //end debugging stuff
1193
1194
1195 }//end of loop on hubs
1196
1197
1198
1199 numNonzVec = m_newColumnNonz;
1200 costVec = m_costVec;
1201 rowIdxVec = m_newColumnRowIdx;
1202 valuesVec = m_newColumnRowValue;
1203 std::cout << "Lower Bound = " << m_lowerBnd << std::endl;
1204
1205
1206 if(testVal != m_totalDemand) {
1207
1208 std::cout << "TOTAL DEMAND = " << m_totalDemand << std::endl;
1209 std::cout << "Test Value = " << testVal << std::endl;
1210 throw ErrorClass( "inconsistent demand calculation" );
1211 }
1212
1213
1214
1215
1216
1217
1218 } catch (const ErrorClass& eclass) {
1219
1220 throw ErrorClass(eclass.errormsg);
1221
1222 }
1223
1224
1225 //set method parameters
1226 numNewColumns = m_numHubs;
1227 lowerBound = m_lowerBnd;
1228
1229 std::cout << "LEAVING GET COLUMNS" << std::endl;
1230 return;
1231
1232}//end getColumns
1233
1234
1235
1530
1531
1532 std::cout << "Executing OSBearcatSolverXkij::getInitialRestrictedMaster2( )" << std::endl;
1533
1534 //this master will have m_numNodes artificial variables
1535 int numVarArt;
1536 //numVarArt = 2*m_numNodes;
1537 numVarArt = m_numNodes;
1538
1539 //numVarArt = 0;
1540
1541 // define the classes
1542 FileUtil *fileUtil = NULL;
1543 OSiLReader *osilreader = NULL;
1544 CoinSolver *solver = NULL;
1545 OSInstance *osinstance = NULL;
1546 // end classes
1547
1548 int i;
1549 int j;
1550 int k;
1551 std::string testFileName;
1552 std::string osil;
1553
1554 std::map<int, std::map<int, std::vector<int> > >::iterator mit;
1555 std::map<int, std::vector<int> >::iterator mit2;
1556 std::vector<int>::iterator vit;
1557
1558 //std::vector< int> indexes;
1559 fileUtil = new FileUtil();
1560
1561 m_osinstanceMaster = NULL;
1562
1563 //add linear constraint coefficients
1564 //number of values will nodes.size() the coefficients in the node constraints
1565 //plus coefficients in convexity constraints which is the number of varaibles
1566 int kountNonz;
1567 int kount;
1569 int numThetaVar = m_numberOfSolutions*m_numHubs;
1570 //the extra is for the artificial variables
1571 double *values = new double[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar + numVarArt];
1572 int *indexes = new int[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar + numVarArt] ;
1573 int *starts = new int[ numThetaVar + 1 + numVarArt];
1574 kount = 0;
1575 starts[ 0] = 0;
1576 int startsIdx;
1577 startsIdx = 0;
1578 startsIdx++;
1579
1580 for(i = 0; i < numVarArt; i++){
1582 m_thetaPnt[ m_numThetaVar++] = 0;
1583
1584 }
1585
1586 std::vector<IndexValuePair*> primalValPair;
1587
1588 //
1589 //stuff for cut generation
1590 //
1591 double* xVar;
1592 int numXVar;
1593 numXVar = m_numNodes*(m_numNodes - 1);
1594 xVar = new double[ numXVar];
1595 //zero this array out
1596 for(i = 0; i < numXVar; i++){
1597
1598 xVar[ i] = 0;
1599
1600 }
1601 //getRows function call return parameters
1602 int numNewRows;
1603 int* numRowNonz = NULL;
1604 int** colIdx = NULL;
1605 double** rowValues = NULL ;
1606 double* rowLB;
1607 double* rowUB;
1608 //end of getRows function call return parameters
1609 //
1610 //end of cut generation stuff
1611 //
1612
1613
1614 try {
1615
1616 if(m_initOSiLFile.size() == 0) throw ErrorClass("OSiL file to generate restricted master missing");
1617 osil = fileUtil->getFileAsString( m_initOSiLFile.c_str());
1618
1619 osilreader = new OSiLReader();
1620 osinstance = osilreader->readOSiL(osil);
1621
1622
1623
1624 //turn on or off fixing the variables
1625
1627 //set kount to the start of the z variables
1628 //go past the x variables
1629 //here is where we fix the z variables
1631 //if we are using SK's heuristic fix the assignment of nodes to hubs
1632 if(m_use1OPTstart == true){
1633 osinstance->bVariablesModified = true;
1634 //get the first solution
1635 mit = m_initSolMap.find( 0);
1636 for ( mit2 = mit->second.begin() ; mit2 != mit->second.end(); mit2++ ){ //we are looping over routes in solution mit
1637
1638
1639
1640 for ( vit = mit2->second.begin() ; vit != mit2->second.end(); vit++ ){
1641
1642
1643 osinstance->instanceData->variables->var[ kount + mit2->first*m_numNodes + *vit]->lb = 1.0;
1644 std::cout << "FIXING LOWER BOUND ON VARIABLE " << osinstance->getVariableNames()[ kount + mit2->first*m_numNodes + *vit ] << std::endl;
1645
1646 }
1647
1648 }
1649 }
1650 //*/
1651
1652
1653
1654 //fill in the cost vector first
1655 //the x vector starts at 2*m_numHubs
1656
1657 int idx1;
1658 int idx2;
1659 idx2 = 0; //zRouteDemand have 0 coefficients in obj
1660 //fill in m_cost from the cost coefficient in osil
1661 for(k = 0; k < m_numHubs; k++){
1662
1663 idx1 = 0;
1664
1665 for(i = 0; i < m_numNodes; i++){
1666
1667 for(j = 0; j < i; j++){
1668
1669 m_cost[k][idx1++ ] = osinstance->instanceData->objectives->obj[0]->coef[ idx2++ ]->value;
1670 }
1671
1672 for(j = i + 1; j < m_numNodes; j++){
1673
1674 m_cost[k][idx1++ ] = osinstance->instanceData->objectives->obj[0]->coef[ idx2++ ]->value;
1675
1676 }
1677 }
1678 }
1679
1680
1681
1682 //get variable names for checking purposes
1683 std::string* varNames;
1684 varNames = osinstance->getVariableNames();
1685
1686
1687 //start building the restricted master here
1689 m_osinstanceMaster->setInstanceDescription("The Initial Restricted Master");
1690
1691 // first the variables
1692 // we have numVarArt] artificial variables
1694
1695 // now add the objective function
1697 //add m_numNodes artificial variables
1698 SparseVector *objcoeff = new SparseVector( m_numberOfSolutions*m_numHubs + numVarArt);
1699
1700
1701 // now the constraints
1703
1704
1705
1706 //add the artificial variables -- they must be first in the model
1707
1708 int varNumber;
1709 varNumber = 0;
1710 std::string masterVarName;
1711 kountNonz = 0;
1712
1713 for(i = 0; i < m_numNodes; i++){
1714
1715
1716 objcoeff->indexes[ varNumber ] = varNumber ;
1717 //if obj too large we get the following error
1718 //Assertion failed: (fabs(obj[i]) < 1.0e25), function createRim,
1719 //file ../../../Clp/src/ClpSimplex.cpp, l
1720
1721 //objcoeff->values[ varNumber ] = OSDBL_MAX;
1722
1723 //objcoeff->values[ varNumber ] = 1.0e24;
1724 objcoeff->values[ varNumber ] = m_osDecompParam.artVarCoeff;
1725
1726 m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt2("AP", i ) ,
1727 0, 1.0, 'C');
1728
1729
1730
1731 values[ kountNonz] = 1;
1732 indexes[ kountNonz++] = i ;
1733 starts[ startsIdx++] = kountNonz;
1734
1735
1736
1737 }
1738
1739 /*
1740 for(i = 0; i < m_numNodes; i++){
1741
1742
1743 objcoeff->indexes[ varNumber ] = varNumber ;
1744
1745 //if obj too large we get the following error
1746 //Assertion failed: (fabs(obj[i]) < 1.0e25), function createRim,
1747 //file ../../../Clp/src/ClpSimplex.cpp, l
1748 //objcoeff->values[ varNumber ] = 1.0e24;
1749 //kipp -- change this number -- even 1.0e24 drives
1750 //clp crazy and gives infeasible when actually feasible
1751 objcoeff->values[ varNumber ] =m_osDecompParam.artVarCoeff;
1752
1753
1754
1755 m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt2("AN", i ) ,
1756 0, OSDBL_MAX, 'C');
1757
1758
1759 values[ kountNonz] = -1;
1760 indexes[ kountNonz++] = i ;
1761 starts[ startsIdx++] = kountNonz;
1762
1763
1764
1765 }
1766 */
1767
1768
1769 //
1770 // now get the primal solution
1771 //solve the model for solution in the osoption object
1772
1773
1774 solver = new CoinSolver();
1775 solver->sSolverName ="cbc";
1776 solver->osinstance = osinstance;
1777 solver->buildSolverInstance();
1778 solver->osoption = m_osoption;
1779 OsiSolverInterface *si = solver->osiSolver;
1780
1781 solver->solve();
1782
1783
1784 //get the solver solution status
1785
1786 std::cout << "Solution Status = " << solver->osresult->getSolutionStatusType( 0 ) << std::endl;
1787
1788
1789
1790
1791
1792
1793 //get the optimal objective function value
1794
1795 primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
1796
1797 //loop over routes again to create master objective coefficients
1798 bool isCutAdded;
1799 isCutAdded = true;
1800 while(isCutAdded == true){
1801
1802 isCutAdded = false;
1803
1804 for(k = 0; k < m_numHubs; k++){
1805 numNewRows = 0;
1806
1807 //lets get the x variables
1808 //the variables for this route should be from 2*numHubs + k*(numNodes - 1*)*(numNodes - 1)
1809 idx1 = 2*m_numHubs + k*m_numNodes*(m_numNodes - 1);
1810 idx2 = idx1 + m_numNodes*(m_numNodes - 1);
1811 //end of x variables
1812
1813 std::cout << "HUB " << k << " VARIABLES" << std::endl;
1814
1815 //generate a cut
1816 //need to do index mapping
1817
1818
1819
1820 for(i = idx1; i < idx2; i++){
1821 if( primalValPair[ i]->value > .01 ){
1822 std::cout << osinstance->getVariableNames()[ primalValPair[ i]->idx ] << std::endl;
1823 std::cout << m_variableNames[ primalValPair[ i]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs ] << " " << primalValPair[ i]->value << std::endl;
1824 //m_thetaIndex[ m_numThetaNonz++ ] = primalValPair[ i]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs;
1825
1826 //scatter operation
1827 xVar[ primalValPair[ i]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs ] = primalValPair[ i]->value;
1828
1829
1830 }
1831 }
1832
1833 //get a cut
1834
1835
1836 getCutsX(xVar, numXVar, numNewRows, numRowNonz,
1837 colIdx,rowValues, rowLB, rowUB);
1838
1839
1840
1841 if(numNewRows >= 1){
1842 isCutAdded = true;
1843 std::cout << "WE HAVE A CUT " << std::endl;
1844 std::cout << "EXPRESS CUT IN X(I, J) SPACE" << std::endl;
1845 for(i = 0; i < numRowNonz[ 0]; i++){
1846
1847 std::cout << m_variableNames[ colIdx[0][ i] ] << std::endl;
1848
1849 }
1850
1851
1852 for(i = 0; i < numNewRows; i++){
1853
1854 //set the variable indexes to the correct values
1855
1856 std::cout << "EXPRESS CUT IN X(K, I, J) SPACE" << std::endl;
1857
1858 for(j = 0; j < numRowNonz[ i]; j++){
1859
1860 colIdx[ i][ j] = colIdx[ i][ j] + k*(m_numNodes - 1)*m_numNodes + 2*m_numHubs ;
1861
1862 std::cout << osinstance->getVariableNames()[ colIdx[ i][ j] ] << std::endl;
1863 }
1864
1865 std::cout << "CUT UPPER BOUND = " << rowUB[ i] << std::endl;
1866
1867
1868 si->addRow(numRowNonz[ i], colIdx[ i], rowValues[ i], rowLB[ i], rowUB[ i] ) ;
1869
1870
1871 }
1872
1873
1874
1875 }
1876
1877 //zero out the vector again
1878 for(i = idx1; i < idx2; i++){
1879 if( primalValPair[ i]->value > .01 ){
1880 xVar[ primalValPair[ i]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs ] = 0;
1881 }
1882
1883 }
1884
1885
1886 //std::cout << osinstance->getVariableNames()[ k ] << std::endl;
1887 //std::cout << osinstance->getVariableNames()[ k + m_numHubs ] << std::endl;
1888 std::cout << "Optimal Objective Value = " << primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value << std::endl;
1889
1890
1891
1892 }//end for on k -- hubs
1893
1894
1895 std::cout << std::endl << std::endl;
1896 if( isCutAdded == true) {
1897
1898 std::cout << "A CUT WAS ADDED, CALL SOLVE AGAIN" << std::endl;
1899 solver->solve();
1900 primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
1901 std::cout << "New Solution Status = " << solver->osresult->getSolutionStatusType( 0 ) << std::endl;
1902 std::cout << "Optimal Objective Value = " << solver->osresult->getObjValue(0, 0) << std::endl;
1903 }
1904
1905
1906 }//end while -- we have an integer solution with no subtours
1907
1908
1909
1910
1911 int i1;
1912 int j1;
1913
1914 m_bestIPValue = 0;
1915
1916 for(k = 0; k < m_numHubs; k++){
1917
1918 //lets get the x variables
1919 //the variables for this route should be from 2*numHubs + k*(numNodes - 1*)*(numNodes - 1)
1920 idx1 = 2*m_numHubs + k*m_numNodes*(m_numNodes - 1);
1921 idx2 = idx1 + m_numNodes*(m_numNodes - 1);
1922
1923 //for(i = idx1; i < idx2; i++){
1924
1925 // if( primalValPair[ i]->value > .1 ){
1926 //std::cout << osinstance->getVariableNames()[ primalValPair[ i]->idx ] << std::endl;
1927 //std::cout << m_variableNames[ primalValPair[ i]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs ] << std::endl;
1928 //m_thetaIndex[ m_numThetaNonz++ ] = primalValPair[ i]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs;
1929 //}
1930
1931 //}
1932
1933
1934
1935 for(i1 = 0; i1 < m_numNodes; i1++){
1936
1937
1938 for(j1 = 0; j1 < i1; j1++){
1939
1940 //std::cout << osinstance->getVariableNames()[ primalValPair[ idx1 ]->idx] << " " << primalValPair[ idx1 ]->value << " " << i1 << " " << j1 << std::endl;
1941
1942 if( primalValPair[ idx1 ]->value > .1 ){
1943
1944 m_thetaIndex[ m_numThetaNonz++ ] = primalValPair[ idx1]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs;
1945
1946 // a positive variable pointing into node j1
1947 if(j1 >= m_numHubs){
1948
1949 values[ kountNonz] = 1;
1950 indexes[ kountNonz++] = j1 - m_numHubs ;
1951
1952
1953 }
1954
1955 }
1956
1957 idx1++;
1958 }//end of for on j1
1959
1960
1961 for(j1 = i1 + 1; j1 < m_numNodes; j1++){
1962
1963 //std::cout << osinstance->getVariableNames()[ primalValPair[ idx1 ]->idx] << " " << primalValPair[ idx1 ]->value << " " << i1 << " " << j1 << std::endl;
1964
1965
1966 if( primalValPair[ idx1 ]->value > .1 ){
1967
1968 m_thetaIndex[ m_numThetaNonz++ ] = primalValPair[ idx1]->idx - k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs;
1969
1970
1971 // a positive variable pointing into node j1
1972 if(j1 >= m_numHubs){
1973
1974 values[ kountNonz] = 1;
1975 indexes[ kountNonz++] = j1 - m_numHubs ;
1976
1977 }
1978
1979 }
1980
1981 idx1++;
1982
1983 }//end of for on j1
1984
1985 }//end of for on i1
1986
1987
1988 //now for convexity row k
1989
1990 values[ kountNonz] = 1;
1991 indexes[ kountNonz++] = m_numNodes - m_numHubs + k ;
1992
1993
1994 std::cout << " m_numThetaVar " << m_numThetaVar << std::endl;
1995
1997 m_thetaCost[ m_numThetaVar++ ] = primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value;
1999
2000 masterVarName = makeStringFromInt2("theta(", k);
2001 masterVarName += makeStringFromInt2(",", 0);
2002 masterVarName += ")";
2003 intVarSet.insert ( std::pair<int,double>(varNumber, 1.0) );
2004 m_osinstanceMaster->addVariable(varNumber++, masterVarName, 0, 1, 'C');
2005
2006 std::cout << "Optimal Objective Value = " << primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value << std::endl;
2007
2008 objcoeff->indexes[ k + numVarArt ] = k + numVarArt ;
2009 objcoeff->values[ k + numVarArt ] = primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value;
2010
2011 m_bestIPValue += primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value;
2012
2013 std::cout << "m_bestIPValue = " << m_bestIPValue << std::endl;
2014 starts[ startsIdx++] = kountNonz;
2015
2016 }//end of index on k
2017
2018 //add the constraints
2019 //add the row saying we must visit each node
2020
2021 for( i = 0; i < m_numNodes - m_numHubs ; i++){
2022
2023 m_osinstanceMaster->addConstraint(i, makeStringFromInt2("visitNode_", i + m_numHubs) , 1.0, 1.0, 0);
2024 }
2025
2026 kount = 0;
2027
2028 //add the convexity row
2029 for( i = m_numNodes - m_numHubs; i < m_numNodes ; i++){
2030
2031 m_osinstanceMaster->addConstraint(i, makeStringFromInt2("convexityRowRoute_", kount++ ) , 1.0, 1.0, 0);
2032 }
2033
2034 m_osinstanceMaster->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
2035
2036
2037 //add the linear constraints coefficients
2039 values, 0, kountNonz - 1, indexes, 0, kountNonz - 1, starts, 0, startsIdx);
2040
2041
2042 primalValPair.clear();
2043 delete solver;
2044 solver = NULL;
2045
2046 delete objcoeff;
2047 objcoeff = NULL;
2048 std::cout << m_osinstanceMaster->printModel( ) << std::endl;
2049 std::cout << "NONZ = " << kountNonz << std::endl;
2050
2051 //exit( 1);
2052 //delete[] values;
2053 //values = NULL;
2054 //delete[] starts;
2055 //starts = NULL;
2056 //delete[] indexes;
2057 //indexes = NULL;
2058 delete osilreader;
2059 osilreader = NULL;
2060
2061
2062
2063
2064
2065 } catch (const ErrorClass& eclass) {
2066 std::cout << std::endl << std::endl << std::endl;
2067 if (osilreader != NULL)
2068 delete osilreader;
2069 if (solver != NULL)
2070 delete solver;
2071
2072
2073 // Problem with the parser
2074 throw ErrorClass(eclass.errormsg);
2075 }
2076
2077 delete fileUtil;
2078 fileUtil = NULL;
2079 delete[] xVar;
2080 xVar = NULL;
2081
2082 return m_osinstanceMaster;
2083}//end generateInitialRestrictedMaster2
2084
2085
2086
2088
2089
2090 std::cout << "Executing getOptions(OSOption *osoption)" << std::endl;
2091 //get any options relevant to OSColGenApp
2092 try{
2093
2094 int i;
2095
2096
2097 std::vector<SolverOption*> solverOptions;
2098 std::vector<SolverOption*>::iterator vit;
2099 std::vector<int>::iterator vit2;
2100 std::vector<int >demand;
2101 std::vector<int >routeCapacity;
2102 std::vector<int >routeMinPickup;
2103
2105 solverOptions = osoption->getSolverOptions("routeSolver");
2106 if (solverOptions.size() == 0) throw ErrorClass( "options for routeSolver not available");
2107 //iterate over the vector
2108
2109 int tmpVal;
2110
2111 std::string routeString; //variable for parsing a category option
2112 std::string solutionString; //variable for parsing a category option
2113 string::size_type pos1; //variable for parsing a category option
2114 string::size_type pos2; //variable for parsing a category option
2115 string::size_type pos3; //variable for parsing a category option
2116
2117
2118 std::map<int, std::map<int, std::vector<int> > >::iterator mit;
2119 std::map<int, std::vector<int> >::iterator mit2;
2120 int solutionNumber;
2121 int routeNumber;
2122
2123
2124 for (vit = solverOptions.begin(); vit != solverOptions.end(); vit++) {
2125
2126
2127 //std::cout << (*vit)->name << std::endl;
2128
2129 //(*vit)->name.compare("initialCol") == 0
2130 //if(rowNames[ i3].find("routeCapacity(1)") == string::npos )
2131
2132 if( (*vit)->name.find("numHubs") != std::string::npos){
2133
2134
2135 std::istringstream hubBuffer( (*vit)->value);
2136 hubBuffer >> m_numHubs;
2137 std::cout << "numHubs = " << m_numHubs << std::endl;
2138
2139 }else{
2140
2141 if((*vit)->name.find("numNodes") != std::string::npos){
2142
2143
2144 std::istringstream numNodesBuffer( (*vit)->value);
2145 numNodesBuffer >> m_numNodes;
2146 std::cout << "numNodes = " << m_numNodes << std::endl;
2147
2148 }else{
2149 if((*vit)->name.find("totalDemand") != std::string::npos){
2150
2151
2152 std::istringstream totalDemandBuffer( (*vit)->value);
2153 totalDemandBuffer >> m_totalDemand;
2154 std::cout << "m_totalDemand = " << m_totalDemand << std::endl;
2155
2156 }else{
2157 if((*vit)->name.find("routeMinPickup") != std::string::npos){
2158
2159
2160 std::istringstream routeMinPickupBuffer( (*vit)->value);
2161 routeMinPickupBuffer >> tmpVal;
2162 routeMinPickup.push_back( tmpVal);
2163 //std::cout << "m_minDemand = " << tmpVal << std::endl;
2164
2165 }else{
2166 if( (*vit)->name.find("demand") != std::string::npos ){
2167
2168
2169 std::istringstream demandBuffer( (*vit)->value);
2170 demandBuffer >> tmpVal;
2171 if(tmpVal <= 0 && demand.size() > m_numHubs) throw ErrorClass("must have strictly positive demand");
2172 if(tmpVal < m_minDemand && demand.size() > m_numHubs ) m_minDemand = tmpVal;
2173 demand.push_back( tmpVal);
2174 //std::cout << "demand = " << tmpVal << std::endl;
2175
2176 }else{
2177 if((*vit)->name.find("routeCapacity") != std::string::npos ){
2178 std::istringstream routeCapacityBuffer( (*vit)->value);
2179 routeCapacityBuffer >> tmpVal;
2180 routeCapacity.push_back( tmpVal);
2181 //std::cout << "m_routeCapacity = " << tmpVal << std::endl;
2182
2183 }else{
2184
2185 if((*vit)->name.find("osilFile") != std::string::npos ){
2186 m_initOSiLFile = (*vit)->value;
2187 std::cout << "m_initOSiLFile = " << m_initOSiLFile << std::endl;
2188
2189 }else{
2190
2191 if( (*vit)->name.find("restrictedMasterSolution") != std::string::npos ){
2192 //std::istringstream buffer( (*vit)->value);
2193 //buffer >> m_numberOfSolutions;
2194
2195 //get the block number and solution number
2196 //first get routeString and soluionString
2197 //parse the category string base on :
2198 pos1 = (*vit)->category.find( ":");
2199 if(pos1 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
2200
2201 //solutionString = (*vit)->category.substr( pos1 + 1, pos2 - pos1 - 1);
2202 solutionString = (*vit)->category.substr( 0, pos1);
2203 routeString = (*vit)->category.substr( pos1 + 1);
2204
2205 pos2 = solutionString.find( "n");
2206 if(pos2 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
2207
2208 std::istringstream solutionBuffer( solutionString.substr( pos2 + 1) );
2209 solutionBuffer >> solutionNumber;
2210 //std::cout << "solution number = " << solutionNumber << std::endl;
2211
2212
2213 pos3 = routeString.find( "e");
2214 if(pos3 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
2215 std::istringstream routeBuffer( routeString.substr( pos3 + 1) );
2216 routeBuffer >> routeNumber;
2217 //std::cout << "route number = " << routeNumber << std::endl;
2218 std::istringstream nodeBuffer( (*vit)->value);
2219 nodeBuffer >> tmpVal;
2220
2221 mit = m_initSolMap.find( solutionNumber );
2222
2223 if( mit != m_initSolMap.end() ){
2224 // we have solution from before
2225 // do we have a new route?
2226
2227 mit2 = mit->second.find( routeNumber);
2228
2229 if(mit2 != mit->second.end() ){
2230 // we have a route from before and solution from before
2231
2232
2233 mit2->second.push_back( tmpVal);
2234
2235
2236 }else{
2237
2238 //we have a new route, but old solution
2239 std::vector<int> tmpVec;
2240 tmpVec.push_back( tmpVal) ;
2241 mit->second.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
2242
2243
2244 }
2245
2246 }else{
2247 // we have a new solution
2248 std::vector<int> tmpVec;
2249 tmpVec.push_back( tmpVal) ;
2250
2251 std::map<int, std::vector<int> > tmpMap;
2252 tmpMap.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
2253 m_initSolMap.insert( std::pair<int, std::map<int, std::vector<int> > >(solutionNumber, tmpMap) ) ;
2254
2255 }
2256 }//if on restricted master solution
2257 else{
2258 if( (*vit)->name.find("maxMasterColumns") != std::string::npos){
2259
2260
2261 std::istringstream maxMasterColumns( (*vit)->value);
2262 maxMasterColumns >> m_maxMasterColumns;
2263 std::cout << "m_maxMasterColumn = " << m_maxMasterColumns << std::endl;
2264
2265 }else{
2266 if( (*vit)->name.find("maxThetaNonz") != std::string::npos){
2267
2268 std::istringstream maxThetaNonz( (*vit)->value);
2269 maxThetaNonz >> m_maxThetaNonz;
2270 std::cout << "m_maxThetaNonz = " << m_maxThetaNonz << std::endl;
2271
2272 }else{
2273 if( (*vit)->name.find("use1OPTstart") != std::string::npos){
2274 m_use1OPTstart = false;
2275 if ( (*vit)->value.find("true") != std::string::npos ) m_use1OPTstart = true;
2276 std::cout << "m_use1OPTstart = " << m_use1OPTstart << std::endl;
2277
2278 }else{
2279 if( (*vit)->name.find("maxBmatrixCon") != std::string::npos ){
2280
2281 std::istringstream maxBmatrixCon( (*vit)->value);
2282 maxBmatrixCon >> m_maxBmatrixCon;
2283 std::cout << "m_maxBmatrixCon = " << m_maxBmatrixCon << std::endl;
2284
2285 }else{
2286 if( (*vit)->name.find("maxBmatrixNonz") != std::string::npos ){
2287
2288 std::istringstream maxBmatrixNonz( (*vit)->value);
2289 maxBmatrixNonz >> m_maxBmatrixNonz;
2290 std::cout << "m_maxBmatrixNonz = " << m_maxBmatrixNonz << std::endl;
2291
2292
2293 }
2294 }
2295 }
2296 }
2297 }
2298 }
2299 }
2300 }
2301 }
2302 }
2303 }
2304 }
2305 }//end if on solver options
2306
2307 }//end for loop on options
2308
2309
2310 //now fill in route capacities
2311 i = 0;
2312 m_routeCapacity = new int[ m_numHubs];
2313 if(m_numHubs != routeCapacity.size( ) ) throw ErrorClass("inconsistent number of HUBS");
2314 for (vit2 = routeCapacity.begin(); vit2 != routeCapacity.end(); vit2++) {
2315
2316 *(m_routeCapacity + i++) = *vit2;
2317
2318 }
2319 routeCapacity.clear();
2320
2321
2322 //now fill in route min pickups
2323 i = 0;
2324 m_routeMinPickup = new int[ m_numHubs];
2325 if(m_numHubs != routeMinPickup.size( ) ) throw ErrorClass("inconsistent number of HUBS");
2326 for (vit2 = routeMinPickup.begin(); vit2 != routeMinPickup.end(); vit2++) {
2327
2328 *(m_routeMinPickup + i++) = *vit2;
2329
2330 }
2331 routeMinPickup.clear();
2332
2333
2334
2335
2336 //now fill in demand
2337 i = 0;
2338 m_demand = new int[ m_numNodes];
2339 if(m_numNodes != demand.size( ) ) throw ErrorClass("inconsistent number of demand nodes");
2340 for (vit2 = demand.begin(); vit2 != demand.end(); vit2++) {
2341
2342 *(m_demand + i++) = *vit2;
2343
2344 }
2345 demand.clear();
2346
2347 //kipp -- fill in numberOfRestricedMasterSolutions from map size
2349
2350
2351 } catch (const ErrorClass& eclass) {
2352
2353 throw ErrorClass(eclass.errormsg);
2354
2355 }
2356
2357}//end getOptions
2358
2359
2360
2361void OSBearcatSolverXkij::getCutsTheta(const double* theta, const int numTheta,
2362 int &numNewRows, int* &numNonz, int** &colIdx,
2363 double** &values, double* &rowLB, double* &rowUB) {
2364 //critical -- the variables that come in the theta variables
2365 //not the x variables, we must convert to x, find a cut in x-space
2366 //and then convert back to theta
2367
2368 int i;
2369 int j;
2370 int k;
2371 int index;
2372 int rowKount;
2373 int tmpKount;
2374 int indexAdjust = m_numNodes - m_numHubs;
2375 double* tmpRhs;
2376 int numSepRows = m_osinstanceSeparation->getConstraintNumber() ;
2377
2378 tmpRhs = new double[ numSepRows ];
2379
2380 for(i = 0; i < numSepRows; i++){
2381
2382 tmpRhs[ i] = 0;
2383 }
2384
2385 try{
2387 //m_numNodes is the number of artificial variables
2388 if(numTheta != m_numThetaVar ) throw
2389 ErrorClass("number of master varibles in OSBearcatSolverXkij::getCuts inconsistent");
2390
2391 //for(i = 0; i < numTheta; i++){
2392
2393 //std::cout << "numTheta = " << numTheta << std::endl;
2394 //std::cout << "m_numThetaVar = " << m_numThetaVar - 1 << std::endl;
2395
2396 //exit( 1);
2397
2398 for(i = 0; i < numTheta; i++){
2399
2400 //get a postive theta
2401 if(theta[ i] > m_osDecompParam.zeroTol){
2402
2403 //get the xij indexes associated with this variable
2404 for(j = m_thetaPnt[ i ]; j < m_thetaPnt[ i + 1 ]; j++ ){
2405
2406 //get the xij index
2407
2408
2409
2410 rowKount = m_separationIndexMap[ m_thetaIndex[ j] ];
2411
2412 //std::cout << "rowKount = " << rowKount <<std::endl;
2413
2414 if(rowKount < OSINT_MAX ){
2415
2416 tmpRhs[ rowKount] -= theta[ i];
2417
2418 }
2419
2420 }
2421 }
2422 }
2423
2424
2425 // don't adjust the kludge row
2426
2427 for(i = indexAdjust; i < numSepRows - 1; i++){
2428
2429 if(-tmpRhs[ i] > 1 + m_osDecompParam.zeroTol ){
2430 // quick and dirty does x_{ij} + x_{ji} <= 1 cut off anything
2431 //std::cout << " tmpRhs[ i] = " << tmpRhs[ i] << std::endl;
2432 //which variable is this
2433 //kipp this an inefficient way of finding i and j -- improve this
2434 int tmpKount = indexAdjust;
2435 for(int i1 = m_numHubs; i1 < m_numNodes; i1++){
2436
2437
2438
2439 for(int j1 = i1+1; j1 < m_numNodes; j1++){
2440
2441 if(tmpKount == i){
2442
2443 //std::cout << "i = " << i1 << std::endl;
2444 //std::cout << "j = " << j1 << std::endl;
2445 //okay generate a cut that says
2446 // x(i1,j1) + x(j1, i1) << 1
2447 //get index for i1,j1
2448 m_Bmatrix[ m_numBmatrixNonz++ ] = i1*(m_numNodes - 1) + j1 - 1 ;
2449 //get index for j1,i1
2450 m_Bmatrix[ m_numBmatrixNonz++ ] = j1*(m_numNodes - 1) + i1 ;
2453
2454 numNewRows = 1;
2455
2456 m_newRowNonz[ 0] = 0;
2457 m_newRowUB[ 0] = 1;
2458 m_newRowLB[ 0] = 0;
2459
2460 //now we have to get the theta column indexes
2461 //scatter the constraint in the x - variables
2462
2463 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2464 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2465
2466
2467 std::cout << " m_Bmatrix[ j] " << m_Bmatrix[ j] << std::endl ;
2468
2469 m_tmpScatterArray[ m_Bmatrix[ j] ] = 1;
2470
2471 }
2472
2473
2474
2475
2476 for(k = 0; k < m_numThetaVar ; k++){
2477
2478 //get the xij indexes in this column
2479 tmpKount = 0;
2480
2481
2482 for(j = m_thetaPnt[k]; j < m_thetaPnt[k + 1] ; j++){
2483
2484 if(m_tmpScatterArray[ m_thetaIndex[ j] ] > 0 ){ //count number of xij for theta_i
2485
2486 std::cout << " Variable " << m_variableNames[ m_thetaIndex[ j] ] << std::endl;
2487
2488 tmpKount++;
2489
2490 }
2491
2492 }//end loop on j
2493
2494 if(tmpKount > 0){
2495 //theta_i has a nonzero coefficient in this row
2496
2497 m_newRowColumnIdx[0][ m_newRowNonz[ 0] ] = k ;
2498 //m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
2499 m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
2500
2501
2502 }
2503
2504 }//end loop on k
2505
2506
2507 //zero out the scatter array again
2508
2509 //::cout << " m_numBmatrixCon - 1 " << m_numBmatrixCon - 1 << std::endl;
2510 //std::cout << " m_pntBmatrix[ m_numBmatrixCon - 1] " << m_pntBmatrix[ m_numBmatrixCon - 1] << std::endl ;
2511 //std::cout << " m_pntBmatrix[ m_numBmatrixCon ] " << m_pntBmatrix[ m_numBmatrixCon ] << std::endl ;
2512 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2513 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2514
2515 m_tmpScatterArray[ m_Bmatrix[ j] ] = 0;
2516
2517 }
2518
2519 numNonz = m_newRowNonz;
2520 colIdx = m_newRowColumnIdx;
2521 values = m_newRowColumnValue;
2522 rowUB = m_newRowUB;
2523 rowLB = m_newRowLB;
2524
2525 delete[] tmpRhs;
2526 tmpRhs = NULL;
2527
2528 //we found a row, add the corresponding artificial variables
2529 //to the transformation matrix
2530 m_numThetaVar++;
2533 //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz;
2534 //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz;
2535
2536 return;
2537
2538 } //end loop on if tmpKount
2539
2540 tmpKount++;
2541
2542 }//loop on j1
2543
2544 }//loop on i1
2545
2546
2547 }//end if on tmpRHS
2548
2549 m_separationClpModel->setRowUpper(i, tmpRhs[ i] );
2550 m_separationClpModel->setRowLower(i, tmpRhs[ i] );
2551
2552 }//end loop on i
2553
2554
2555 //std::cout << m_osinstanceSeparation->printModel() << std::endl;
2556
2557
2558 std::vector<int> dualIdx;
2559 std::vector<int>::iterator vit1;
2560 std::vector<int>::iterator vit2;
2561
2562 //if the objective function value is greater than zero we have a cut
2563 //the cut is based on the nodes with dual value - 1
2564
2565 for(k = 0; k < indexAdjust; k++){
2566 //std::cout << std::endl << std::endl;
2567
2568
2569 m_separationClpModel->setRowUpper(k, 0.0);
2570 //we don't need output
2571
2572 m_separationClpModel->setLogLevel( 0);
2573
2574 m_separationClpModel->primal();
2575
2576 if(m_separationClpModel->getObjValue() > m_osDecompParam.zeroTol){
2577 std::cout << "DOING SEPARATION FOR NODE " << k + m_numHubs << std::endl;
2578 std::cout << "SEPERATION OBJ VALUE = " << m_separationClpModel->getObjValue() << std::endl;
2579 numNewRows = 1;
2580
2581 for(i = 0; i < m_numNodes - m_numHubs ; i++){
2582 //std::cout << m_osinstanceSeparation->getConstraintNames()[ i] << " = " << m_separationClpModel->getRowPrice()[ i] << std::endl;
2583 if( m_separationClpModel->getRowPrice()[ i] - m_osDecompParam.zeroTol <= -1) dualIdx.push_back( i) ;
2584 }
2585
2586 for (vit1 = dualIdx.begin(); vit1 != dualIdx.end(); vit1++) {
2587
2588 i = *vit1 + m_numHubs;
2589
2590 for (vit2 = dualIdx.begin(); vit2 != dualIdx.end(); vit2++) {
2591
2592 j = *vit2 + m_numHubs;
2593
2594 if( i > j ){
2595
2596 index = i*(m_numNodes -1) + j;
2597 std::cout << "CUT VARIABLE = " << m_variableNames[ index ] <<std::endl;
2598 m_Bmatrix[ m_numBmatrixNonz++ ] = index ;
2599
2600 }else{
2601
2602 if( i < j ){
2603
2604 index = i*(m_numNodes -1) + j - 1;
2605 std::cout << "CUT VARIABLE = " << m_variableNames[ index ] <<std::endl;
2606 m_Bmatrix[ m_numBmatrixNonz++ ] = index ;
2607
2608 }
2609 }
2610
2611 }//end for on vit2
2612 }//end for on vit1
2613
2614 //add the tour-breaking cut
2617
2618 // multiply the transformation matrix times this cut to get the cut in theta space
2619 // do the usual trick and scatter m_Bmatrix into a dense vector
2620
2621 //reset
2622 // don't adjust the kludge row
2623 for(i = indexAdjust; i < numSepRows - 1; i++){
2624
2625 m_separationClpModel->setRowUpper(i, 0.0 );
2626 m_separationClpModel->setRowLower(i, 0.0 );
2627
2628
2629 }
2630 m_separationClpModel->setRowUpper(k, 1.0);
2631 delete[] tmpRhs;
2632 tmpRhs = NULL;
2633
2634
2635 m_newRowNonz[ 0] = 0;
2636 m_newRowUB[ 0] = dualIdx.size() - 1;
2637 m_newRowLB[ 0] = 0;
2638
2639 dualIdx.clear();
2640
2641 //now we have to get the theata column indexes
2642 //scatter the constraint in the x - variables
2643
2644 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2645 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2646
2647 m_tmpScatterArray[ m_Bmatrix[ j] ] = 1;
2648
2649 }
2650
2651
2652
2653
2654 for(i = 0; i < m_numThetaVar ; i++){
2655
2656 //get the xij indexes in this column
2657 tmpKount = 0;
2658 for(j = m_thetaPnt[i]; j < m_thetaPnt[i + 1] ; j++){
2659
2660 if(m_tmpScatterArray[ m_thetaIndex[ j] ] > 0 ){ //count number of xij for theta_i
2661
2662 tmpKount++;
2663
2664 }
2665
2666 }//end loop on j
2667
2668 if(tmpKount > 0){
2669 //theta_i has a nonzero coefficient in this row
2670
2671 m_newRowColumnIdx[0][ m_newRowNonz[ 0] ] = i ;
2672
2673 m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
2674
2675
2676 }
2677
2678 }//end loop on i
2679
2680
2681 //zero out the scatter array again
2682
2683 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
2684 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
2685
2686 m_tmpScatterArray[ m_Bmatrix[ j] ] = 0;
2687
2688 }
2689
2690
2691
2692 numNonz = m_newRowNonz;
2693 colIdx = m_newRowColumnIdx;
2694 values = m_newRowColumnValue;
2695 rowUB = m_newRowUB;
2696 rowLB = m_newRowLB;
2697 m_numThetaVar++;
2700 //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz; //first artificial variable
2701 //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz; //second artificial varaible
2702
2703 return;
2704
2705
2706
2707 }//end if on obj value
2708 m_separationClpModel->setRowUpper(k, 1.0);
2709 dualIdx.clear();
2710
2711 }//end loop on k
2712
2713 //if we are here there was no cut
2714 //reset
2715 // don't adjust the kludge row
2716 for(i = indexAdjust; i < numSepRows - 1; i++){
2717
2718 m_separationClpModel->setRowUpper(i, 0.0 );
2719 m_separationClpModel->setRowLower(i, 0.0 );
2720
2721
2722 }
2723
2724 delete[] tmpRhs;
2725 tmpRhs = NULL;
2726
2727 } catch (const ErrorClass& eclass) {
2728
2729 throw ErrorClass(eclass.errormsg);
2730
2731 }
2732
2733
2734
2735}//end getCutsTheta
2736
2737
2738
2739
2740
2741
2742
2743void OSBearcatSolverXkij::getCutsX(const double* x, const int numX,
2744 int &numNewRows, int* &numNonz, int** &colIdx,
2745 double** &values, double* &rowLB, double* &rowUB) {
2746 //critical -- we are assuming that the size of x is going to be
2747 // m_numNodes*(m_numNodes - 1)
2748
2749 int i;
2750 int j;
2751 int k;
2752 int index;
2753 int rowKount;
2754
2755
2756 int indexAdjust = m_numNodes - m_numHubs;
2757 double* tmpRhs;
2758 int numSepRows = m_osinstanceSeparation->getConstraintNumber() ;
2759
2760 tmpRhs = new double[ numSepRows ];
2761
2762 for(i = 0; i < numSepRows; i++){
2763
2764 tmpRhs[ i] = 0;
2765 }
2766
2767 try{
2769
2770 for(i = 0; i < numX; i++){
2771
2772 //get a postive theta
2773 if(x[ i] > m_osDecompParam.zeroTol){
2774
2775 //the row index for x_{ij}
2776 rowKount = m_separationIndexMap[ i ];
2777
2778 if(rowKount < OSINT_MAX ){
2779
2780 tmpRhs[ rowKount] -= x[ i];
2781
2782 }
2783
2784 }
2785 }// end i loop
2786
2787 for(i = indexAdjust; i < numSepRows - 1; i++){
2788
2789 if(-tmpRhs[ i] > 1 + m_osDecompParam.zeroTol ){
2790
2791 // quick and dirty does x_{ij} + x_{ji} <= 1 cut off anything
2792 //std::cout << " tmpRhs[ i] = " << tmpRhs[ i] << std::endl;
2793 //which variable is this
2794 //kipp this an inefficient way of finding i and j -- improve this
2795 int tmpKount = indexAdjust;
2796 for(int i1 = m_numHubs; i1 < m_numNodes; i1++){
2797
2798 for(int j1 = i1+1; j1 < m_numNodes; j1++){
2799
2800 if(tmpKount == i){
2801
2802 numNewRows = 1;
2803
2804 m_newRowNonz[ 0] = 2;
2805 m_newRowUB[ 0] = 1;
2806 m_newRowLB[ 0] = 0;
2807
2808 m_newRowColumnIdx[ 0][ 0 ] = i1*(m_numNodes - 1) + j1 - 1;
2809 m_newRowColumnIdx[ 0][ 1 ] = j1*(m_numNodes - 1) + i1;
2810 m_newRowColumnValue[ 0][ 0] = 1;
2811 m_newRowColumnValue[ 0][ 1] = 1;
2812
2813 numNonz = m_newRowNonz;
2814 colIdx = m_newRowColumnIdx;
2815 values = m_newRowColumnValue;
2816 rowUB = m_newRowUB;
2817 rowLB = m_newRowLB;
2818
2819 delete[] tmpRhs;
2820 tmpRhs = NULL;
2821 return;
2822
2823
2824
2825 }
2826
2827 tmpKount++;
2828
2829 }// end loop on j1
2830
2831 }//end loop on i1
2832
2833
2834 }//end if on tmpRHS
2835
2836 m_separationClpModel->setRowUpper(i, tmpRhs[ i] );
2837 m_separationClpModel->setRowLower(i, tmpRhs[ i] );
2838
2839 }//end loop on i
2840
2841
2842 //std::cout << m_osinstanceSeparation->printModel() << std::endl;
2843
2844
2845 std::vector<int> dualIdx;
2846 std::vector<int>::iterator vit1;
2847 std::vector<int>::iterator vit2;
2848
2849 //if the objective function value is greater than zero we have a cut
2850 //the cut is based on the nodes with dual value - 1
2851
2852 for(k = 0; k < indexAdjust; k++){
2853 std::cout << std::endl << std::endl;
2854
2855
2856 m_separationClpModel->setRowUpper(k, 0.0);
2857
2858
2859 m_separationClpModel->primal();
2860
2861 if(m_separationClpModel->getObjValue() > m_osDecompParam.zeroTol){
2862 std::cout << "DOING SEPARATION FOR NODE " << k + m_numHubs << std::endl;
2863 std::cout << "SEPERATION OBJ = " << m_separationClpModel->getObjValue() << std::endl;
2864 numNewRows = 1;
2865 m_newRowNonz[ 0] = 0;
2866 m_newRowLB[ 0] = 0;
2867
2868 for(i = 0; i < m_numNodes - m_numHubs ; i++){
2869 //std::cout << m_osinstanceSeparation->getConstraintNames()[ i] << " = " << m_separationClpModel->getRowPrice()[ i] << std::endl;
2870 if( m_separationClpModel->getRowPrice()[ i] - m_osDecompParam.zeroTol <= -1) dualIdx.push_back( i) ;
2871 }
2872
2873 for (vit1 = dualIdx.begin(); vit1 != dualIdx.end(); vit1++) {
2874
2875 i = *vit1 + m_numHubs;
2876
2877 for (vit2 = dualIdx.begin(); vit2 != dualIdx.end(); vit2++) {
2878
2879 j = *vit2 + m_numHubs;
2880
2881 if( i > j ){
2882
2883 index = i*(m_numNodes -1) + j;
2884 std::cout << "CUT VARIABLE = " << m_variableNames[ index] <<std::endl;
2885 m_newRowColumnValue[ 0][ m_newRowNonz[ 0] ] = 1.0;
2886 m_newRowColumnIdx[ 0][ m_newRowNonz[ 0]++ ] = index;
2887
2888 }else{
2889
2890 if( i < j ){
2891
2892 index = i*(m_numNodes -1) + j - 1;
2893 std::cout << "CUT VARIABLE = " << m_variableNames[ index] <<std::endl;
2894 m_newRowColumnValue[ 0][ m_newRowNonz[ 0] ] = 1.0;
2895 m_newRowColumnIdx[ 0][ m_newRowNonz[ 0]++ ] = index;
2896
2897 }
2898 }
2899
2900 }//end for on vit2
2901 }//end for on vit1
2902
2903
2904 m_newRowUB[ 0] = dualIdx.size() - 1;
2905
2906 dualIdx.clear();
2907 //reset
2908 // don't adjust the kludge row
2909 for(i = indexAdjust; i < numSepRows - 1; i++){
2910
2911 m_separationClpModel->setRowUpper(i, 0.0 );
2912 m_separationClpModel->setRowLower(i, 0.0 );
2913
2914
2915 }
2916 m_separationClpModel->setRowUpper(k, 1.0);
2917 delete[] tmpRhs;
2918 tmpRhs = NULL;
2919
2920
2921 numNonz = m_newRowNonz;
2922 colIdx = m_newRowColumnIdx;
2923 values = m_newRowColumnValue;
2924 rowUB = m_newRowUB;
2925 rowLB = m_newRowLB;
2926
2927 return;
2928
2929
2930
2931 }//end if on obj value
2932 m_separationClpModel->setRowUpper(k, 1.0);
2933 dualIdx.clear();
2934
2935 }//end loop on k
2936
2937 //if we are here there was no cut
2938 //reset
2939 // don't adjust the kludge row
2940 for(i = indexAdjust; i < numSepRows - 1; i++){
2941
2942 m_separationClpModel->setRowUpper(i, 0.0 );
2943 m_separationClpModel->setRowLower(i, 0.0 );
2944
2945
2946 }
2947
2948 delete[] tmpRhs;
2949 tmpRhs = NULL;
2950
2951 } catch (const ErrorClass& eclass) {
2952
2953 throw ErrorClass(eclass.errormsg);
2954
2955 }
2956
2957
2958}//end getCutsX
2959
2960
2961void OSBearcatSolverXkij::calcReducedCost( const double* yA, const double* yB){
2962
2963 int k;
2964 int i;
2965 int j;
2966 int l;
2967 int kount;
2968
2969 int tmpVal;
2970 tmpVal = m_numNodes - 1;
2971
2972 for(k = 0; k < m_numHubs; k++){
2973 kount = 0;
2974
2975 for(l = 1; l <= m_upperBoundL[ k]; l++){
2976
2977
2978 for(i = 0; i< m_numNodes; i++){
2979
2980
2981
2982 for(j = 0; j < i; j++){
2983
2984 if(j < m_numHubs){
2985
2986 m_rc[k][ kount++] = l*m_cost[k][ i*tmpVal + j ] ;
2987
2988 }else{
2989
2990 m_rc[k][ kount++] = l*m_cost[k][ i*tmpVal + j ] - yA[ j - m_numHubs] ;
2991 }
2992
2993
2994 }
2995
2996
2997
2998 for(j = i + 1; j < m_numNodes; j++){
2999
3000
3001 if(j < m_numHubs){
3002
3003 m_rc[k][ kount++] = l*m_cost[k][ i*tmpVal + j - 1 ];
3004
3005 } else {
3006
3007
3008 m_rc[k][ kount++] = l*m_cost[k][ i*tmpVal + j - 1 ] - yA[ j - m_numHubs ];
3009
3010 }
3011
3012 }
3013
3014
3015 }
3016
3017
3018 }//end l loop
3019
3020
3021 }//end hub loop
3022
3023 //now adjust for tour breaking constraints
3024
3025
3026 int startPnt ;
3027
3028 for(i = 0; i < m_numBmatrixCon; i++){
3029
3030 //get the xij
3031
3032 for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
3033
3034
3035
3036 for(k = 0; k < m_numHubs; k++){
3037
3038
3039 for(l = 1; l <= m_upperBoundL[ k]; l++){
3040
3041 //startPnt = k*m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes) + (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
3042 startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
3043
3044 m_rc[ k][ startPnt + m_Bmatrix[ j] ] -= yB[ i];
3045
3046 }
3047
3048 }
3049
3050
3051 }
3052
3053 }
3054
3055}//end calcReducedCost
3056
3057
3059
3060 int i;
3061 int j;
3062 int kount;
3063
3064 kount = 0;
3065
3066 for(i = 0; i< m_numNodes; i++){
3067
3068 //if we have (i, j) where j is hub then do not subtract off phi[ j]
3069 for(j = 0; j < i; j++){
3070
3071 m_variableNames[ kount] = makeStringFromInt2("x(" , i);
3072 m_variableNames[ kount] += makeStringFromInt2( "," , j);
3073 m_variableNames[ kount] += ")";
3074 //std::cout << "GAIL VARIABLE NAME " << m_variableNames[ kount] << std::endl;
3075
3076 kount++;
3077
3078 }
3079
3080 for(j = i + 1; j < m_numNodes; j++){
3081
3082 m_variableNames[ kount] = makeStringFromInt2("x(" , i);
3083 m_variableNames[ kount] += makeStringFromInt2( "," , j);
3084 m_variableNames[ kount] += ")";
3085
3086 //std::cout << "GAIL VARIABLE NAME " << m_variableNames[ kount] << std::endl;
3087 kount++;
3088
3089 }
3090
3091
3092 }
3093}//end createVariableNames
3094
3096
3097 //arrays for the coupling constraint matrix
3098 //this is in the x variable space, not theta
3099 //int* m_pntAmatrix;
3100 //int* m_Amatrix;
3101
3102
3103 int i;
3104 int j;
3105 int numNonz;
3106
3107 //loop over nodes
3108 m_pntAmatrix[ 0] = 0;
3109 numNonz = 0;
3110
3111 for(j = m_numHubs; j < m_numNodes; j++){
3112
3113
3114 for(i = 0; i < j; i++){
3115
3116 m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j - 1 ;
3117
3118 }
3119
3120 for(i = j + 1; i < m_numNodes; i++){
3121
3122 m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j ;
3123
3124 }
3125
3126 m_pntAmatrix[ j - m_numHubs + 1] = numNonz;
3127
3128 }
3129
3130 /*
3131 for(i = 0; i < m_numNodes - m_numHubs; i++){
3132
3133 for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
3134
3135 std::cout << m_variableNames[ m_Amatrix[ j ] ] << std::endl;
3136
3137 }
3138
3139 }
3140 */
3141
3142}//end createAmatrix
3143
3144void OSBearcatSolverXkij::pauHana( std::vector<int> &m_zOptIndexes, int numNodes, int numColsGen){
3145
3146 std::cout << std::endl;
3147 std::cout << " PAU HANA TIME! " << std::endl;
3148 double cost;
3149 cost = 0;
3150 std::vector<int>::iterator vit;
3151 try{
3152 int i;
3153 int j;
3154
3155
3156 //we better NOT have any artifical variables positive
3157 //for(i = 0; i < numVarArt ; i++){
3158 //
3159 // if(theta[ i] > m_osDecompParam.zeroTol) throw ErrorClass("we have a positive artificial variable");
3160 //}
3161
3162 //for(i = 0; i < m_numThetaVar ; i++){
3163
3164 //cost += theta[ i]*m_thetaCost[ i ];
3165 //std::cout << "COLUMN VALUE = " << theta[ i] << std::endl;
3166
3167 //}
3168
3169
3170 for(vit = m_zOptIndexes.begin() ; vit != m_zOptIndexes.end(); vit++){
3171
3172 i = *vit;
3173 std::cout << "x variables for column " << i << std::endl;
3174
3175 cost += m_thetaCost[ i ];
3176
3177 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
3178
3179 std::cout << "INDEX = " << m_thetaIndex[ j] << std::endl;
3180 std::cout << m_variableNames[ m_thetaIndex[ j] ] << " = " << 1 << std::endl;
3181
3182 }
3183
3184 }
3185
3186
3187 std::cout << "cost = " << cost << std::endl << std::endl;
3188
3189 std::cout << std::endl << std::endl;
3190 std::cout << "LOWER BOUND VALUE = " << m_bestLPValue << std::endl;
3191 std::cout << "FINAL BEST IP SOLUTION VALUE = " << m_bestIPValue << std::endl;
3192 std::cout << "NUMBER OF COLUMNS IN FINAL MASTER = " << m_numThetaVar << std::endl;
3193 //std::cout << "NUMBER OF GENERATED COLUMNS = " << m_numThetaVar - 2*m_numNodes - 2*m_numBmatrixCon << std::endl;
3194 //the original master has m_numHubs + m_numNodes columns
3195 std::cout << "NUMBER OF GENERATED COLUMNS = " << numColsGen << std::endl;
3196 std::cout << "NUMBER OF GENERATED CUTS = " << m_numBmatrixCon << std::endl;
3197 std::cout << "NUMBER OF NODES = " << numNodes << std::endl;
3198 std::cout << " PAU!!!" << std::endl;
3199
3200 std::cout << std::endl << std::endl;
3201
3202
3203
3204
3205 std::cout << std::endl << std::endl;
3206 }catch (const ErrorClass& eclass) {
3207
3208 throw ErrorClass(eclass.errormsg);
3209
3210 }
3211
3212}//end pauHana -- no pun intended
3213
3214
3216
3217
3218
3219
3221
3222 //add linear constraint coefficients
3223 //number of values will nodes.size() the coefficients in the node constraints
3224 //plus coefficients in convexity constraints which is the number of varaibles
3225 int kountNonz;
3226 int kount;
3227 int startsIdx;
3228 //we build these on nodes that do not include the hubs
3229 int numYvar = (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1);
3230 int numVvar = m_numNodes - m_numHubs;
3231 //the plus 1 is for the kludge row
3232 int numCon = (m_numNodes - m_numHubs) + (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1)/2 + 1;
3233 double *values = new double[ 2*numYvar + 2*numVvar];
3234 int *indexes = new int[ 2*numYvar + 2*numVvar];
3235 int *starts = new int[ numYvar + numVvar + 1];
3236 starts[ 0] = 0;
3237 startsIdx = 0;
3238 startsIdx++;
3239 kountNonz = 0;
3240 int i;
3241 int j;
3242
3243
3244 std::string separationVarName;
3245 std::string separationConName;
3246
3247 try {
3248
3250
3251 //start building the separation instance
3252
3253 m_osinstanceSeparation->setInstanceDescription("The Tour Breaking Separation Problem");
3254
3255
3256 // now the constraints
3258
3259
3260 //add the node rows
3261 for( i = 0; i < m_numNodes - m_numHubs ; i++){
3262
3263 m_osinstanceSeparation->addConstraint(i, makeStringFromInt2("nodeRow_", i+ m_numHubs ) , 0.0, 1.0, 0);
3264
3265 }
3266
3267 //add the variable rows rows
3268
3269 int rowKounter;
3270 rowKounter = m_numNodes - m_numHubs;
3271
3272 for(i = m_numHubs; i < m_numNodes; i++){
3273
3274
3275
3276 for(j = i+1; j < m_numNodes; j++){
3277 separationConName = makeStringFromInt2("Row_(", i);
3278 separationConName += makeStringFromInt2(",", j);
3279 separationConName += ")";
3280
3281 m_osinstanceSeparation->addConstraint(rowKounter++, separationConName , 0, 0, 0);
3282 }
3283
3284 }
3285
3286 // the klude row so we have +/-1 in every column
3287 m_osinstanceSeparation->addConstraint(rowKounter++, "kludgeRow" , 0, m_numNodes, 0);
3288
3289 // the variables
3290 m_osinstanceSeparation->setVariableNumber( numYvar + numVvar);
3291
3292
3293
3294 std::cout << "NUMBER OF VARIABLES SET = " << numYvar + numVvar << std::endl;
3295 //add the v variables
3296 for(i = 0; i < numVvar; i++){
3297
3298 separationVarName = makeStringFromInt2("v", i + m_numHubs);
3299
3300 m_osinstanceSeparation->addVariable(i, separationVarName, 0, 1, 'C');
3301
3302 values[ kountNonz ] = -1.0;
3303 indexes[ kountNonz ] = i;
3304 kountNonz++;
3305
3306 values[ kountNonz ] = 1.0;
3307 indexes[ kountNonz ] = rowKounter - 1;
3308 kountNonz++;
3309
3310
3311
3312 starts[ startsIdx++ ] = kountNonz;
3313
3314
3315 }
3316 //add the y variables
3317 kount = numVvar;
3318
3319 int i1;
3320 int j1;
3321 int kountCon;
3322 kountCon = m_numNodes - m_numHubs;
3323
3324 for(i1 = 0; i1 < m_numNodes - m_numHubs; i1++){
3325
3326
3327 i = i1 + m_numHubs;
3328
3329 for(j1 = i1 + 1; j1 < m_numNodes - m_numHubs; j1++){
3330
3331
3332 j = j1 + m_numHubs;
3333
3334 separationVarName = makeStringFromInt2("y(", i);
3335 separationVarName += makeStringFromInt2(",", j);
3336 separationVarName += ")";
3337 m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
3338
3339 //map the variable to row kountCon
3340
3341 // i < j case
3342 m_separationIndexMap[ i*(m_numNodes - 1) + (j - 1) ] = kountCon;
3343
3344 values[ kountNonz ] = 1.0;
3345 indexes[ kountNonz ] = i1;
3346 kountNonz++;
3347
3348 values[ kountNonz ] = -1.0;
3349 indexes[ kountNonz ] = kountCon ;
3350 kountNonz++;
3351
3352 starts[ startsIdx++ ] = kountNonz;
3353
3354
3355
3356
3357 separationVarName = makeStringFromInt2("y(", j );
3358 separationVarName += makeStringFromInt2(",", i);
3359 separationVarName += ")";
3360 m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
3361
3362 values[ kountNonz ] = 1.0;
3363 indexes[ kountNonz ] = j1;
3364 kountNonz++;
3365
3366 // i < j case
3367 m_separationIndexMap[ j*(m_numNodes - 1) + i ] = kountCon;
3368
3369 values[ kountNonz ] = -1.0;
3370 indexes[ kountNonz ] = kountCon ;
3371 kountNonz++;
3372
3373 starts[ startsIdx++ ] = kountNonz;
3374
3375
3376 kountCon++;
3377
3378
3379 }
3380
3381 }
3382
3383 std::cout << "NUMBER OF VARIABLES ADDED = " << kount << std::endl;
3384
3385 // now add the objective function
3387 SparseVector *objcoeff = new SparseVector( numVvar);
3388
3389
3390 for(i = 0; i < numVvar; i++){
3391
3392 objcoeff->indexes[ i] = i;
3393 objcoeff->values[ i] = 1.0;
3394
3395 }
3396
3397
3398
3399
3400
3401 m_osinstanceSeparation->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
3402 //now for the nonzeros
3403 //add the linear constraints coefficients
3405 values, 0, kountNonz - 1, indexes, 0, kountNonz - 1, starts, 0, startsIdx);
3406
3407
3408
3409 //std::cout << m_osinstanceSeparation->printModel( ) << std::endl;
3410 //
3411 delete objcoeff;
3412
3413
3414 // now create the Clp model
3415
3416
3417 //below is temporty see if we can setup as a Clp network problem
3418 CoinPackedMatrix* matrix;
3419 bool columnMajor = true;
3420 double maxGap = 0;
3421 matrix = new CoinPackedMatrix(
3422 columnMajor, //Column or Row Major
3427 columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->indexes : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->indexes, //Pointer to start of minor dimension indexes -- change to allow for row storage
3429 0, 0, maxGap );
3430
3431 ClpNetworkMatrix network( *matrix);
3432
3433 m_separationClpModel = new ClpSimplex();
3434
3435 //if( m_osinstanceSeparation->getObjectiveMaxOrMins()[0] == "min") osiSolver->setObjSense(1.0);
3436 m_separationClpModel->setOptimizationDirection( 1);
3441 );
3442
3443 m_separationClpModel->factorization()->maximumPivots(200 + m_separationClpModel->numberRows() / 100);
3444
3445
3446 delete matrix;
3447
3448 }catch (const ErrorClass& eclass) {
3449
3450 throw ErrorClass(eclass.errormsg);
3451
3452 }
3453
3454 return NULL;
3455}//end getSeparationInstance
3456
3457
3458
3459int OSBearcatSolverXkij::getBranchingVar(const double* theta, const int numThetaVar ) {
3460
3461 int varIdx;
3462 varIdx = -1;
3463 int i;
3464 int j;
3465 int numVar = m_numNodes*m_numNodes - m_numHubs ;
3466
3467 double from1Distance;
3468 double from0Distance;
3469 double fraction;
3470 double minFraction;
3471
3472 double *xvalues;
3473
3474
3475 xvalues = new double[ numVar];
3476 for(i = 0; i < numVar; i++){
3477 xvalues[ i] = 0;
3478 }
3479
3480 try{
3481 if(numThetaVar != m_numThetaVar) throw ErrorClass("inconsistent number of variables in getBranchingVar");
3482 //loop over the fractional thetas
3483 for(i = 0; i < m_numThetaVar; i++){
3484
3485 if( ( theta[ i ] > m_osDecompParam.zeroTol ) && ( theta[ i ] < 1 - m_osDecompParam.zeroTol ) ){
3486
3487 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
3488
3489 xvalues[ m_thetaIndex[ j] ] += theta[ i ] ;
3490
3491 }
3492
3493 }
3494
3495
3496 }
3497
3498 //let's branch on a variable in and out of hub first
3499 minFraction = 1.0;
3500 //ideally we find minFraction very close to .5
3501
3502 for(i = 0; i < m_numHubs; i++){
3503
3504 for( j = 0; j < i; j++){
3505
3506 //j < i so the index is i*(m_numNodes - 1) + j
3507 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
3508 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
3509 fraction = std::max(from1Distance, from0Distance);
3510 //try to find fractional variable that is the closest to .5
3511 if(fraction < minFraction){
3512
3513 minFraction = fraction;
3514 varIdx = i*(m_numNodes - 1) + j;
3515 }
3516
3517 }
3518
3519 for(j = i + 1; j < m_numNodes; j++){
3520
3521 //j < i so the index is i*(m_numNodes - 1) + j - 1
3522 //j < i so the index is i*(m_numNodes - 1) + j
3523 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
3524 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
3525 fraction = std::max(from1Distance, from0Distance);
3526 //try to find fractional variable that is the closest to .5
3527 if(fraction < minFraction) {
3528
3529 minFraction = fraction;
3530 varIdx = i*(m_numNodes - 1) + j - 1;
3531 }
3532
3533
3534 }
3535
3536 }
3537
3538 //if we have a candidate among arcs in/out of hubs, take it
3539
3540 if(minFraction > 1 - m_osDecompParam.zeroTol){
3541
3542 for(i = m_numHubs; i < m_numNodes; i++){
3543
3544 for( j = 0; j < i; j++){
3545
3546 //j < i so the index is i*(m_numNodes - 1) + j
3547 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
3548 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
3549 fraction = std::max(from1Distance, from0Distance);
3550 //try to find fractional variable that is the closest to .5
3551 if(fraction < minFraction) {
3552
3553 minFraction = fraction;
3554 varIdx = i*(m_numNodes - 1) + j ;
3555 }
3556
3557 }
3558
3559 for(j = i + 1; j < m_numNodes; j++){
3560
3561 //j < i so the index is i*(m_numNodes - 1) + j - 1
3562 //j < i so the index is i*(m_numNodes - 1) + j
3563 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
3564 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
3565 fraction = std::max(from1Distance, from0Distance);
3566 //try to find fractional variable that is the closest to .5
3567 if(fraction < minFraction) {
3568
3569 minFraction = fraction;
3570 varIdx = i*(m_numNodes - 1) + j - 1;
3571 }
3572
3573 }
3574
3575 }
3576
3577 }//end of if on minFraction
3578 std::cout << " HERE IS GAIL 1" << std::endl;
3579
3580 //zero out the scatter array
3581
3582 delete[] xvalues;
3583 std::cout << " HERE IS GAIL 2" << std::endl;
3584 xvalues = NULL;
3585
3586 return varIdx;
3587
3588 }catch (const ErrorClass& eclass) {
3589
3590 delete[] xvalues;
3591 xvalues = NULL;
3592
3593 throw ErrorClass(eclass.errormsg);
3594
3595 }
3596
3597
3598}//end getBranchingVar Dense
3599
3600
3601
3602int OSBearcatSolverXkij::getBranchingVar(const int* thetaIdx, const double* theta,
3603 const int numThetaVar) {
3604
3605 int varIdx;
3606 varIdx = -1;
3607 int i;
3608 int j;
3609 int numVar = m_numNodes*m_numNodes - m_numHubs ;
3610 double from1Distance;
3611 double from0Distance;
3612 double fraction;
3613 double minFraction;
3614
3615 double *xvalues;
3616
3617
3618 xvalues = new double[ numVar];
3619 for(i = 0; i < numVar; i++){
3620 xvalues[ i] = 0;
3621 }
3622
3623 try{
3624 //loop over the fractional thetas
3625 //working with a sparse matrix
3626 for(i = 0; i < numThetaVar; i++){
3627
3628 if( ( theta[ i ] > m_osDecompParam.zeroTol ) && ( theta[ i ] < 1 - m_osDecompParam.zeroTol ) ){
3629
3630 for(j = m_thetaPnt[ thetaIdx[ i] ]; j < m_thetaPnt[ thetaIdx[ i] + 1] ; j++){
3631
3632 xvalues[ m_thetaIndex[ j] ] += theta[ i ] ;
3633
3634 }
3635
3636 }
3637
3638
3639 }
3640
3641
3642 //let's branch on a variable in and out of hub first
3643 minFraction = 1.0;
3644 //ideally we find minFraction very close to .5
3645
3646 for(i = 0; i < m_numHubs; i++){
3647
3648 for( j = 0; j < i; j++){
3649
3650 //j < i so the index is i*(m_numNodes - 1) + j
3651 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
3652 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
3653 fraction = std::max(from1Distance, from0Distance);
3654 //try to find fractional variable that is the closest to .5
3655 if(fraction < minFraction){
3656
3657 minFraction = fraction;
3658 varIdx = i*(m_numNodes - 1) + j;
3659 }
3660
3661 }
3662
3663 for(j = i + 1; j < m_numNodes; j++){
3664
3665 //j < i so the index is i*(m_numNodes - 1) + j - 1
3666 //j < i so the index is i*(m_numNodes - 1) + j
3667 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
3668 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
3669 fraction = std::max(from1Distance, from0Distance);
3670 //try to find fractional variable that is the closest to .5
3671 if(fraction < minFraction) {
3672
3673 minFraction = fraction;
3674 varIdx = i*(m_numNodes - 1) + j - 1;
3675 }
3676
3677
3678 }
3679
3680 }
3681
3682 //if we have a candidate among arcs in/out of hubs, take it
3683
3684 std::cout << "MIN FRACTION = " << minFraction << std::endl;
3685
3686 if(minFraction > 1 - m_osDecompParam.zeroTol){
3687
3688 for(i = m_numHubs; i < m_numNodes; i++){
3689
3690
3691
3692 for( j = 0; j < i; j++){
3693
3694 //j < i so the index is i*(m_numNodes - 1) + j
3695 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
3696 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
3697 fraction = std::max(from1Distance, from0Distance);
3698 //try to find fractional variable that is the closest to .5
3699 if(fraction < minFraction) {
3700
3701 minFraction = fraction;
3702 varIdx = i*(m_numNodes - 1) + j ;
3703 }
3704
3705 }
3706
3707 for(j = i + 1; j < m_numNodes; j++){
3708
3709 //j < i so the index is i*(m_numNodes - 1) + j - 1
3710 //j < i so the index is i*(m_numNodes - 1) + j
3711 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
3712 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
3713 fraction = std::max(from1Distance, from0Distance);
3714 //try to find fractional variable that is the closest to .5
3715 if(fraction < minFraction) {
3716
3717 minFraction = fraction;
3718 varIdx = i*(m_numNodes - 1) + j - 1;
3719 }
3720
3721 }
3722
3723 }
3724
3725 }//end of if on minFraction
3726
3727 //zero out the scatter array
3728
3729 delete[] xvalues;
3730 xvalues = NULL;
3731
3732 return varIdx;
3733
3734 }catch (const ErrorClass& eclass) {
3735
3736 delete[] xvalues;
3737 xvalues = NULL;
3738
3739 throw ErrorClass(eclass.errormsg);
3740
3741 }
3742
3743
3744}//end getBranchingVar Sparse
3745
3746
3747void OSBearcatSolverXkij::getBranchingCut(const double* thetaVar, const int numThetaVar,
3748 const std::map<int, int> &varConMap, int &varIdx, int &numNonz,
3749 int* &indexes, double* &values) {
3750
3751 //get a branching variable
3752 int i;
3753 int j;
3754 int kount;
3755 numNonz = 0;
3756 //keep numNonz at zero if there is no cut
3757 //there will be no cut if the xij is in conVarMap
3758
3759 try{
3760
3761 if(numThetaVar != m_numThetaVar) throw ErrorClass("inconsistent number of variables in getBranchingCut");
3762
3763
3764 varIdx = getBranchingVar(thetaVar, numThetaVar );
3765
3766 std::cout << "Branching on Variable: " << m_variableNames[ varIdx] << std::endl;
3767
3768 //if this variable is in the map, then we just return with the index,
3769 //if this variable is NOT in the map then we add a cut
3770
3771 if( varConMap.find( varIdx) == varConMap.end() ){
3772
3773 for(i = 0; i < m_numThetaVar; i++){
3774
3775 kount = 0;
3776
3777 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
3778
3779 if ( m_thetaIndex[ j] == varIdx) kount++ ;
3780
3781 }
3782
3783 //count is the number times variable i appears in the constraint
3784
3785 if(kount > 0){
3786
3787 branchCutIndexes[ numNonz] = i;
3788 branchCutValues[ numNonz++] = kount ;
3789
3790 }
3791
3792 }
3793
3794
3795 //add varIdx cut to B matrix
3796 m_Bmatrix[ m_numBmatrixNonz++] = varIdx;
3799
3800 //make sure to add artificial variables
3801 //of course they have no nonzero elements in
3802 //the transformation matrix
3803 m_numThetaVar++;
3806 //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz; //first artificial variable
3807 //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz; //second artificial variable
3808
3809
3810 }//end of if on checking for map membership
3811
3812 //set return arguments
3813
3814 indexes = branchCutIndexes;
3815 values = branchCutValues;
3816
3817 return;
3818
3819 }catch (const ErrorClass& eclass) {
3820
3821 throw ErrorClass(eclass.errormsg);
3822
3823 }
3824
3825}//end getBranchingCut dense
3826
3827
3828void OSBearcatSolverXkij::getBranchingCut(const int* thetaIdx, const double* thetaVar,
3829 const int numThetaVar, const std::map<int, int> &varConMap,
3830 int &varIdx, int &numNonz, int* &indexes, double* &values) {
3831
3832 //get a branching variable
3833 int i;
3834 int j;
3835 int kount;
3836 numNonz = 0;
3837 //keep numNonz at zero if there is no cut
3838 //there will be no cut if the xij is in conVarMap
3839
3840 try{
3841
3842
3843
3844 varIdx = getBranchingVar(thetaIdx, thetaVar, numThetaVar );
3845
3846 std::cout << "Branching on Variable: " << m_variableNames[ varIdx] << std::endl;
3847
3848 //if this variable is in the map, then we just return with the index,
3849 //if this variable is NOT in the map then we add a cut
3850
3851 if( varConMap.find( varIdx) == varConMap.end() ){
3852
3853
3854
3855
3856
3857
3858 for(i = 0; i < m_numThetaVar; i++){
3859
3860 kount = 0;
3861
3862 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
3863
3864 if ( m_thetaIndex[ j] == varIdx) kount++ ;
3865
3866 }
3867
3868 //count is the number times variable i appears in the constraint
3869
3870 if(kount > 0){
3871
3872 branchCutIndexes[ numNonz] = i;
3873 branchCutValues[ numNonz++] = kount ;
3874
3875 }
3876
3877 }
3878
3879
3880 //add varIdx cut to B matrix
3881 m_Bmatrix[ m_numBmatrixNonz++] = varIdx;
3884
3885 //make sure to add artificial variables
3886 //of course they have no nonzero elements in
3887 //the transformation matrix
3888 m_numThetaVar++;
3891 //m_thetaPnt[ m_numThetaVar++] = m_numThetaNonz; //first artificial variable
3892 //m_thetaPnt[ m_numThetaVar] = m_numThetaNonz; // second artificial variable
3893
3894
3895 }//end of if on checking for map membership
3896
3897 //set return arguments
3898
3899 indexes = branchCutIndexes;
3900 values = branchCutValues;
3901
3902 return;
3903
3904 }catch (const ErrorClass& eclass) {
3905
3906 throw ErrorClass(eclass.errormsg);
3907
3908 }
3909
3910}//end getBranchingCut sparse
3911
3912
3914
3915 try{
3916 //kipp -- stil not done we depend on SKs solution
3917 //let's get the initial assignment of nodes to hubs
3918 //this is in m_initSolMap which is calculated when we
3919 //call getOptions( OSOption *osoption);
3920
3921 if(m_initSolMap.size() == 0) getOptions( m_osoption);
3922
3923 //get cost vector
3924
3925 //get demand vector
3926
3927 //now improve on m_initSolMap
3928
3929 }catch (const ErrorClass& eclass) {
3930
3931 throw ErrorClass(eclass.errormsg);
3932
3933 }
3934
3935
3936}//end getInitialSolution
3937
3938
3939void OSBearcatSolverXkij::resetMaster( std::map<int, int> &inVars, OsiSolverInterface *si){
3940
3941 int i;
3942 int j;
3943
3944 int kount;
3945 int numNonz;
3946
3947 std::map<int, int>::iterator mit;
3948 //temporarily create memory for the columns we keep
3949 int numVars = inVars.size();
3950 int numVarArt;
3951 //there 2*m_numNodes in the A matrix
3952 //there are m_numBmatrixCon B matrix constraints
3953 //numVarArt = 2*m_numNodes + 2*m_numBmatrixCon;
3954 numVarArt = m_numNodes + m_numBmatrixCon;
3955
3956 //arrays for the new osinstance
3957 std::vector<double> valuesVec;
3958 double *values = NULL;
3959
3960 std::vector<int> indexesVec;
3961 int *indexes = NULL ;
3962
3963 int *starts = new int[ numVars + 1 + numVarArt];
3964
3965 int startsIdx;
3966
3967
3968
3969 //temporay holders
3970 int* thetaPntTmp;
3971 int* thetaIndexTmp;
3972 int* tmpConvexity = new int[ m_numThetaVar];
3973
3974 //get the number of nonzeros that we need
3975 numNonz = 0;
3976
3977 for(mit = inVars.begin(); mit != inVars.end(); mit++){
3978
3979 numNonz += m_thetaPnt[mit->first + 1 ] - m_thetaPnt[ mit->first ];
3980 }
3981
3982 //allocate the memory
3983 thetaPntTmp = new int[ numVars + 1];
3984 thetaIndexTmp = new int[ numNonz];
3985
3986
3987 //error check
3988 for(mit = inVars.begin(); mit != inVars.end(); mit++){
3989
3990
3991 //std::cout << "VARIABLE INDEX = " << mit->first << " OBJ COEF = " << si->getObjCoefficients()[ mit->first ] << std::endl;
3992 if( convexityRowIndex[ mit->first] == -1) throw ErrorClass( "we have an artificial variable in reset master");
3993
3994
3995 }
3996
3997 //fill in the temporary arrays
3998 kount = 0;
3999 numNonz = 0;
4000 thetaPntTmp[ kount] = 0;
4001
4002 for(mit = inVars.begin(); mit != inVars.end(); mit++){
4003
4004 //std::cout << "mit->first " << mit->first << " mit->second " << mit->second << std::endl;
4005
4006 kount++;
4007
4008 for(i = m_thetaPnt[ mit->first ]; i < m_thetaPnt[mit->first + 1 ]; i++){
4009
4010 thetaIndexTmp[ numNonz++] = m_thetaIndex[ i];
4011
4012 //std::cout << "Column = " << mit->first << " Variable " << m_variableNames[ m_thetaIndex[ i] ] << std::endl;
4013
4014 }
4015
4016 thetaPntTmp[ kount] = numNonz;
4017
4018 //std::cout << "kount = " << kount << " thetaPntTmp[ kount] = " << thetaPntTmp[ kount] << std::endl;
4019 //readjust numbering to take into account artificial variables
4020 //mit->second += numVarArt;
4021 //kipp -- double check calculation below
4022 inVars[ mit->first] = numVarArt + kount - 1 ;
4023
4024 }
4025
4026 //std::cout << "kount = " << kount << std::endl;
4027 //std::cout << "numVars = " << numVars << std::endl;
4028
4029
4030
4031 //copy old values of convexityRowIndex
4032 for(i = 0; i < m_numThetaVar; i++) tmpConvexity[ i] = convexityRowIndex[ i];
4033
4034 //reset the theta pointers
4035 //first the artificial variables
4036 m_numThetaVar = 0;
4037 m_numThetaNonz = 0;
4038 for(i = 0; i < numVarArt; i++){
4039
4041 m_thetaPnt[ m_numThetaVar++] = 0;
4042
4043
4044 }
4045 //now fill in the other pointers from the temp arrarys
4046 //std::cout << "Number of artificial variables = " << numVarArt << std::endl;
4047 intVarSet.clear();
4048 for(mit = inVars.begin(); mit != inVars.end(); mit++){
4049
4050
4051 intVarSet.insert ( std::pair<int,double>(mit->second, 1.0) );
4052
4053 //std::cout << " m_numThetaVar = " << m_numThetaVar << " m_numThetaNonz = " << m_numThetaNonz << std::endl;
4054 //std::cout << "Variable number " << mit->first << " OBJ coefficient = " << si->getObjCoefficients()[ mit->first] << std::endl;
4055
4056 convexityRowIndex[ m_numThetaVar] = tmpConvexity[ mit->first];
4057
4059
4060 for(j = thetaPntTmp[ mit->second - numVarArt]; j < thetaPntTmp[ mit->second - numVarArt + 1 ]; j++){
4061
4062
4063 m_thetaIndex[ m_numThetaNonz ] = thetaIndexTmp[ m_numThetaNonz] ;
4065
4066 }
4067
4068 }
4069
4071 //std::cout << " number of art vars = " << numVarArt << std::endl;
4072 //std::cout << " m_numThetaVar = " << m_numThetaVar << std::endl;
4073 //std::cout << " m_numThetaNonz = " << m_numThetaNonz << std::endl;
4074 //done with the transformation matrix
4075
4076
4077
4078 //
4079 //write old master --- just for testing
4080 si->writeLp( "gailTest" );
4081
4082 //now create the formulation
4083
4084 //first get each column of the new master
4085 //first take care of the artificial variables
4086 numNonz = 0;
4087 startsIdx = 0;
4088 starts[ startsIdx++] = numNonz;
4089
4090 for(i = 0; i < numVarArt; i++){
4091 numNonz++;
4092 starts[ startsIdx++] = numNonz;
4093 valuesVec.push_back( 1.0);
4094 indexesVec.push_back( i);
4095
4096 }
4097
4098
4099 int rowCount;
4100
4101 int numAmatrixRows;
4102 numAmatrixRows = m_numNodes - m_numHubs;
4103
4104 for(mit = inVars.begin(); mit != inVars.end(); mit++){
4105
4106 //std::cout << "CONVEXITY ROW = " << convexityRowIndex[ mit->second] << std::endl;
4107 valuesVec.push_back( 1.0);
4108 indexesVec.push_back( numAmatrixRows + convexityRowIndex[ mit->second] );
4109 //increment numNonz by 1 for the convexity row
4110 numNonz++;
4111
4112 for(j = m_thetaPnt[ mit->second ]; j < m_thetaPnt[ mit->second + 1 ]; j++){
4113
4115
4116 //std::cout << "Column = " << mit->second << " Variable " << m_variableNames[ m_thetaIndex[ j] ] << std::endl;
4117
4118 }
4119
4120
4121
4122 //multiply the sparse array by each A matrix constraint
4123 for(i = 0; i < numAmatrixRows; i++){
4124
4125 rowCount = 0;
4126
4127 for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
4128
4129 //m_Amatrix[ j] is a variable index -- this logic works
4130 //since the Amatrix coefficient is 1 -- we don't need a value
4131 //it indexes variable that points into the node
4132 rowCount += m_tmpScatterArray[ m_Amatrix[ j] ];
4133
4134
4135 }
4136
4137 if(rowCount > 0){
4138
4139 numNonz++;
4140
4141 //std::cout << "Column = " << mit->second << " Nonzero in A marix row " << i << " Value = " << rowCount << std::endl;
4142 valuesVec.push_back( rowCount);
4143 indexesVec.push_back( i);
4144
4145
4146 }
4147
4148
4149 }//end loop on coupling constraints
4150
4151
4152 //multiply the sparse array by each B matrix constraint
4153 for(i = 0; i < m_numBmatrixCon; i++){
4154
4155 rowCount = 0;
4156
4157 for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
4158
4159 //m_Amatrix[ j] is a variable index -- this logic works
4160 //since the Amatrix coefficient is 1 -- we don't need a value
4161 //it indexes variable that points into the node
4162 rowCount += m_tmpScatterArray[ m_Bmatrix[ j] ];
4163
4164
4165 }
4166
4167 if(rowCount > 0){
4168 numNonz++;
4169
4170 //std::cout << "Column = " << mit->first << " Nonzero in B matrix row " << i + m_numNodes<< " Value = " << rowCount << std::endl;
4171
4172 valuesVec.push_back( rowCount);
4173 indexesVec.push_back( i + m_numNodes);
4174 }
4175
4176
4177 }//end loop on B matrix constraints
4178
4179
4180 //zero out the scatter array
4181
4182 for(j = m_thetaPnt[ mit->second ]; j < m_thetaPnt[ mit->second + 1 ]; j++){
4183
4185
4186 }
4187
4188 starts[ startsIdx++] = numNonz;
4189
4190 }
4191
4192
4193 //for(i = 0; i < startsIdx; i++) std::cout << "starts[ i] = " << starts[ i] << std::endl;
4194 values = new double[ numNonz];
4195 indexes = new int[ numNonz];
4196
4197 if(numNonz != valuesVec.size() ) throw ErrorClass("dimension problem in reset");
4198 if(numNonz != indexesVec.size() ) throw ErrorClass("dimension problem in reset");
4199
4200 for(i = 0; i < numNonz; i++){
4201
4202 values[ i] = valuesVec[i];
4203 indexes[ i] = indexesVec[i];
4204
4205 }
4206
4207
4208
4209 //construct the new master
4210 //create an OSInstance from the tmp arrays
4211 // delete the old m_osinstanceMaster
4212
4213 delete m_osinstanceMaster;
4214 m_osinstanceMaster = NULL;
4215
4216 //start building the restricted master here
4218 m_osinstanceMaster->setInstanceDescription("The Restricted Master");
4219
4220 // first the variables
4221 // we have numVarArt] artificial variables
4222 m_osinstanceMaster->setVariableNumber( numVars + numVarArt );
4223
4224 // now add the objective function
4225 //m_osinstanceMaster->setObjectiveNumber( 1);
4226 //add m_numNodes artificial variables
4227 SparseVector *objcoeff = new SparseVector( numVars + numVarArt);
4228
4229
4230
4231 // now the constraints
4233
4234
4235 //add the artificial variables first
4236
4237 int varNumber;
4238 varNumber = 0;
4239
4240
4241 //define the artificial variables
4242 for(i = 0; i < numVarArt; i++){
4243
4244 objcoeff->indexes[ varNumber ] = varNumber ;
4245
4246 objcoeff->values[ varNumber ] = m_osDecompParam.artVarCoeff;
4247
4249
4250 m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt2("x", i ) ,
4251 0, 1.0, 'C');
4252
4253
4254 }
4255
4256 // now the theta variables
4257 kount = 0;
4258 for(mit = inVars.begin(); mit != inVars.end(); mit++){
4259
4260 objcoeff->indexes[ varNumber ] = varNumber ;
4261
4262 objcoeff->values[ varNumber ] = si->getObjCoefficients()[ mit->first] ;
4263
4264 m_thetaCost[ varNumber] = si->getObjCoefficients()[ mit->first];
4265
4266 m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt2("x", kount + numVarArt ) ,
4267 0, 1.0, 'C');
4268
4269 kount++;
4270
4271
4272
4273 }
4274
4275
4276
4277 for(i = 0; i < m_numNodes; i++){
4278
4280 1.0, 1.0, 0);
4281
4282 }
4283
4284
4285 for(i = m_numNodes; i < m_numBmatrixCon + m_numNodes; i++){
4286
4288 si->getRowLower()[ i], si->getRowUpper()[ i], 0);
4289
4290
4291 }
4292
4293
4294 // now add the objective function
4296 m_osinstanceMaster->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
4297
4298 //add the linear constraints coefficients
4300 values, 0, numNonz - 1, indexes, 0, numNonz - 1, starts, 0, startsIdx);
4301
4302
4303 //std::cout << m_osinstanceMaster->printModel( ) << std::endl;
4304
4305 //garbage collection
4306 delete[] tmpConvexity;
4307 tmpConvexity = NULL;
4308 delete[] thetaPntTmp;
4309 thetaPntTmp = NULL;
4310 delete[] thetaIndexTmp;
4311 thetaIndexTmp = NULL;
4312 delete objcoeff;
4313 objcoeff = NULL;
4314}//end resetMaster
4315
4316
4317
4318std::string makeStringFromInt2(std::string theString, int theInt){
4319 ostringstream outStr;
4320 outStr << theString;
4321 outStr << theInt;
4322 return outStr.str();
4323}//end makeStringFromInt
4324
4325
4326
4327
4328
std::string makeStringFromInt2(std::string theString, int theInt)
OSOption * osoption
#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
class used to make it easy to read and write files.
Definition: OSFileUtil.h:38
std::string getFileAsString(const char *fname)
read a file and return contents as a string.
Definition: OSFileUtil.cpp:35
Variables * variables
variables is a pointer to a Variables object
Definition: OSInstance.h:2185
Objectives * objectives
objectives is a pointer to a Objectives object
Definition: OSInstance.h:2188
OSInstance * getSeparationInstance()
virtual OSInstance * getInitialRestrictedMaster()
OSInstance* OSBearcatSolverXkij::getInitialRestrictedMaster( ){.
int * m_routeCapacity
the route capacity – bus seating limit this can vary with the route/hub
virtual void resetMaster(std::map< int, int > &inVars, OsiSolverInterface *si)
INPUT:
double qrouteCost(const int &k, const int &l, const double *c, int *kountVar)
kipp – document
OSInstance * m_osinstanceSeparation
int * m_upperBoundL
largest possible L-value on a route – this will be the minimum of m_routeCapacity and total demand
double ** m_rc
the reduced cost vector for each k, we asssume order is (l, i, j)
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...
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...
int m_minDemand
m_minDemand is the value of the minimum demand node – it is not the minimum demand that must be carri...
int * convexityRowIndex
conconvexityRowIndex holds the index of the convexity row that the theta columns are in.
void getOptions(OSOption *osoption)
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...
~OSBearcatSolverXkij()
Default Destructor.
int * m_routeMinPickup
the minimum number of students that we pickup on a route this can vary with the route/hub
OSBearcatSolverXkij()
Default Constructor.
void getInitialSolution()
generate an intitial feasible solution in theta space for the initial master
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...
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...
int * m_separationIndexMap
m_separationIndexMap maps the variable index into the appropriate row in the separation problem for t...
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...
virtual void pauHana(std::vector< int > &m_zOptIndexes, int numNodes, int numColsGen)
int * m_lowerBoundL
smallest possible L-value on a route for now this will equal
ClpSimplex * m_separationClpModel
int * m_demand
m_demand is the vector of node demands
virtual void initializeDataStructures()
allocate memory and initialize arrays
int m_maxThetaNonz
m_maxMasterNonz is the maximumn number of nonzero elements we allow in the transformation matrix betw...
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 ...
int getBranchingVar(const double *theta, const int numThetaVar)
RETURN VALUES: return the integer index of a fractional x_{ij} variable.
int m_upperBoundLMax
largest possible L-value over all routes
double ** m_cost
the distance/cost vectors
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_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_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 bVariablesModified
bVariablesModified is true if the variables data has been modified.
Definition: OSInstance.h:2288
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.
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
InstanceData * instanceData
A pointer to an InstanceData object.
Definition: OSInstance.h:2278
double * getVariableLowerBounds()
Get variable lower bounds.
int getVariableNumber()
Get number of variables.
bool setInstanceDescription(std::string description)
set the instance description.
std::string * getVariableNames()
Get variable names.
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
std::string getSolutionStatusType(int solIdx)
Get the [i]th optimization solution status type, where i equals the given solution index.
Definition: OSResult.cpp:2051
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
Used to read an OSiL string.
Definition: OSiLReader.h:38
OSInstance * readOSiL(const std::string &osil)
parse the OSiL model instance.
Definition: OSiLReader.cpp:53
double value
value is the value of the objective function coefficient corresponding to the variable with index idx
Definition: OSInstance.h:128
ObjCoef ** coef
coef is pointer to an array of ObjCoef object pointers
Definition: OSInstance.h:176
Objective ** obj
coef is pointer to an array of ObjCoef object pointers
Definition: OSInstance.h:205
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
double lb
lb corresponds to the optional attribute that holds the variable lower bound.
Definition: OSInstance.h:56
Variable ** var
Here we define a pointer to an array of var pointers.
Definition: OSInstance.h:97
const int OSINT_MAX
Definition: OSParameters.h:94
const double OSDBL_MAX
Definition: OSParameters.h:93