source: branches/devel/Cbc/examples/OsiBranchLink.cpp @ 529

Last change on this file since 529 was 476, checked in by forrest, 13 years ago

osi stuff

File size: 44.3 KB
Line 
1// Copyright (C) 2006, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "OsiSolverInterface.hpp"
13#include "OsiBranchLink.hpp"
14#include "CoinError.hpp"
15#include "CoinHelperFunctions.hpp"
16#include "CoinPackedMatrix.hpp"
17#include "CoinWarmStartBasis.hpp"
18
19// Default Constructor
20OsiOldLink::OsiOldLink ()
21  : OsiSOS(),
22    numberLinks_(0)
23{
24}
25
26// Useful constructor (which are indices)
27OsiOldLink::OsiOldLink (const OsiSolverInterface * solver,  int numberMembers,
28           int numberLinks, int first , const double * weights, int identifier)
29  : OsiSOS(),
30    numberLinks_(numberLinks)
31{
32  numberMembers_ = numberMembers;
33  members_ = NULL;
34  sosType_ = 1;
35  if (numberMembers_) {
36    weights_ = new double[numberMembers_];
37    members_ = new int[numberMembers_*numberLinks_];
38    if (weights) {
39      memcpy(weights_,weights,numberMembers_*sizeof(double));
40    } else {
41      for (int i=0;i<numberMembers_;i++)
42        weights_[i]=i;
43    }
44    // weights must be increasing
45    int i;
46    double last=-COIN_DBL_MAX;
47    for (i=0;i<numberMembers_;i++) {
48      assert (weights_[i]>last+1.0e-12);
49      last=weights_[i];
50    }
51    for (i=0;i<numberMembers_*numberLinks_;i++) {
52      members_[i]=first+i;
53    }
54  } else {
55    weights_ = NULL;
56  }
57}
58
59// Useful constructor (which are indices)
60OsiOldLink::OsiOldLink (const OsiSolverInterface * solver,  int numberMembers,
61           int numberLinks, int sosType, const int * which , const double * weights, int identifier)
62  : OsiSOS(),
63    numberLinks_(numberLinks)
64{
65  numberMembers_ = numberMembers;
66  members_ = NULL;
67  sosType_ = 1;
68  if (numberMembers_) {
69    weights_ = new double[numberMembers_];
70    members_ = new int[numberMembers_*numberLinks_];
71    if (weights) {
72      memcpy(weights_,weights,numberMembers_*sizeof(double));
73    } else {
74      for (int i=0;i<numberMembers_;i++)
75        weights_[i]=i;
76    }
77    // weights must be increasing
78    int i;
79    double last=-COIN_DBL_MAX;
80    for (i=0;i<numberMembers_;i++) {
81      assert (weights_[i]>last+1.0e-12);
82      last=weights_[i];
83    }
84    for (i=0;i<numberMembers_*numberLinks_;i++) {
85      members_[i]= which[i];
86    }
87  } else {
88    weights_ = NULL;
89  }
90}
91
92// Copy constructor
93OsiOldLink::OsiOldLink ( const OsiOldLink & rhs)
94  :OsiSOS(rhs)
95{
96  numberLinks_ = rhs.numberLinks_;
97  if (numberMembers_) {
98    delete [] members_;
99    members_ = CoinCopyOfArray(rhs.members_,numberMembers_*numberLinks_);
100  }
101}
102
103// Clone
104OsiObject *
105OsiOldLink::clone() const
106{
107  return new OsiOldLink(*this);
108}
109
110// Assignment operator
111OsiOldLink & 
112OsiOldLink::operator=( const OsiOldLink& rhs)
113{
114  if (this!=&rhs) {
115    OsiSOS::operator=(rhs);
116    delete [] members_;
117    numberLinks_ = rhs.numberLinks_;
118    if (numberMembers_) {
119      members_ = CoinCopyOfArray(rhs.members_,numberMembers_*numberLinks_);
120    } else {
121      members_ = NULL;
122    }
123  }
124  return *this;
125}
126
127// Destructor
128OsiOldLink::~OsiOldLink ()
129{
130}
131
132// Infeasibility - large is 0.5
133double 
134OsiOldLink::infeasibility(const OsiBranchingInformation * info,int & whichWay) const
135{
136  int j;
137  int firstNonZero=-1;
138  int lastNonZero = -1;
139  const double * solution = info->solution_;
140  //const double * lower = info->lower_;
141  const double * upper = info->upper_;
142  double integerTolerance = info->integerTolerance_;
143  double weight = 0.0;
144  double sum =0.0;
145
146  // check bounds etc
147  double lastWeight=-1.0e100;
148  int base=0;
149  for (j=0;j<numberMembers_;j++) {
150    for (int k=0;k<numberLinks_;k++) {
151      int iColumn = members_[base+k];
152      if (lastWeight>=weights_[j]-1.0e-7)
153        throw CoinError("Weights too close together in OsiLink","infeasibility","OsiLink");
154      lastWeight = weights_[j];
155      double value = CoinMax(0.0,solution[iColumn]);
156      sum += value;
157      if (value>integerTolerance&&upper[iColumn]) {
158        // Possibly due to scaling a fixed variable might slip through
159        if (value>upper[iColumn]+1.0e-8) {
160#ifdef OSI_DEBUG
161          printf("** Variable %d (%d) has value %g and upper bound of %g\n",
162                 iColumn,j,value,upper[iColumn]);
163#endif
164        } 
165        value = CoinMin(value,upper[iColumn]);
166        weight += weights_[j]*value;
167        if (firstNonZero<0)
168          firstNonZero=j;
169        lastNonZero=j;
170      }
171    }
172    base += numberLinks_;
173  }
174  double valueInfeasibility;
175  whichWay=1;
176  whichWay_=1;
177  if (lastNonZero-firstNonZero>=sosType_) {
178    // find where to branch
179    assert (sum>0.0);
180    weight /= sum;
181    valueInfeasibility = lastNonZero-firstNonZero+1;
182    valueInfeasibility *= 0.5/((double) numberMembers_);
183    //#define DISTANCE
184#ifdef DISTANCE
185    assert (sosType_==1); // code up
186    /* may still be satisfied.
187       For LOS type 2 we might wish to move coding around
188       and keep initial info in model_ for speed
189    */
190    int iWhere;
191    bool possible=false;
192    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
193      if (fabs(weight-weights_[iWhere])<1.0e-8) {
194        possible=true;
195        break;
196      }
197    }
198    if (possible) {
199      // One could move some of this (+ arrays) into model_
200      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
201      const double * element = matrix->getMutableElements();
202      const int * row = matrix->getIndices();
203      const CoinBigIndex * columnStart = matrix->getVectorStarts();
204      const int * columnLength = matrix->getVectorLengths();
205      const double * rowSolution = solver->getRowActivity();
206      const double * rowLower = solver->getRowLower();
207      const double * rowUpper = solver->getRowUpper();
208      int numberRows = matrix->getNumRows();
209      double * array = new double [numberRows];
210      CoinZeroN(array,numberRows);
211      int * which = new int [numberRows];
212      int n=0;
213      int base=numberLinks_*firstNonZero;
214      for (j=firstNonZero;j<=lastNonZero;j++) {
215        for (int k=0;k<numberLinks_;k++) {
216          int iColumn = members_[base+k];
217          double value = CoinMax(0.0,solution[iColumn]);
218          if (value>integerTolerance&&upper[iColumn]) {
219            value = CoinMin(value,upper[iColumn]);
220            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
221              int iRow = row[j];
222              double a = array[iRow];
223              if (a) {
224                a += value*element[j];
225                if (!a)
226                  a = 1.0e-100;
227              } else {
228                which[n++]=iRow;
229                a=value*element[j];
230                assert (a);
231              }
232              array[iRow]=a;
233            }
234          }
235        }
236        base += numberLinks_;
237      }
238      base=numberLinks_*iWhere;
239      for (int k=0;k<numberLinks_;k++) {
240        int iColumn = members_[base+k];
241        const double value = 1.0;
242        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
243          int iRow = row[j];
244          double a = array[iRow];
245          if (a) {
246            a -= value*element[j];
247            if (!a)
248              a = 1.0e-100;
249          } else {
250            which[n++]=iRow;
251            a=-value*element[j];
252            assert (a);
253          }
254          array[iRow]=a;
255        }
256      }
257      for (j=0;j<n;j++) {
258        int iRow = which[j];
259        // moving to point will increase row solution by this
260        double distance = array[iRow];
261        if (distance>1.0e-8) {
262          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
263            possible=false;
264            break;
265          }
266        } else if (distance<-1.0e-8) {
267          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
268            possible=false;
269            break;
270          } 
271        }
272      }
273      for (j=0;j<n;j++)
274        array[which[j]]=0.0;
275      delete [] array;
276      delete [] which;
277      if (possible) {
278        valueInfeasibility=0.0;
279        printf("possible %d %d %d\n",firstNonZero,lastNonZero,iWhere);
280      }
281    }
282#endif
283  } else {
284    valueInfeasibility = 0.0; // satisfied
285  }
286  infeasibility_=valueInfeasibility;
287  otherInfeasibility_=1.0-valueInfeasibility;
288  return valueInfeasibility;
289}
290
291// This looks at solution and sets bounds to contain solution
292double
293OsiOldLink::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
294{
295  int j;
296  int firstNonZero=-1;
297  int lastNonZero = -1;
298  const double * solution = info->solution_;
299  const double * upper = info->upper_;
300  double integerTolerance = info->integerTolerance_;
301  double weight = 0.0;
302  double sum =0.0;
303
304  int base=0;
305  for (j=0;j<numberMembers_;j++) {
306    for (int k=0;k<numberLinks_;k++) {
307      int iColumn = members_[base+k];
308      double value = CoinMax(0.0,solution[iColumn]);
309      sum += value;
310      if (value>integerTolerance&&upper[iColumn]) {
311        weight += weights_[j]*value;
312        if (firstNonZero<0)
313          firstNonZero=j;
314        lastNonZero=j;
315      }
316    }
317    base += numberLinks_;
318  }
319#ifdef DISTANCE
320  if (lastNonZero-firstNonZero>sosType_-1) {
321    /* may still be satisfied.
322       For LOS type 2 we might wish to move coding around
323       and keep initial info in model_ for speed
324    */
325    int iWhere;
326    bool possible=false;
327    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
328      if (fabs(weight-weights_[iWhere])<1.0e-8) {
329        possible=true;
330        break;
331      }
332    }
333    if (possible) {
334      // One could move some of this (+ arrays) into model_
335      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
336      const double * element = matrix->getMutableElements();
337      const int * row = matrix->getIndices();
338      const CoinBigIndex * columnStart = matrix->getVectorStarts();
339      const int * columnLength = matrix->getVectorLengths();
340      const double * rowSolution = solver->getRowActivity();
341      const double * rowLower = solver->getRowLower();
342      const double * rowUpper = solver->getRowUpper();
343      int numberRows = matrix->getNumRows();
344      double * array = new double [numberRows];
345      CoinZeroN(array,numberRows);
346      int * which = new int [numberRows];
347      int n=0;
348      int base=numberLinks_*firstNonZero;
349      for (j=firstNonZero;j<=lastNonZero;j++) {
350        for (int k=0;k<numberLinks_;k++) {
351          int iColumn = members_[base+k];
352          double value = CoinMax(0.0,solution[iColumn]);
353          if (value>integerTolerance&&upper[iColumn]) {
354            value = CoinMin(value,upper[iColumn]);
355            for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
356              int iRow = row[j];
357              double a = array[iRow];
358              if (a) {
359                a += value*element[j];
360                if (!a)
361                  a = 1.0e-100;
362              } else {
363                which[n++]=iRow;
364                a=value*element[j];
365                assert (a);
366              }
367              array[iRow]=a;
368            }
369          }
370        }
371        base += numberLinks_;
372      }
373      base=numberLinks_*iWhere;
374      for (int k=0;k<numberLinks_;k++) {
375        int iColumn = members_[base+k];
376        const double value = 1.0;
377        for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
378          int iRow = row[j];
379          double a = array[iRow];
380          if (a) {
381            a -= value*element[j];
382            if (!a)
383              a = 1.0e-100;
384          } else {
385            which[n++]=iRow;
386            a=-value*element[j];
387            assert (a);
388          }
389          array[iRow]=a;
390        }
391      }
392      for (j=0;j<n;j++) {
393        int iRow = which[j];
394        // moving to point will increase row solution by this
395        double distance = array[iRow];
396        if (distance>1.0e-8) {
397          if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) {
398            possible=false;
399            break;
400          }
401        } else if (distance<-1.0e-8) {
402          if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) {
403            possible=false;
404            break;
405          } 
406        }
407      }
408      for (j=0;j<n;j++)
409        array[which[j]]=0.0;
410      delete [] array;
411      delete [] which;
412      if (possible) {
413        printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere);
414        firstNonZero=iWhere;
415        lastNonZero=iWhere;
416      }
417    }
418  }
419#else
420  assert (lastNonZero-firstNonZero<sosType_) ;
421#endif
422  base=0;
423  for (j=0;j<firstNonZero;j++) {
424    for (int k=0;k<numberLinks_;k++) {
425      int iColumn = members_[base+k];
426      solver->setColUpper(iColumn,0.0);
427    }
428    base += numberLinks_;
429  }
430  // skip
431  base += numberLinks_;
432  for (j=lastNonZero+1;j<numberMembers_;j++) {
433    for (int k=0;k<numberLinks_;k++) {
434      int iColumn = members_[base+k];
435      solver->setColUpper(iColumn,0.0);
436    }
437    base += numberLinks_;
438  }
439  // go to coding as in OsiSOS
440  abort();
441  return -1.0;
442}
443
444// Redoes data when sequence numbers change
445void 
446OsiOldLink::resetSequenceEtc(int numberColumns, const int * originalColumns)
447{
448  int n2=0;
449  for (int j=0;j<numberMembers_*numberLinks_;j++) {
450    int iColumn = members_[j];
451    int i;
452    for (i=0;i<numberColumns;i++) {
453      if (originalColumns[i]==iColumn)
454        break;
455    }
456    if (i<numberColumns) {
457      members_[n2]=i;
458      weights_[n2++]=weights_[j];
459    }
460  }
461  if (n2<numberMembers_) {
462    printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2/numberLinks_);
463    numberMembers_=n2/numberLinks_;
464  }
465}
466
467// Creates a branching object
468OsiBranchingObject * 
469OsiOldLink::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
470{
471  int j;
472  const double * solution = info->solution_;
473  double tolerance = info->primalTolerance_;
474  const double * upper = info->upper_;
475  int firstNonFixed=-1;
476  int lastNonFixed=-1;
477  int firstNonZero=-1;
478  int lastNonZero = -1;
479  double weight = 0.0;
480  double sum =0.0;
481  int base=0;
482  for (j=0;j<numberMembers_;j++) {
483    for (int k=0;k<numberLinks_;k++) {
484      int iColumn = members_[base+k];
485      if (upper[iColumn]) {
486        double value = CoinMax(0.0,solution[iColumn]);
487        sum += value;
488        if (firstNonFixed<0)
489          firstNonFixed=j;
490        lastNonFixed=j;
491        if (value>tolerance) {
492          weight += weights_[j]*value;
493          if (firstNonZero<0)
494            firstNonZero=j;
495          lastNonZero=j;
496        }
497      }
498    }
499    base += numberLinks_;
500  }
501  assert (lastNonZero-firstNonZero>=sosType_) ;
502  // find where to branch
503  assert (sum>0.0);
504  weight /= sum;
505  int iWhere;
506  double separator=0.0;
507  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
508    if (weight<weights_[iWhere+1])
509      break;
510  if (sosType_==1) {
511    // SOS 1
512    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
513  } else {
514    // SOS 2
515    if (iWhere==firstNonFixed)
516      iWhere++;;
517    if (iWhere==lastNonFixed-1)
518      iWhere = lastNonFixed-2;
519    separator = weights_[iWhere+1];
520  }
521  // create object
522  OsiBranchingObject * branch;
523  branch = new OsiOldLinkBranchingObject(solver,this,way,separator);
524  return branch;
525}
526OsiOldLinkBranchingObject::OsiOldLinkBranchingObject()
527  :OsiSOSBranchingObject()
528{
529}
530
531// Useful constructor
532OsiOldLinkBranchingObject::OsiOldLinkBranchingObject (OsiSolverInterface * solver,
533                                              const OsiOldLink * set,
534                                              int way ,
535                                              double separator)
536  :OsiSOSBranchingObject(solver,set,way,separator)
537{
538}
539
540// Copy constructor
541OsiOldLinkBranchingObject::OsiOldLinkBranchingObject ( const OsiOldLinkBranchingObject & rhs) :OsiSOSBranchingObject(rhs)
542{
543}
544
545// Assignment operator
546OsiOldLinkBranchingObject & 
547OsiOldLinkBranchingObject::operator=( const OsiOldLinkBranchingObject& rhs)
548{
549  if (this != &rhs) {
550    OsiSOSBranchingObject::operator=(rhs);
551  }
552  return *this;
553}
554OsiBranchingObject * 
555OsiOldLinkBranchingObject::clone() const
556{ 
557  return (new OsiOldLinkBranchingObject(*this));
558}
559
560
561// Destructor
562OsiOldLinkBranchingObject::~OsiOldLinkBranchingObject ()
563{
564}
565double
566OsiOldLinkBranchingObject::branch(OsiSolverInterface * solver)
567{
568  const OsiOldLink * set =
569    dynamic_cast <const OsiOldLink *>(originalObject_) ;
570  assert (set);
571  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
572  branchIndex_++;
573  int numberMembers = set->numberMembers();
574  const int * which = set->members();
575  const double * weights = set->weights();
576  int numberLinks = set->numberLinks();
577  //const double * lower = info->lower_;
578  //const double * upper = solver->getColUpper();
579  // *** for way - up means fix all those in down section
580  if (way<0) {
581    int i;
582    for ( i=0;i<numberMembers;i++) {
583      if (weights[i] > value_)
584        break;
585    }
586    assert (i<numberMembers);
587    int base=i*numberLinks;;
588    for (;i<numberMembers;i++) {
589      for (int k=0;k<numberLinks;k++) {
590        int iColumn = which[base+k];
591        solver->setColUpper(iColumn,0.0);
592      }
593      base += numberLinks;
594    }
595  } else {
596    int i;
597    int base=0;
598    for ( i=0;i<numberMembers;i++) { 
599      if (weights[i] >= value_) {
600        break;
601      } else {
602        for (int k=0;k<numberLinks;k++) {
603          int iColumn = which[base+k];
604          solver->setColUpper(iColumn,0.0);
605        }
606        base += numberLinks;
607      }
608    }
609    assert (i<numberMembers);
610  }
611  return 0.0;
612}
613// Print what would happen 
614void
615OsiOldLinkBranchingObject::print(const OsiSolverInterface * solver)
616{
617  const OsiOldLink * set =
618    dynamic_cast <const OsiOldLink *>(originalObject_) ;
619  assert (set);
620  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
621  int numberMembers = set->numberMembers();
622  int numberLinks = set->numberLinks();
623  const double * weights = set->weights();
624  const int * which = set->members();
625  const double * upper = solver->getColUpper();
626  int first=numberMembers;
627  int last=-1;
628  int numberFixed=0;
629  int numberOther=0;
630  int i;
631  int base=0;
632  for ( i=0;i<numberMembers;i++) {
633    for (int k=0;k<numberLinks;k++) {
634      int iColumn = which[base+k];
635      double bound = upper[iColumn];
636      if (bound) {
637        first = CoinMin(first,i);
638        last = CoinMax(last,i);
639      }
640    }
641    base += numberLinks;
642  }
643  // *** for way - up means fix all those in down section
644  base=0;
645  if (way<0) {
646    printf("SOS Down");
647    for ( i=0;i<numberMembers;i++) {
648      if (weights[i] > value_) 
649        break;
650      for (int k=0;k<numberLinks;k++) {
651        int iColumn = which[base+k];
652        double bound = upper[iColumn];
653        if (bound)
654          numberOther++;
655      }
656      base += numberLinks;
657    }
658    assert (i<numberMembers);
659    for (;i<numberMembers;i++) {
660      for (int k=0;k<numberLinks;k++) {
661        int iColumn = which[base+k];
662        double bound = upper[iColumn];
663        if (bound)
664          numberFixed++;
665      }
666      base += numberLinks;
667    }
668  } else {
669    printf("SOS Up");
670    for ( i=0;i<numberMembers;i++) {
671      if (weights[i] >= value_)
672        break;
673      for (int k=0;k<numberLinks;k++) {
674        int iColumn = which[base+k];
675        double bound = upper[iColumn];
676        if (bound)
677          numberFixed++;
678      }
679      base += numberLinks;
680    }
681    assert (i<numberMembers);
682    for (;i<numberMembers;i++) {
683      for (int k=0;k<numberLinks;k++) {
684        int iColumn = which[base+k];
685        double bound = upper[iColumn];
686        if (bound)
687          numberOther++;
688      }
689      base += numberLinks;
690    }
691  }
692  assert ((numberFixed%numberLinks)==0);
693  assert ((numberOther%numberLinks)==0);
694  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
695         value_,first,weights[first],last,weights[last],numberFixed/numberLinks,
696         numberOther/numberLinks);
697}
698// Default Constructor
699OsiBiLinear::OsiBiLinear ()
700  : OsiObject2(),
701    coefficient_(0.0),
702    xMeshSize_(0.0),
703    yMeshSize_(0.0),
704    xSatisfied_(1.0e-6),
705    ySatisfied_(1.0e-6),
706    xySatisfied_(1.0e-6),
707    xyBranchValue_(0.0),
708    xColumn_(-1),
709    yColumn_(-1),
710    firstLambda_(-1),
711    branchingStrategy_(0),
712    xRow_(-1),
713    yRow_(-1),
714    xyRow_(-1),
715    convexity_(-1),
716    chosen_(-1)
717{
718}
719
720// Useful constructor
721OsiBiLinear::OsiBiLinear (OsiSolverInterface * solver, int xColumn,
722                          int yColumn, int xyRow, double coefficient,
723                          double xMesh, double yMesh,
724                          int numberExistingObjects,const OsiObject ** objects )
725  : OsiObject2(),
726    coefficient_(coefficient),
727    xMeshSize_(xMesh),
728    yMeshSize_(yMesh),
729    xSatisfied_(1.0e-6),
730    ySatisfied_(1.0e-6),
731    xySatisfied_(1.0e-6),
732    xyBranchValue_(0.0),
733    xColumn_(xColumn),
734    yColumn_(yColumn),
735    firstLambda_(-1),
736    branchingStrategy_(0),
737    xRow_(-1),
738    yRow_(-1),
739    xyRow_(xyRow),
740    convexity_(-1),
741    chosen_(-1)
742{
743  double columnLower[4];
744  double columnUpper[4];
745  double objective[4];
746  double rowLower[3];
747  double rowUpper[3];
748  CoinBigIndex starts[5];
749  int index[16];
750  double element[16];
751  int i;
752  starts[0]=0;
753  // rows
754  int numberRows = solver->getNumRows();
755  // convexity
756  rowLower[0]=1.0;
757  rowUpper[0]=1.0;
758  convexity_ = numberRows;
759  starts[1]=0;
760  // x
761  rowLower[1]=0.0;
762  rowUpper[1]=0.0;
763  index[0]=xColumn_;
764  element[0]=-1.0;
765  xRow_ = numberRows+1;
766  starts[2]=1;
767  int nAdd=2;
768  if (xColumn!=yColumn) {
769    rowLower[2]=0.0;
770    rowUpper[2]=0.0;
771    index[1]=yColumn;
772    element[1]=-1.0;
773    nAdd=3;
774    yRow_ = numberRows+2;
775    starts[3]=2;
776  } else {
777    yRow_=-1;
778    branchingStrategy_=1;
779  }
780  solver->addRows(nAdd,starts,index,element,rowLower,rowUpper);
781  int n=0;
782  // order is LxLy, LxUy, UxLy and UxUy
783  firstLambda_ = solver->getNumCols();
784  // bit sloppy as theoretically could be infeasible but otherwise need to do more work
785  double xB[2];
786  double yB[2];
787  const double * lower = solver->getColLower();
788  const double * upper = solver->getColUpper();
789  xB[0]=lower[xColumn_];
790  xB[1]=upper[xColumn_];
791  // adjust
792  double distance;
793  int steps;
794  distance = xB[1]-xB[0];
795  steps = (int) ((distance+0.5*xMeshSize_)/xMeshSize_);
796  distance = xB[0]+xMeshSize_*steps;
797  if (fabs(xB[1]-distance)>1.0e-9) {
798    printf("bad x mesh %g %g %g -> %g\n",xB[0],xMeshSize_,xB[1],distance);
799    xB[1]=distance;
800    solver->setColUpper(xColumn_,distance);
801  }
802  yB[0]=lower[yColumn_];
803  yB[1]=upper[yColumn_];
804  distance = yB[1]-yB[0];
805  steps = (int) ((distance+0.5*yMeshSize_)/yMeshSize_);
806  distance = yB[0]+yMeshSize_*steps;
807  if (fabs(yB[1]-distance)>1.0e-9) {
808    printf("bad y mesh %g %g %g -> %g\n",yB[0],yMeshSize_,yB[1],distance);
809    yB[1]=distance;
810    solver->setColUpper(yColumn_,distance);
811  }
812  for (i=0;i<4;i++) {
813    double x = (i<2) ? xB[0] : xB[1];
814    double y = ((i&1)==0) ? yB[0] : yB[1];
815    columnLower[i]=0.0;
816    columnUpper[i]=COIN_DBL_MAX;
817    objective[i]=0.0;
818    double value;
819    // xy
820    value=coefficient_*x*y;
821    if (fabs(value)<1.0e-12)
822      value = 1.0e-12;
823    element[n]=value;
824    index[n++]=xyRow_;
825    // convexity
826    value=1.0;
827    element[n]=value;
828    index[n++]=0+numberRows;
829    // x
830    value=coefficient_*x;
831    if (fabs(value)<1.0e-12)
832      value = 1.0e-12;
833    element[n]=value;
834    index[n++]=1+numberRows;
835    if (xColumn_!=yColumn_) {
836      // y
837      value=coefficient_*y;
838      if (fabs(value)<1.0e-12)
839      value = 1.0e-12;
840      element[n]=value;
841      index[n++]=2+numberRows;
842    }
843    starts[i+1]=n;
844  }
845  solver->addCols(4,starts,index,element,columnLower,columnUpper,objective);
846  // At least one has to have a mesh
847  if (!xMeshSize_&&(!yMeshSize_||yRow_<0)) {
848    printf("one of x and y must have a mesh size\n");
849    abort();
850  } else if (yRow_>=0) {
851    if (!xMeshSize_)
852      branchingStrategy_ = 2;
853    else if (!yMeshSize_)
854      branchingStrategy_ = 1;
855  }
856  // Now add constraints to link in x and or y to existing ones.
857  bool xDone=false;
858  bool yDone=false;
859  // order is LxLy, LxUy, UxLy and UxUy
860  for (i=numberExistingObjects-1;i>=0;i--) {
861    const OsiObject * obj = objects[i];
862    const OsiBiLinear * obj2 =
863      dynamic_cast <const OsiBiLinear *>(obj) ;
864    if (obj2) {
865      if (xColumn_==obj2->xColumn_&&!xDone) {
866        // make sure y equal
867        double rhs=0.0;
868        CoinBigIndex starts[2];
869        int index[4];
870        double element[4]= {1.0,1.0,-1.0,-1.0};
871        starts[0]=0;
872        starts[1]=4;
873        index[0]=firstLambda_+0;
874        index[1]=firstLambda_+1;
875        index[2]=obj2->firstLambda_+0;
876        index[3]=obj2->firstLambda_+1;
877        solver->addRows(1,starts,index,element,&rhs,&rhs);
878        xDone=true;
879      }
880      if (yColumn_==obj2->yColumn_&&yRow_>=0&&!yDone) {
881        // make sure x equal
882        double rhs=0.0;
883        CoinBigIndex starts[2];
884        int index[4];
885        double element[4]= {1.0,1.0,-1.0,-1.0};
886        starts[0]=0;
887        starts[1]=4;
888        index[0]=firstLambda_+0;
889        index[1]=firstLambda_+2;
890        index[2]=obj2->firstLambda_+0;
891        index[3]=obj2->firstLambda_+2;
892        solver->addRows(1,starts,index,element,&rhs,&rhs);
893        yDone=true;
894      }
895    }
896  }
897}
898
899// Copy constructor
900OsiBiLinear::OsiBiLinear ( const OsiBiLinear & rhs)
901  :OsiObject2(rhs),
902   coefficient_(rhs.coefficient_),
903   xMeshSize_(rhs.xMeshSize_),
904   yMeshSize_(rhs.yMeshSize_),
905   xSatisfied_(rhs.xSatisfied_),
906   ySatisfied_(rhs.ySatisfied_),
907   xySatisfied_(rhs.xySatisfied_),
908   xyBranchValue_(rhs.xyBranchValue_),
909   xColumn_(rhs.xColumn_),
910   yColumn_(rhs.yColumn_),
911   firstLambda_(rhs.firstLambda_),
912   branchingStrategy_(rhs.branchingStrategy_),
913   xRow_(rhs.xRow_),
914   yRow_(rhs.yRow_),
915   xyRow_(rhs.xyRow_),
916   convexity_(rhs.convexity_),
917   chosen_(rhs.chosen_)
918{
919}
920
921// Clone
922OsiObject *
923OsiBiLinear::clone() const
924{
925  return new OsiBiLinear(*this);
926}
927
928// Assignment operator
929OsiBiLinear & 
930OsiBiLinear::operator=( const OsiBiLinear& rhs)
931{
932  if (this!=&rhs) {
933    OsiObject2::operator=(rhs);
934    coefficient_ = rhs.coefficient_;
935    xMeshSize_ = rhs.xMeshSize_;
936    yMeshSize_ = rhs.yMeshSize_;
937    xSatisfied_ = rhs.xSatisfied_;
938    ySatisfied_ = rhs.ySatisfied_;
939    xySatisfied_ = rhs.xySatisfied_;
940    xyBranchValue_ = rhs.xyBranchValue_;
941    xColumn_ = rhs.xColumn_;
942    yColumn_ = rhs.yColumn_;
943    firstLambda_ = rhs.firstLambda_;
944    branchingStrategy_ = rhs.branchingStrategy_;
945    xRow_ = rhs.xRow_;
946    yRow_ = rhs.yRow_;
947    xyRow_ = rhs.xyRow_;
948    convexity_ = rhs.convexity_;
949    chosen_ = rhs.chosen_;
950  }
951  return *this;
952}
953
954// Destructor
955OsiBiLinear::~OsiBiLinear ()
956{
957}
958
959// Infeasibility - large is 0.5
960double 
961OsiBiLinear::infeasibility(const OsiBranchingInformation * info,int & whichWay) const
962{
963  // order is LxLy, LxUy, UxLy and UxUy
964  double xB[2];
965  double yB[2];
966  xB[0]=info->lower_[xColumn_];
967  xB[1]=info->upper_[xColumn_];
968  yB[0]=info->lower_[yColumn_];
969  yB[1]=info->upper_[yColumn_];
970  double x = info->solution_[xColumn_];
971  x = CoinMax(x,xB[0]);
972  x = CoinMin(x,xB[1]);
973  double y = info->solution_[yColumn_];
974  y = CoinMax(y,yB[0]);
975  y = CoinMin(y,yB[1]);
976  int j;
977#ifndef NDEBUG
978  double xLambda = 0.0;
979  double yLambda = 0.0;
980  for (j=0;j<4;j++) {
981    int iX = j>>1;
982    int iY = j&1;
983    xLambda += xB[iX]*info->solution_[firstLambda_+j];
984    yLambda += yB[iY]*info->solution_[firstLambda_+j];
985  }
986  assert (fabs(x-xLambda)<1.0e-4);
987  assert (fabs(y-yLambda)<1.0e-4);
988#endif
989  // If x or y not satisfied then branch on that
990  double distance;
991  int steps;
992  bool xSatisfied;
993  double xNew;
994  distance = x-xB[0];
995  if (xMeshSize_) {
996    steps = (int) ((distance+0.5*xMeshSize_)/xMeshSize_);
997    xNew = xB[0]+steps*xMeshSize_;
998    assert (xNew<=xB[1]+1.0e-5);
999    xSatisfied =  (fabs(xNew-x)<xSatisfied_);
1000  } else {
1001    xSatisfied=true;
1002  }
1003  bool ySatisfied;
1004  double yNew;
1005  distance = y-yB[0];
1006  if (yMeshSize_) {
1007    steps = (int) ((distance+0.5*yMeshSize_)/yMeshSize_);
1008    yNew = yB[0]+steps*yMeshSize_;
1009    assert (yNew<=yB[1]+1.0e-5);
1010    ySatisfied =  (fabs(yNew-y)<ySatisfied_)||!y;
1011  } else {
1012    ySatisfied=true;
1013  }
1014  /* There are several possibilities
1015     1 - one or both are unsatisfied and branching strategy tells us what to do
1016     2 - both are unsatisfied and branching strategy is 0
1017     3 - both are satisfied but xy is not
1018         3a one has bounds within satisfied_ - other does not
1019         (or neither have but branching strategy tells us what to do)
1020         3b neither do - and branching strategy does not tell us
1021         3c both do - treat as feasible knowing another copy of object will fix
1022     4 - both are satisfied and xy is satisfied - as 3c
1023  */
1024  chosen_=-1;
1025  xyBranchValue_=COIN_DBL_MAX;
1026  whichWay_=0;
1027  if ( !xSatisfied) {
1028    if (!ySatisfied) {
1029      if (branchingStrategy_==0) {
1030        // If pseudo shadow prices then see what would happen
1031        if (info->defaultDual_>=0.0) {
1032          // need coding here
1033          if (fabs(x-xNew)>fabs(y-yNew)) {
1034            chosen_=0;
1035            xyBranchValue_=x;
1036          } else {
1037            chosen_=1;
1038            xyBranchValue_=y;
1039          }
1040        } else {
1041          if (fabs(x-xNew)>fabs(y-yNew)) {
1042            chosen_=0;
1043            xyBranchValue_=x;
1044          } else {
1045            chosen_=1;
1046            xyBranchValue_=y;
1047          }
1048        }
1049      } else if (branchingStrategy_==1) {
1050        chosen_=0;
1051        xyBranchValue_=x;
1052      } else {
1053        chosen_=1;
1054        xyBranchValue_=y;
1055      }
1056    } else {
1057      // y satisfied
1058      chosen_=0;
1059      xyBranchValue_=x;
1060    }
1061  } else {
1062    // x satisfied
1063    if (!ySatisfied) {
1064      chosen_=1;
1065      xyBranchValue_=y;
1066    } else {
1067      /*
1068        3 - both are satisfied but xy is not
1069         3a one has bounds within satisfied_ - other does not
1070         (or neither have but branching strategy tells us what to do)
1071         3b neither do - and branching strategy does not tell us
1072         3c both do - treat as feasible knowing another copy of object will fix
1073        4 - both are satisfied and xy is satisfied - as 3c
1074      */
1075      double xyTrue = x*y;
1076      double xyLambda = 0.0;
1077      for (j=0;j<4;j++) {
1078        int iX = j>>1;
1079        int iY = j&1;
1080        xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j];
1081      }
1082      if (fabs(xyLambda-xyTrue)<xySatisfied_) {
1083        // satisfied
1084      } else {
1085        if (xB[1]-xB[0]>=xSatisfied_&&xMeshSize_) {
1086          if (yB[1]-yB[0]>=ySatisfied_&&yMeshSize_) {
1087            if (branchingStrategy_==0) {
1088              // If pseudo shadow prices then see what would happen
1089              if (info->defaultDual_>=0.0) {
1090                // need coding here
1091                if (xB[1]-xB[0]>yB[1]-yB[0]) {
1092                  chosen_=0;
1093                  xyBranchValue_=0.5*(xB[0]+xB[1]);
1094                } else {
1095                  chosen_=1;
1096                  xyBranchValue_=0.5*(yB[0]+yB[1]);
1097                }
1098              } else {
1099                if (xB[1]-xB[0]>yB[1]-yB[0]) {
1100                  chosen_=0;
1101                  xyBranchValue_=0.5*(xB[0]+xB[1]);
1102                } else {
1103                  chosen_=1;
1104                  xyBranchValue_=0.5*(yB[0]+yB[1]);
1105                }
1106              }
1107            } else if (branchingStrategy_==1) {
1108              chosen_=0;
1109              xyBranchValue_=0.5*(xB[0]+xB[1]);
1110            } else {
1111              chosen_=1;
1112              xyBranchValue_=0.5*(yB[0]+yB[1]);
1113            }
1114          } else {
1115            // y satisfied
1116            chosen_=0;
1117            xyBranchValue_=0.5*(xB[0]+xB[1]);
1118          }
1119        } else if (yB[1]-yB[0]>=ySatisfied_&&yMeshSize_) {
1120          chosen_=1;
1121          xyBranchValue_=0.5*(yB[0]+yB[1]);
1122        } else {
1123          // treat as satisfied
1124        }
1125      }
1126    }
1127  }
1128  if (chosen_==-1) {
1129    infeasibility_=0.0;
1130  } else if (chosen_==0) {
1131    infeasibility_ = CoinMax(fabs(xyBranchValue_-x),1.0e-12);
1132    assert (xyBranchValue_>=info->lower_[xColumn_]&&xyBranchValue_<=info->upper_[xColumn_]);
1133  } else {
1134    infeasibility_ = CoinMax(fabs(xyBranchValue_-y),1.0e-12);
1135    assert (xyBranchValue_>=info->lower_[yColumn_]&&xyBranchValue_<=info->upper_[yColumn_]);
1136  }
1137  if (info->defaultDual_<0.0) {
1138    // not using pseudo shadow prices
1139    otherInfeasibility_ = 1.0-infeasibility_;
1140  } else {
1141    abort();
1142  }
1143  whichWay=whichWay_;
1144  return infeasibility_;
1145}
1146
1147// This looks at solution and sets bounds to contain solution
1148double
1149OsiBiLinear::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
1150{
1151  // order is LxLy, LxUy, UxLy and UxUy
1152  double xB[2];
1153  double yB[2];
1154  xB[0]=info->lower_[xColumn_];
1155  xB[1]=info->upper_[xColumn_];
1156  yB[0]=info->lower_[yColumn_];
1157  yB[1]=info->upper_[yColumn_];
1158  double x = info->solution_[xColumn_];
1159  double y = info->solution_[yColumn_];
1160  int j;
1161#ifndef NDEBUG
1162  double xLambda = 0.0;
1163  double yLambda = 0.0;
1164  for (j=0;j<4;j++) {
1165    int iX = j>>1;
1166    int iY = j&1;
1167    xLambda += xB[iX]*info->solution_[firstLambda_+j];
1168    yLambda += yB[iY]*info->solution_[firstLambda_+j];
1169  }
1170  if (fabs(x-xLambda)>1.0e-4||
1171      fabs(y-yLambda)>1.0e-4)
1172    printf("feasibleregion x %d given %g lambda %g y %d given %g lambda %g\n",
1173           xColumn_,x,xLambda,
1174           yColumn_,y,yLambda);
1175#endif
1176  double infeasibility=0.0;
1177  double distance;
1178  int steps;
1179  double xNew=x;
1180  distance = x-xB[0];
1181  if (xMeshSize_) {
1182    steps = (int) ((distance+0.5*xMeshSize_)/xMeshSize_);
1183    xNew = xB[0]+steps*xMeshSize_;
1184    assert (xNew<=xB[1]+1.0e-5);
1185    infeasibility +=  fabs(xNew-x);
1186    solver->setColLower(xColumn_,xNew);
1187    solver->setColUpper(xColumn_,xNew);
1188  }
1189  double yNew=y;
1190  distance = y-yB[0];
1191  if (yMeshSize_) {
1192    steps = (int) ((distance+0.5*yMeshSize_)/yMeshSize_);
1193    yNew = yB[0]+steps*yMeshSize_;
1194    assert (yNew<=yB[1]+1.0e-5);
1195    infeasibility +=  fabs(yNew-y);
1196    solver->setColLower(yColumn_,yNew);
1197    solver->setColUpper(yColumn_,yNew);
1198  }
1199  double xyTrue = xNew*yNew;
1200  double xyLambda = 0.0;
1201  for (j=0;j<4;j++) {
1202    int iX = j>>1;
1203    int iY = j&1;
1204    xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j];
1205  }
1206  infeasibility += fabs(xyTrue-xyLambda);
1207  return infeasibility;
1208}
1209
1210// Redoes data when sequence numbers change
1211void 
1212OsiBiLinear::resetSequenceEtc(int numberColumns, const int * originalColumns)
1213{
1214  int i;
1215  for (i=0;i<numberColumns;i++) {
1216    if (originalColumns[i]==firstLambda_)
1217      break;
1218  }
1219  if (i<numberColumns) {
1220    firstLambda_ = i;
1221    for (int j=0;j<4;j++) {
1222      assert (originalColumns[j+i]-firstLambda_==j);
1223    }
1224  } else {
1225    printf("lost set\n");
1226    abort();
1227  }
1228  // rows will be out anyway
1229  abort();
1230}
1231
1232// Creates a branching object
1233OsiBranchingObject * 
1234OsiBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
1235{
1236  // create object
1237  OsiBranchingObject * branch;
1238  assert (chosen_==0||chosen_==1);
1239  if (chosen_==0) 
1240    assert (xyBranchValue_>=info->lower_[xColumn_]&&xyBranchValue_<=info->upper_[xColumn_]);
1241  else
1242    assert (xyBranchValue_>=info->lower_[yColumn_]&&xyBranchValue_<=info->upper_[yColumn_]);
1243  branch = new OsiBiLinearBranchingObject(solver,this,way,xyBranchValue_,chosen_);
1244  return branch;
1245}
1246// Does work of branching
1247void 
1248OsiBiLinear::newBounds(OsiSolverInterface * solver, int way, short xOrY, double separator) const
1249{
1250  int iColumn;
1251  double mesh;
1252  if (xOrY==0) {
1253    iColumn=xColumn_;
1254    mesh=xMeshSize_;
1255  } else {
1256    iColumn=yColumn_;
1257    mesh=yMeshSize_;
1258  }
1259  double lower = solver->getColLower()[iColumn];
1260  double distance;
1261  int steps;
1262  double zNew=separator;
1263  distance = separator-lower;
1264  assert (mesh);
1265  steps = (int) ((distance+0.5*mesh)/mesh);
1266  zNew = lower+steps*mesh;
1267  assert (zNew<=solver->getColUpper()[iColumn]+1.0e-5);
1268#ifndef NDEBUG
1269    double oldUpper = solver->getColUpper()[iColumn] ;
1270    double oldLower = solver->getColLower()[iColumn] ;
1271#endif
1272  if (way<0) {
1273    if (zNew>separator)
1274      zNew -= mesh;
1275#ifndef NDEBUG
1276    double oldUpper = solver->getColUpper()[iColumn] ;
1277    assert (oldUpper>zNew-1.0e-8);
1278    if (oldUpper<zNew+1.0e-8)
1279      printf("null change on columnUpper %d - bounds %g,%g\n",iColumn,oldLower,oldUpper);
1280#endif
1281    solver->setColUpper(iColumn,zNew);
1282  } else {
1283    if (zNew<separator)
1284      zNew += mesh;
1285#ifndef NDEBUG
1286    double oldLower = solver->getColLower()[iColumn] ;
1287    assert (oldLower<zNew+1.0e-8);
1288    if (oldLower>zNew-1.0e-8)
1289      printf("null change on columnLower %d - bounds %g,%g\n",iColumn,oldLower,oldUpper);
1290#endif
1291    solver->setColLower(iColumn,zNew);
1292  }
1293#if 0
1294  // always free up lambda
1295  for (int i=firstLambda_;i<firstLambda_+4;i++) {
1296    solver->setColLower(i,0.0);
1297    solver->setColUpper(i,COIN_DBL_MAX);
1298  }
1299#endif
1300}
1301// Updates coefficients
1302void 
1303OsiBiLinear::updateCoefficients(const double * lower, const double * upper,
1304                                CoinPackedMatrix * matrix, CoinWarmStartBasis * basis) const
1305{
1306  double * element = matrix->getMutableElements();
1307  const int * row = matrix->getIndices();
1308  const CoinBigIndex * columnStart = matrix->getVectorStarts();
1309  //const int * columnLength = matrix->getVectorLengths();
1310  // order is LxLy, LxUy, UxLy and UxUy
1311  double xB[2];
1312  double yB[2];
1313  xB[0]=lower[xColumn_];
1314  xB[1]=upper[xColumn_];
1315  yB[0]=lower[yColumn_];
1316  yB[1]=upper[yColumn_];
1317  //printf("x %d (%g,%g) y %d (%g,%g)\n",
1318  // xColumn_,xB[0],xB[1],
1319  // yColumn_,yB[0],yB[1]);
1320  CoinWarmStartBasis::Status status[4];
1321  int numStruct = basis->getNumStructural()-firstLambda_;
1322  for (int j=0;j<4;j++) {
1323    status[j]=(j<numStruct) ? basis->getStructStatus(j+firstLambda_) : CoinWarmStartBasis::atLowerBound;
1324    int iX = j>>1;
1325    double x = xB[iX];
1326    int iY = j&1;
1327    double y = yB[iY];
1328    CoinBigIndex k = columnStart[j+firstLambda_];
1329    double value;
1330    // xy
1331    value=coefficient_*x*y;
1332    assert (row[k]==xyRow_);
1333#if BI_PRINT > 1
1334    printf("j %d xy (%d,%d) coeff from %g to %g\n",j,xColumn_,yColumn_,element[k],value);
1335#endif
1336    element[k++]=value;
1337    // convexity
1338    assert (row[k]==convexity_);
1339    k++;
1340    // x
1341    value=coefficient_*x;
1342#if BI_PRINT > 1
1343    printf("j %d x (%d) coeff from %g to %g\n",j,xColumn_,element[k],value);
1344#endif
1345    assert (row[k]==xRow_);
1346    element[k++]=value;
1347    if (yRow_>=0) {
1348      // y
1349      value=coefficient_*y;
1350#if BI_PRINT > 1
1351      printf("j %d y (%d) coeff from %g to %g\n",j,yColumn_,element[k],value);
1352#endif
1353      assert (row[k]==yRow_);
1354      element[k++]=value;
1355    }
1356  }
1357 
1358  if (xB[0]==xB[1]) {
1359    if (yB[0]==yB[1]) {
1360      // only one basic
1361      bool first=true;
1362      for (int j=0;j<4;j++) {
1363        if (status[j]==CoinWarmStartBasis::basic) {
1364          if (first) {
1365            first=false;
1366          } else {
1367            basis->setStructStatus(j+firstLambda_,CoinWarmStartBasis::atLowerBound);
1368#if BI_PRINT
1369            printf("zapping %d (x=%d,y=%d)\n",j,xColumn_,yColumn_);
1370#endif
1371          }
1372        }
1373      }
1374    } else {
1375      if (status[0]==CoinWarmStartBasis::basic&&
1376          status[2]==CoinWarmStartBasis::basic) {
1377        basis->setStructStatus(2+firstLambda_,CoinWarmStartBasis::atLowerBound);
1378#if BI_PRINT
1379        printf("zapping %d (x=%d,y=%d)\n",2,xColumn_,yColumn_);
1380#endif
1381      }
1382      if (status[1]==CoinWarmStartBasis::basic&&
1383          status[3]==CoinWarmStartBasis::basic) {
1384        basis->setStructStatus(3+firstLambda_,CoinWarmStartBasis::atLowerBound);
1385#if BI_PRINT
1386        printf("zapping %d (x=%d,y=%d)\n",3,xColumn_,yColumn_);
1387#endif
1388      }
1389    }
1390  } else if (yB[0]==yB[1]) {
1391    if (status[0]==CoinWarmStartBasis::basic&&
1392        status[1]==CoinWarmStartBasis::basic) {
1393      basis->setStructStatus(1+firstLambda_,CoinWarmStartBasis::atLowerBound);
1394#if BI_PRINT
1395      printf("zapping %d (x=%d,y=%d)\n",1,xColumn_,yColumn_);
1396#endif
1397    }
1398    if (status[2]==CoinWarmStartBasis::basic&&
1399        status[3]==CoinWarmStartBasis::basic) {
1400      basis->setStructStatus(3+firstLambda_,CoinWarmStartBasis::atLowerBound);
1401#if BI_PRINT
1402      printf("zapping %d (x=%d,y=%d)\n",3,xColumn_,yColumn_);
1403#endif
1404    }
1405  }
1406}
1407// This does NOT set mutable stuff
1408double 
1409OsiBiLinear::checkInfeasibility(const OsiBranchingInformation * info) const
1410{
1411  int way;
1412  double saveInfeasibility = infeasibility_;
1413  int saveWhichWay = whichWay_;
1414  double saveXyBranchValue = xyBranchValue_;
1415  short saveChosen = chosen_;
1416  double value = infeasibility(info,way);
1417  infeasibility_ = saveInfeasibility;
1418  whichWay_ = saveWhichWay;
1419  xyBranchValue_ = saveXyBranchValue;
1420  chosen_ = saveChosen;
1421  return value;
1422}
1423OsiBiLinearBranchingObject::OsiBiLinearBranchingObject()
1424  :OsiTwoWayBranchingObject(),
1425   chosen_(0)
1426{
1427}
1428
1429// Useful constructor
1430OsiBiLinearBranchingObject::OsiBiLinearBranchingObject (OsiSolverInterface * solver,
1431                                                        const OsiBiLinear * set,
1432                                                        int way ,
1433                                                        double separator,
1434                                                        int chosen)
1435  :OsiTwoWayBranchingObject(solver,set,way,separator),
1436   chosen_(chosen)
1437{
1438  assert (chosen_>=0&&chosen_<2);
1439}
1440
1441// Copy constructor
1442OsiBiLinearBranchingObject::OsiBiLinearBranchingObject ( const OsiBiLinearBranchingObject & rhs) 
1443  :OsiTwoWayBranchingObject(rhs),
1444   chosen_(rhs.chosen_)
1445{
1446}
1447
1448// Assignment operator
1449OsiBiLinearBranchingObject & 
1450OsiBiLinearBranchingObject::operator=( const OsiBiLinearBranchingObject& rhs)
1451{
1452  if (this != &rhs) {
1453    OsiTwoWayBranchingObject::operator=(rhs);
1454    chosen_ = rhs.chosen_;
1455  }
1456  return *this;
1457}
1458OsiBranchingObject * 
1459OsiBiLinearBranchingObject::clone() const
1460{ 
1461  return (new OsiBiLinearBranchingObject(*this));
1462}
1463
1464
1465// Destructor
1466OsiBiLinearBranchingObject::~OsiBiLinearBranchingObject ()
1467{
1468}
1469double
1470OsiBiLinearBranchingObject::branch(OsiSolverInterface * solver)
1471{
1472  const OsiBiLinear * set =
1473    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
1474  assert (set);
1475  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
1476  branchIndex_++;
1477  set->newBounds(solver, way, chosen_, value_);
1478  return 0.0;
1479}
1480// Print what would happen 
1481void
1482OsiBiLinearBranchingObject::print(const OsiSolverInterface * solver)
1483{
1484  const OsiBiLinear * set =
1485    dynamic_cast <const OsiBiLinear *>(originalObject_) ;
1486  assert (set);
1487  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
1488  int iColumn = (chosen_==1) ? set->xColumn() : set->yColumn();
1489  printf("OsiBiLinear would branch %s on %c variable %d from value %g\n",
1490         (way<0) ? "down" : "up",
1491         (chosen_==0) ? 'X' : 'Y', iColumn, value_);
1492}
1493/** Default Constructor
1494
1495  Equivalent to an unspecified binary variable.
1496*/
1497OsiSimpleFixedInteger::OsiSimpleFixedInteger ()
1498  : OsiSimpleInteger()
1499{
1500}
1501
1502/** Useful constructor
1503
1504  Loads actual upper & lower bounds for the specified variable.
1505*/
1506OsiSimpleFixedInteger::OsiSimpleFixedInteger (const OsiSolverInterface * solver, int iColumn)
1507  : OsiSimpleInteger(solver,iColumn)
1508{
1509}
1510
1511 
1512// Useful constructor - passed solver index and original bounds
1513OsiSimpleFixedInteger::OsiSimpleFixedInteger ( int iColumn, double lower, double upper)
1514  : OsiSimpleInteger(iColumn,lower,upper)
1515{
1516}
1517
1518// Useful constructor - passed simple integer
1519OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleInteger &rhs)
1520  : OsiSimpleInteger(rhs)
1521{
1522}
1523
1524// Copy constructor
1525OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleFixedInteger & rhs)
1526  :OsiSimpleInteger(rhs)
1527
1528{
1529}
1530
1531// Clone
1532OsiObject *
1533OsiSimpleFixedInteger::clone() const
1534{
1535  return new OsiSimpleFixedInteger(*this);
1536}
1537
1538// Assignment operator
1539OsiSimpleFixedInteger & 
1540OsiSimpleFixedInteger::operator=( const OsiSimpleFixedInteger& rhs)
1541{
1542  if (this!=&rhs) {
1543    OsiSimpleInteger::operator=(rhs);
1544  }
1545  return *this;
1546}
1547
1548// Destructor
1549OsiSimpleFixedInteger::~OsiSimpleFixedInteger ()
1550{
1551}
1552// Infeasibility - large is 0.5
1553double 
1554OsiSimpleFixedInteger::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
1555{
1556  double value = info->solution_[columnNumber_];
1557  value = CoinMax(value, info->lower_[columnNumber_]);
1558  value = CoinMin(value, info->upper_[columnNumber_]);
1559  double nearest = floor(value+(1.0-0.5));
1560  if (nearest>value) { 
1561    whichWay=1;
1562  } else {
1563    whichWay=0;
1564  }
1565  infeasibility_ = fabs(value-nearest);
1566  bool satisfied=false;
1567  if (infeasibility_<=info->integerTolerance_) {
1568    otherInfeasibility_ = 1.0;
1569    satisfied=true;
1570    if (info->lower_[columnNumber_]!=info->upper_[columnNumber_])
1571      infeasibility_ = 1.0e-5;
1572    else
1573      infeasibility_ = 0.0;
1574  } else if (info->defaultDual_<0.0) {
1575    otherInfeasibility_ = 1.0-infeasibility_;
1576  } else {
1577    const double * pi = info->pi_;
1578    const double * activity = info->rowActivity_;
1579    const double * lower = info->rowLower_;
1580    const double * upper = info->rowUpper_;
1581    const double * element = info->elementByColumn_;
1582    const int * row = info->row_;
1583    const CoinBigIndex * columnStart = info->columnStart_;
1584    const int * columnLength = info->columnLength_;
1585    double direction = info->direction_;
1586    double downMovement = value - floor(value);
1587    double upMovement = 1.0-downMovement;
1588    double valueP = info->objective_[columnNumber_]*direction;
1589    CoinBigIndex start = columnStart[columnNumber_];
1590    CoinBigIndex end = start + columnLength[columnNumber_];
1591    double upEstimate = 0.0;
1592    double downEstimate = 0.0;
1593    if (valueP>0.0)
1594      upEstimate = valueP*upMovement;
1595    else
1596      downEstimate -= valueP*downMovement;
1597    double tolerance = info->primalTolerance_;
1598    for (CoinBigIndex j=start;j<end;j++) {
1599      int iRow = row[j];
1600      if (lower[iRow]<-1.0e20) 
1601        assert (pi[iRow]<=1.0e-4);
1602      if (upper[iRow]>1.0e20) 
1603        assert (pi[iRow]>=-1.0e-4);
1604      valueP = pi[iRow]*direction;
1605      double el2 = element[j];
1606      double value2 = valueP*el2;
1607      double u=0.0;
1608      double d=0.0;
1609      if (value2>0.0)
1610        u = value2;
1611      else
1612        d = -value2;
1613      // if up makes infeasible then make at least default
1614      double newUp = activity[iRow] + upMovement*el2;
1615      if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance)
1616        u = CoinMax(u,info->defaultDual_);
1617      upEstimate += u*upMovement;
1618      // if down makes infeasible then make at least default
1619      double newDown = activity[iRow] - downMovement*el2;
1620      if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance)
1621        d = CoinMax(d,info->defaultDual_);
1622      downEstimate += d*downMovement;
1623    }
1624    if (downEstimate>=upEstimate) {
1625      infeasibility_ = CoinMax(1.0e-12,upEstimate);
1626      otherInfeasibility_ = CoinMax(1.0e-12,downEstimate);
1627      whichWay = 1;
1628    } else {
1629      infeasibility_ = CoinMax(1.0e-12,downEstimate);
1630      otherInfeasibility_ = CoinMax(1.0e-12,upEstimate);
1631      whichWay = 0;
1632    }
1633  }
1634  if (preferredWay_>=0&&!satisfied)
1635    whichWay = preferredWay_;
1636  whichWay_=whichWay;
1637  return infeasibility_;
1638}
1639// Creates a branching object
1640OsiBranchingObject * 
1641OsiSimpleFixedInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const 
1642{
1643  double value = info->solution_[columnNumber_];
1644  value = CoinMax(value, info->lower_[columnNumber_]);
1645  value = CoinMin(value, info->upper_[columnNumber_]);
1646  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
1647  double nearest = floor(value+0.5);
1648  double integerTolerance = info->integerTolerance_;
1649  if (fabs(value-nearest)<integerTolerance) {
1650    // adjust value
1651    if (nearest!=info->upper_[columnNumber_])
1652      value = nearest+2.0*integerTolerance;
1653    else
1654      value = nearest-2.0*integerTolerance;
1655  }
1656  OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way,
1657                                             value);
1658  return branch;
1659}
1660
Note: See TracBrowser for help on using the repository browser.