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