1 /* $Id$
2 *
3 * Name: CouenneSOS.cpp
4 * Author: Pietro Belotti
5 * Purpose: find SOS constraints in problem and add them to list of
6 * branching objects
7 *
8 * (C) Carnegie-Mellon University, 2008-10.
9 * This file is licensed under the Eclipse Public License (EPL)
10 */
11
12 #include <vector>
13
14 #include "CouenneExprGroup.hpp"
15 #include "CouenneExprAux.hpp"
16
17 #include "CbcModel.hpp"
18 #include "CbcBranchActual.hpp"
19 #include "CbcCutGenerator.hpp"
20 #include "CbcCompareActual.hpp"
21
22 #include "CouenneProblem.hpp"
23
24 //using namespace Osi;
25 using namespace Couenne;
26
27 /// find SOS objects
findSOS(CbcModel * CbcModelPtr,OsiSolverInterface * solver,OsiObject ** objects)28 int CouenneProblem::findSOS (CbcModel *CbcModelPtr,
29 OsiSolverInterface *solver,
30 OsiObject **objects) {
31
32 // check auxiliaries defined as
33 // x_i binary. Disable it and add relative SOS to array "objects"
34
35 int nSOS = 0;
36
37 for (std::vector <exprVar *>::const_iterator v = variables_.begin ();
38 v != variables_.end (); ++v)
39
40 if (((*v) -> Type () == AUX) &&
41 ((*v) -> Multiplicity () > 0) &&
42 ((*v) -> sign () == expression::AUX_EQ) &&
43 ((*v) -> Image () -> code () == COU_EXPRGROUP)) {
44
45 expression *img = (*v) -> Image ();
46
47 exprGroup *group = dynamic_cast <exprGroup *> (img -> isaCopy () ?
48 img -> Copy () :
49 img);
50 if (!group)
51 continue;
52
53 int wind = (*v) -> Index ();
54 CouNumber cterm = group -> getc0 ();
55 bool
56 defVar = true,
57 invertSOS = false;
58
59 // now check if this is
60 //
61 // defvar==true:
62 // 1) an auxiliary fixed to one ==> it's a SOS if its image is x1+x2+...+xk
63 // 1a) an auxiliary fixed to -1 ==> it's a SOS if its image is -x1-x2-...-xk (see minlplib/csched2)
64 //
65 // defvar==false:
66 // 2) an auxiliary fixed to zero ==> it's a SOS if its image is -x1-x2-...-xk+1
67 // 2a) or if its image is x1+x2+...+xk-1
68
69 if (fabs (cterm - 1.) < COUENNE_EPS) {defVar = false;}
70 else if (fabs (cterm + 1.) < COUENNE_EPS) {defVar = false; invertSOS = true;}
71 else if (fabs (cterm) > COUENNE_EPS) continue; // and defVar is true
72
73 if (defVar) { // implies cterm == 0
74 if ((fabs (Lb (wind) + 1.) < COUENNE_EPS) && (fabs (Ub (wind) + 1.) < COUENNE_EPS)) invertSOS = true;
75 else if (!((fabs (Lb (wind) - 1.) < COUENNE_EPS) && (fabs (Ub (wind) - 1.) < COUENNE_EPS))) continue;
76 } else
77
78 if ((fabs (Lb (wind)) > COUENNE_EPS) ||
79 (fabs (Ub (wind)) > COUENNE_EPS))
80 continue;
81
82 size_t lsz = group -> lcoeff (). size ();
83
84 if (((lsz <= 2) && defVar) ||
85 ((lsz <= 1) && !defVar))
86 continue;
87
88 // there are two possibilities:
89 //
90 // 1) w is defined as w = 1 - \sum_{i=0}^n x_i -- defvar = false
91 // 2) w is defined as \sum_{i=0}^n x_i and w \in [1,1] -- defvar = true
92
93 bool
94 intSOS = (*v) -> isInteger (),
95 isSOS = true,
96 onlyOrigVars = true; // if SOS constraint only contains
97 // original variables, it has been
98 // spotted by Cbc already
99
100 exprGroup::lincoeff &lcoe = group -> lcoeff ();
101 exprGroup::lincoeff::iterator l = lcoe. begin ();
102
103 for (;l != lcoe. end (); ++l) {
104
105 if ((fabs (l -> second - (invertSOS ? -1. : 1.)) > COUENNE_EPS) || // wrong coefficient?
106 (fabs (Lb (l -> first -> Index ())) > COUENNE_EPS)) { // positive lower bound?
107
108 isSOS = false;
109 break;
110
111 } else
112 if (!(l -> first -> isInteger ()))
113 intSOS = false;
114
115 if (l -> first -> Index () >= nOrigVars_) //
116 onlyOrigVars = false;
117 }
118
119 if (!isSOS || !intSOS)// || onlyOrigVars)
120 continue;
121
122 // printf ("----- found SOS: ");
123 // (*v) -> print (); printf (" := ");
124 // (*v) -> Image () -> print (); printf ("\n");
125
126 // it is a SOS -- if intSOS==true, it's also integer
127
128 int
129 indStart = defVar ? 0 : 1,
130 nelem = indStart + lcoe. size (),
131 *indices = new int [nelem];
132
133 if (!defVar)
134 indices [0] = (*v) -> Index ();
135
136 for (int i=indStart, j=0; i<nelem; i++)
137 indices [i] = lcoe [j++]. first -> Index ();
138
139 // TODO: if use Cbc, add CbcSOSBranchingObject
140
141 //CouenneSOSObject *newsos = new CouenneSOSObject (solver, nelem, indices, NULL, 1, jnlst_, true, true);
142 //OsiSOS *newsos = new OsiSOS (solver, nelem, indices, NULL, 1);
143 CbcSOS *newsos = new CbcSOS (CbcModelPtr, nelem, indices, NULL, nSOS, 1);
144
145 objects [nSOS] = newsos;
146 // as in BonBabSetupBase.cpp:675
147 newsos -> setPriority (10);
148 newsos -> setIntegerValued (intSOS);
149
150 nSOS++;
151 }
152
153 if (nSOS)
154 jnlst_ -> Printf (Ipopt::J_ERROR, J_COUENNE, "%d SOS constraint%s found\n", nSOS, nSOS == 1 ? "" : "s");
155
156 return nSOS;
157 }
158