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.
Definition: OSCoinSolver.h:38
virtual void buildSolverInstance()
The implementation of the corresponding virtual function.
OsiSolverInterface * osiSolver
osiSolver is the osi solver object – in this case clp, glpk, cbc, cplex, symphony or dylp
Definition: OSCoinSolver.h:93
virtual void solve()
The implementation of the corresponding virtual function.
std::string sSolverName
sSolverName is the name of the Coin solver used, e.g.
OSInstance * osinstance
osinstance holds the problem instance in-memory as an OSInstance object
OSOption * osoption
osoption holds the solver options in-memory as an OSOption object
OSResult * osresult
osresult holds the solution or result of the model in-memory as an OSResult object
used for throwing exceptions.
Definition: OSErrorClass.h:32
std::string errormsg
errormsg is the error that is causing the exception to be thrown
Definition: OSErrorClass.h:42
int * m_demand
m_demand is the vector of node demands
std::map< int, std::string > m_tmpVarMap
double getRouteDistance(int numNodes, int hubIndex, CoinSolver *solver, std::vector< int > zk, double *xVar)
call this method to get a minimum TSP tour for a given assignment of nodes to routes INPUT: int numNo...
double * m_cost
the distance/cost vectors
double calcNonlinearRelax(std::vector< double > &m_zRootLPx_vals)
calculate the nonlinear relaxation value for an LP solution
std::string * m_nodeName
m_nodeName is the vector of node names
int * m_upperBoundL
largest possible L-value on a route – this will be the minimum of m_routeCapacity and total demand
void getCutsMultiCommod(const double *thetaVar, const int numThetaVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)
This is the routine that generates the multi-item cuts.
std::string m_initOSiLFile
int * m_BmatrixRowIndex
m_BmatrixRowIndex holds the index of the convexity row that the constraint corresponds to,...
CoinSolver * getTSP(int numNodes, double *cost)
call this method to get a TSP instance
int * m_lowerBoundL
smallest possible L-value on a route for now this will equal
int * m_hubPoint
m_hubPoint[ k] points to the the k'th hub that we use in getOptL
ClpSimplex * m_separationClpModel
int * m_separationIndexMap
m_separationIndexMap maps the variable index into the appropriate row in the separation problem for t...
double qrouteCost(const int &k, const int &l, const double *c, int *kountVar)
kipp – document
void getVariableIndexMap()
this method will populate: std::map<std::pair<int, int>,int>m_xVarIndexMap this gives us
int * m_routeMinPickup
the minimum number of students that we pickup on a route this can vary with the route/hub
double ** m_rc
the reduced cost vector for each k, we asssume order is (l, i, j)
int m_maxThetaNonz
m_maxMasterNonz is the maximumn number of nonzero elements we allow in the transformation matrix betw...
virtual OSInstance * getInitialRestrictedMaster()
OSInstance* OSBearcatSolverXij::getInitialRestrictedMaster( ){.
~OSBearcatSolverXij()
Default Destructor.
void getInitialSolution()
generate an intitial feasible solution in theta space for the initial master
std::map< std::pair< int, int >, int > m_xVarIndexMap
m_xVarIndexMap takes a variable indexed by (i,j) and returns the index of the variable in one dimensi...
void getOptions(OSOption *osoption)
virtual void resetMaster(std::map< int, int > &inVars, OsiSolverInterface *si)
INPUT:
int m_upperBoundLMax
largest possible L-value over all routes
virtual void initializeDataStructures()
allocate memory and initialize arrays
virtual void pauHana(std::vector< int > &m_zOptIndexes, std::vector< double > &m_zRootLPx_vals, int numNodes, int numColsGen, std::string message)
OSInstance * m_osinstanceSeparation
std::map< int, std::map< int, std::vector< int > > > m_initSolMap
the index on the outer map is on the solution number, the index on the inner map indexes the route nu...
void permuteHubs()
this method will calculate a permuation of the hubs so that they are in ascending order,...
virtual void getBranchingCut(const double *thetaVar, const int numThetaVar, const std::map< int, int > &varConMap, int &varIdx, int &numNonz, int *&indexes, double *&values)
RETURN VALUES: varIdx – the variable number x_{ij} for branching numNonz – number of theta indexes in...
int getBranchingVar(const double *theta, const int numThetaVar)
RETURN VALUES: return the integer index of a fractional x_{ij} variable.
int m_minDemand
m_minDemand is the value of the minimum demand node – it is not the minimum demand that must be carri...
std::vector< CoinSolver * > m_multCommodCutSolvers
m_multCommodCutSolvers is a vector of solvers, one solver for each hub, used to find multicommodity f...
OSInstance * getSeparationInstance()
void getOptL(double **c)
bool m_costSetInOption
m_costSetInOption is true if the costs are set using the OSOption file
void getFeasibleSolution()
call this method to get generate an instance of the generalized assignment problem and find a feasibl...
int * m_routeCapacity
the route capacity – bus seating limit this can vary with the route/hub
OSBearcatSolverXij()
Default Constructor.
bool OneOPT()
try and find a feasible solution, return false if solution not feasible
bool m_use1OPTstart
if m_use1OPTstart is true we use the option file to fix the nodes to hubs found by SK's 1OPT heuristi...
int * m_convexityRowIndex
m_convexityRowIndex holds the index of the convexity row that the theta columns are in.
virtual void getColumns(const double *yA, const int numARows, const double *yB, const int numBRows, int &numNewColumns, int *&numNonz, double *&cost, int **&rowIdx, double **&values, double &lowerBound)
RETURN VALUES: int numNewColumns – number of new columns generated int* numNonz – number of nonzeros ...
void getCutsX(const double *xVar, const int numXVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)
RETURN VALUES: int numNewRows – number of new rows generated int* numNonz – number of nonzeros in eac...
void getCutsTheta(const double *thetaVar, const int numThetaVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)
RETURN VALUES: int numNewRows – number of new rows generated int* numNonz – number of nonzeros in eac...
CoinSolver * getMultiCommodInstance(int hubIndex)
call this method to get an instance that is used to generate a multicommodity cut
void calcReducedCost(const double *yA, const double *yB)
calculate the reduced costs c – input of the objective function costs yA – dual values on node assign...
double artVarCoeff
artVarCoeff is the coefficient of the artificial variable in the objective function
Definition: OSDecompParam.h:55
double zeroTol
we terminate column generation when the reduced costs are not smaller than zeroTol
Definition: OSDecompParam.h:50
OSDecompParam m_osDecompParam
share the parameters with the decomposition solver
int m_maxBmatrixCon
m_maxBmatrixCon is the maximum number of B matrix constraints it is the number of tour breaking const...
double m_rootLPValue
double m_bestIPValue
int m_numNodes
m_numNodes is the number of nodes (both pickup and hub) in the model
int m_maxMasterRows
m_maxMasterColumns is the maximumn number of rows we allow in the master, in this application it is e...
int m_numBmatrixCon
m_numBmatrixCon is the number of constraints in B - 1, we have the -1 because: m_pntBmatrix[ k] point...
int m_numHubs
m_numHubs is the number of hubs/routes
std::string * m_variableNames
std::set< std::pair< int, double > > intVarSet
intVarSet holds and std::pair where the first element is the index of an integer variable and the sec...
OSInstance * m_osinstanceMaster
OSOption * m_osoption
int m_maxMasterColumns
m_maxMasterColumns is the maximumn number of columns we allow in the master
int m_maxBmatrixNonz
m_maxBmatrixNonz is the maximum number of nonzero elements in the B matrix constraints
double * m_BmatrixVal
double m_bestLPValue
The in-memory representation of an OSiL instance..
Definition: OSInstance.h:2263
double * getConstraintLowerBounds()
Get constraint lower bounds.
double * getVariableUpperBounds()
Get variable upper bounds.
bool setConstraintNumber(int number)
set the number of constraints.
bool addVariable(int index, std::string name, double lowerBound, double upperBound, char type)
add a variable.
std::string printModel()
Print the infix representation of the problem.
bool addConstraint(int index, std::string name, double lowerBound, double upperBound, double constant)
add a constraint.
int getConstraintNumber()
Get number of constraints.
bool bConstraintsModified
bConstraintsModified is true if the constraints data has been modified.
Definition: OSInstance.h:2298
int getLinearConstraintCoefficientNumber()
Get number of specified (usually nonzero) linear constraint coefficient values.
SparseMatrix * getLinearConstraintCoefficientsInRowMajor()
Get linear constraint coefficients in row major.
SparseMatrix * getLinearConstraintCoefficientsInColumnMajor()
Get linear constraint coefficients in column major.
bool setInstanceSource(std::string source)
set the instance source.
double ** getDenseObjectiveCoefficients()
getDenseObjectiveCoefficients.
bool setLinearConstraintCoefficients(int numberOfValues, bool isColumnMajor, double *values, int valuesBegin, int valuesEnd, int *indexes, int indexesBegin, int indexesEnd, int *starts, int startsBegin, int startsEnd)
set linear constraint coefficients
double * getVariableLowerBounds()
Get variable lower bounds.
int getVariableNumber()
Get number of variables.
bool setInstanceDescription(std::string description)
set the instance description.
bool addObjective(int index, std::string name, std::string maxOrMin, double constant, double weight, SparseVector *objectiveCoefficients)
add an objective.
bool setObjectiveNumber(int number)
set the number of objectives.
bool setVariableNumber(int number)
set the number of variables.
double * getConstraintUpperBounds()
Get constraint upper bounds.
The Option Class.
Definition: OSOption.h:3565
std::vector< SolverOption * > getSolverOptions(std::string solver_name)
Get the options associated with a given solver.
Definition: OSOption.cpp:4508
The Result Class.
Definition: OSResult.h:2549
std::string getSolutionStatusType(int solIdx)
Get the [i]th optimization solution status type, where i equals the given solution index.
Definition: OSResult.cpp:2051
double getOptimalObjValue(int objIdx, int solIdx)
Get one solution's optimal objective value.
Definition: OSResult.cpp:3065
std::vector< IndexValuePair * > getOptimalPrimalVariableValues(int solIdx)
Get one solution of optimal primal variable values.
Definition: OSResult.cpp:2215
double getObjValue(int solIdx, int objIdx)
Definition: OSResult.cpp:3050
int * indexes
indexes holds an integer array of rowIdx (or colIdx) elements in coefMatrix (AMatrix).
Definition: OSGeneral.h:258
int * starts
starts holds an integer array of start elements in coefMatrix (AMatrix), which points to the start of...
Definition: OSGeneral.h:252
double * values
values holds a double array of value elements in coefMatrix (AMatrix), which contains nonzero element...
Definition: OSGeneral.h:264
a sparse vector data structure
Definition: OSGeneral.h:123
double * values
values holds a double array of nonzero values.
Definition: OSGeneral.h:164
int * indexes
indexes holds an integer array of indexes whose corresponding values are nonzero.
Definition: OSGeneral.h:159
bool bDeleteArrays
bDeleteArrays is true if we delete the arrays in garbage collection set to true by default
Definition: OSGeneral.h:149
#define OSINT_MAX
Definition: OSParameters.h:122
#define OSDBL_MAX
Definition: OSParameters.h:115
OSResult * osresult