# Changeset 520

Ignore:
Timestamp:
Apr 27, 2007 3:39:58 PM (13 years ago)
Message:

limit power upper/lower bounds to get reasonable coefficients

Location:
branches/Couenne/Couenne/src
Files:
17 edited

### Legend:

Unmodified
Removed

 r506 CouNumber x, CouNumber l, CouNumber u, bool is_global) const { CouNumber opp_slope = - fprime (x);
• ## branches/Couenne/Couenne/src/convex/generateCuts.cpp

 r515 for (int i=0, j = problem_ -> nAuxs (); j--; i++) { expression * image = problem_ -> Aux (i) -> Image (); // TODO: check if list contains all and only aux's to cut /*expression * image = problem_ -> Aux (i) -> Image (); if (   (image -> Linearity () > LINEAR)          // if linear, no need to cut twice && (image -> dependsOn (changed, nchanged))  // if expression does not depend on )                                            // changed variables, do not cut )*/                                            // changed variables, do not cut problem_ -> Aux (i) -> generateCuts (si, cs, this); }
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprExp.cpp

 r480 #include #include #include #include argument_ -> getBounds (le, ue); CouNumber x = (*argument_) (), CouNumber x = (cg -> isFirst ()) ? 0 : powNewton ((*argument_) (), (*aux) (), exp, exp, exp), l = (*le) (), u = (*ue) (); int ns = cg -> nSamples (); // change bounds to get finite coefficients //  CouNumber fact = 2 * ns; /*if (x > 0) { if (l < log (COUENNE_EPS))      l = x / fact - 1; if (u > log (COUENNE_INFINITY)) u = x * fact + 1; } else { if (l < log (COUENNE_EPS))      l = x * fact - 1; if (u > log (COUENNE_INFINITY)) u = x / fact + 1; }*/ // add tangents with finite coefficients if (l < log (COUENNE_EPS))      l = x - ns;
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprGroup.cpp

 r460 OsiRowCut *cut   = new OsiRowCut; CouNumber rhs = c0_; CouNumber rhs = - c0_; // first, make room for aux variable if (curr -> Type () == CONST) // constant term in sum rhs += curr -> Value (); rhs -= curr -> Value (); else {                        // variable coeff [++nterms] = 1.; cut -> setRow (nterms+1, index, coeff); rhs = - rhs; cut -> setUb (rhs);
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprInv.cpp

 r476 // derivative of inv (x) inline CouNumber oppInvSqr (register CouNumber x) {return (- inv (x*x));} #define MIN_DENOMINATOR 1e-10 aux -> getBounds (wle, wue); expression *xe = argument_; CouNumber x = (*xe)  (); CouNumber x; int w_ind = aux       -> Index (), } int sign = (l > 0) ? +1 : -1; // bound cg -> addEnvelope (cs, sign, inv, oppInvSqr, w_ind, x_ind, x, l, u); cg -> addEnvelope (cs, (l > 0) ? +1 : -1, inv, oppInvSqr, w_ind, x_ind, (cg -> isFirst ()) ? // place it somewhere in the interval (we won't care) ((l > COUENNE_EPS) ? l : u) : // not first call, gotta replace it where it gives deepest cut powNewton ((*argument_) (), (*aux) (), inv, oppInvSqr, inv_dblprime), l, u); delete xle;
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprLog.cpp

 r480 #include #include #include #include #include argument_ -> getBounds (le, ue); CouNumber x = (*argument_) (), CouNumber x = (cg -> isFirst ()) ? 1 : powNewton ((*argument_) (), (*aux) (), log, inv, oppInvSqr), l = (*le) (), u = (*ue) ();
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprPow-envelope.cpp

 r506 // orthogonal with line through current- and tangent point) //  x = powNewton (x, y, power_k, power_k_prime, power_k_dblprime); if (!(cg -> isFirst ())) x = powNewton (x, y, power_k, power_k_prime, power_k_dblprime); if      (x nSamples ())); CouNumber step     = 1 + log (1. + (double) (cg -> nSamples ())), powThres = pow (COU_MAX_COEFF, 1./k); if (l < - COUENNE_INFINITY + 1) { if (l < - powThres + 1) { l = x - step; if (u > COUENNE_INFINITY - 1) if (u > powThres - 1) u = x + step; } else if (u > COUENNE_INFINITY - 1) if (u > powThres - 1) u = x + step;
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprPow-getBounds.cpp

 r476 int rndexp; bool isInt    =  fabs (expon - (rndexp = FELINE_round (expon))) < COUENNE_EPS; bool isInt    =  fabs (expon - (rndexp = COUENNE_round (expon))) < COUENNE_EPS; bool isInvInt = !isInt && ((fabs (expon) > COUENNE_EPS) && (fabs (1/expon - (rndexp = FELINE_round (1/expon))) < COUENNE_EPS)); (fabs (1/expon - (rndexp = COUENNE_round (1/expon))) < COUENNE_EPS)); if ((isInt || isInvInt) && (rndexp % 2) && (rndexp > 0)) {
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprPow.cpp

 r507 int x_ind = xe  -> Index (); CouNumber w = (*aux) (), x = (*xe)  (), CouNumber w, x, l = (*xle) (), u = (*xue) (); } bool isInt    =            fabs (k    - (double) (intk = FELINE_round (k)))    < COUENNE_EPS, isInvInt = !isInt && (fabs (1./k - (double) (intk = FELINE_round (1./k))) < COUENNE_EPS); bool isInt    =            fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS, isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS); // two macro-cases: if (isInvInt) { if (cg -> isFirst ()) { w = (l>0) ? 1 : (u<0) ? -1 : 0; x = 0; } q = safe_pow (q, k); sign = -1; } else sign = 1; else { if (cg -> isFirst ()) { x = (l>0) ? l : (u<0) ? u : 0; w = 0; } sign = 1; } CouNumber powThres = pow (COU_MAX_COEFF, 1./k); // don't want big coefficients // lower envelope if (l > -COUENNE_INFINITY) { if (l > -powThres) { if (u > q * l) { // upper x is after "turning point", add lower envelope addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, q*l, u, sign); // upper envelope if (u < COUENNE_INFINITY) { if (u < powThres) { if (l < q * u) { // lower x is before "turning point", add upper envelope addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, q*u, -sign); && (!isInvInt || !(intk % 2)) && (l < - COUENNE_EPS) && (u < (l=0)))        // CAUTION! l is updated here, if negative && (u < (l=0)))        // CAUTION! l updated if negative (k real) return; // upper envelope -- when k negative, add only if CouNumber powThreshold = pow (1e20, 1./k); CouNumber powThres = pow (COU_MAX_COEFF, 1./k), // don't want big coefficients powStep  = 1; if ((  (k > COUENNE_EPS) || (l > COUENNE_EPS)             // bounds do not contain 0 || (l > COUENNE_EPS)       // bounds do not contain 0 || (u < - COUENNE_EPS)) && (l > - powThreshold) &&  // and are finite (u <   powThreshold)) (l > - powThres) &&    // and are finite (u <   powThres)) cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), -sign); // similarly, pay attention not to add infinite slopes #define POWER_RANGE 1e2; if (cg -> isFirst()) { x = (k<0) ? ((u<0) ? u : (l>0) ? l : 0) : 0; w = 0; } if (k > COUENNE_EPS) { if (u >   powThreshold) u = x + POWER_RANGE; if (l < - powThreshold) l = x - POWER_RANGE; if (u >   powThres) u = x + powStep; if (l < - powThres) l = x - powStep; } else { if (fabs (l) < COUENNE_EPS) l =  1. / POWER_RANGE; // l --> 0+ if (fabs (u) < COUENNE_EPS) u = -1. / POWER_RANGE; // u --> 0- if (fabs (l) < COUENNE_EPS) l =  1. / powThres; // l --> 0+ if (fabs (u) < COUENNE_EPS) u = -1. / powThres; // u --> 0- }
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprSub.cpp

 r476 expression *y = arglist_ [1]; int wind = w -> Index (); int xind = x -> Index (); int yind = y -> Index (); int wi = w -> Index (); int xi = x -> Index (); int yi = y -> Index (); if (x->Type () == CONST) { // (c - y) or (c - d) CouNumber c = x -> Value (); if (y->Type() == CONST) cg -> createCut (cs, c - y->Value(), 0, wind, 1.,   -1, 0., -1, 0., true); else                    cg -> createCut (cs, c,              0, wind, 1., yind, 1., -1, 0., true); } if (x->Type () == CONST) // (c - y) or (c - d) if (y->Type() == CONST) cg -> createCut (cs, x->Value()-y->Value(), 0, wi, 1, -1, 0, -1, 0, true); else                    cg -> createCut (cs, x->Value(),            0, wi, 1, yi, 1, -1, 0, true); else // (x - y) or (x - d) if (y->Type() == CONST) cg -> createCut (cs, y->Value(),     0, wind, -1., xind, 1., -1, 0.,true); else                    cg -> createCut (cs, 0.,          0, wind, -1., xind, 1., yind, -1.,true); if (y->Type() == CONST) cg -> createCut (cs, y->Value(),     0, wi, -1., xi, 1., -1, 0.,true); else                    cg -> createCut (cs, 0.,          0, wi, -1., xi, 1., yi, -1.,true); }
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprSum.cpp

 r450 // generate convexification cut for constraint w = this /// generate convexification cut for constraint w = this void exprSum::generateCuts (exprAux *w, const OsiSolverInterface &si, CouNumber rhs = 0; // first, make room for aux variable /// first, make room for aux variable coeff [0] = -1.; index [0] = w -> Index (); int nv = 1; // scan arglist for (aux) variables and constants /// scan arglist for (aux) variables and constants for (int i=0; i Type () == CONST) rhs += arglist_ [i] -> Value (); rhs -= arglist_ [i] -> Value (); else { coeff [nv]   = 1.; cut -> setRow (nv, index, coeff); cut -> setUb (-rhs); cut -> setLb (-rhs); cut -> setUb (rhs); cut -> setLb (rhs); delete [] index; delete [] coeff; // added only once, it is global /// added only once, it is global cut -> setGloballyValid ();
• ## branches/Couenne/Couenne/src/expression/exprAux.cpp

 r515 print (std::cout);  printf (" := "); image_ -> print (std::cout); printf (" [%.3e,%.3e] ---> ", expression::Lbound (varIndex_), expression::Ubound (varIndex_)); int index; if ((image_ -> Argument ()) && ((index = image_ -> Argument () -> Index ()) >= 0)) printf ("[%.3e,%.3e] ", expression::Lbound (index), expression::Ubound (index)); else if (image_ -> ArgList ()) for (int i=0; i nArgs (); i++) if ((index = image_ -> ArgList () [i] -> Index ()) >= 0) printf ("[%.3e,%.3e] ", expression::Lbound (index), expression::Ubound (index)); printf("\n"); for (int jj=j; jj < cs.sizeRowCuts ();jj++)
• ## branches/Couenne/Couenne/src/expression/operators/exprPow.cpp

 r498 double exponent = arglist_ [1] -> Value (); if (fabs (exponent - FELINE_round (exponent)) > COUENNE_EPS) if (fabs (exponent - COUENNE_round (exponent)) > COUENNE_EPS) return NONLINEAR; if (arglist_ [1] -> Type () == CONST) { int expInt = (int) FELINE_round (exponent); int expInt = (int) COUENNE_round (exponent); if (arglist_ [0] -> Linearity () == LINEAR) { int intk; // integer (or integer inverse of) exponent bool isint    =           (k    - (intk = FELINE_round (k))    < COUENNE_EPS), // k   is integer isinvint = !isint && (1./k - (intk = FELINE_round (1./k)) < COUENNE_EPS); // 1/k is integer bool isint    =           (k    - (intk = COUENNE_round (k))    < COUENNE_EPS), // k   is integer isinvint = !isint && (1./k - (intk = COUENNE_round (1./k)) < COUENNE_EPS); // 1/k is integer CouNumber wl = l [wind], // lower w
• ## branches/Couenne/Couenne/src/include/CouennePrecisions.h

 r427 #include #define COUENNE_EPS       1e-7 /* keep it at least 1e-7 or strange things happen */ /* keep it at least 1e-7 or strange things happen */ #define COUENNE_EPS       1e-7 #define COUENNE_EPS_SIMPL 1e-20 #define COUENNE_INFINITY  1e+50 #define COU_MAX_COEFF     1e6 #define FELINE_round(x) ((int) (floor ((x)+0.5))) #define COUENNE_round(x) ((int) (floor ((x)+0.5))) #endif
• ## branches/Couenne/Couenne/src/include/operators/exprInv.h

 r395 // the operator itself /// the operator itself inline CouNumber inv (register CouNumber arg) {return 1.0 / arg;} /// derivative of inv (x) inline CouNumber oppInvSqr (register CouNumber x) {return (- inv (x*x));} // class inverse (1/f(x)) /// inv_dblprime, second derivative of inv (x) inline CouNumber inv_dblprime (register CouNumber x) {return (2 * inv (x*x*x));} /// class inverse (1/f(x)) class exprInv: public exprUnary { public: // Constructors, destructor /// Constructors, destructor exprInv  (expression *al): exprUnary (al) {} //< non-leaf expression, with argument list // cloning method /// cloning method expression *clone () const {return new exprInv (argument_ -> clone ());} inline unary_function F () {return inv;} // output "1/argument" /// output "1/argument" void print (std::ostream&) const; // differentiation /// differentiation expression *differentiate (int index); // get a measure of "how linear" the expression is (see CouenneTypes.h) /// get a measure of "how linear" the expression is (see CouenneTypes.h) virtual inline int Linearity () { if (argument_ -> Type () == CONST) return CONSTANT; } // Get lower and upper bound of an expression (if any) /// Get lower and upper bound of an expression (if any) void getBounds (expression *&, expression *&); // generate equality between *this and *w /// generate equality between *this and *w void generateCuts (exprAux *w, const OsiSolverInterface &si, OsiCuts &cs, const CouenneCutGenerator *cg); /// /// code for comparisons virtual enum expr_type code () {return COU_EXPRINV;}
• ## branches/Couenne/Couenne/src/include/operators/exprPow.h

 r506 if ((base < 0) && ((fabs (exponent - (rndexp = FELINE_round (exponent))) < COUENNE_EPS) || ((fabs (exponent - (rndexp = COUENNE_round (exponent))) < COUENNE_EPS) || ((fabs (exponent) > COUENNE_EPS) && (fabs (1. / exponent - (rndexp = FELINE_round (1. / exponent))) < COUENNE_EPS))) (fabs (1. / exponent - (rndexp = COUENNE_round (1. / exponent))) < COUENNE_EPS))) && (rndexp % 2)) return (- pow (- base, exponent));