source: trunk/Cbc/src/CbcSymmetry.cpp @ 2092

Last change on this file since 2092 was 2092, checked in by forrest, 4 years ago

symmetry and diving

File size: 37.8 KB
Line 
1/* $Id: CbcSymmetry.cpp 925 2012-11-27 19:11:04Z stefan $
2 *
3 * hacked from CouenneSymmetry.cpp
4 * Name:    CbcSymmetry.cpp
5 * Author:  Jim Ostrowski (the good bits - rest JJHF)
6 * Purpose: methods for exploiting symmetry
7 * Date:    October 13, 2010
8 *
9 * This file is licensed under the Eclipse Public License (EPL)
10 */
11//#define PRINT_MORE 1
12#include <stdio.h>
13
14#ifdef COIN_HAS_NTY
15
16#include <cassert>
17#include <vector>
18#include <algorithm>
19#include <ostream>
20#include <iterator>
21
22#include "CbcSymmetry.hpp"
23#include "CbcBranchingObject.hpp"
24#include "CoinTime.hpp"
25#define NAUTY_MAX_LEVEL 2000
26#if NAUTY_MAX_LEVEL
27extern int nauty_maxalllevel;
28#endif
29/* Deliberately not threadsafe to save effort
30   Just for statistics
31   and not worth gathering across threads
32   can redo later
33 */
34static int nautyBranchCalls_ = 0;
35static int lastNautyBranchCalls_ = 0;
36static int nautyBranchSucceeded_ = 0;
37static int nautyFixCalls_ = 0;
38static int lastNautyFixCalls_ = 0;
39static int nautyFixSucceeded_ = 0;
40static double nautyTime_ = 0.0;
41static double nautyFixes_= 0.0; 
42static double nautyOtherBranches_ = 0.0;
43
44void Node::node(int i, double c , double l, double u, int cod, int s){
45  index = i;
46  coeff = c;
47  lb = l;
48  ub = u;
49  color = -1;
50  code = cod;
51  sign = s;
52}
53
54inline bool CbcSymmetry::compare (register Node &a, register Node &b) const {
55  if(a.get_code() == b.get_code() )
56    if(a.get_coeff() == b.get_coeff() )
57      if(a.get_sign() == b.get_sign() )
58        if( fabs ( a.get_lb() - b.get_lb() ) <= COUENNE_HACKED_EPS )
59          if( fabs ( a.get_ub() - b.get_ub() ) <= COUENNE_HACKED_EPS )
60            return 1; 
61  return 0;
62}
63
64void CbcSymmetry::Compute_Symmetry() const{
65
66
67  std::sort(node_info_. begin (), node_info_. end (), node_sort);
68
69  for (std::vector <Node>:: iterator i = node_info_. begin (); i != node_info_. end (); ++i) 
70    (*i).color_vertex(-1);
71 
72  int color = 1;
73  for (std::vector <Node>:: iterator i = node_info_. begin (); i != node_info_. end (); ++i) {
74    if( (*i).get_color() == -1){
75      (*i).color_vertex(color);
76#ifdef PRINT_MORE
77      printf ("Graph vertex %d is given color %d\n", (*i).get_index(), color);
78#endif
79      nauty_info_ -> color_node((*i).get_index(), color);
80      for (std::vector <Node>:: iterator j = i+1; j != node_info_. end (); ++j)
81        if( compare( (*i) , (*j) ) ==1){
82          (*j).color_vertex(color);
83          nauty_info_ -> color_node((*j).get_index(),color);
84#ifdef PRINT_MORE
85          printf ("Graph vertex %d is given color %d, the same as vertex %d\n", (*j).get_index(), color, (*i).get_index());
86#endif
87        }
88      //       else
89      // j = node_info_. end();
90      color++;
91    }
92  }
93
94  //Print_Orbits ();
95  nauty_info_ -> computeAuto();
96  //Print_Orbits ();
97}
98
99int
100CbcSymmetry::statsOrbits(CbcModel * model, int type) const
101{
102  char general[200];
103  int returnCode=0;
104  bool printSomething=true;
105  if (type) {
106    double branchSuccess=0.0;
107    if (nautyBranchSucceeded_) 
108      branchSuccess = nautyOtherBranches_/nautyBranchSucceeded_;
109    double fixSuccess=0.0;
110    if (nautyFixSucceeded_) 
111      fixSuccess = nautyFixes_/nautyFixSucceeded_;
112    if (nautyBranchCalls_>lastNautyBranchCalls_||
113        nautyFixCalls_>lastNautyFixCalls_) {
114      sprintf(general,"Orbital branching tried %d times, succeeded %d times - average extra %7.3f, fixing %d times (%d, %7.3f)",
115              nautyBranchCalls_,nautyBranchSucceeded_,branchSuccess,
116              nautyFixCalls_,nautyFixSucceeded_,fixSuccess);
117      lastNautyBranchCalls_=nautyBranchCalls_;
118      lastNautyFixCalls_=nautyFixCalls_;
119    } else {
120      printSomething=false;
121    }
122  } else {
123    returnCode = nauty_info_->getNumGenerators();
124    if (!nauty_info_->errorStatus()) {
125      if (returnCode && numberUsefulOrbits_) {
126        sprintf (general,"Nauty: %d orbits (%d useful covering %d variables), %d generators, group size: %g - dense size %d, sparse %d - took %g seconds",
127                 nauty_info_->getNumOrbits(),numberUsefulOrbits_,numberUsefulObjects_,
128                 nauty_info_ -> getNumGenerators () ,
129                 nauty_info_ -> getGroupSize (),
130                 whichOrbit_[0],whichOrbit_[1],nautyTime_);
131      } else {
132        if ((model->moreSpecialOptions2()&(128|256))!=(128|256)) 
133          sprintf(general,"Nauty did not find any useful orbits in time %g",nautyTime_);
134        else
135          sprintf(general,"Nauty did not find any useful orbits - but keeping Nauty on");
136      }
137    } else {
138      // error
139      sprintf(general,"Nauty failed with error code %d (%g seconds)",
140              nauty_info_->errorStatus(),nautyTime_);
141      model->setMoreSpecialOptions2(model->moreSpecialOptions2()&~(128|256));
142    }
143  }
144  if (printSomething)
145    model->messageHandler()->message(CBC_GENERAL,
146                                     model->messages())
147      << general << CoinMessageEol ;
148  return returnCode;
149}
150 
151void CbcSymmetry::Print_Orbits () const {
152
153  //printf ("num gens = %d, num orbits = %d \n", nauty_info_ -> getNumGenerators(), nauty_info_ -> getNumOrbits() );
154
155  std::vector<std::vector<int> > *new_orbits = nauty_info_->getOrbits();
156
157  printf ("Nauty: %d generators, group size: %.0g",
158          //  nauty_info_->getNumOrbits(),
159          nauty_info_ -> getNumGenerators () ,
160          nauty_info_ -> getGroupSize ());
161
162  int nNonTrivialOrbits = 0;
163
164  for (unsigned int i = 0; i < new_orbits -> size(); i++) {
165    if ((*new_orbits)[i].size() > 1) 
166      nNonTrivialOrbits++;
167    else continue;
168
169    // int orbsize = (*new_orbits)[i].size();
170    // printf( "Orbit %d [size: %d] [", i, orbsize);
171    // copy ((*new_orbits)[i].begin(), (*new_orbits)[i].end(),
172    //    std::ostream_iterator<int>(std::cout, " "));
173    // printf("] \n");
174  }
175
176  printf (" (%d non-trivial orbits).\n", nNonTrivialOrbits);
177
178#if 1
179  if (nNonTrivialOrbits) {
180
181    int orbCnt = 0;
182
183    std::vector<std::vector<int> > *orbits = nauty_info_ -> getOrbits ();
184
185    for   (std::vector <std::vector<int> >::iterator i = orbits -> begin ();  i != orbits -> end (); ++i) {
186
187      printf ("Orbit %d: ", orbCnt++);
188
189      for (std::vector<int>::iterator j = i -> begin (); j != i -> end (); ++j)
190        printf (" %d", *j);
191
192      printf ("\n");
193    }
194  }
195#endif
196
197
198#if 0
199  if (nNonTrivialOrbits)
200    for (int i=0; i< nVars (); i++) {
201
202      std::vector< int > *branch_orbit = Find_Orbit (i);
203
204      if (branch_orbit -> size () > 1) {
205        printf ("x%04d: ", i);
206
207        for (std::vector<int>::iterator it = branch_orbit -> begin (); it != branch_orbit -> end (); ++it)
208          printf ("%d ", *it);
209        printf ("\n");
210      }
211    }
212#endif
213  delete new_orbits;
214}
215void 
216CbcSymmetry::fillOrbits()
217{
218  for (int i=0;i<numberColumns_;i++)
219    whichOrbit_[i]=-1;
220  numberUsefulOrbits_=0;
221  numberUsefulObjects_=0;
222
223  std::vector<std::vector<int> > *orbits = nauty_info_ -> getOrbits ();
224
225  for   (std::vector <std::vector<int> >::iterator i = orbits -> begin ();  i != orbits -> end (); ++i) {
226    int nUseful=0;
227    int jColumn=-2;
228    for (std::vector<int>::iterator j = i -> begin (); j != i -> end (); ++j) {
229      int iColumn=*j;
230      if (iColumn<numberColumns_) {
231        whichOrbit_[iColumn]=numberUsefulOrbits_;
232        nUseful++;
233        jColumn=iColumn;
234      }
235    }
236    if (nUseful>1) {
237      numberUsefulOrbits_++;
238      numberUsefulObjects_ += nUseful;
239    } else if (jColumn>=0) {
240      assert (nUseful);
241      whichOrbit_[jColumn]=-2;
242    }
243  }
244  delete orbits;
245}
246int 
247CbcSymmetry::largestOrbit(const double * lower, const double * upper) const
248{
249  int * counts = new int[numberUsefulOrbits_];
250  memset(counts,0,numberUsefulOrbits_*sizeof(int));
251  for (int i=0;i<numberColumns_;i++) {
252    int iOrbit=whichOrbit_[i];
253    if (iOrbit>=0) {
254      if (lower[i]==0.0&&upper[i]==1.0)
255        counts[iOrbit]++;
256    }
257  }
258  int iOrbit=-1;
259  int maxOrbit=0;
260  for (int i=0;i<numberUsefulOrbits_;i++) {
261    if (counts[i]>maxOrbit) {
262      maxOrbit=counts[i];
263      iOrbit=i;
264    }
265  }
266  delete [] counts;
267  return iOrbit;
268}
269
270std::vector<int> *CbcSymmetry::Find_Orbit(int index) const{
271
272  std::vector<int> *orbit = new std::vector <int>;
273  int which_orbit = -1;
274  std::vector<std::vector<int> > *new_orbits = nauty_info_->getOrbits();
275
276  for (unsigned int i = 0; i < new_orbits -> size(); i++) {
277    for (unsigned int j = 0; j < (*new_orbits)[i].size(); j++) {
278      //   for (std::vector <int>:: iterator j = new_orbits[i].begin(); new_orbits[i].end(); ++j){
279      if( (*new_orbits)[i][j] ==  index)
280        which_orbit = i;
281    }
282  }
283 
284  //  for (std::vector <int>:: iterator j = new_orbits[which_orbit].begin(); new_orbits[which_orbit].end(), ++j)
285  for (unsigned int j = 0; j < (*new_orbits)[which_orbit].size(); j++) 
286    orbit -> push_back ((*new_orbits)[which_orbit][j]);
287
288  delete new_orbits;
289
290  return orbit;
291}
292
293
294void CbcSymmetry::ChangeBounds (const double * new_lb, const double * new_ub, 
295                                int num_cols, bool justFixedAtOne) const {
296  if (justFixedAtOne)
297    nautyFixCalls_++;
298  else
299    nautyBranchCalls_++;
300  std::sort(node_info_. begin (), node_info_. end (), index_sort);
301
302  for (int  i = 0; i < num_cols; i++) {
303    //   printf("Var %d  lower bound: %f   upper bound %f \n", i, new_lb[i], new_ub[i]);
304   
305    assert (node_info_[i].get_index () == i);
306    double newLower = new_lb[i];
307    double newUpper = new_ub[i];
308    if (justFixedAtOne) {
309      // free up all fixed at zero
310      if (!newLower) 
311        newUpper = 1.0;
312    }
313    node_info_[i ].bounds ( newLower , newUpper );
314    //printf("Var %d  INPUT lower bound: %f   upper bound %f \n", i, node_info_[i].get_lb(), node_info_[i].get_ub());
315  }
316}
317void CbcSymmetry::setupSymmetry (const OsiSolverInterface & solver) {
318  double startCPU = CoinCpuTime ();
319  const double *objective = solver.getObjCoefficients() ;
320  const double *columnLower = solver.getColLower() ;
321  const double *columnUpper = solver.getColUpper() ;
322  int numberColumns = solver.getNumCols() ;
323  int numberRows = solver.getNumRows();
324  int iRow, iColumn;
325 
326  // Row copy
327  CoinPackedMatrix matrixByRow(*solver.getMatrixByRow());
328  const double * elementByRow = matrixByRow.getElements();
329  const int * column = matrixByRow.getIndices();
330  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
331  const int * rowLength = matrixByRow.getVectorLengths();
332
333  const double * rowLower = solver.getRowLower();
334  const double * rowUpper = solver.getRowUpper();
335  //  // Find Coefficients
336
337  /// initialize nauty
338
339  int num_affine = 0;
340
341  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
342    if (objective[iColumn] && objective[iColumn]!=1.0)
343      num_affine++;
344  }
345  for (iRow = 0; iRow < numberRows; iRow++) {
346    for (CoinBigIndex j = rowStart[iRow]; 
347         j < rowStart[iRow] + rowLength[iRow]; j++) {
348      int jColumn = column[j];
349      double value = elementByRow[j];
350      if (value!=1.0)
351        num_affine++;
352    }
353  }
354
355  // Create Nauty object
356
357  int coef_count= numberRows + numberColumns +1;
358  int nc = num_affine + coef_count;
359  // create graph (part 1)
360
361  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
362    Node var_vertex;
363    var_vertex.node(iColumn,0.0,columnLower[iColumn],columnUpper[iColumn],-1,-1 );
364    node_info_.push_back(var_vertex);
365  }
366  // objective
367  int index = numberColumns;
368  {
369    Node vertex;
370    vertex.node( index , 0.0 , -COIN_DBL_MAX, COIN_DBL_MAX, 
371                 COUENNE_HACKED_EXPRGROUP, 0);
372    node_info_.push_back( vertex);
373  }
374
375  // compute space for sparse
376  size_t * v = NULL;
377  int * d = NULL;
378  int * e = NULL;
379  bool sparse=false;
380  double spaceDense = nc+WORDSIZE-1;
381  spaceDense *= nc+WORDSIZE-1;
382  spaceDense /= WORDSIZE;
383  int spaceSparse = 0;
384  {
385    size_t numberElements = 0;
386    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
387      double value = objective[iColumn];
388      if (value) {
389        if (value==1.0) {
390          numberElements+=2;
391        } else {
392          numberElements+=4;
393          coef_count ++;
394        }
395      }
396    }
397    for (iRow = 0; iRow < numberRows; iRow++) {
398      for (CoinBigIndex j = rowStart[iRow]; 
399           j < rowStart[iRow] + rowLength[iRow]; j++) {
400        int jColumn = column[j];
401        double value = elementByRow[j];
402        if (value==1.0) {
403          numberElements+=2;
404        } else {
405          numberElements+=4;
406          coef_count ++;
407        }
408      }
409    }
410    spaceSparse = 2*nc+numberElements;
411    //printf("Space for sparse is %d for dense %g\n",
412    //     spaceSparse,spaceDense);
413#ifdef NTY_TRACES
414    bool goSparse = true;
415#else
416    bool goSparse = (spaceSparse<spaceDense);
417#endif
418    // for now always sparse
419    goSparse=true;
420    if (goSparse) {
421      sparse=true;
422      v = new size_t [nc+1];
423      d = new int [nc];
424      e = new int [numberElements];
425      size_t * counts = new size_t [coef_count+1];
426      memset(counts,0,coef_count*sizeof(size_t));
427      coef_count= numberRows + numberColumns +1;
428      // create graph (part 2)
429      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
430        double value = objective[iColumn];
431        if (value) {
432          if (value==1.0) {
433            counts[index]++;
434            counts[iColumn]++;
435          } else {
436            counts[index]++;
437            counts[coef_count]+=2;
438            counts[iColumn]++;
439            coef_count ++;
440          }
441        }
442      }
443      index++;
444      for (iRow = 0; iRow < numberRows; iRow++) {
445        for (CoinBigIndex j = rowStart[iRow]; 
446             j < rowStart[iRow] + rowLength[iRow]; j++) {
447          int jColumn = column[j];
448          double value = elementByRow[j];
449          if (value==1.0) {
450            counts[index]++;
451            counts[jColumn]++;
452          } else {
453            counts[index]++;
454            counts[coef_count]+=2;
455            counts[jColumn]++;
456            coef_count ++;
457          }
458        }
459        index++;
460      }
461      // create graph (part 3)
462      assert (nc==coef_count);
463      numberElements=0;
464      v[0]=0;
465      for (int i=0;i<nc;i++) {
466        int count=counts[i];
467        d[i]=count;
468        counts[i]=v[i];
469        numberElements+=count;;
470        v[i+1]=numberElements;
471      }
472      index = numberColumns;
473      coef_count= numberRows + numberColumns +1;
474      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
475        double value = objective[iColumn];
476        if (value) {
477          int where;
478          if (value==1.0) {
479            where = counts[index];
480            counts[index]++;
481            e[where]=iColumn;
482            where = counts[iColumn];
483            counts[iColumn]++;
484            e[where]=index;
485          } else {
486            Node coef_vertex;
487            coef_vertex.node( coef_count, value, value, value, -2, 0 );
488            node_info_.push_back(coef_vertex);
489            where = counts[index];
490            counts[index]++;
491            e[where]=coef_count;
492            where = counts[coef_count];
493            counts[coef_count]++;
494            e[where]=index;
495            where = counts[coef_count];
496            counts[coef_count]++;
497            e[where]=iColumn;
498            where = counts[iColumn];
499            counts[iColumn]++;
500            e[where]=coef_count;
501            coef_count ++;
502          }
503      }
504      }
505      index++;
506      for (iRow = 0; iRow < numberRows; iRow++) {
507        Node vertex;
508        vertex.node( index , 0.0 , rowLower[iRow], rowUpper[iRow],
509                     COUENNE_HACKED_EXPRGROUP, 0);
510        node_info_.push_back( vertex);
511        for (CoinBigIndex j = rowStart[iRow]; 
512             j < rowStart[iRow] + rowLength[iRow]; j++) {
513          int jColumn = column[j];
514          double value = elementByRow[j];
515          int where;
516          if (value==1.0) {
517            where = counts[index];
518            counts[index]++;
519            e[where]=jColumn;
520            where = counts[jColumn];
521            counts[jColumn]++;
522            e[where]=index;
523          } else {
524            Node coef_vertex;
525            coef_vertex.node( coef_count, value, value, value, -2, 0 );
526            node_info_.push_back(coef_vertex);
527            where = counts[index];
528            counts[index]++;
529            e[where]=coef_count;
530            where = counts[coef_count];
531            counts[coef_count]++;
532            e[where]=index;
533            where = counts[coef_count];
534            counts[coef_count]++;
535            e[where]=jColumn;
536            where = counts[jColumn];
537            counts[jColumn]++;
538            e[where]=coef_count;
539            coef_count ++;
540          }
541        }
542        index++;
543      }
544      delete [] counts;
545    }
546  }
547
548  nauty_info_ = new CbcNauty(nc,v,d,e); 
549  delete [] v;
550  delete [] d;
551  delete [] e;
552  if (!sparse) {
553    // create graph (part 2)
554    coef_count= numberRows + numberColumns +1;
555    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
556      double value = objective[iColumn];
557      if (value) {
558        if (value==1.0) {
559          nauty_info_->addElement(index,iColumn);
560          nauty_info_->addElement(iColumn,index);
561        } else {
562          Node coef_vertex;
563          coef_vertex.node( coef_count, value, value, value, -2, 0 );
564          node_info_.push_back(coef_vertex);
565          nauty_info_->addElement(index,  coef_count);
566          nauty_info_->addElement(coef_count, index);
567          nauty_info_->addElement(coef_count,  iColumn);
568          nauty_info_->addElement(iColumn, coef_count);
569          coef_count ++;
570        }
571      }
572    }
573    index++;
574    for (iRow = 0; iRow < numberRows; iRow++) {
575      Node vertex;
576      vertex.node( index , 0.0 , rowLower[iRow], rowUpper[iRow],
577                   COUENNE_HACKED_EXPRGROUP, 0);
578      node_info_.push_back( vertex);
579      for (CoinBigIndex j = rowStart[iRow]; 
580           j < rowStart[iRow] + rowLength[iRow]; j++) {
581        int jColumn = column[j];
582        double value = elementByRow[j];
583        if (value==1.0) {
584          nauty_info_->addElement(index,jColumn);
585          nauty_info_->addElement(jColumn,index);
586        } else {
587          Node coef_vertex;
588          coef_vertex.node( coef_count, value, value, value, -2, 0 );
589          node_info_.push_back(coef_vertex);
590          nauty_info_->addElement(index,  coef_count);
591          nauty_info_->addElement(coef_count, index);
592          nauty_info_->addElement(coef_count,  jColumn);
593          nauty_info_->addElement(jColumn, coef_count);
594          coef_count ++;
595        }
596      }
597      index++;
598    }
599  }
600  numberColumns_ = numberColumns;
601  whichOrbit_ = new int [2*numberColumns_];
602  nautyBranchCalls_ = 0;
603  nautyBranchSucceeded_ = 0;
604  nautyFixCalls_ = 0;
605  nautyFixSucceeded_ = 0;
606  nautyTime_ = 0.0;
607  nautyFixes_= 0.0; 
608  nautyOtherBranches_ = 0.0;
609  Compute_Symmetry ();
610  fillOrbits();
611  //whichOrbit_[2]=numberUsefulOrbits_;
612  //Print_Orbits ();
613  // stats in array
614  if (spaceDense<COIN_INT_MAX)
615    whichOrbit_[0]=spaceDense;
616  else
617    whichOrbit_[0]=COIN_INT_MAX;
618  whichOrbit_[1]=spaceSparse;
619  double endCPU = CoinCpuTime ();
620  nautyTime_ = endCPU-startCPU;
621}
622// Fixes variables using orbits (returns number fixed)
623int 
624CbcSymmetry::orbitalFixing(OsiSolverInterface * solver)
625{
626  int numberColumns = solver->getNumCols();
627  char * status = new char [numberColumns];
628  ChangeBounds(solver->getColLower(),
629               solver->getColUpper(),
630               solver->getNumCols(),true);
631  Compute_Symmetry();
632  fillOrbits();
633  int n=0;
634  //#define PRINT_MORE 1
635  const int * alternativeOrbits = whichOrbit();
636  if (alternativeOrbits) {
637    for (int i=0;i<numberColumns;i++) {
638      char type='0';
639      if (solver->getColUpper()[i]) {
640        if (solver->getColLower()[i]) {
641          type='1';
642        } else {
643          double value=solver->getColSolution()[i];
644          if (value<0.0001)
645            type='L';
646          else if (value>0.9999)
647            type='U';
648          else
649            type='X';
650        }
651      }
652      status[i]=type;
653    }
654    n=0;
655    for (int i=0;i<numberColumns;i++) {
656      if (status[i]!='0'&&status[i]!='1') {
657        int iOrbit=alternativeOrbits[i];
658        if (iOrbit<0)
659          continue;
660        for (int j=i+1;j<numberColumns;j++) {
661          if (status[j]=='0'&&alternativeOrbits[j]==iOrbit) {
662#if PRINT_MORE>1
663            printf("In alternative orbit %d - %d free (%c), %d fixed to 0\n",
664                   iOrbit,i,status[i],j);
665#endif
666            status[i]='0'; // can fix on both branches
667            solver->setColUpper(i,0.0);
668            n++;
669            break;
670          }
671        }
672      }
673    }
674  }
675  delete [] status;
676  if (n) {
677    nautyFixSucceeded_++;
678    nautyFixes_ += n;
679#if PRINT_MORE
680    printf("%d orbital fixes\n",n);
681#endif
682  }
683  return n;
684}
685// Default Constructor
686CbcSymmetry::CbcSymmetry ()
687  : nauty_info_(NULL),
688    numberColumns_(0),
689    numberUsefulOrbits_(0),
690    numberUsefulObjects_(0),
691    whichOrbit_(NULL)
692{
693}
694// Copy constructor
695CbcSymmetry::CbcSymmetry ( const CbcSymmetry & rhs)
696{
697  node_info_ = rhs.node_info_;
698  nauty_info_ = new CbcNauty(*rhs.nauty_info_);
699  numberUsefulOrbits_ = rhs.numberUsefulOrbits_;
700  numberUsefulObjects_ = rhs.numberUsefulObjects_;
701  numberColumns_ = rhs.numberColumns_;
702  if (rhs.whichOrbit_) 
703    whichOrbit_=CoinCopyOfArray(rhs.whichOrbit_,numberColumns_);
704  else
705    whichOrbit_ = NULL;
706}
707
708// Assignment operator
709CbcSymmetry &
710CbcSymmetry::operator=( const CbcSymmetry & rhs)
711{
712  if (this != &rhs) {
713    delete nauty_info_;
714    node_info_ = rhs.node_info_;
715    nauty_info_ = new CbcNauty(*rhs.nauty_info_);
716    delete [] whichOrbit_;
717    numberColumns_ = rhs.numberColumns_;
718    numberUsefulOrbits_ = rhs.numberUsefulOrbits_;
719    numberUsefulObjects_ = rhs.numberUsefulObjects_;
720    if (rhs.whichOrbit_) 
721      whichOrbit_=CoinCopyOfArray(rhs.whichOrbit_,numberColumns_);
722    else
723      whichOrbit_ = NULL;
724  }
725  return *this;
726}
727
728// Destructor
729CbcSymmetry::~CbcSymmetry ()
730{
731  delete nauty_info_; 
732  delete [] whichOrbit_;
733}
734
735CbcNauty::CbcNauty(int vertices, const size_t * v, const int * d, const int * e)
736{
737  //printf("Need sparse nauty - wordsize %d\n",WORDSIZE);
738  n_ = vertices;
739  m_ = (n_ + WORDSIZE - 1)/WORDSIZE;
740  if (v)
741    nel_ = v[n_];
742  else
743    nel_ = 0;
744
745  //printf ("size of long = %d (%d)\nwordsize = %d\nn,m = %d,%d\n",
746  //          SIZEOF_LONG, sizeof (long), WORDSIZE, n_, m_);
747
748  nauty_check (WORDSIZE, m_, n_, NAUTYVERSIONID);
749
750  /// Apparently sizes are skewed on 64bit machines
751
752#define MULTIPLIER 1
753
754  if (!nel_) {
755    G_ = (graph *) malloc(MULTIPLIER * m_ * n_ * sizeof(int));
756    GSparse_ = NULL;
757  } else {
758    G_ = NULL;
759    GSparse_ = (sparsegraph *) malloc(sizeof(sparsegraph));
760    SG_INIT(*GSparse_);
761    SG_ALLOC(*GSparse_,n_,nel_,"malloc");
762    GSparse_->nv = n_; /* Number of vertices */
763    GSparse_->nde = nel_;
764  }
765  lab_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int)); 
766  ptn_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int));
767  active_ = NULL;
768  orbits_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int));
769#ifndef NTY_TRACES
770  options_ = (optionblk *) malloc(MULTIPLIER * sizeof(optionblk));
771  stats_ = (statsblk *) malloc(MULTIPLIER * sizeof(statsblk));
772#else
773  options_ = (TracesOptions *) malloc(MULTIPLIER * sizeof(TracesOptions));
774  stats_ = (TracesStats *) malloc(MULTIPLIER * sizeof(TracesStats));
775#endif
776  worksize_ = 100*m_;
777  workspace_ = (setword *) malloc(MULTIPLIER * worksize_*sizeof(setword));
778  canonG_ = NULL;
779  if ((G_ == 0&&GSparse_ == 0) || lab_ == 0 || ptn_ == 0 || 
780      orbits_ == 0 || options_ == 0 || stats_ == 0 ||
781      workspace_ == 0) assert(0);
782
783  // Zero allocated memory
784  if (G_) {
785    memset(G_, 0, m_*n_*sizeof(int));
786  } else {
787    //for (int i=0;i<n_;i++) {
788    //GSparse_->v[i]=v[i];
789    //}
790    memcpy(GSparse_->v,v,n_*sizeof(size_t));
791    memcpy(GSparse_->d,d,n_*sizeof(int));
792    memcpy(GSparse_->e,e,nel_*sizeof(int));
793  }
794  memset(lab_, 0, n_*sizeof(int));
795  memset(ptn_, 0, n_*sizeof(int));
796  memset(orbits_, 0, n_*sizeof(int));
797  memset(workspace_, 0, worksize_*sizeof(setword));
798#ifndef NTY_TRACES
799  memset(options_, 0,MULTIPLIER * sizeof(optionblk));
800#else
801  memset(options_, 0,MULTIPLIER * sizeof(TracesOptions));
802#endif
803
804  // Set the options you want
805#ifndef NTY_TRACES
806  options_->getcanon = FALSE;
807  options_->digraph = FALSE;
808  options_->writeautoms = FALSE;
809  options_->writemarkers = FALSE;
810  options_->defaultptn = TRUE;
811  options_->cartesian = FALSE;
812  options_->linelength = 78;
813  options_->outfile = NULL;
814  options_->userrefproc = NULL;
815  options_->userautomproc = NULL;
816  options_->userlevelproc = NULL;
817  options_->usernodeproc = NULL;
818  //  options_->usertcellproc = NULL;
819  options_->invarproc = NULL;
820  options_->tc_level = 100;
821  options_->mininvarlevel = 0;
822  options_->maxinvarlevel = 1;
823  options_->invararg = 0;
824  options_->dispatch = &dispatch_graph;
825#else
826  options_->getcanon = FALSE;
827  options_->writeautoms = FALSE;
828  options_->cartesian = FALSE;
829  options_->digraph = FALSE;
830  options_->defaultptn = TRUE;
831  options_->linelength = 78;
832#endif
833  if (G_) {
834    // Make an empty graph
835    for (int j = 0; j < n_; j++) {
836      set *gv = GRAPHROW(G_, j, m_);
837      EMPTYSET(gv, m_);
838    }
839  }
840
841  vstat_ = new int[n_];
842  clearPartitions();
843  afp_ = NULL;
844}
845
846CbcNauty::~CbcNauty()
847{
848  if (G_) free(G_);
849  if (GSparse_) {
850    SG_FREE(*GSparse_);
851    free(GSparse_);
852  }
853  if (lab_) free(lab_);
854  if (ptn_) free(ptn_);
855  if (active_) free(active_);
856  if (orbits_) free(orbits_);
857  if (options_) free(options_);
858  if (stats_) free(stats_);
859  if (workspace_) free(workspace_);
860  if (canonG_) free(canonG_);
861  if (vstat_) delete [] vstat_;
862}
863// Copy constructor
864CbcNauty::CbcNauty ( const CbcNauty & rhs)
865{
866  n_ = rhs.n_;
867  m_ = rhs.m_;
868  nel_ = rhs.nel_;
869  G_ = NULL;
870  GSparse_ = NULL;
871  if (!nel_) {
872    G_ = (graph *) malloc(MULTIPLIER * m_ * n_ * sizeof(int));
873  } else {
874    GSparse_ = (sparsegraph *) malloc(sizeof(sparsegraph));
875    SG_INIT(*GSparse_);
876    SG_ALLOC(*GSparse_,n_,nel_,"malloc");
877    GSparse_->nv = n_; /* Number of vertices */
878    GSparse_->nde = nel_;
879  }
880  lab_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int)); 
881  ptn_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int));
882  orbits_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int));
883#ifndef NTY_TRACES
884  options_ = (optionblk *) malloc(MULTIPLIER * sizeof(optionblk));
885  stats_ = (statsblk *) malloc(MULTIPLIER * sizeof(statsblk));
886#else
887  options_ = (TracesOptions *) malloc(MULTIPLIER * sizeof(TracesOptions));
888  stats_ = (TracesStats *) malloc(MULTIPLIER * sizeof(TracesStats));
889#endif
890  worksize_ = 100*m_;
891  workspace_ = (setword *) malloc(MULTIPLIER * worksize_*sizeof(setword));
892  vstat_ = new int[n_];
893  canonG_ = NULL;
894  if ((G_ == 0 && GSparse_ == 0) || lab_ == 0 || ptn_ == 0 || 
895      orbits_ == 0 || options_ == 0 || stats_ == 0 ||
896      workspace_ == 0) assert(0);
897
898  // Copy allocated memory
899  if (G_) {
900    memcpy(G_, rhs.G_, m_*n_*sizeof(int));
901  } else {
902    memcpy(GSparse_->v,rhs.GSparse_->v,n_*sizeof(size_t));
903    memcpy(GSparse_->d,rhs.GSparse_->d,n_*sizeof(int));
904    memcpy(GSparse_->e,rhs.GSparse_->e,nel_*sizeof(int));
905  }
906  memcpy(lab_, rhs.lab_, n_*sizeof(int));
907  memcpy(ptn_, rhs.ptn_, n_*sizeof(int));
908  memcpy(orbits_, rhs.orbits_, n_*sizeof(int));
909  memcpy(workspace_, rhs.workspace_, worksize_*sizeof(setword));
910#ifndef NTY_TRACES
911  memcpy(options_, rhs.options_, MULTIPLIER * sizeof(optionblk));
912  memcpy(stats_, rhs.stats_, MULTIPLIER * sizeof(statsblk));
913#else
914  memcpy(options_, rhs.options_,MULTIPLIER * sizeof(TracesOptions));
915  memcpy(stats_, rhs.stats_, MULTIPLIER * sizeof(TracesStats));
916#endif
917  memcpy(vstat_,rhs.vstat_,n_*sizeof(int));
918
919  // ? clearPartitions();
920  active_ = NULL;
921  afp_ = rhs.afp_; // ? no copy ?
922}
923
924// Assignment operator
925CbcNauty &
926CbcNauty::operator=( const CbcNauty & rhs)
927{
928  if (this != &rhs) {
929    if (G_) free(G_);
930    if (GSparse_) {
931      SG_FREE(*GSparse_);
932      free(GSparse_);
933    }
934    if (lab_) free(lab_);
935    if (ptn_) free(ptn_);
936    if (active_) free(active_);
937    if (orbits_) free(orbits_);
938    if (options_) free(options_);
939    if (stats_) free(stats_);
940    if (workspace_) free(workspace_);
941    if (canonG_) free(canonG_);
942    if (vstat_) delete [] vstat_;
943    {
944      n_ = rhs.n_;
945      m_ = rhs.m_;
946      nel_ = rhs.nel_;
947      G_ = NULL;
948      GSparse_ = NULL;
949      if (!nel_) {
950        G_ = (graph *) malloc(MULTIPLIER * m_ * n_ * sizeof(int));
951      } else {
952        GSparse_ = (sparsegraph *) malloc(sizeof(sparsegraph));
953        SG_INIT(*GSparse_);
954        SG_ALLOC(*GSparse_,n_,nel_,"malloc");
955        GSparse_->nv = n_; /* Number of vertices */
956        GSparse_->nde = nel_;
957      }
958      lab_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int)); 
959      ptn_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int));
960      orbits_ = (int *) malloc(MULTIPLIER * n_ * sizeof(int));
961#ifndef NTY_TRACES
962      options_ = (optionblk *) malloc(MULTIPLIER * sizeof(optionblk));
963      stats_ = (statsblk *) malloc(MULTIPLIER * sizeof(statsblk));
964#else
965      options_ = (TracesOptions *) malloc(MULTIPLIER * sizeof(TracesOptions));
966      stats_ = (TracesStats *) malloc(MULTIPLIER * sizeof(TracesStats));
967#endif
968      worksize_ = 100*m_;
969      workspace_ = (setword *) malloc(MULTIPLIER * worksize_*sizeof(setword));
970      vstat_ = new int[n_];
971      canonG_ = NULL;
972      if ((G_ == 0 && GSparse_ == 0) || lab_ == 0 || ptn_ == 0 || 
973          orbits_ == 0 || options_ == 0 || stats_ == 0 ||
974          workspace_ == 0) assert(0);
975     
976      // Copy allocated memory
977      if (!nel_) {
978        memcpy(G_, rhs.G_, m_*n_*sizeof(int));
979      } else {
980        memcpy(GSparse_->v,rhs.GSparse_->v,n_*sizeof(size_t));
981        memcpy(GSparse_->d,rhs.GSparse_->d,n_*sizeof(int));
982        memcpy(GSparse_->e,rhs.GSparse_->e,nel_*sizeof(int));
983      }
984      memcpy(lab_, rhs.lab_, n_*sizeof(int));
985      memcpy(ptn_, rhs.ptn_, n_*sizeof(int));
986      memcpy(orbits_, rhs.orbits_, n_*sizeof(int));
987      memcpy(workspace_, rhs.workspace_, worksize_*sizeof(setword));
988#ifndef NTY_TRACES
989      memcpy(options_, rhs.options_, MULTIPLIER * sizeof(optionblk));
990      memcpy(stats_, rhs.stats_, MULTIPLIER * sizeof(statsblk));
991#else
992      memcpy(options_, rhs.options_,MULTIPLIER * sizeof(TracesOptions));
993      memcpy(stats_, rhs.stats_, MULTIPLIER * sizeof(TracesStats));
994#endif
995      memcpy(vstat_,rhs.vstat_,n_*sizeof(int));
996     
997      // ? clearPartitions();
998      active_ = NULL;
999      afp_ = rhs.afp_; // ? no copy ?
1000    }
1001  }
1002  return *this;
1003}
1004
1005void
1006CbcNauty::addElement(int ix, int jx)
1007{
1008  // Right now die if bad index.  Can throw exception later
1009  //printf("addelement %d %d \n", ix, jx);
1010  assert(ix < n_ && jx < n_);
1011  if(ix != jx){  //No Loops
1012    set *gv = GRAPHROW(G_, ix, m_);
1013    ADDELEMENT(gv, jx);
1014    set *gv2 = GRAPHROW(G_, jx, m_);
1015    ADDELEMENT(gv2, ix);
1016    autoComputed_ = false;
1017  }
1018}
1019
1020void 
1021CbcNauty::clearPartitions()
1022{
1023  for (int j = 0; j < n_; j++) {
1024    vstat_[j] = 1;
1025    //printf("vstat %d = %d", j, vstat_[j]);
1026  }
1027
1028  autoComputed_ = false;
1029}
1030
1031void
1032CbcNauty::computeAuto()
1033{
1034
1035  //  if (autoComputed_) return;
1036
1037  //double startCPU = CoinCpuTime ();
1038
1039  options_->defaultptn = FALSE;
1040
1041  // Here we only implement the partitions
1042  // [ fix1 | fix0 (union) free | constraints ]
1043  int ix = 0;
1044 
1045  for( int color = 1; color <= n_; color++){
1046    for (int j = 0; j < n_; j++) {
1047      if (vstat_[j] == color) {
1048        lab_[ix] = j;
1049        ptn_[ix] = color;
1050        ix++;
1051      }
1052    }
1053     if (ix > 0) ptn_[ix-1] = 0;
1054  }
1055 
1056  /*
1057  for (int j = 0; j < n_; j++)
1058    printf("ptn %d = %d      lab = %d \n", j, ptn_[j], lab_[j]);
1059  */
1060
1061  // Should be number of columns
1062  assert(ix == n_);
1063  // Now the constraints if needed
1064 
1065  // Compute Partition
1066   
1067  if (G_) {
1068#ifndef NTY_TRACES
1069    nauty(G_, lab_, ptn_, active_, orbits_, options_, 
1070          stats_, workspace_, worksize_, m_, n_, canonG_);
1071#else
1072    abort();
1073#endif
1074  } else {
1075#if NAUTY_MAX_LEVEL
1076    nauty_maxalllevel=NAUTY_MAX_LEVEL;
1077#endif
1078#ifndef NTY_TRACES
1079    options_->dispatch = &dispatch_sparse;
1080    sparsenauty(GSparse_, lab_, ptn_, orbits_, options_, 
1081          stats_, NULL);
1082#else
1083    //options_->dispatch = &dispatch_sparse;
1084    Traces(GSparse_, lab_, ptn_, orbits_, options_, 
1085          stats_, NULL);
1086#endif
1087  }
1088  autoComputed_ = true;
1089
1090  //double endCPU = CoinCpuTime ();
1091
1092  //nautyTime_ += endCPU - startCPU;
1093  // Need to make sure all generators are written
1094  if (afp_) fflush(afp_);   
1095}
1096
1097void
1098CbcNauty::deleteElement(int ix, int jx)
1099{
1100  // Right now die if bad index.  Can throw exception later
1101  assert(ix < n_ && jx < n_);
1102  set *gv = GRAPHROW(G_, ix, m_);
1103  if (ISELEMENT(gv, jx)) {
1104    DELELEMENT(gv, jx);
1105  } 
1106  autoComputed_ = false;
1107}
1108
1109double
1110CbcNauty::getGroupSize() const
1111{
1112  if (!autoComputed_) return -1.0;
1113  return( stats_->grpsize1 * pow(10.0, (double) stats_->grpsize2) );
1114}
1115
1116int
1117CbcNauty::getNumGenerators() const
1118{
1119  if (!autoComputed_) return -1;
1120  return(stats_->numgenerators);
1121}
1122
1123int
1124CbcNauty::getNumOrbits() const
1125{
1126  if (!autoComputed_) return -1;
1127  return(stats_->numorbits);
1128}
1129
1130std::vector<std::vector<int> >
1131*CbcNauty::getOrbits() const
1132{
1133  std::vector<std::vector<int> > *orb = new std::vector<std::vector<int> >;
1134  if (!autoComputed_) return orb;
1135  orb -> resize(getNumOrbits());
1136  std::multimap<int, int> orbmap;
1137  std::set<int> orbkeys;
1138  for (int j = 0; j < n_; j++) {
1139    orbkeys.insert(orbits_[j]);
1140    orbmap.insert(std::make_pair(orbits_[j], j));
1141  }
1142
1143  int orbix = 0;
1144  for (std::set<int>::iterator it = orbkeys.begin();
1145       it != orbkeys.end(); ++it) {
1146    std::multimap<int, int>::iterator pos;
1147    for (pos = orbmap.lower_bound(*it);
1148         pos != orbmap.upper_bound(*it); ++pos) {
1149      (*orb)[orbix].push_back(pos->second);
1150    }
1151    orbix++;
1152  }
1153
1154  assert(orbix == getNumOrbits());
1155  return orb;
1156}
1157
1158void
1159CbcNauty::getVstat(double *v, int nv)
1160{
1161  assert(nv == n_);
1162  memcpy(v, vstat_, nv * sizeof(VarStatus));
1163}
1164
1165/*
1166bool
1167CbcNauty::isAllFixOneOrbit(const std::vector<int> &orbit) const
1168{
1169
1170  for(std::vector<int>::const_iterator it = orbit.begin();
1171      it != orbit.end(); ++it) {
1172    if (*it >= n_) return false;
1173    if (vstat_[*it] != FIX_AT_ONE) return false;
1174  }
1175  return true;
1176}
1177
1178bool
1179CbcNauty::isAllFreeOrbit(const std::vector<int> &orbit) const
1180{
1181  for(std::vector<int>::const_iterator it = orbit.begin();
1182      it != orbit.end(); ++it) {
1183    if (*it >= n_) return false;
1184    if (vstat_[*it] != FREE) return false;
1185  }
1186  return true;
1187}
1188
1189bool
1190CbcNauty::isConstraintOrbit(const std::vector<int> &orbit) const
1191{
1192  for(std::vector<int>::const_iterator it = orbit.begin();
1193      it != orbit.end(); ++it) {
1194    if (*it >= n_) return true;
1195  }
1196  return false;
1197 
1198}
1199
1200bool
1201CbcNauty::isMixedFreeZeroOrbit(const std::vector<int> &orbit) const
1202{
1203  bool containsFree = false;
1204  bool containsZero = false;
1205
1206  for(std::vector<int>::const_iterator it = orbit.begin();
1207      it != orbit.end(); ++it) {
1208    if (*it >= n_) return false;   
1209    if (vstat_[*it] == FREE) containsFree = true;
1210    if (vstat_[*it] == FIX_AT_ZERO) containsZero = true;   
1211    if (containsFree && containsZero) break;
1212  } 
1213  return (containsFree && containsZero);
1214}
1215*/
1216
1217void 
1218CbcNauty::setWriteAutoms(const std::string &fname)
1219{
1220  afp_ = fopen(fname.c_str(), "w");
1221  options_->writeautoms = TRUE;
1222#ifndef NTY_TRACES
1223  options_->writemarkers = FALSE;
1224#endif
1225  options_->outfile = afp_;
1226
1227}
1228
1229void 
1230CbcNauty::unsetWriteAutoms()
1231{
1232  fclose(afp_);
1233  options_->writeautoms = FALSE;
1234}
1235
1236// Default Constructor
1237CbcOrbitalBranchingObject::CbcOrbitalBranchingObject()
1238        : CbcBranchingObject(),
1239          column_(-1),
1240          numberOther_(0),
1241          numberExtra_(0),
1242          fixToZero_(NULL)
1243{
1244}
1245
1246// Useful constructor
1247CbcOrbitalBranchingObject::CbcOrbitalBranchingObject (CbcModel * model, int column,
1248                                                      int way ,
1249                                                      int numberExtra,
1250                                                      const int * extraToZero)
1251  : CbcBranchingObject(model, -1, way, 0.5),
1252    column_(column),
1253    numberOther_(0),
1254    numberExtra_(0),
1255    fixToZero_(NULL)
1256{
1257  CbcSymmetry * symmetryInfo = model->symmetryInfo();
1258  assert (symmetryInfo);
1259  // Filled in (hopefully)
1260  const int * orbit = symmetryInfo->whichOrbit();
1261  int iOrbit=orbit[column];
1262  assert (iOrbit>=0);
1263  int numberColumns = model->getNumCols();
1264  numberOther_=-1;
1265  for (int i=0;i<numberColumns;i++) {
1266    if (orbit[i]==iOrbit)
1267      numberOther_++;
1268  }
1269  assert (numberOther_>0);
1270  nautyBranchSucceeded_++;
1271  nautyOtherBranches_ += numberOther_;
1272  numberExtra_ = numberExtra;
1273  fixToZero_ = new int [numberOther_+numberExtra_];
1274  int n=0;
1275  for (int i=0;i<numberColumns;i++) {
1276    if (orbit[i]==iOrbit && i != column)
1277      fixToZero_[n++]=i;
1278  }
1279  for (int i=0;i<numberExtra;i++) {
1280      fixToZero_[n++]=extraToZero[i];
1281  }
1282}
1283
1284// Copy constructor
1285CbcOrbitalBranchingObject::CbcOrbitalBranchingObject (const CbcOrbitalBranchingObject & rhs)
1286        : CbcBranchingObject(rhs),
1287          column_(rhs.column_),
1288          numberOther_(rhs.numberOther_),
1289          numberExtra_(rhs.numberExtra_)
1290{
1291  fixToZero_ = CoinCopyOfArray(rhs.fixToZero_,numberOther_+numberExtra_);
1292}
1293
1294// Assignment operator
1295CbcOrbitalBranchingObject &
1296CbcOrbitalBranchingObject::operator=( const CbcOrbitalBranchingObject & rhs)
1297{
1298    if (this != &rhs) {
1299        CbcBranchingObject::operator=(rhs);
1300        delete [] fixToZero_;
1301        column_=rhs.column_;
1302        numberOther_=rhs.numberOther_;
1303        numberExtra_=rhs.numberExtra_;
1304        fixToZero_ = CoinCopyOfArray(rhs.fixToZero_,numberOther_+numberExtra_);
1305    }
1306    return *this;
1307}
1308CbcBranchingObject *
1309CbcOrbitalBranchingObject::clone() const
1310{
1311    return (new CbcOrbitalBranchingObject(*this));
1312}
1313
1314
1315// Destructor
1316CbcOrbitalBranchingObject::~CbcOrbitalBranchingObject ()
1317{
1318  delete [] fixToZero_;
1319}
1320
1321double
1322CbcOrbitalBranchingObject::branch()
1323{
1324  decrementNumberBranchesLeft();
1325  if (model_->logLevel()>1)
1326    print();
1327  OsiSolverInterface * solver = model_->solver();
1328  if (way_ < 0) {
1329    solver->setColUpper(column_,0.0);
1330    for ( int i = 0; i < numberOther_+numberExtra_; i++) {
1331      solver->setColUpper(fixToZero_[i],0.0);
1332    }
1333    way_ = 1;     // Swap direction
1334 } else {
1335    solver->setColLower(column_,1.0);
1336    for (int i = numberOther_; i < numberOther_+numberExtra_; i++) {
1337      solver->setColUpper(fixToZero_[i],0.0);
1338    }
1339    way_ = -1;    // Swap direction
1340  }
1341  return 0.0;
1342}
1343/* Update bounds in solver as in 'branch' and update given bounds.
1344   branchState is -1 for 'down' +1 for 'up' */
1345void
1346CbcOrbitalBranchingObject::fix(OsiSolverInterface * solver,
1347                           double * lower, double * upper,
1348                           int branchState) const
1349{
1350  if (branchState < 0) {
1351    upper[column_]=0.0;
1352    for ( int i = 0; i < numberOther_+numberExtra_; i++) {
1353      upper[fixToZero_[i]]=0.0;;
1354    }
1355  } else {
1356    lower[column_]=1.0;
1357    for (int i = numberOther_; i < numberOther_+numberExtra_; i++) {
1358      upper[fixToZero_[i]]=0.0;;
1359    }
1360  }
1361}
1362// Print what would happen
1363void
1364CbcOrbitalBranchingObject::print()
1365{
1366    if (way_ < 0) {
1367      printf("Orbital Down - to zero %d",column_);
1368      for ( int i = 0; i < numberOther_+numberExtra_; i++) {
1369        printf(" %d",fixToZero_[i]);
1370      }
1371    } else {
1372      printf("Orbital Up - to one %d, to zero",column_);
1373      for (int i = numberOther_; i < numberOther_+numberExtra_; i++) {
1374        printf(" %d",fixToZero_[i]);
1375      }
1376    }
1377    printf("\n");
1378}
1379
1380/** Compare the original object of \c this with the original object of \c
1381    brObj. Assumes that there is an ordering of the original objects.
1382    This method should be invoked only if \c this and brObj are of the same
1383    type.
1384    Return negative/0/positive depending on whether \c this is
1385    smaller/same/larger than the argument.
1386*/
1387int
1388CbcOrbitalBranchingObject::compareOriginalObject
1389(const CbcBranchingObject* brObj) const
1390{
1391  const CbcOrbitalBranchingObject* br =
1392    dynamic_cast<const CbcOrbitalBranchingObject*>(brObj);
1393  assert(!br);
1394  abort();
1395  return 0;
1396}
1397
1398/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1399    same type and must have the same original object, but they may have
1400    different feasible regions.
1401    Return the appropriate CbcRangeCompare value (first argument being the
1402    sub/superset if that's the case). In case of overlap (and if \c
1403    replaceIfOverlap is true) replace the current branching object with one
1404    whose feasible region is the overlap.
1405*/
1406CbcRangeCompare
1407CbcOrbitalBranchingObject::compareBranchingObject
1408(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1409{
1410  const CbcOrbitalBranchingObject* br =
1411    dynamic_cast<const CbcOrbitalBranchingObject*>(brObj);
1412  assert(!br);
1413  abort();
1414  return CbcRangeDisjoint;
1415}
1416#endif
1417
Note: See TracBrowser for help on using the repository browser.