Changeset 1049


Ignore:
Timestamp:
Jan 26, 2014 7:59:21 AM (6 years ago)
Author:
pbelotti
Message:

fix in exprTrilinear when one variable is fixed

Location:
trunk/Couenne/src/convex/operators
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Couenne/src/convex/operators/conv-exprMul-genCuts.cpp

    r490 r1049  
    126126  // out of the hyperbola's belly
    127127
    128   CouNumber x0 = (*(arglist_ [0])) (),
    129             y0 = (*(arglist_ [1])) ();
    130 
    131   unifiedProdCuts (cg, cs,
    132                    xi, x0,      xl, xu,
    133                    yi, y0,      yl, yu,
    134                    wi, (*w) (), wl, wu,
     128  unifiedProdCuts (cg, cs,
     129                   xi, (*(arglist_ [0])) (), xl, xu,
     130                   yi, (*(arglist_ [1])) (), yl, yu,
     131                   wi, (*w) (),              wl, wu,
    135132                   //sign == expression::AUX_LEQ ? -COIN_DBL_MAX : wl,
    136133                   //sign == expression::AUX_GEQ ?  COIN_DBL_MAX : wu,
  • trunk/Couenne/src/convex/operators/conv-exprTrilinear-gencuts.cpp

    r1047 r1049  
    793793#endif
    794794    } else {
    795     // compute the 6 permutations of the 3 variables
    796     ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
    797     permutation3(ind,ibnd);
    798     int i, flagg=0, idx=0;
    799     i = 0;
    800     while(i < 6 && flagg == 0) {
    801       if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
    802          <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
    803          vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
    804          <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
    805         {
    806           idx = i;   // store the index of the permutation satisfying the condition
    807           flagg = 1;  // condition is satisfied
    808         }
    809       i++;
    810     }
    811     if (flagg==0) {
    812       std::cout << "ERROR!!!" << std::endl; exit(0);
    813     }
    814     }
    815     v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
     795      // compute the 6 permutations of the 3 variables
     796      ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
     797      permutation3(ind,ibnd);
     798      int i, flagg=0, idx=0;
     799      i = 0;
     800      while(i < 6 && flagg == 0) {
     801        if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
     802           <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
     803           vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
     804           <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
     805          {
     806            idx = i;   // store the index of the permutation satisfying the condition
     807            flagg = 1;  // condition is satisfied
     808          }
     809        i++;
     810      }
     811      if (flagg==0) {
     812        std::cout << "ERROR!!!" << std::endl; exit(0);
     813      }
     814   
     815      v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
     816    }
    816817
    817818    double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
     
    11521153  int varInd [4];
    11531154
     1155  enum auxSign sign = cg -> Problem () -> Var (w -> Index ()) -> sign ();
     1156
    11541157  for (int i=0; i<3; i++)
    11551158    varInd [i] = args [i] -> Index ();
     
    11601163  std::vector <std::vector <double> > cutCoeff;
    11611164  std::vector <double>                cutLb, cutUb;
     1165
     1166  // These cuts rely on three non-fixed variables: if l_i=u_i for i in
     1167  // 1,2,3, there is a NaN due to division by (u_i - l_i)
     1168
     1169  int n_var_fixed = 0, isFixed [3] = {0,0,0};
     1170  double fixed_prod = 1.;
     1171
     1172  for (int i=0; i<3; ++i) {
     1173
     1174    double lb, ub;
     1175
     1176    ArgList () [i] -> getBounds (lb, ub);
     1177
     1178    if (fabs (ub - lb) < COUENNE_EPS) {
     1179      isFixed [i] = 1;
     1180      ++ n_var_fixed;
     1181      fixed_prod *= (.5*(lb+ub));
     1182    }
     1183  }
     1184
     1185  if (n_var_fixed) {
     1186
     1187    if ((fixed_prod < COUENNE_EPS) || (n_var_fixed == 3)) {
     1188
     1189      if (!(cg->createCut (cs,
     1190                           (sign == expression::AUX_LEQ) ? - COIN_DBL_MAX : fixed_prod,
     1191                           (sign == expression::AUX_GEQ) ?   COIN_DBL_MAX : fixed_prod,
     1192                           w -> Index (), 1))) {
     1193
     1194        cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "exprTriLin: variable should be fixed but cut can't be added: ");
     1195        w -> print ();
     1196        cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "\n");
     1197      }
     1198
     1199    } else {
     1200
     1201      // We have either one or two fixed variables
     1202
     1203      switch (n_var_fixed) {
     1204
     1205      case 1: // this is w = f1 * v2 * v3, where one variable is
     1206              // fixed. Revert to exprMul
     1207
     1208        {
     1209          int
     1210            xi = (!isFixed [0]) ? 0 : (!isFixed [1]) ? 1 : 2,
     1211            yi = (!isFixed [2]) ? 2 : (!isFixed [1]) ? 1 : 0;
     1212
     1213          assert (xi != 2);
     1214          assert (yi != 0);
     1215
     1216          xi = varInd [xi];
     1217          yi = varInd [yi];
     1218
     1219          double
     1220            *lb = cg -> Problem () -> Lb (),
     1221            *ub = cg -> Problem () -> Ub (),
     1222            xl = lb [xi],
     1223            xu = ub [xi],
     1224            yl = lb [yi],
     1225            yu = ub [yi],
     1226            wl = lb [varInd [3]],
     1227            wu = ub [varInd [3]];
     1228
     1229          unifiedProdCuts (cg, cs,
     1230                           xi,         (*(arglist_ [0])) (), xl, xu,
     1231                           yi,         (*(arglist_ [1])) (), yl, yu,
     1232                           varInd [3], (*w) (),              wl, wu,
     1233                           chg, sign);
     1234        }
     1235
     1236        break;
     1237
     1238      case 2:
     1239
     1240        { // easy... cut is w = f1 * f2 * v3, where f* are fixed and v* is remaining
     1241
     1242          int varInd = (!isFixed [0]) ? 0 : (!isFixed [1]) ? 1 : 2;
     1243
     1244          double
     1245            lb = ((sign == expression::AUX_LEQ) ? - COIN_DBL_MAX : 0),
     1246            ub = ((sign == expression::AUX_GEQ) ?   COIN_DBL_MAX : 0);
     1247
     1248          if (!(cg->createCut (cs, lb, ub, w -> Index (), 1, w -> ArgList () [varInd] -> Index (), -fixed_prod))) {
     1249            cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "exprTriLin: variable should be fixed but cut can't be added: "); w -> print ();
     1250            cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "\n");
     1251          }
     1252        }
     1253
     1254        break;
     1255
     1256      default:
     1257        cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "exprTriLin: Error, there should be one or two fixed variables");
     1258        break;
     1259      }
     1260    }
     1261
     1262    return;
     1263  }
    11621264
    11631265  TriLinCuts (cg -> Problem () -> Lb (),
     
    11921294    }
    11931295
    1194     std::copy (cutIndices [i].begin (), cutIndices [i].end (), ind);
    1195     std::copy (cutCoeff   [i].begin (), cutCoeff   [i].end (), coe);
    1196 
    1197     OsiRowCut cut (cutLb [i], cutUb [i], 4, 4, ind, coe);
    1198     //cut.print ();
    1199 
    1200     delete [] ind;
    1201     delete [] coe;
    1202 
    1203     if (cg -> Problem () -> bestSol ()) {
    1204 
    1205       // check validity of cuts by verifying they don't cut the
    1206       // optimum in the node
    1207 
    1208       double *sol = cg -> Problem () -> bestSol ();
    1209       const double
    1210         *lb = cg -> Problem () -> Lb (),
    1211         *ub = cg -> Problem () -> Ub ();
    1212 
    1213       int nVars = cg -> Problem () -> nVars ();
    1214 
    1215       bool optIn = true;
    1216 
    1217       for (int i=0; i< nVars; i++)
    1218         if ((sol [i] < lb [i] - COUENNE_EPS) ||
    1219             (sol [i] > ub [i] + COUENNE_EPS)) {
    1220           optIn = false;
    1221           break;
    1222         }
    1223 
    1224       if (optIn) {
    1225 
    1226         for (unsigned int i=0; i<cutIndices.size (); i++) {
    1227 
    1228           double chs = 0.;
    1229 
    1230           for (unsigned int j=0; j<cutIndices[i].size(); j++)
    1231             chs += cutCoeff [i] [j] * sol [cutIndices [i] [j]];
    1232 
    1233           if ((chs < cutLb [i] - COUENNE_EPS) ||
    1234               (chs > cutUb [i] + COUENNE_EPS)) {
    1235 
    1236             printf ("cut %d violates optimum:\n", i);
    1237 
    1238             if (cutLb [i] > -COUENNE_INFINITY) printf ("%g <= ", cutLb [i]);
    1239             for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g x%d ", cutCoeff [i] [j],       cutIndices [i] [j]);  printf ("\n = ");
    1240             for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g *%g ", cutCoeff [i] [j],  sol [cutIndices [i] [j]]); printf ("\n = ");
    1241             for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g ",     cutCoeff [i] [j] * sol [cutIndices [i] [j]]); printf (" = %g", chs);
    1242             if (cutUb [i] <  COUENNE_INFINITY) printf (" <= %g", cutUb [i]);
    1243             printf ("\n");
    1244 
     1296    if ((cutLb [i] > - COUENNE_INFINITY) ||
     1297        (cutUb [i] <   COUENNE_INFINITY)) {
     1298
     1299      std::copy (cutIndices [i].begin (), cutIndices [i].end (), ind);
     1300      std::copy (cutCoeff   [i].begin (), cutCoeff   [i].end (), coe);
     1301
     1302      OsiRowCut cut (cutLb [i], cutUb [i], 4, 4, ind, coe);
     1303      //cut.print ();
     1304
     1305      delete [] ind;
     1306      delete [] coe;
     1307
     1308      if (cg -> Problem () -> bestSol ()) {
     1309
     1310        // check validity of cuts by verifying they don't cut the
     1311        // optimum in the node
     1312
     1313        double *sol = cg -> Problem () -> bestSol ();
     1314        const double
     1315          *lb = cg -> Problem () -> Lb (),
     1316          *ub = cg -> Problem () -> Ub ();
     1317
     1318        int nVars = cg -> Problem () -> nVars ();
     1319
     1320        bool optIn = true;
     1321
     1322        for (int i=0; i< nVars; i++)
     1323          if ((sol [i] < lb [i] - COUENNE_EPS) ||
     1324              (sol [i] > ub [i] + COUENNE_EPS)) {
     1325            optIn = false;
     1326            break;
     1327          }
     1328
     1329        if (optIn) {
     1330
     1331          for (unsigned int i=0; i<cutIndices.size (); i++) {
     1332
     1333            double chs = 0.;
     1334
     1335            for (unsigned int j=0; j<cutIndices[i].size(); j++)
     1336              chs += cutCoeff [i] [j] * sol [cutIndices [i] [j]];
     1337
     1338            if ((chs < cutLb [i] - COUENNE_EPS) ||
     1339                (chs > cutUb [i] + COUENNE_EPS)) {
     1340
     1341              printf ("cut %d violates optimum:\n", i);
     1342
     1343              if (cutLb [i] > -COUENNE_INFINITY) printf ("%g <= ", cutLb [i]);
     1344              for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g x%d ", cutCoeff [i] [j],       cutIndices [i] [j]);  printf ("\n = ");
     1345              for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g *%g ", cutCoeff [i] [j],  sol [cutIndices [i] [j]]); printf ("\n = ");
     1346              for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g ",     cutCoeff [i] [j] * sol [cutIndices [i] [j]]); printf (" = %g", chs);
     1347              if (cutUb [i] <  COUENNE_INFINITY) printf (" <= %g", cutUb [i]);
     1348              printf ("\n");
     1349
     1350            }
    12451351          }
    12461352        }
    12471353      }
    1248     }
    1249 
    1250     cs.insert (cut);
     1354
     1355      cs.insert (cut);
     1356    }
    12511357  }
    12521358  //delete [] varInd;
Note: See TracChangeset for help on using the changeset viewer.