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