# source:trunk/Couenne/src/problem/CouenneSymmetry.cpp@792

Last change on this file since 792 was 792, checked in by pbelotti, 2 years ago

fixing orbital branching in non-strong branching. Setting default priorities of continuous/integer variables to <100 to get ahead of CbcObjects?. Fixing some memory leaks

• Property svn:keywords set to `Author Date Id Revision`
File size: 13.9 KB
Line
1/* \$Id\$
2 *
3 * Name:    CouenneSymmetry.cpp
4 * Author:  Jim Ostrowski
5 * Purpose: methods for exploiting symmetry
6 * Date:    October 13, 2010
7 *
9 */
10
11#include <cassert>
12#include <vector>
13#include <algorithm>
14#include <ostream>
15#include <iterator>
16#include <stdio.h>
17
18#include "CouenneExprVar.hpp"
19#include "CouenneExprGroup.hpp"
20
21#include "CouenneProblem.hpp"
22
23using namespace Couenne;
24
25#ifdef COIN_HAS_NTY
26
27#include "Nauty.h"
28
29void Node::node(int i, double c , double l, double u, int cod, int s){
30  index = i;
31  coeff = c;
32  lb = l;
33  ub = u;
34  color = -1;
35  code = cod;
36  sign = s;
37}
38
39inline bool CouenneProblem::compare (register Node &a, register Node &b) const {
40  if(a.get_code() == b.get_code() )
41    if(a.get_coeff() == b.get_coeff() )
42      if(a.get_sign() == b.get_sign() )
43        if( fabs ( a.get_lb() - b.get_lb() ) <= COUENNE_EPS )
44          if( fabs ( a.get_ub() - b.get_ub() ) <= COUENNE_EPS )
45            return 1;
46  return 0;
47}
48
49/*
50bool CouenneProblem::node_sort (Node a, Node  b){
51  bool is_less = 0;
52
53  if(a.get_code() < b.get_code() )
54    is_less = 1;
55  else {
56    if(a.get_code() == b.get_code() )
57      if(a.get_coeff() < b.get_coeff() )
58        is_less = 1;
59      else{
60        if(a.get_coeff() ==  b.get_coeff() )
61          if(a.get_lb() < b.get_lb())
62            is_less = 1;
63          else{
64            if(a.get_lb() == b.get_lb())
65              if(a.get_ub() < b.get_ub())
66                is_less = 1;
67              else{
68                if(a.get_index() < b.get_index())
69                  is_less = 1;
70              }
71          }
72      }
73  }
74  return is_less;
75}
76bool CouenneProblem::index_sort (Node a, Node b){
77  return (a.get_index() < b.get_index() );
78}
79*/
80
81void CouenneProblem::sym_setup (){
82
83  //  // Find Coefficients
84
85  /// initialize nauty
86
87  int num_affine = 0;
88
89  for (std::vector <exprVar *>:: iterator i = Variables (). begin ();
90       i != Variables (). end (); ++i) {
91
92    if ((*i) -> Type () == AUX) {
93      if ((*i) -> Image () -> code () == COU_EXPRDIV) {
94              num_affine ++;
95      }
96      if ((*i) -> Image () -> code () != COU_EXPRGROUP) {
97        if ((*i) -> Image () -> Type () == N_ARY) {
98          for (int a=0; a < (*i) -> Image () -> nArgs(); a++) {
99            expression *arg = (*i) -> Image () -> ArgList () [a];
100
101            if (arg -> Type () == CONST) {
102              num_affine ++;
103
104            }
105          }
106        }
107      }
108      if ((*i) -> Image () -> code () == COU_EXPRGROUP) {
109
110        exprGroup *e = dynamic_cast <exprGroup *> ((*i) -> Image ());
111
112        // add a node for e -> getC0 ();
113        if (e -> getc0 () != 0 ){
114          num_affine ++;
115        }
116
117        // for each term add nodes for their non-one coefficients and their variable
118
119        for (exprGroup::lincoeff::iterator el = e ->lcoeff().begin (); el != e -> lcoeff ().end (); ++el) {
120          if ( el -> second !=1){
121            num_affine ++;
122          }
123        }
124      }
125    }
126  }
127
128  // Create global Nauty object
129
130  int nc = num_affine + nVars ();
131  // printf (" There are   %d  coefficient vertices in the graph \n", num_affine);
132  //printf (" Graph has    %d  vertices \n", nc);
133
134  nauty_info = new Nauty(nc);
135  // create graph
136
137  int coef_count= nVars ();
138  for (std::vector <exprVar *>:: iterator i =  Variables (). begin ();
139       i != Variables (). end (); ++i) {
140
141    //    printf ("I have code %d \n",  (*i) ->  Image() -> code() );
142
143    if ((*i) -> Type () == AUX) {
144      //printf ("aux is %d with code %d \n", (*i) -> Index (), (*i) -> Image () -> code() );
145      // this is an auxiliary variable
146
147      Node vertex;
148      vertex.node( (*i) -> Index () , 0.0 , (*i) -> lb () , (*i) -> ub () ,  (*i) -> Image () -> code(), (*i)-> sign() );
149      //printf(" sign of aux %d \n", (*i) -> sign () );
150      node_info.push_back( vertex);
151
152      // add node in nauty graph for its index, (*i) -> Index ()
153
154      if ((*i) -> Image () -> Type () == N_ARY) {
155
156        if ((*i) -> Image () -> code () == COU_EXPRDIV) {
157          expression *arg = (*i) -> Image () -> ArgList () [0];
158          nauty_info->addElement((*i) -> Index (),  arg -> Index ());
159          expression *arg2 = (*i) -> Image () -> ArgList () [1];
160          nauty_info->addElement((*i) -> Index (),  coef_count);
161          nauty_info->addElement( coef_count,  arg2 -> Index ());
162          Node coef_vertex;
163          coef_vertex.node( coef_count, -1, -1 ,-1, -2 , 0);
164          node_info.push_back(coef_vertex);
165          coef_count ++;
166        }
167
168        if ((*i) -> Image () -> code () != COU_EXPRGROUP) {
169
170          for (int a=0; a < (*i) -> Image () -> nArgs(); a++) {
171            expression *arg = (*i) -> Image () -> ArgList () [a];
172
173            if (arg -> Type () != CONST) {
174              //printf (" add edge  %d , %d\n", (*i) -> Index (),  arg -> Index ());
175              nauty_info->addElement((*i) -> Index (),  arg -> Index ());
176              nauty_info->addElement( arg -> Index (), (*i) -> Index ());
177            }
178
179            else{
180
181              assert (arg -> Type () == CONST);
182
183              // printf (" add new vertex to graph, coef # %d, value %g \n", coef_count, arg -> Value() );
184              // printf (" add edge aux index %d ,  coef index %d\n", (*i) -> Index (),  coef_count);
185              nauty_info->addElement((*i) -> Index (),  coef_count);
186              nauty_info->addElement( coef_count, (*i) -> Index ());
187
188
189              Node coef_vertex;
190              coef_vertex.node( coef_count, arg -> Value(), arg -> Value() , arg -> Value(), -2 , 0);
191              node_info.push_back(coef_vertex);
192              coef_count ++;
193            }
194
195          }
196        }
197
198
199        if ((*i) -> Image () -> code () == COU_EXPRGROUP) {
200
201          // dynamic_cast it to an exprGroup
202          exprGroup *e = dynamic_cast <exprGroup *> ((*i) -> Image ());
203
204          // add a node for e -> getC0 ();
205          if (e -> getc0 () != 0 ){
206            Node coef_vertex;
207            coef_vertex.node( coef_count, e -> getc0(), e -> getc0() , e -> getc0(), -2, 0 );
208            node_info.push_back(coef_vertex);
209
210            //printf ("Add coef vertex to graph (coef value   %f) \n", e -> getc0 () );
211            //printf (" add edge aux index %d ,  coef index %d\n", (*i) -> Index (), coef_count);
212            nauty_info->addElement((*i) -> Index (),  coef_count);
213            nauty_info->addElement( coef_count, (*i) -> Index ());
214
215
216            coef_count ++;
217          }
218
219          // for each term add nodes for their non-one coefficients and their variable
220
221          for (exprGroup::lincoeff::iterator el = e ->lcoeff().begin (); el != e -> lcoeff ().end (); ++el) {
222
223            if ( el -> second ==1){
224              //printf (" add edge index %d ,  index %d\n", (*i) -> Index (), el -> first -> Index()    );
225              nauty_info->addElement((*i) -> Index (),  el -> first -> Index());
226              nauty_info->addElement( el -> first -> Index (), (*i) -> Index ());
227            }
228            else{
229              //printf (" add new vertex to graph, coef # %d with coef %f \n", coef_count, el -> second);
230              Node coef_vertex;
231              coef_vertex.node( coef_count, el -> second, el -> second, el -> second, -2, 0 );
232              node_info.push_back(coef_vertex);
233
234              //printf (" add edge aux index %d ,  coef index %d\n", (*i) -> Index (), coef_count);
235              nauty_info->addElement((*i) -> Index (),  coef_count);
236              nauty_info->addElement( coef_count, (*i) -> Index ());
237
238              // printf (" add edge coef index %d ,  2nd index %d\n", coef_count,  el -> first -> Index()  );
239              nauty_info->addElement(coef_count,  el -> first -> Index());
240              nauty_info->addElement( el -> first -> Index (), coef_count);
241              coef_count ++;
242            }
243            // coefficient = el -> second
244
245            // variable index is el -> first -> Index ()
246          }
247
248        }
249
250      }
251      else if ((*i) -> Image () -> Type () == UNARY) {
252        //      printf ("variable is unary  %d\n", (*i) -> Index ());
253        expression *arg = (*i) -> Image () -> Argument () ;
254        nauty_info->addElement( arg-> Index(), (*i) -> Index() );
255        //printf (" add edge aux index %d ,  coef index %d\n", (*i) -> Index (), arg-> Index());
256      }
257      else if ((*i) -> Image () -> Type () == AUX) {
258        //printf ("variable is AUX  %d\n", (*i) -> Index ());
259        nauty_info->addElement((*i) -> Index (), (*i) -> Image() -> Index());
260        nauty_info->addElement( (*i) -> Image() -> Index(), (*i) -> Index() );
261        //printf (" add edge aux index %d ,  coef index %d\n", (*i) -> Index (), (*i) -> Image() -> Index());
262      }
263      else if ((*i) -> Image () -> Type () == VAR) {
264        //printf ("variable is VAR  %d, image %d \n", (*i) -> Index (), (*i) -> Image() -> Index());
265        nauty_info->addElement((*i) -> Index (), (*i) -> Image() -> Index());
266        nauty_info->addElement( (*i) -> Image() -> Index(), (*i) -> Index() );
267
268        //printf (" add edge aux index %d ,  coef index %d\n", (*i) -> Index (), (*i) -> Image() -> Index());
269      }
270    }
271    else {
272      // printf ("variable is %d\n", (*i) -> Index ());
273      Node var_vertex;
274
275      // Bounds of +- infinity make the compare function likely to return a false negative. Rather than add inf as a boud, I use lb-1 (or ub +1
276      if( (*i) -> ub() >= COUENNE_INFINITY && (*i) -> lb() <= - COUENNE_INFINITY){
277        var_vertex.node( (*i) -> Index () , 0 , 1 , 0 ,  -1, -1 );
278        node_info.push_back(var_vertex);
279        //printf( "var info index %d, lb %f, ub %f \n",(*i) -> Index () , 1 , 0 ) ;
280      }
281      else  if( (*i) -> ub() >= COUENNE_INFINITY ){
282        var_vertex.node( (*i) -> Index () , 0 , (*i) -> lb () , (*i) -> lb() -1 ,  -1, -1 );
283        node_info.push_back(var_vertex);
284        //printf( "var info index %d, lb %f, ub %f \n",(*i) -> Index () , (*i) -> lb () , (*i) -> lb () -1 ) ;
285      }
286      else  if( (*i) -> lb() <= - COUENNE_INFINITY){
287        var_vertex.node( (*i) -> Index () , 0 , (*i) -> ub () +1 , (*i) -> ub () ,  -1, -1 );
288        node_info.push_back(var_vertex);
289        //printf( "var info index %d, lb %f, ub %f \n",(*i) -> Index () , (*i) -> ub () +1 , (*i) -> ub () ) ;
290      }
291      else{
292        var_vertex.node( (*i) -> Index () , 0 , (*i) -> lb () , (*i) -> ub () ,  -1, -1 );
293        //printf( "var info index %d, lb %f, ub %f \n",(*i) -> Index () , (*i) -> lb () , (*i) -> ub () ) ;
294        // var_vertex.get_index() , var_vertex.get_coeff() , var_vertex.get_lb() , var_vertex.get_ub() ,  var_vertex.get_code() );
295        node_info.push_back(var_vertex);
296        // this is an original variable
297      }
298    }
299  }
300
301}
302
303
304void CouenneProblem::Compute_Symmetry() const{
305
306   //  ChangeBounds (Lb (), Ub (), nVars ());
307
308  // jnlst_ -> Printf(Ipopt::J_VECTOR, J_BRANCHING,"== Computing Symmetry\n");
309  // for (int i = 0; i < nVars (); i++)
310  //   if (Var (i) -> Multiplicity () > 0)
311  //     jnlst_->Printf(Ipopt::J_VECTOR, J_BRANCHING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
312  //                 X  (i), Lb (i), Ub (i));
313  // jnlst_->Printf(Ipopt::J_VECTOR, J_BRANCHING,"=============================\n");
314
315  std::sort(node_info. begin (), node_info. end (), node_sort);
316
317  for (std::vector <Node>:: iterator i = node_info. begin (); i != node_info. end (); ++i)
318    (*i).color_vertex(-1);
319
320  int color = 1;
321  for (std::vector <Node>:: iterator i = node_info. begin (); i != node_info. end (); ++i) {
322    if( (*i).get_color() == -1){
323      (*i).color_vertex(color);
324      //printf ("Graph vertex %d is given color %d\n", (*i).get_index(), color);
325      nauty_info -> color_node((*i).get_index(), color);
326      for (std::vector <Node>:: iterator j = i+1; j != node_info. end (); ++j)
327        if( compare( (*i) , (*j) ) ==1){
328          (*j).color_vertex(color);
329          nauty_info -> color_node((*j).get_index(),color);
330          //printf ("Graph vertex %d is given color %d, the same as vertex %d\n", (*j).get_index(), color, (*i).get_index());
331        }
332      //       else
333      // j = node_info. end();
334      color++;
335    }
336  }
337
338  //Print_Orbits ();
339
340  nauty_info -> computeAuto();
341}
342
343
344void CouenneProblem::Print_Orbits () const {
345
346  //printf ("num gens = %d, num orbits = %d \n", nauty_info -> getNumGenerators(), nauty_info -> getNumOrbits() );
347
348  std::vector<std::vector<int> > *new_orbits = nauty_info->getOrbits();
349
350  printf("There are %d orbits and %d generators",
351  nauty_info->getNumOrbits(),
352  nauty_info->getNumGenerators());
353
354  int nNonTrivialOrbits = 0;
355
356  for (unsigned int i = 0; i < new_orbits -> size(); i++) {
357    if ((*new_orbits)[i].size() > 1)
358      nNonTrivialOrbits++;
359    else continue;
360
361    // int orbsize = (*new_orbits)[i].size();
362    // printf( "Orbit %d [size: %d] [", i, orbsize);
363    // copy ((*new_orbits)[i].begin(), (*new_orbits)[i].end(),
364    //    std::ostream_iterator<int>(std::cout, " "));
365    // printf("] \n");
366  }
367
368  printf (" (%d non-trivial orbits).\n", nNonTrivialOrbits);
369
370#if 0
371  if (nNonTrivialOrbits)
372    for (int i=0; i< nVars (); i++) {
373
374      std::vector< int > *branch_orbit = Find_Orbit (i);
375
376      if (branch_orbit -> size () > 1) {
377        printf ("x%04d: ", i);
378
379        for (std::vector<int>::iterator it = branch_orbit -> begin (); it != branch_orbit -> end (); ++it)
380          printf ("%d ", *it);
381        printf ("\n");
382      }
383    }
384#endif
385  delete new_orbits;
386}
387
388std::vector<int> *CouenneProblem::Find_Orbit(int index) const{
389
390  std::vector<int> *orbit = new std::vector <int>;
391  int which_orbit = -1;
392  std::vector<std::vector<int> > *new_orbits = nauty_info->getOrbits();
393
394  for (unsigned int i = 0; i < new_orbits -> size(); i++) {
395    for (unsigned int j = 0; j < (*new_orbits)[i].size(); j++) {
396      //   for (std::vector <int>:: iterator j = new_orbits[i].begin(); new_orbits[i].end(); ++j){
397      if( (*new_orbits)[i][j] ==  index)
398        which_orbit = i;
399    }
400  }
401
402  //  for (std::vector <int>:: iterator j = new_orbits[which_orbit].begin(); new_orbits[which_orbit].end(), ++j)
403  for (unsigned int j = 0; j < (*new_orbits)[which_orbit].size(); j++)
404    orbit -> push_back ((*new_orbits)[which_orbit][j]);
405
406  delete new_orbits;
407
408  return orbit;
409}
410
411
412void CouenneProblem::ChangeBounds (const double * new_lb, const double * new_ub, int num_cols) const {
413  assert (num_cols == nVars ()); // replaced Variables () . size () as Variables () is not a const method
414  std::sort(node_info. begin (), node_info. end (), index_sort);
415
416  for (int  i = 0; i < num_cols; i++) {
417    //   printf("Var %d  lower bound: %f   upper bound %f \n", i, new_lb[i], new_ub[i]);
418
419    assert (node_info[i].get_index () == i);
420    node_info[i ].bounds ( new_lb[i] , new_ub[i] );
421    // printf("Var %d  INPUT lower bound: %f   upper bound %f \n", i, node_info[i].get_lb(), node_info[i].get_ub());
422  }
423}
424#endif
425
426void CouenneProblem::setupSymmetry () {
427
428#ifdef COIN_HAS_NTY
429  sym_setup ();
430  Compute_Symmetry ();
431  Print_Orbits ();
432#else
433  if (orbitalBranching_) {
434    printf ("\
435Couenne: Warning, you have set orbital_branching but Nauty is not available.\n\
436Reconfigure with appropriate options --with-nauty-lib=/path/to/libnauty.* and --with-nauty-incdir=/path/to/nauty/include/files/\n");
437  }
438#endif
439}
Note: See TracBrowser for help on using the repository browser.