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.
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
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.
std::string errormsg
errormsg is the error that is causing the exception to be thrown
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.
Variables * variables
variables is a pointer to a Variables object
Objectives * objectives
objectives is a pointer to a Objectives object
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
double zeroTol
we terminate column generation when the reduced costs are not smaller than zeroTol
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...
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
The in-memory representation of an OSiL instance..
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.
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.
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.
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.
std::string getSolutionStatusType(int solIdx)
Get the [i]th optimization solution status type, where i equals the given solution index.
std::vector< IndexValuePair * > getOptimalPrimalVariableValues(int solIdx)
Get one solution of optimal primal variable values.
double getObjValue(int solIdx, int objIdx)
Used to read an OSiL string.
Definition OSiLReader.h:38
OSInstance * readOSiL(const std::string &osil)
parse the OSiL model instance.
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
#define OSINT_MAX
#define OSDBL_MAX