source: trunk/Clp/test/OsiClpSolverInterfaceTest.cpp @ 1687

Last change on this file since 1687 was 1663, checked in by lou, 9 years ago

Add EPL license notice in test.

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