source: branches/split/Clp/test/OsiClpSolverInterfaceTest.cpp @ 1527

Last change on this file since 1527 was 1527, checked in by stefan, 10 years ago

add OsiClp? including OsiClp? unittest

File size: 50.1 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include "OsiConfig.h"
7
8#include <cassert>
9//#include <cstdlib>
10//#include <cstdio>
11//#include <iostream>
12
13#include "OsiClpSolverInterface.hpp"
14#include "OsiCuts.hpp"
15#include "OsiRowCut.hpp"
16#include "OsiColCut.hpp"
17#include "CoinMessage.hpp"
18#include "ClpMessage.hpp"
19#include "ClpFactorization.hpp"
20#include "CoinModel.hpp"
21#include "CoinIndexedVector.hpp"
22
23//#############################################################################
24
25class OsiClpMessageTest :
26   public CoinMessageHandler {
27
28public:
29  virtual int print() ;
30  OsiClpMessageTest();
31};
32
33OsiClpMessageTest::OsiClpMessageTest() : CoinMessageHandler()
34{
35}
36int
37OsiClpMessageTest::print()
38{
39  if (currentMessage().externalNumber()==0&&currentSource()=="Clp") 
40    std::cout<<"This is not actually an advertisement by Dash Associates - just my feeble attempt to test message handling and language - JJHF"<<std::endl;
41  else if (currentMessage().externalNumber()==5&&currentSource()=="Osi") 
42    std::cout<<"End of search trapped"<<std::endl;
43  return CoinMessageHandler::print();
44}
45
46//--------------------------------------------------------------------------
47// test EKKsolution methods.
48int
49OsiClpSolverInterfaceUnitTest(const std::string & mpsDir, const std::string & netlibDir)
50{
51 
52  // Test default constructor
53  {
54    OsiClpSolverInterface m;
55    assert( m.rowsense_==NULL );
56    assert( m.rhs_==NULL );
57    assert( m.rowrange_==NULL );
58    assert( m.matrixByRow_==NULL );
59    assert( m.ws_==NULL);
60    assert( m.itlimOrig_==9999999);
61    assert( m.lastAlgorithm_==0);
62    assert( m.integerInformation_==NULL);
63  }
64 
65 
66  {   
67    CoinRelFltEq eq;
68    OsiClpSolverInterface m;
69    std::string fn = mpsDir+"exmip1";
70    m.readMps(fn.c_str(),"mps");
71   
72    {
73      OsiClpSolverInterface im;   
74     
75      assert( im.getNumCols() == 0 ); 
76     
77      assert( im.getModelPtr()!=NULL );
78      // Test reset
79      im.reset();
80      assert( im.rowsense_==NULL );
81      assert( im.rhs_==NULL );
82      assert( im.rowrange_==NULL );
83      assert( im.matrixByRow_==NULL );
84      assert( im.ws_==NULL);
85      assert( im.itlimOrig_==9999999);
86      assert( im.lastAlgorithm_==0);
87      assert( im.integerInformation_==NULL);
88    }
89   
90    // Test copy constructor and assignment operator
91    {
92      OsiClpSolverInterface lhs;
93      {     
94        OsiClpSolverInterface im(m);       
95       
96        OsiClpSolverInterface imC1(im);
97        assert( imC1.getModelPtr()!=im.getModelPtr() );
98        assert( imC1.getNumCols() == im.getNumCols() );
99        assert( imC1.getNumRows() == im.getNumRows() );   
100       
101        OsiClpSolverInterface imC2(im);
102        assert( imC2.getModelPtr()!=im.getModelPtr() );
103        assert( imC2.getNumCols() == im.getNumCols() );
104        assert( imC2.getNumRows() == im.getNumRows() ); 
105       
106        assert( imC2.getModelPtr()!=imC1.getModelPtr() );
107       
108        lhs=imC2;
109      }
110      // Test that lhs has correct values even though rhs has gone out of scope
111     
112      assert( lhs.getModelPtr() != m.getModelPtr() );
113      assert( lhs.getNumCols() == m.getNumCols() );
114      assert( lhs.getNumRows() == m.getNumRows() );     
115    }
116    // Test clone
117    {
118      OsiClpSolverInterface oslSi(m);
119      OsiSolverInterface * siPtr = &oslSi;
120      OsiSolverInterface * siClone = siPtr->clone();
121      OsiClpSolverInterface * oslClone = dynamic_cast<OsiClpSolverInterface*>(siClone);
122      assert( oslClone != NULL );
123      assert( oslClone->getModelPtr() != oslSi.getModelPtr() );
124      assert( oslClone->getNumRows() == oslSi.getNumRows() );
125      assert( oslClone->getNumCols() == m.getNumCols() );
126     
127      delete siClone;
128    }
129 
130    // test infinity
131    {
132      OsiClpSolverInterface si;
133      assert( eq(si.getInfinity(),OsiClpInfinity));
134    }     
135   
136    // Test setting solution
137    {
138      OsiClpSolverInterface m1(m);
139      int i;
140
141      double * cs = new double[m1.getNumCols()];
142      for ( i = 0;  i < m1.getNumCols();  i++ ) 
143        cs[i] = i + .5;
144      m1.setColSolution(cs);
145      for ( i = 0;  i < m1.getNumCols();  i++ ) 
146        assert(m1.getColSolution()[i] == i + .5);
147     
148      double * rs = new double[m1.getNumRows()];
149      for ( i = 0;  i < m1.getNumRows();  i++ ) 
150        rs[i] = i - .5;
151      m1.setRowPrice(rs);
152      for ( i = 0;  i < m1.getNumRows();  i++ ) 
153        assert(m1.getRowPrice()[i] == i - .5);
154
155      delete [] cs;
156      delete [] rs;
157    }
158   
159   
160    // Test fraction Indices
161    {
162      OsiClpSolverInterface fim;
163      std::string fn = mpsDir+"exmip1";
164      fim.readMps(fn.c_str(),"mps");
165      // exmip1.mps has 2 integer variables with index 2 & 3
166      assert(  fim.isContinuous(0) );
167      assert(  fim.isContinuous(1) );
168      assert( !fim.isContinuous(2) );
169      assert( !fim.isContinuous(3) );
170      assert(  fim.isContinuous(4) );
171     
172      assert( !fim.isInteger(0) );
173      assert( !fim.isInteger(1) );
174      assert(  fim.isInteger(2) );
175      assert(  fim.isInteger(3) );
176      assert( !fim.isInteger(4) );
177     
178      assert( !fim.isBinary(0) );
179      assert( !fim.isBinary(1) );
180      assert(  fim.isBinary(2) );
181      assert(  fim.isBinary(3) );
182      assert( !fim.isBinary(4) );
183     
184      assert( !fim.isIntegerNonBinary(0) );
185      assert( !fim.isIntegerNonBinary(1) );
186      assert( !fim.isIntegerNonBinary(2) );
187      assert( !fim.isIntegerNonBinary(3) );
188      assert( !fim.isIntegerNonBinary(4) );
189      // Test fractionalIndices
190      {
191        double sol[]={2.9,3.0};
192        memcpy(fim.modelPtr_->primalColumnSolution()+2,sol,2*sizeof(double));
193        OsiVectorInt fi = fim.getFractionalIndices(1e-5);
194        assert( fi.size() == 1 );
195        assert( fi[0]==2 );
196       
197        // Set integer variables very close to integer values
198        sol[0]=5 + .00001/2.;
199        sol[1]=8 - .00001/2.;
200        memcpy(fim.modelPtr_->primalColumnSolution()+2,sol,2*sizeof(double));
201        fi = fim.getFractionalIndices(1e-5);
202        assert( fi.size() == 0 );
203       
204        // Set integer variables close, but beyond tolerances
205        sol[0]=5 + .00001*2.;
206        sol[1]=8 - .00001*2.;
207        memcpy(fim.modelPtr_->primalColumnSolution()+2,sol,2*sizeof(double));
208        fi = fim.getFractionalIndices(1e-5);
209        assert( fi.size() == 2 );
210        assert( fi[0]==2 );
211        assert( fi[1]==3 );
212      }
213
214
215     
216      // Change date so column 2 & 3 are integerNonBinary
217      double ub[]={5.0,6.0};
218      memcpy(fim.modelPtr_->columnUpper()+2,ub,2*sizeof(double));
219      assert( !fim.isBinary(0) );
220      assert( !fim.isBinary(1) );
221      assert( !fim.isBinary(2) );
222      assert( !fim.isBinary(3) );
223      assert( !fim.isBinary(4) );
224     
225      assert( !fim.isIntegerNonBinary(0) );
226      assert( !fim.isIntegerNonBinary(1) );
227      assert(  fim.isIntegerNonBinary(2) );
228      assert(  fim.isIntegerNonBinary(3) );
229      assert( !fim.isIntegerNonBinary(4) );
230    }
231    // Test some catches
232      if (!OsiClpHasNDEBUG())
233    {
234      OsiClpSolverInterface solver;
235#ifndef NDEBUG  /* in optimized mode, the following code crashes */
236      try {
237        solver.setObjCoeff(0,0.0);
238      }
239      catch (CoinError e) {
240        std::cout<<"Correct throw"<<std::endl;
241      }
242#endif
243      std::string fn = mpsDir+"exmip1";
244      solver.readMps(fn.c_str(),"mps");
245      try {
246        solver.setObjCoeff(0,0.0);
247      }
248      catch (CoinError e) {
249        std::cout<<"** Incorrect throw"<<std::endl;
250        abort();
251      }
252#ifndef NDEBUG
253      try {
254        int index[]={0,20};
255        double value[]={0.0,0.0,0.0,0.0};
256        solver.setColSetBounds(index,index+2,value);
257      }
258      catch (CoinError e) {
259        std::cout<<"Correct throw"<<std::endl;
260      }
261#endif
262    }
263    // Test apply cuts method
264    {     
265      OsiClpSolverInterface im(m);
266      OsiCuts cuts;
267     
268      // Generate some cuts
269      {
270        // Get number of rows and columns in model
271        int nr=im.getNumRows();
272        int nc=im.getNumCols();
273        assert( nr == 5 );
274        assert( nc == 8 );
275       
276        // Generate a valid row cut from thin air
277        int c;
278        {
279          int *inx = new int[nc];
280          for (c=0;c<nc;c++) inx[c]=c;
281          double *el = new double[nc];
282          for (c=0;c<nc;c++) el[c]=((double)c)*((double)c);
283         
284          OsiRowCut rc;
285          rc.setRow(nc,inx,el);
286          rc.setLb(-100.);
287          rc.setUb(100.);
288          rc.setEffectiveness(22);
289         
290          cuts.insert(rc);
291          delete[]el;
292          delete[]inx;
293        }
294       
295        // Generate valid col cut from thin air
296        {
297          const double * oslColLB = im.getColLower();
298          const double * oslColUB = im.getColUpper();
299          int *inx = new int[nc];
300          for (c=0;c<nc;c++) inx[c]=c;
301          double *lb = new double[nc];
302          double *ub = new double[nc];
303          for (c=0;c<nc;c++) lb[c]=oslColLB[c]+0.001;
304          for (c=0;c<nc;c++) ub[c]=oslColUB[c]-0.001;
305         
306          OsiColCut cc;
307          cc.setLbs(nc,inx,lb);
308          cc.setUbs(nc,inx,ub);
309         
310          cuts.insert(cc);
311          delete [] ub;
312          delete [] lb;
313          delete [] inx;
314        }
315       
316        {
317          // Generate a row and column cut which are ineffective
318          OsiRowCut * rcP= new OsiRowCut;
319          rcP->setEffectiveness(-1.);
320          cuts.insert(rcP);
321          assert(rcP==NULL);
322         
323          OsiColCut * ccP= new OsiColCut;
324          ccP->setEffectiveness(-12.);
325          cuts.insert(ccP);
326          assert(ccP==NULL);
327        }
328        {
329          //Generate inconsistent Row cut
330          OsiRowCut rc;
331          const int ne=1;
332          int inx[ne]={-10};
333          double el[ne]={2.5};
334          rc.setRow(ne,inx,el);
335          rc.setLb(3.);
336          rc.setUb(4.);
337          assert(!rc.consistent());
338          cuts.insert(rc);
339        }
340        {
341          //Generate inconsistent col cut
342          OsiColCut cc;
343          const int ne=1;
344          int inx[ne]={-10};
345          double el[ne]={2.5};
346          cc.setUbs(ne,inx,el);
347          assert(!cc.consistent());
348          cuts.insert(cc);
349        }
350        {
351          // Generate row cut which is inconsistent for model m
352          OsiRowCut rc;
353          const int ne=1;
354          int inx[ne]={10};
355          double el[ne]={2.5};
356          rc.setRow(ne,inx,el);
357          assert(rc.consistent());
358          assert(!rc.consistent(im));
359          cuts.insert(rc);
360        }
361        {
362          // Generate col cut which is inconsistent for model m
363          OsiColCut cc;
364          const int ne=1;
365          int inx[ne]={30};
366          double el[ne]={2.0};
367          cc.setLbs(ne,inx,el);
368          assert(cc.consistent());
369          assert(!cc.consistent(im));
370          cuts.insert(cc);
371        }
372        {
373          // Generate col cut which is infeasible
374          OsiColCut cc;
375          const int ne=1;
376          int inx[ne]={0};
377          double el[ne]={2.0};
378          cc.setUbs(ne,inx,el);
379          cc.setEffectiveness(1000.);
380          assert(cc.consistent());
381          assert(cc.consistent(im));
382          assert(cc.infeasible(im));
383          cuts.insert(cc);
384        }
385      }
386      assert(cuts.sizeRowCuts()==4);
387      assert(cuts.sizeColCuts()==5);
388     
389      OsiSolverInterface::ApplyCutsReturnCode rc = im.applyCuts(cuts);
390      assert( rc.getNumIneffective() == 2 );
391      assert( rc.getNumApplied() == 2 );
392      assert( rc.getNumInfeasible() == 1 );
393      assert( rc.getNumInconsistentWrtIntegerModel() == 2 );
394      assert( rc.getNumInconsistent() == 2 );
395      assert( cuts.sizeCuts() == rc.getNumIneffective() +
396        rc.getNumApplied() +
397        rc.getNumInfeasible() +
398        rc.getNumInconsistentWrtIntegerModel() +
399        rc.getNumInconsistent() );
400    }
401   
402    {   
403      OsiClpSolverInterface oslSi(m);
404      int nc = oslSi.getNumCols();
405      int nr = oslSi.getNumRows();
406      const double * cl = oslSi.getColLower();
407      const double * cu = oslSi.getColUpper();
408      const double * rl = oslSi.getRowLower();
409      const double * ru = oslSi.getRowUpper();
410      assert( nc == 8 );
411      assert( nr == 5 );
412      assert( eq(cl[0],2.5) );
413      assert( eq(cl[1],0.0) );
414      assert( eq(cu[1],4.1) );
415      assert( eq(cu[2],1.0) );
416      assert( eq(rl[0],2.5) );
417      assert( eq(rl[4],3.0) );
418      assert( eq(ru[1],2.1) );
419      assert( eq(ru[4],15.0) );
420     
421      const double * cs = oslSi.getColSolution();
422      assert( eq(cs[0],2.5) );
423      assert( eq(cs[7],0.0) );
424     
425      assert( !eq(cl[3],1.2345) );
426      oslSi.setColLower( 3, 1.2345 );
427      assert( eq(oslSi.getColLower()[3],1.2345) );
428     
429      assert( !eq(cu[4],10.2345) );
430      oslSi.setColUpper( 4, 10.2345 );
431      assert( eq(oslSi.getColUpper()[4],10.2345) );
432
433      double objValue = oslSi.getObjValue();
434      assert( eq(objValue,3.5) );
435
436      assert( eq( oslSi.getObjCoefficients()[0],  1.0) );
437      assert( eq( oslSi.getObjCoefficients()[1],  0.0) );
438      assert( eq( oslSi.getObjCoefficients()[2],  0.0) );
439      assert( eq( oslSi.getObjCoefficients()[3],  0.0) );
440      assert( eq( oslSi.getObjCoefficients()[4],  2.0) );
441      assert( eq( oslSi.getObjCoefficients()[5],  0.0) );
442      assert( eq( oslSi.getObjCoefficients()[6],  0.0) );
443      assert( eq( oslSi.getObjCoefficients()[7], -1.0) );
444    }
445   
446    // Test matrixByRow method
447    { 
448      const OsiClpSolverInterface si(m);
449      const CoinPackedMatrix * smP = si.getMatrixByRow();
450      // LL:      const OsiClpPackedMatrix * osmP = dynamic_cast<const OsiClpPackedMatrix*>(smP);
451      // LL: assert( osmP!=NULL );
452     
453      CoinRelFltEq eq;
454      const double * ev = smP->getElements();
455      assert( eq(ev[0],   3.0) );
456      assert( eq(ev[1],   1.0) );
457      assert( eq(ev[2],  -2.0) );
458      assert( eq(ev[3],  -1.0) );
459      assert( eq(ev[4],  -1.0) );
460      assert( eq(ev[5],   2.0) );
461      assert( eq(ev[6],   1.1) );
462      assert( eq(ev[7],   1.0) );
463      assert( eq(ev[8],   1.0) );
464      assert( eq(ev[9],   2.8) );
465      assert( eq(ev[10], -1.2) );
466      assert( eq(ev[11],  5.6) );
467      assert( eq(ev[12],  1.0) );
468      assert( eq(ev[13],  1.9) );
469     
470      const CoinBigIndex * mi = smP->getVectorStarts();
471      assert( mi[0]==0 );
472      assert( mi[1]==5 );
473      assert( mi[2]==7 );
474      assert( mi[3]==9 );
475      assert( mi[4]==11 );
476      assert( mi[5]==14 );
477     
478      const int * ei = smP->getIndices();
479      assert( ei[0]  ==  0 );
480      assert( ei[1]  ==  1 );
481      assert( ei[2]  ==  3 );
482      assert( ei[3]  ==  4 );
483      assert( ei[4]  ==  7 );
484      assert( ei[5]  ==  1 );
485      assert( ei[6]  ==  2 );
486      assert( ei[7]  ==  2 );
487      assert( ei[8]  ==  5 );
488      assert( ei[9]  ==  3 );
489      assert( ei[10] ==  6 );
490      assert( ei[11] ==  0 );
491      assert( ei[12] ==  4 );
492      assert( ei[13] ==  7 );   
493     
494      assert( smP->getMajorDim() == 5 ); 
495      assert( smP->getNumElements() == 14 );
496     
497    }
498    // Test adding several cuts
499    {
500      OsiClpSolverInterface fim;
501      std::string fn = mpsDir+"exmip1";
502      fim.readMps(fn.c_str(),"mps");
503      // exmip1.mps has 2 integer variables with index 2 & 3
504      fim.initialSolve();
505      OsiRowCut cuts[3];
506     
507     
508      // Generate one ineffective cut plus two trivial cuts
509      int c;
510      int nc = fim.getNumCols();
511      int *inx = new int[nc];
512      for (c=0;c<nc;c++) inx[c]=c;
513      double *el = new double[nc];
514      for (c=0;c<nc;c++) el[c]=1.0e-50+((double)c)*((double)c);
515     
516      cuts[0].setRow(nc,inx,el);
517      cuts[0].setLb(-100.);
518      cuts[0].setUb(500.);
519      cuts[0].setEffectiveness(22);
520      el[4]=0.0; // to get inf later
521     
522      for (c=2;c<4;c++) {
523        el[0]=1.0;
524        inx[0]=c;
525        cuts[c-1].setRow(1,inx,el);
526        cuts[c-1].setLb(1.);
527        cuts[c-1].setUb(100.);
528        cuts[c-1].setEffectiveness(c);
529      }
530      fim.writeMps("x1.mps");
531      fim.applyRowCuts(3,cuts);
532      fim.writeMps("x2.mps");
533      // resolve - should get message about zero elements
534      fim.resolve();
535      fim.writeMps("x3.mps");
536      // check integer solution
537      const double * cs = fim.getColSolution();
538      CoinRelFltEq eq;
539      assert( eq(cs[2],   1.0) );
540      assert( eq(cs[3],   1.0) );
541      // check will find invalid matrix
542      el[0]=1.0/el[4];
543      inx[0]=0;
544      cuts[0].setRow(nc,inx,el);
545      cuts[0].setLb(-100.);
546      cuts[0].setUb(500.);
547      cuts[0].setEffectiveness(22);
548      fim.applyRowCut(cuts[0]);
549      // resolve - should get message about zero elements
550      fim.resolve();
551      assert (fim.isAbandoned());
552      delete[]el;
553      delete[]inx;
554    }
555    // Test using bad basis
556    {
557      OsiClpSolverInterface fim;
558      std::string fn = mpsDir+"exmip1";
559      fim.readMps(fn.c_str(),"mps");
560      fim.initialSolve();
561      int numberRows = fim.getModelPtr()->numberRows();
562      int * rowStatus = new int[numberRows];
563      int numberColumns = fim.getModelPtr()->numberColumns();
564      int * columnStatus = new int[numberColumns];
565      fim.getBasisStatus(columnStatus,rowStatus);
566      int i;
567      for (i=0;i<numberRows;i++) {
568        rowStatus[i]=2;
569      }       
570      fim.setBasisStatus(columnStatus,rowStatus);
571      fim.resolve();
572      delete [] rowStatus;
573      delete [] columnStatus;
574    }
575        // Test matrixByCol method
576    {
577 
578      const OsiClpSolverInterface si(m);
579      const CoinPackedMatrix * smP = si.getMatrixByCol();
580      // LL:      const OsiClpPackedMatrix * osmP = dynamic_cast<const OsiClpPackedMatrix*>(smP);
581      // LL: assert( osmP!=NULL );
582     
583      CoinRelFltEq eq;
584      const double * ev = smP->getElements();
585      assert( eq(ev[0],   3.0) );
586      assert( eq(ev[1],   5.6) );
587      assert( eq(ev[2],   1.0) );
588      assert( eq(ev[3],   2.0) );
589      assert( eq(ev[4],   1.1) );
590      assert( eq(ev[5],   1.0) );
591      assert( eq(ev[6],  -2.0) );
592      assert( eq(ev[7],   2.8) );
593      assert( eq(ev[8],  -1.0) );
594      assert( eq(ev[9],   1.0) );
595      assert( eq(ev[10],  1.0) );
596      assert( eq(ev[11], -1.2) );
597      assert( eq(ev[12], -1.0) );
598      assert( eq(ev[13],  1.9) );
599     
600      const CoinBigIndex * mi = smP->getVectorStarts();
601      assert( mi[0]==0 );
602      assert( mi[1]==2 );
603      assert( mi[2]==4 );
604      assert( mi[3]==6 );
605      assert( mi[4]==8 );
606      assert( mi[5]==10 );
607      assert( mi[6]==11 );
608      assert( mi[7]==12 );
609      assert( mi[8]==14 );
610     
611      const int * ei = smP->getIndices();
612      assert( ei[0]  ==  0 );
613      assert( ei[1]  ==  4 );
614      assert( ei[2]  ==  0 );
615      assert( ei[3]  ==  1 );
616      assert( ei[4]  ==  1 );
617      assert( ei[5]  ==  2 );
618      assert( ei[6]  ==  0 );
619      assert( ei[7]  ==  3 );
620      assert( ei[8]  ==  0 );
621      assert( ei[9]  ==  4 );
622      assert( ei[10] ==  2 );
623      assert( ei[11] ==  3 );
624      assert( ei[12] ==  0 );
625      assert( ei[13] ==  4 );   
626     
627      assert( smP->getMajorDim() == 8 ); 
628      assert( smP->getNumElements() == 14 );
629
630      assert( smP->getSizeVectorStarts()==9 );
631      assert( smP->getMinorDim() == 5 );
632     
633    }
634    //--------------
635    // Test rowsense, rhs, rowrange, matrixByRow
636    {
637      OsiClpSolverInterface lhs;
638      {     
639        assert( m.rowrange_==NULL );
640        assert( m.rowsense_==NULL );
641        assert( m.rhs_==NULL );
642        assert( m.matrixByRow_==NULL );
643       
644        OsiClpSolverInterface siC1(m);     
645        assert( siC1.rowrange_==NULL );
646        assert( siC1.rowsense_==NULL );
647        assert( siC1.rhs_==NULL );
648        assert( siC1.matrixByRow_==NULL );
649
650        const char   * siC1rs  = siC1.getRowSense();
651        assert( siC1rs[0]=='G' );
652        assert( siC1rs[1]=='L' );
653        assert( siC1rs[2]=='E' );
654        assert( siC1rs[3]=='R' );
655        assert( siC1rs[4]=='R' );
656       
657        const double * siC1rhs = siC1.getRightHandSide();
658        assert( eq(siC1rhs[0],2.5) );
659        assert( eq(siC1rhs[1],2.1) );
660        assert( eq(siC1rhs[2],4.0) );
661        assert( eq(siC1rhs[3],5.0) );
662        assert( eq(siC1rhs[4],15.) ); 
663       
664        const double * siC1rr  = siC1.getRowRange();
665        assert( eq(siC1rr[0],0.0) );
666        assert( eq(siC1rr[1],0.0) );
667        assert( eq(siC1rr[2],0.0) );
668        assert( eq(siC1rr[3],5.0-1.8) );
669        assert( eq(siC1rr[4],15.0-3.0) );
670       
671        const CoinPackedMatrix * siC1mbr = siC1.getMatrixByRow();
672        assert( siC1mbr != NULL );
673       
674        const double * ev = siC1mbr->getElements();
675        assert( eq(ev[0],   3.0) );
676        assert( eq(ev[1],   1.0) );
677        assert( eq(ev[2],  -2.0) );
678        assert( eq(ev[3],  -1.0) );
679        assert( eq(ev[4],  -1.0) );
680        assert( eq(ev[5],   2.0) );
681        assert( eq(ev[6],   1.1) );
682        assert( eq(ev[7],   1.0) );
683        assert( eq(ev[8],   1.0) );
684        assert( eq(ev[9],   2.8) );
685        assert( eq(ev[10], -1.2) );
686        assert( eq(ev[11],  5.6) );
687        assert( eq(ev[12],  1.0) );
688        assert( eq(ev[13],  1.9) );
689       
690        const CoinBigIndex * mi = siC1mbr->getVectorStarts();
691        assert( mi[0]==0 );
692        assert( mi[1]==5 );
693        assert( mi[2]==7 );
694        assert( mi[3]==9 );
695        assert( mi[4]==11 );
696        assert( mi[5]==14 );
697       
698        const int * ei = siC1mbr->getIndices();
699        assert( ei[0]  ==  0 );
700        assert( ei[1]  ==  1 );
701        assert( ei[2]  ==  3 );
702        assert( ei[3]  ==  4 );
703        assert( ei[4]  ==  7 );
704        assert( ei[5]  ==  1 );
705        assert( ei[6]  ==  2 );
706        assert( ei[7]  ==  2 );
707        assert( ei[8]  ==  5 );
708        assert( ei[9]  ==  3 );
709        assert( ei[10] ==  6 );
710        assert( ei[11] ==  0 );
711        assert( ei[12] ==  4 );
712        assert( ei[13] ==  7 );   
713       
714        assert( siC1mbr->getMajorDim() == 5 ); 
715        assert( siC1mbr->getNumElements() == 14 );
716       
717
718        assert( siC1rs  == siC1.getRowSense() );
719        assert( siC1rhs == siC1.getRightHandSide() );
720        assert( siC1rr  == siC1.getRowRange() );
721
722        // Change OSL Model by adding free row
723        OsiRowCut rc;
724        rc.setLb(-DBL_MAX);
725        rc.setUb( DBL_MAX);
726        OsiCuts cuts;
727        cuts.insert(rc);
728        siC1.applyCuts(cuts);
729             
730        // Since model was changed, test that cached
731        // data is now freed.
732        assert( siC1.rowrange_==NULL );
733        assert( siC1.rowsense_==NULL );
734        assert( siC1.rhs_==NULL );
735        // now updated! assert( siC1.matrixByRow_==NULL );
736       
737        siC1rs  = siC1.getRowSense();
738        assert( siC1rs[0]=='G' );
739        assert( siC1rs[1]=='L' );
740        assert( siC1rs[2]=='E' );
741        assert( siC1rs[3]=='R' );
742        assert( siC1rs[4]=='R' );
743        assert( siC1rs[5]=='N' );
744
745        siC1rhs = siC1.getRightHandSide();
746        assert( eq(siC1rhs[0],2.5) );
747        assert( eq(siC1rhs[1],2.1) );
748        assert( eq(siC1rhs[2],4.0) );
749        assert( eq(siC1rhs[3],5.0) );
750        assert( eq(siC1rhs[4],15.) ); 
751        assert( eq(siC1rhs[5],0.0 ) ); 
752
753        siC1rr  = siC1.getRowRange();
754        assert( eq(siC1rr[0],0.0) );
755        assert( eq(siC1rr[1],0.0) );
756        assert( eq(siC1rr[2],0.0) );
757        assert( eq(siC1rr[3],5.0-1.8) );
758        assert( eq(siC1rr[4],15.0-3.0) );
759        assert( eq(siC1rr[5],0.0) );
760   
761        lhs=siC1;
762      }
763      // Test that lhs has correct values even though siC1 has gone out of scope   
764      assert( lhs.rowrange_==NULL );
765      assert( lhs.rowsense_==NULL );
766      assert( lhs.rhs_==NULL ); 
767      assert( lhs.matrixByRow_==NULL ); 
768     
769      const char * lhsrs  = lhs.getRowSense();
770      assert( lhsrs[0]=='G' );
771      assert( lhsrs[1]=='L' );
772      assert( lhsrs[2]=='E' );
773      assert( lhsrs[3]=='R' );
774      assert( lhsrs[4]=='R' );
775      assert( lhsrs[5]=='N' );
776     
777      const double * lhsrhs = lhs.getRightHandSide();
778      assert( eq(lhsrhs[0],2.5) );
779      assert( eq(lhsrhs[1],2.1) );
780      assert( eq(lhsrhs[2],4.0) );
781      assert( eq(lhsrhs[3],5.0) );
782      assert( eq(lhsrhs[4],15.) ); 
783      assert( eq(lhsrhs[5],0.0) ); 
784     
785      const double *lhsrr  = lhs.getRowRange();
786      assert( eq(lhsrr[0],0.0) );
787      assert( eq(lhsrr[1],0.0) );
788      assert( eq(lhsrr[2],0.0) );
789      assert( eq(lhsrr[3],5.0-1.8) );
790      assert( eq(lhsrr[4],15.0-3.0) );
791      assert( eq(lhsrr[5],0.0) );     
792     
793      const CoinPackedMatrix * lhsmbr = lhs.getMatrixByRow();
794      assert( lhsmbr != NULL );       
795      const double * ev = lhsmbr->getElements();
796      assert( eq(ev[0],   3.0) );
797      assert( eq(ev[1],   1.0) );
798      assert( eq(ev[2],  -2.0) );
799      assert( eq(ev[3],  -1.0) );
800      assert( eq(ev[4],  -1.0) );
801      assert( eq(ev[5],   2.0) );
802      assert( eq(ev[6],   1.1) );
803      assert( eq(ev[7],   1.0) );
804      assert( eq(ev[8],   1.0) );
805      assert( eq(ev[9],   2.8) );
806      assert( eq(ev[10], -1.2) );
807      assert( eq(ev[11],  5.6) );
808      assert( eq(ev[12],  1.0) );
809      assert( eq(ev[13],  1.9) );
810     
811      const CoinBigIndex * mi = lhsmbr->getVectorStarts();
812      assert( mi[0]==0 );
813      assert( mi[1]==5 );
814      assert( mi[2]==7 );
815      assert( mi[3]==9 );
816      assert( mi[4]==11 );
817      assert( mi[5]==14 );
818     
819      const int * ei = lhsmbr->getIndices();
820      assert( ei[0]  ==  0 );
821      assert( ei[1]  ==  1 );
822      assert( ei[2]  ==  3 );
823      assert( ei[3]  ==  4 );
824      assert( ei[4]  ==  7 );
825      assert( ei[5]  ==  1 );
826      assert( ei[6]  ==  2 );
827      assert( ei[7]  ==  2 );
828      assert( ei[8]  ==  5 );
829      assert( ei[9]  ==  3 );
830      assert( ei[10] ==  6 );
831      assert( ei[11] ==  0 );
832      assert( ei[12] ==  4 );
833      assert( ei[13] ==  7 );   
834     
835      int md = lhsmbr->getMajorDim();
836      assert(  md == 6 ); 
837      assert( lhsmbr->getNumElements() == 14 );
838    }
839   
840  }
841
842  // Test add/delete columns
843  {   
844    OsiClpSolverInterface m;
845    std::string fn = mpsDir+"p0033";
846    m.readMps(fn.c_str(),"mps");
847    double inf = m.getInfinity();
848
849    CoinPackedVector c0;
850    c0.insert(0, 4);
851    c0.insert(1, 1);
852    m.addCol(c0, 0, inf, 3);
853    m.initialSolve();
854    double objValue = m.getObjValue();
855    CoinRelFltEq eq(1.0e-2);
856    assert( eq(objValue,2520.57) );
857    // Try deleting first column
858    int * d = new int[1];
859    d[0]=0;
860    m.deleteCols(1,d);
861    delete [] d;
862    d=NULL;
863    m.resolve();
864    objValue = m.getObjValue();
865    assert( eq(objValue,2520.57) );
866    // Try deleting column we added
867    int iCol = m.getNumCols()-1;
868    m.deleteCols(1,&iCol);
869    m.resolve();
870    objValue = m.getObjValue();
871    assert( eq(objValue,2520.57) );
872
873  }
874  // Test branch and bound
875  {   
876    OsiClpSolverInterface m;
877    std::string fn = mpsDir+"p0033";
878    m.readMps(fn.c_str(),"mps");
879    // reduce printout
880    m.setHintParam(OsiDoReducePrint,true,OsiHintTry);
881    // test maximization
882    int n  = m.getNumCols();
883    int i;
884    double * obj2 = new double[n];
885    const double * obj = m.getObjCoefficients();
886    for (i=0;i<n;i++) {
887      obj2[i]=-obj[i];
888    }
889    m.setObjective(obj2);
890    delete [] obj2;
891    m.setObjSense(-1.0);
892    // Save bounds
893    double * saveUpper = CoinCopyOfArray(m.getColUpper(),n);
894    double * saveLower = CoinCopyOfArray(m.getColLower(),n);
895    m.branchAndBound();
896    // reset bounds
897    m.setColUpper(saveUpper);
898    m.setColLower(saveLower);
899    delete [] saveUpper;
900    delete [] saveLower;
901    // reset cutoff - otherwise we won't get solution
902    m.setDblParam(OsiDualObjectiveLimit,-COIN_DBL_MAX);
903    m.branchAndBound();
904    double objValue = m.getObjValue();
905    CoinRelFltEq eq(1.0e-2);
906    assert( eq(objValue,-3089) );
907    const double * cs = m.getColSolution();
908    for ( i=0;i<n;i++) {
909      if (cs[i]>1.0e-7)
910        printf("%d has value %g\n",i,cs[i]);
911    }
912  }
913  // Test matt
914  if (fopen("../Mps/Infeas/galenet.mps","r")) {   
915    OsiClpSolverInterface m;
916    m.readMps("../Mps/Infeas/galenet","mps");
917    m.setHintParam(OsiDoPresolveInResolve, true, OsiHintDo);
918    m.resolve();
919   
920    std::vector<double *> rays = m.getDualRays(1);
921    std::cout << "Dual Ray: " << std::endl;
922    for(int i = 0; i < m.getNumRows(); i++){
923      if(fabs(rays[0][i]) > 0.00001)
924        std::cout << i << " : " << rays[0][i] << std::endl;
925    }
926   
927    std::cout << "isProvenOptimal = " << m.isProvenOptimal() << std::endl;
928    std::cout << "isProvenPrimalInfeasible = " << m.isProvenPrimalInfeasible()
929         << std::endl;
930   
931    delete [] rays[0];
932   
933  }
934  // Test infeasible bounds
935  {
936    OsiClpSolverInterface solver;
937    std::string fn = mpsDir+"exmip1";
938    solver.readMps(fn.c_str(),"mps");
939    int index[]={0};
940    double value[]={1.0,0.0};
941    solver.setColSetBounds(index,index+1,value);
942    solver.setHintParam(OsiDoPresolveInInitial, false, OsiHintDo);
943    solver.initialSolve();
944    assert (!solver.isProvenOptimal());
945  }
946
947  // Build a model
948  {   
949    OsiClpSolverInterface model;
950    std::string fn = mpsDir+"p0033";
951    model.readMps(fn.c_str(),"mps");
952    // Point to data
953    int numberRows = model.getNumRows();
954    const double * rowLower = model.getRowLower();
955    const double * rowUpper = model.getRowUpper();
956    int numberColumns = model.getNumCols();
957    const double * columnLower = model.getColLower();
958    const double * columnUpper = model.getColUpper();
959    const double * columnObjective = model.getObjCoefficients();
960    // get row copy
961    CoinPackedMatrix rowCopy = *model.getMatrixByRow();
962    const int * column = rowCopy.getIndices();
963    const int * rowLength = rowCopy.getVectorLengths();
964    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
965    const double * element = rowCopy.getElements();
966   
967    // solve
968    model.initialSolve();
969    // Now build new model
970    CoinModel build;
971    // Row bounds
972    int iRow;
973    for (iRow=0;iRow<numberRows;iRow++) {
974      build.setRowBounds(iRow,rowLower[iRow],rowUpper[iRow]);
975      build.setRowName(iRow,model.getRowName(iRow).c_str());
976    }
977    // Column bounds and objective
978    int iColumn;
979    for (iColumn=0;iColumn<numberColumns;iColumn++) {
980      build.setColumnLower(iColumn,columnLower[iColumn]);
981      build.setColumnUpper(iColumn,columnUpper[iColumn]);
982      build.setObjective(iColumn,columnObjective[iColumn]);
983      build.setColumnName(iColumn,model.getColName(iColumn).c_str());
984    }
985    // Adds elements one by one by row (backwards by row)
986    for (iRow=numberRows-1;iRow>=0;iRow--) {
987      int start = rowStart[iRow];
988      for (int j=start;j<start+rowLength[iRow];j++) 
989        build(iRow,column[j],element[j]);
990    }
991    // Now create Model
992    OsiClpSolverInterface model2;
993    model2.loadFromCoinModel(build);
994    model2.initialSolve();
995    // Save - should be continuous
996    model2.writeMps("continuous");
997    int * whichInteger = new int[numberColumns];
998    for (iColumn=0;iColumn<numberColumns;iColumn++) 
999      whichInteger[iColumn]=iColumn;
1000    // mark as integer
1001    model2.setInteger(whichInteger,numberColumns);
1002    delete [] whichInteger;
1003    // save - should be integer
1004    model2.writeMps("integer");
1005    // check names are there
1006    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1007      printf("%d name %s\n",iColumn,model2.getColName(iColumn).c_str());
1008    }
1009   
1010    // Now do with strings attached
1011    // Save build to show how to go over rows
1012    CoinModel saveBuild = build;
1013    build = CoinModel();
1014    // Column bounds
1015    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1016      build.setColumnLower(iColumn,columnLower[iColumn]);
1017      build.setColumnUpper(iColumn,columnUpper[iColumn]);
1018    }
1019    // Objective - half the columns as is and half with multiplier of "1.0+multiplier"
1020    // Pick up from saveBuild (for no reason at all)
1021    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1022      double value = saveBuild.objective(iColumn);
1023      if (iColumn*2<numberColumns) {
1024        build.setObjective(iColumn,columnObjective[iColumn]);
1025      } else {
1026        // create as string
1027        char temp[100];
1028        sprintf(temp,"%g + abs(%g*multiplier)",value,value);
1029        build.setObjective(iColumn,temp);
1030      }
1031    }
1032    // It then adds rows one by one but for half the rows sets their values
1033    //      with multiplier of "1.0+1.5*multiplier"
1034    for (iRow=0;iRow<numberRows;iRow++) {
1035      if (iRow*2<numberRows) {
1036        // add row in simple way
1037        int start = rowStart[iRow];
1038        build.addRow(rowLength[iRow],column+start,element+start,
1039                     rowLower[iRow],rowUpper[iRow]);
1040      } else {
1041        // As we have to add one by one let's get from saveBuild
1042        CoinModelLink triple=saveBuild.firstInRow(iRow);
1043        while (triple.column()>=0) {
1044          int iColumn = triple.column();
1045          if (iColumn*2<numberColumns) {
1046            // just value as normal
1047            build(iRow,triple.column(),triple.value());
1048          } else {
1049            // create as string
1050            char temp[100];
1051            sprintf(temp,"%g + (1.5*%g*multiplier)",triple.value(), triple.value());
1052            build(iRow,iColumn,temp);
1053          }
1054          triple=saveBuild.next(triple);
1055        }
1056        // but remember to do rhs
1057        build.setRowLower(iRow,rowLower[iRow]);
1058        build.setRowUpper(iRow,rowUpper[iRow]);
1059      }
1060    }
1061    // If small switch on error printing
1062    if (numberColumns<50)
1063      build.setLogLevel(1);
1064    int numberErrors=model2.loadFromCoinModel(build);
1065    // should fail as we never set multiplier
1066    assert (numberErrors);
1067    build.associateElement("multiplier",0.0);
1068    numberErrors=model2.loadFromCoinModel(build);
1069    assert (!numberErrors);
1070    model2.initialSolve();
1071    // It then loops with multiplier going from 0.0 to 2.0 in increments of 0.1
1072    for (double multiplier=0.0;multiplier<2.0;multiplier+= 0.1) {
1073      build.associateElement("multiplier",multiplier);
1074      numberErrors=model2.loadFromCoinModel(build,true);
1075      assert (!numberErrors);
1076      model2.resolve();
1077    }
1078  }
1079  // Solve an lp by hand
1080  {   
1081    OsiClpSolverInterface m;
1082    std::string fn = mpsDir+"p0033";
1083    m.readMps(fn.c_str(),"mps");
1084    m.setObjSense(-1.0);
1085    m.getModelPtr()->messageHandler()->setLogLevel(4);
1086    m.initialSolve();
1087    m.getModelPtr()->factorization()->maximumPivots(5);
1088    m.setObjSense(1.0);
1089    // clone to test pivot as well as primalPivot
1090    OsiSolverInterface * mm =m.clone();
1091    // enable special mode
1092    m.enableSimplexInterface(true);
1093    // need to do after clone
1094    mm->enableSimplexInterface(true);
1095    // we happen to know that variables are 0-1 and rows are L
1096    int numberIterations=0;
1097    int numberColumns = m.getNumCols();
1098    int numberRows = m.getNumRows();
1099    double * fakeCost = new double[numberColumns];
1100    double * duals = new double [numberRows];
1101    double * djs = new double [numberColumns];
1102    const double * solution = m.getColSolution();
1103    memcpy(fakeCost,m.getObjCoefficients(),numberColumns*sizeof(double));
1104    while (1) {
1105      const double * dj;
1106      const double * dual;
1107      if ((numberIterations&1)==0) {
1108        // use given ones
1109        dj = m.getReducedCost();
1110        dual = m.getRowPrice();
1111      } else {
1112        // create
1113        dj = djs;
1114        dual = duals;
1115        m.getReducedGradient(djs,duals,fakeCost);
1116      }
1117      int i;
1118      int colIn=9999;
1119      int direction=1;
1120      double best=1.0e-6;
1121      // find most negative reduced cost
1122      // Should check basic - but should be okay on this problem
1123      for (i=0;i<numberRows;i++) {
1124        double value=dual[i];
1125        if (value>best) {
1126          direction=-1;
1127          best=value;
1128          colIn=-i-1;
1129        }
1130      }
1131      for (i=0;i<numberColumns;i++) {
1132        double value=dj[i];
1133        if (value<-best&&solution[i]<1.0e-6) {
1134          direction=1;
1135          best=-value;
1136          colIn=i;
1137        } else if (value>best&&solution[i]>1.0-1.0e-6) {
1138          direction=-1;
1139          best=value;
1140          colIn=i;
1141        }
1142      }
1143      if (colIn==9999)
1144        break; // should be optimal
1145      int colOut;
1146      int outStatus;
1147      double theta;
1148      int returnCode=m.primalPivotResult(colIn,direction,
1149                                         colOut,outStatus,theta,NULL);
1150      assert (!returnCode);
1151      printf("out %d, direction %d theta %g\n",
1152             colOut,outStatus,theta);
1153      if (colIn!=colOut) {
1154        returnCode=mm->pivot(colIn,colOut,outStatus);
1155        assert (returnCode>=0);
1156      } else {
1157        // bound flip (so pivot does not make sense)
1158        returnCode=mm->primalPivotResult(colIn,direction,
1159                                         colOut,outStatus,theta,NULL);
1160        assert (!returnCode);
1161      }
1162      numberIterations++;
1163    }
1164    delete mm;
1165    delete [] fakeCost;
1166    delete [] duals;
1167    delete [] djs;
1168    // exit special mode
1169    m.disableSimplexInterface();
1170    m.getModelPtr()->messageHandler()->setLogLevel(4);
1171    m.messageHandler()->setLogLevel(0);
1172    m.resolve();
1173    assert (!m.getIterationCount());
1174    m.setObjSense(-1.0);
1175    m.initialSolve();
1176  }
1177  // Do parametrics on the objective by hand
1178  {   
1179    OsiClpSolverInterface m;
1180    std::string fn = mpsDir+"p0033";
1181    m.readMps(fn.c_str(),"mps");
1182    ClpSimplex * simplex = m.getModelPtr();
1183    simplex->messageHandler()->setLogLevel(4);
1184    m.initialSolve();
1185    simplex->factorization()->maximumPivots(5);
1186    simplex->messageHandler()->setLogLevel(63);
1187    m.setObjSense(1.0);
1188    // enable special mode
1189    m.enableSimplexInterface(true);
1190    int numberIterations=0;
1191    int numberColumns = m.getNumCols();
1192    int numberRows = m.getNumRows();
1193    double * changeCost = new double[numberColumns];
1194    double * duals = new double [numberRows];
1195    double * djs = new double [numberColumns];
1196    // As getReducedGradient mucks about with innards of Clp get arrays to save
1197    double * dualsNow = new double [numberRows];
1198    double * djsNow = new double [numberColumns];
1199   
1200    const double * solution = m.getColSolution();
1201    int i;
1202    // Set up change cost
1203    for (i=0;i<numberColumns;i++)
1204      changeCost[i]=1.0+0.1*i;;
1205    // Range of investigation
1206    double totalChange=100.0;
1207    double totalDone=0.0;
1208    while (true) {
1209      // Save current
1210      // (would be more accurate to start from scratch)
1211      memcpy(djsNow, m.getReducedCost(),numberColumns*sizeof(double));
1212      memcpy(dualsNow,m.getRowPrice(),numberRows*sizeof(double));
1213      // Get reduced gradient of changeCost
1214      m.getReducedGradient(djs,duals,changeCost);
1215      int colIn=9999;
1216      int direction=1;
1217      // We are going up to totalChange but we have done some
1218      double best=totalChange-totalDone;
1219      // find best ratio
1220      // Should check basic - but should be okay on this problem
1221      // We are cheating here as we know L rows
1222      // Really should be using L and U status but I modified from previous example
1223      for (i=0;i<numberRows;i++) {
1224        if (simplex->getRowStatus(i)==ClpSimplex::basic) {
1225          assert (fabs(dualsNow[i])<1.0e-4&&fabs(duals[i])<1.0e-4);
1226        } else {
1227          assert (dualsNow[i]<1.0e-4);
1228          if (duals[i]>1.0e-8) {
1229            if (dualsNow[i]+best*duals[i]>0.0) {
1230              best = CoinMax(-dualsNow[i]/duals[i],0.0);
1231              direction=-1;
1232              colIn=-i-1;
1233            }
1234          }
1235        }
1236      }
1237      for (i=0;i<numberColumns;i++) {
1238        if (simplex->getColumnStatus(i)==ClpSimplex::basic) {
1239          assert (fabs(djsNow[i])<1.0e-4&&fabs(djs[i])<1.0e-4);
1240        } else {
1241          if (solution[i]<1.0e-6) {
1242            assert (djsNow[i]>-1.0e-4);
1243            if (djs[i]<-1.0e-8) {
1244              if (djsNow[i]+best*djs[i]<0.0) {
1245                best = CoinMax(-djsNow[i]/djs[i],0.0);
1246                direction=1;
1247                colIn=i;
1248              }
1249            }
1250          } else if (solution[i]>1.0-1.0e-6) {
1251            assert (djsNow[i]<1.0e-4);
1252            if (djs[i]>1.0e-8) {
1253              if (djsNow[i]+best*djs[i]>0.0) {
1254                best = CoinMax(-djsNow[i]/djs[i],0.0);
1255                direction=-1;
1256                colIn=i;
1257              }
1258            }
1259          }
1260        }
1261      }
1262      if (colIn==9999)
1263        break; // should be optimal
1264      // update objective - djs is spare array
1265      const double * obj = m.getObjCoefficients();
1266      for (i=0;i<numberColumns;i++) {
1267        djs[i]=obj[i]+best*changeCost[i];
1268      }
1269      totalDone += best;
1270      printf("Best change %g, total %g\n",best,totalDone);
1271      m.setObjectiveAndRefresh(djs);
1272      int colOut;
1273      int outStatus;
1274      double theta;
1275      int returnCode=m.primalPivotResult(colIn,direction,colOut,outStatus,theta,NULL);
1276      assert (!returnCode);
1277      double objValue = m.getObjValue(); // printed one may not be accurate
1278      printf("in %d out %d, direction %d theta %g, objvalue %g\n",
1279             colIn,colOut,outStatus,theta,objValue);
1280      numberIterations++;
1281    }
1282    // update objective to totalChange- djs is spare array
1283    const double * obj = m.getObjCoefficients();
1284    double best = totalChange-totalDone;
1285    for (i=0;i<numberColumns;i++) {
1286      djs[i]=obj[i]+best*changeCost[i];
1287    }
1288    delete [] changeCost;
1289    delete [] duals;
1290    delete [] djs;
1291    delete [] dualsNow;
1292    delete [] djsNow;
1293    // exit special mode
1294    m.disableSimplexInterface();
1295    simplex->messageHandler()->setLogLevel(4);
1296    m.resolve();
1297    assert (!m.getIterationCount());
1298  }
1299  // Solve an lp when interface is on
1300  {   
1301    OsiClpSolverInterface m;
1302    std::string fn = mpsDir+"p0033";
1303    m.readMps(fn.c_str(),"mps");
1304    // enable special mode
1305    m.setHintParam(OsiDoScale,false,OsiHintDo);
1306    m.setHintParam(OsiDoPresolveInInitial,false,OsiHintDo);
1307    m.setHintParam(OsiDoDualInInitial,false,OsiHintDo);
1308    m.setHintParam(OsiDoPresolveInResolve,false,OsiHintDo);
1309    m.setHintParam(OsiDoDualInResolve,false,OsiHintDo);
1310    m.enableSimplexInterface(true);
1311    m.initialSolve();
1312  }
1313  // Check tableau stuff when simplex interface is on
1314  {   
1315    OsiClpSolverInterface m;
1316    /*
1317       Wolsey : Page 130
1318       max 4x1 -  x2
1319       7x1 - 2x2    <= 14
1320       x2    <= 3
1321       2x1 - 2x2    <= 3
1322       x1 in Z+, x2 >= 0
1323    */
1324   
1325    double inf_ = m.getInfinity();
1326    int n_cols = 2;
1327    int n_rows = 3;
1328   
1329    double obj[2] = {-4.0, 1.0};
1330    double collb[2] = {0.0, 0.0};
1331    double colub[2] = {inf_, inf_};
1332    double rowlb[3] = {-inf_, -inf_, -inf_};
1333    double rowub[3] = {14.0, 3.0, 3.0};
1334   
1335    int rowIndices[5] =  {0,     2,    0,    1,    2};
1336    int colIndices[5] =  {0,     0,    1,    1,    1};
1337    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
1338    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
1339   
1340    m.loadProblem(M, collb, colub, obj, rowlb, rowub);
1341    m.enableSimplexInterface(true);
1342   
1343    m.initialSolve();
1344   
1345    //check that the tableau matches wolsey (B-1 A)
1346    // slacks in second part of binvA
1347    double * binvA = (double*) malloc((n_cols+n_rows) * sizeof(double));
1348   
1349    printf("B-1 A");
1350    int i;
1351    for( i = 0; i < n_rows; i++){
1352      m.getBInvARow(i, binvA,binvA+n_cols);
1353      printf("\nrow: %d -> ",i);
1354      for(int j=0; j < n_cols+n_rows; j++){
1355        printf("%g, ", binvA[j]);
1356      }
1357    }
1358    printf("\n");
1359    printf("And by column");
1360    for( i = 0; i < n_cols+n_rows; i++){
1361      m.getBInvACol(i, binvA);
1362      printf("\ncolumn: %d -> ",i);
1363      for(int j=0; j < n_rows; j++){
1364        printf("%g, ", binvA[j]);
1365      }
1366    }
1367    printf("\n");
1368    // and when doing as expert
1369    m.setSpecialOptions(m.specialOptions()|512);
1370    ClpSimplex * clp = m.getModelPtr();
1371    /* Do twice -
1372       first time with enableSimplexInterface still set
1373       then without and with scaling
1374    */
1375    for (int iPass=0;iPass<2;iPass++) {
1376      const double * rowScale = clp->rowScale();
1377      const double * columnScale = clp->columnScale();
1378      if (!iPass)
1379        assert (!rowScale);
1380      else
1381        assert (rowScale); // only true for this example
1382      /* has to be exactly correct as in OsiClpsolverInterface.cpp
1383         (also redo each pass as may change
1384      */
1385      CoinIndexedVector * rowArray = clp->rowArray(1);
1386      CoinIndexedVector * columnArray = clp->columnArray(0);
1387      int n;
1388      int * which;
1389      double * array;
1390      printf("B-1 A");
1391      for( i = 0; i < n_rows; i++){
1392        m.getBInvARow(i, binvA,binvA+n_cols);
1393        printf("\nrow: %d -> ",i);
1394        int j;
1395        // First columns
1396        n = columnArray->getNumElements();
1397        which = columnArray->getIndices();
1398        array = columnArray->denseVector();
1399        for(j=0; j < n; j++){
1400          int k=which[j];
1401          if (!columnScale) {
1402            printf("(%d %g), ", k, array[k]);
1403          } else {
1404            printf("(%d %g), ", k, array[k]/columnScale[k]);
1405          }
1406          // zero out
1407          array[k]=0.0;
1408        }
1409        // say empty
1410        columnArray->setNumElements(0);
1411        // and check (would not be in any production code)
1412        columnArray->checkClear();
1413        // now rows
1414        n = rowArray->getNumElements();
1415        which = rowArray->getIndices();
1416        array = rowArray->denseVector();
1417        for(j=0; j < n; j++){
1418          int k=which[j];
1419          if (!rowScale) {
1420            printf("(%d %g), ", k+n_cols, array[k]);
1421          } else {
1422            printf("(%d %g), ", k+n_cols, array[k]*rowScale[k]);
1423          }
1424          // zero out
1425          array[k]=0.0;
1426        }
1427        // say empty
1428        rowArray->setNumElements(0);
1429        // and check (would not be in any production code)
1430        rowArray->checkClear();
1431      }
1432      printf("\n");
1433      printf("And by column (trickier)");
1434      const int * pivotVariable = clp->pivotVariable();
1435      for( i = 0; i < n_cols+n_rows; i++){
1436        m.getBInvACol(i, binvA);
1437        printf("\ncolumn: %d -> ",i);
1438        n = rowArray->getNumElements();
1439        which = rowArray->getIndices();
1440        array = rowArray->denseVector();
1441        for(int j=0; j < n; j++){
1442          int k=which[j];
1443          // need to know pivot variable for +1/-1 (slack) and row/column scaling
1444          int pivot = pivotVariable[k];
1445          if (pivot<n_cols) {
1446            // scaled coding is in just in case
1447            if (!columnScale) {
1448              printf("(%d %g), ", k, array[k]);
1449            } else {
1450              printf("(%d %g), ", k, array[k]*columnScale[pivot]);
1451            }
1452          } else {
1453            if (!rowScale) {
1454              printf("(%d %g), ", k, -array[k]);
1455            } else {
1456              printf("(%d %g), ", k, -array[k]/rowScale[pivot-n_cols]);
1457            }
1458          }
1459          // zero out
1460          array[k]=0.0;
1461        }
1462        // say empty
1463        rowArray->setNumElements(0);
1464        // and check (would not be in any production code)
1465        rowArray->checkClear();
1466      }
1467      printf("\n");
1468      // now deal with next pass
1469      if (!iPass) {
1470        m.disableSimplexInterface();
1471        // see if we can get scaling for testing
1472        clp->scaling(1);
1473        m.enableFactorization();
1474      } else {
1475        // may not be needed - but cleaner
1476        m.disableFactorization();
1477      }
1478    }
1479    m.setSpecialOptions(m.specialOptions()&~512);
1480    free(binvA);
1481    // Do using newer interface
1482    m.enableFactorization();
1483    {
1484      CoinIndexedVector * rowArray = new CoinIndexedVector(n_rows);
1485      CoinIndexedVector * columnArray = new CoinIndexedVector(n_cols);
1486      int n;
1487      int * which;
1488      double * array;
1489      printf("B-1 A");
1490      for( i = 0; i < n_rows; i++){
1491        m.getBInvARow(i, columnArray,rowArray);
1492        printf("\nrow: %d -> ",i);
1493        int j;
1494        // First columns
1495        n = columnArray->getNumElements();
1496        which = columnArray->getIndices();
1497        array = columnArray->denseVector();
1498        for(j=0; j < n; j++){
1499          int k=which[j];
1500          printf("(%d %g), ", k, array[k]);
1501          // zero out
1502          array[k]=0.0;
1503        }
1504        // say empty (if I had not zeroed array[k] I would use ->clear())
1505        columnArray->setNumElements(0);
1506        // and check (would not be in any production code)
1507        columnArray->checkClear();
1508        // now rows
1509        n = rowArray->getNumElements();
1510        which = rowArray->getIndices();
1511        array = rowArray->denseVector();
1512        for(j=0; j < n; j++){
1513          int k=which[j];
1514          printf("(%d %g), ", k+n_cols, array[k]);
1515          // zero out
1516          array[k]=0.0;
1517        }
1518        // say empty
1519        rowArray->setNumElements(0);
1520        // and check (would not be in any production code)
1521        rowArray->checkClear();
1522      }
1523      printf("\n");
1524      printf("And by column");
1525      const int * pivotVariable = clp->pivotVariable();
1526      for( i = 0; i < n_cols+n_rows; i++){
1527        m.getBInvACol(i, rowArray);
1528        printf("\ncolumn: %d -> ",i);
1529        n = rowArray->getNumElements();
1530        which = rowArray->getIndices();
1531        array = rowArray->denseVector();
1532        for(int j=0; j < n; j++){
1533          int k=which[j];
1534          printf("(%d %g), ", k, array[k]);
1535          // zero out
1536          array[k]=0.0;
1537        }
1538        // say empty
1539        rowArray->setNumElements(0);
1540        // and check (would not be in any production code)
1541        rowArray->checkClear();
1542      }
1543      printf("\n");
1544      delete rowArray;
1545      delete columnArray;
1546    }
1547    // may not be needed - but cleaner
1548    m.disableFactorization();
1549  }
1550  // Check tableau stuff when simplex interface is off
1551  {   
1552    OsiClpSolverInterface m;
1553    /*
1554       Wolsey : Page 130
1555       max 4x1 -  x2
1556       7x1 - 2x2    <= 14
1557       x2    <= 3
1558       2x1 - 2x2    <= 3
1559       x1 in Z+, x2 >= 0
1560    */
1561   
1562    double inf_ = m.getInfinity();
1563    int n_cols = 2;
1564    int n_rows = 3;
1565   
1566    double obj[2] = {-4.0, 1.0};
1567    double collb[2] = {0.0, 0.0};
1568    double colub[2] = {inf_, inf_};
1569    double rowlb[3] = {-inf_, -inf_, -inf_};
1570    double rowub[3] = {14.0, 3.0, 3.0};
1571   
1572    int rowIndices[5] =  {0,     2,    0,    1,    2};
1573    int colIndices[5] =  {0,     0,    1,    1,    1};
1574    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
1575    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
1576   
1577    m.loadProblem(M, collb, colub, obj, rowlb, rowub);
1578    // Test with scaling
1579    m.setHintParam(OsiDoScale,true);
1580    // Tell code to keep factorization
1581    m.setupForRepeatedUse(3,0);
1582   
1583    m.initialSolve();
1584    // Repeated stuff only in resolve
1585    // m.resolve(); (now in both (just for (3,0))
1586   
1587    //check that the tableau matches wolsey (B-1 A)
1588    // slacks in second part of binvA
1589    double * binvA = (double*) malloc((n_cols+n_rows) * sizeof(double));
1590    // and getBasics
1591    int * pivots = new int[n_rows];
1592    m.getBasics(pivots);
1593   
1594    printf("B-1 A");
1595    int i;
1596    for( i = 0; i < n_rows; i++){
1597      m.getBInvARow(i, binvA,binvA+n_cols);
1598      printf("\nrow: %d (pivot %d) -> ",i,pivots[i]);
1599      for(int j=0; j < n_cols+n_rows; j++){
1600        printf("%g, ", binvA[j]);
1601      }
1602    }
1603    printf("\n");
1604    printf("And by column");
1605    for( i = 0; i < n_cols+n_rows; i++){
1606      m.getBInvACol(i, binvA);
1607      printf("\ncolumn: %d -> ",i);
1608      for(int j=0; j < n_rows; j++){
1609        printf("%g, ", binvA[j]);
1610      }
1611    }
1612    printf("\n");
1613    free(binvA);
1614    delete [] pivots;
1615  }
1616  // Do common solverInterface testing
1617  {
1618    OsiClpSolverInterface m;
1619    return OsiSolverInterfaceCommonUnitTest(&m, mpsDir,netlibDir);
1620  }
1621}
Note: See TracBrowser for help on using the repository browser.