source: stable/2.9/Cbc/src/CbcSymmetry.cpp @ 2208

Last change on this file since 2208 was 2208, checked in by stefan, 3 years ago

merge r2207 from trunk

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