00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <stdio.h>
00012
00013 #include "CouenneExprQuad.hpp"
00014 #include "CouenneJournalist.hpp"
00015
00016 #include "CouenneProblem.hpp"
00017 #include "CouenneExprAux.hpp"
00018 #include "CouenneExprMul.hpp"
00019 #include "CouenneExprPow.hpp"
00020 #include "CouenneLQelems.hpp"
00021
00022 using namespace Couenne;
00023
00034
00035 void CouenneProblem::decomposeTerm (expression *term,
00036 CouNumber initCoe,
00037 CouNumber &c0,
00038 LinMap &lmap,
00039 QuadMap &qmap) {
00040
00041 switch (term -> code ()) {
00042
00043
00044
00045 case COU_EXPRCONST:
00046 c0 += initCoe * term -> Value ();
00047 break;
00048
00049 case COU_EXPRVAR:
00050 lmap.insert (term -> Index (), initCoe);
00051 break;
00052
00053 case COU_EXPROPP:
00054 decomposeTerm (term -> Argument (), -initCoe, c0, lmap, qmap);
00055 break;
00056
00057 case COU_EXPRSUB:
00058 decomposeTerm (term -> ArgList () [0], initCoe, c0, lmap, qmap);
00059 decomposeTerm (term -> ArgList () [1], -initCoe, c0, lmap, qmap);
00060 break;
00061
00062 case COU_EXPRQUAD: {
00063
00064 exprQuad *t = dynamic_cast <exprQuad *> (term -> isaCopy () ?
00065 term -> Copy () :
00066 term);
00067 exprQuad::sparseQ &M = t -> getQ ();
00068
00069 for (exprQuad::sparseQ::iterator row = M.begin ();
00070 row != M.end (); ++row) {
00071
00072 int xind = row -> first -> Index ();
00073
00074 for (exprQuad::sparseQcol::iterator col = row -> second.begin ();
00075 col != row -> second.end (); ++col) {
00076 qmap.insert (xind, col -> first -> Index (), initCoe * col -> second);
00077 }
00078 }
00079 }
00080
00081 case COU_EXPRGROUP: {
00082
00083 exprGroup *t = dynamic_cast <exprGroup *> (term -> isaCopy () ?
00084 term -> Copy () :
00085 term);
00086 exprGroup::lincoeff &lcoe = t -> lcoeff ();
00087
00088
00089 for (int n = lcoe.size (), i=0; n--; i++)
00090 lmap.insert (lcoe [i].first -> Index (), initCoe * lcoe [i].second);
00091
00092 c0 += initCoe * t -> getc0 ();
00093 }
00094
00095 case COU_EXPRSUM: {
00096
00097 expression **al = term -> ArgList ();
00098 for (int i = term -> nArgs (); i--;)
00099 decomposeTerm (*al++, initCoe, c0, lmap, qmap);
00100
00101 } break;
00102
00103
00104
00105
00106
00107 case COU_EXPRMUL: {
00108
00109 std::map <int, CouNumber> indices;
00110 CouNumber coe = initCoe;
00111
00112
00113 flattenMul (term, coe, indices);
00114
00115 if (jnlst_ -> ProduceOutput (Ipopt::J_ALL, J_REFORMULATE)) {
00116 printf ("from flattenmul: [%g] ", coe);
00117 for (std::map <int, CouNumber>::iterator itt = indices.begin ();
00118 itt != indices.end(); ++itt)
00119 printf (" %d,%g",
00120 itt -> first,
00121 itt -> second);
00122 printf ("\n");
00123 }
00124
00125
00126 switch (indices.size ()) {
00127
00128 case 0:
00129 c0 += coe;
00130 break;
00131
00132 case 1: {
00133
00134 std::map <int, CouNumber>::iterator one = indices.begin ();
00135 int index = one -> first;
00136 CouNumber expon = one -> second;
00137
00138 if (fabs (expon - 1) < COUENNE_EPS) lmap.insert (index, coe);
00139 else if (fabs (expon - 2) < COUENNE_EPS) qmap.insert (index, index, coe);
00140 else {
00141 exprAux *aux = addAuxiliary
00142 (new exprPow (new exprClone (Var (index)),
00143 new exprConst (expon)));
00144
00145
00146
00147 lmap.insert (aux -> Index (), coe);
00148 }
00149 } break;
00150
00151 case 2: {
00152
00153 int ind0, ind1;
00154
00155 std::map <int, CouNumber>::iterator one = indices.begin (),
00156 two = one;
00157 ++two;
00158
00159
00160 if (fabs (one -> second - 1) > COUENNE_EPS) {
00161 exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (one -> first)),
00162 new exprConst (one -> second)));
00163 ind0 = aux -> Index ();
00164 } else ind0 = one -> first;
00165
00166
00167 if (fabs (two -> second - 1) > COUENNE_EPS) {
00168 exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (two -> first)),
00169 new exprConst (two -> second)));
00170 ind1 = aux -> Index ();
00171 } else ind1 = two -> first;
00172
00173 qmap.insert (ind0, ind1, coe);
00174 } break;
00175
00176 default: {
00177
00178
00179
00180 expression **al = new expression * [indices.size ()];
00181 std::map <int, CouNumber>::iterator one = indices.begin ();
00182
00183 for (int i=0; one != indices.end (); ++one, i++)
00184 if (fabs (one -> second - 1) > COUENNE_EPS) {
00185 exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (one -> first)),
00186 new exprConst (one -> second)));
00187 al [i] = new exprClone (aux);
00188 } else al [i] = new exprClone (Var (one -> first));
00189
00190
00191
00192
00193 exprMul *mul = new exprMul (al, indices.size ());
00194 exprAux *aux = mul -> standardize (this);
00195 lmap.insert (aux -> Index (), coe);
00196
00197 } break;
00198 }
00199
00200 } break;
00201
00202 case COU_EXPRPOW: {
00203
00204 expression **al = term -> ArgList ();
00205
00206 if (al [1] -> Type () != CONST) {
00207
00208
00209
00210
00211 expression *aux = term -> standardize (this);
00212 if (!aux) aux = term;
00213 lmap.insert (aux -> Index (), initCoe);
00214
00215 } else {
00216
00217
00218 expression *aux = (*al) -> standardize (this);
00219
00220 if (!aux)
00221 aux = *al;
00222
00223 CouNumber expon = al [1] -> Value ();
00224 int ind = aux -> Index ();
00225
00226 if (fabs (expon - 1) < COUENNE_EPS) lmap.insert (ind, initCoe);
00227 else if (fabs (expon - 2) < COUENNE_EPS) qmap.insert (ind, ind, initCoe);
00228 else {
00229 exprAux *aux2 = addAuxiliary
00230 (new exprPow (new exprClone (aux), new exprConst (expon)));
00231 lmap.insert (aux2 -> Index (), initCoe);
00232 }
00233 }
00234 } break;
00235
00236 default: {
00237
00238 expression *aux = term -> standardize (this);
00239 if (!aux)
00240 aux = term;
00241 lmap.insert (aux -> Index (), initCoe);
00242 } break;
00243 }
00244 }