source: trunk/Cbc/src/CbcSymmetry.cpp

Last change on this file was 2491, checked in by forrest, 6 months ago

better SOS in mipstart, ctrl-c back, improve symmetric

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