/home/coin/SVN-release/OS-2.4.0/Couenne/src/cut/sdpcuts/linquad_cuts.cpp

Go to the documentation of this file.
00001 /* $Id: linquad_cuts.cpp 508 2011-02-15 21:52:44Z pbelotti $
00002  *
00003  * Name:    linquad_cuts.cpp
00004  * Author:  Andrea Qualizza
00005  * Purpose: 
00006  *
00007  * This file is licensed under the Eclipse Public License (EPL)
00008  */
00009 
00010 #include <linquad_cuts.hpp>
00011 #include <stdio.h>
00012 #include <math.h>
00013 
00014 #include <CglCutGenerator.hpp>
00015 #include <tracer.hpp>
00016 #include <misc_util.hpp>
00017 
00018 // finds the value xk. (xk,xk^2) is the point on the parabola y=x^2
00019 // closest to the point (xc,yc)
00020 // 
00021 // xk is the zero of the function 2 x^3 - 2 yc x + x - xc = 0
00022 
00023 
00024 double f_  (double x) {return x*x;}
00025 double fp_ (double x) {return 2*x;}
00026 double fpp_(double x) {return 2;}
00027 
00028 double powNewton(double xc, double yc, double (*f)(double),double (*fp)(double),double (*fpp)(double)) {
00029         // Find a zero to the function
00030         //
00031         // F(x) = x - xc + f'(x) (f(x) - yc)
00032         //
00033         // where f(x) is either x^k, exp(x), or log(x).
00034         // The derivative of F(x) is
00035         //
00036         // F'(x) = 1 + f''(x) (f(x) - yc) + (f'(x))^2
00037         //
00038         // Apply usual update:
00039         //
00040         // x(k+1) = x(k) - F(x(k))/F'(x(k))
00041         double xk       = xc;
00042         double fk       = f(xk);
00043         double fpk      = fp(xk);
00044         double F        = fpk * fk;
00045         double Fp       = 1+ fpp(xk) * fk + fpk * fpk;
00046         
00047 
00048         for (int k = NEWTON_MAX_ITER; k--;) {
00049                 xk -= F / Fp;
00050                 
00051                 fk  = f(xk) - yc;
00052                 fpk = fp(xk);
00053                 F   = xk - xc + fpk * fk;
00054                 
00055                 //    printf ("xk = %g; F = %g, fk = %g, fpk = %g\n", xk, F, fk, fpk);
00056                 
00057                 if (fabs(F) < NEWTON_POW_TOLERANCE) 
00058                         break;
00059                 Fp  = 1 + fpp(xk) * fk + fpk * fpk;
00060         }
00061         
00062         return xk;
00063 }
00064 
00065 
00066 void linQuadCutGen(const double *sol, int n, OsiCuts &cs,Tracer *tracer) {
00067         // quadratic diagonal cuts
00068         //      deepest cut: use tangent to quadratic cut passing throgh the 
00069         //      closest point to sol on the quadratic function
00070         //      
00071         //
00072         int origcuts = cs.sizeRowCuts();
00073         Timer linquad_timer;
00074         linquad_timer.start();
00075         for(int i=0;i<n;i++) {
00076                 double xc,yc;
00077 
00078                 //tangent line of region X_ii >= x_i^2 through the closest point 
00079                 //to (sol[x_i],sol[indexQ(i,i,n)]. (Here Xii is our y)
00080                 xc = powNewton(sol[i],sol[indexQ(i,i,n)],&f_,&fp_,&fpp_);
00081                 yc = f_(xc);
00082                 generateTangentDiagonalEntryCut(n,i,cs,xc,yc,sol,true);
00083 
00084                 //tangent line of region X_ii >= x_i^2 through the point (sol[x_i],sol[x_i]*sol[x_i])
00085                 xc = sol[i];
00086                 yc = f_(xc);
00087                 generateTangentDiagonalEntryCut(n,i,cs,xc,yc,sol,true);
00088 
00089                 //tangent line of region X_ii >= x_i^2 through the point
00090                 //(sqrt(sol[indexQ(i,i,n)]),sol[x_i]*sol[x_i])
00091                 if(sol[indexQ(i,i,n)] > 0) {
00092                         xc = sqrt(sol[indexQ(i,i,n)]);
00093                         yc = sol[indexQ(i,i,n)];
00094                         generateTangentDiagonalEntryCut(n,i,cs,xc,yc,sol,true);
00095                 }
00096         }
00097         tracer->setLinquadTime(linquad_timer.time());
00098         tracer->setLinquadTotalCuts(cs.sizeRowCuts() - origcuts);
00099 }
00100 
00101 
00102 
00103 void linQuadCutGenOriginalBounds(const double *xlb, const double *xub, int n, OsiCuts &cs,Tracer *tracer) {
00104 #define LINQUAD_BOUNDS_CUTS_PARTS       3
00105         int origcuts = cs.sizeRowCuts();
00106         Timer linquad_timer;
00107         linquad_timer.start();
00108         for (int i=0;i<n;i++) {
00109                 double xc,yc;
00110                 double interval = xub[i] - xlb[i];
00111                 double step = interval / LINQUAD_BOUNDS_CUTS_PARTS;
00112                 int cnt=0;
00113                 xc = xub[i];
00114                 while (cnt<=LINQUAD_BOUNDS_CUTS_PARTS) {
00115                         yc = f_(xc);
00116                         generateTangentDiagonalEntryCut(n,i,cs,xc,yc,NULL,false);       
00117                         cnt++;
00118                         xc += step;
00119                 }       
00120         }
00121         tracer->setLinquadTime(linquad_timer.time());
00122         tracer->setLinquadTotalCuts(cs.sizeRowCuts() - origcuts);
00123 }
00124 
00125 
00126 void generateTangentDiagonalEntryCut(int n,int i,OsiCuts &cs,double xc,double yc,const double* sol, bool ifViolated) {
00127 // (xc,yc) is a point that satisfies X_ii = x_i^2
00128 // now we compute the tangent line of the region at this point which will be
00129 // (X_ii - yc) = 2 xc (x_i - xc) -->  - 2*xc*x_i + X_ii = yc - 2*xc^2
00130 // the cut is simply : - 2 x_i + X_ii >= yc - 2*xc^2
00131 // sol can be NULL if ifViolated=false
00132                 double *coeff;
00133                 int *ind;
00134                 double rhs;
00135                 OsiRowCut *cut;
00136 
00137                 cut   = new OsiRowCut;
00138                 coeff = new double [2];
00139                 ind   = new int    [2];
00140                 coeff[0]        = -2.0 *xc;
00141                 coeff[1]        = 1.0 ;
00142                 ind[0]          = i ;
00143                 ind[1]          = indexQ(i,i,n);
00144                 rhs     = yc - 2 * xc * xc;
00145                 cut -> setRow (2, ind, coeff);
00146                 cut -> setLb (rhs);
00147 
00148                 if ( !( ifViolated ) || ( cut->violated(sol) <= 0) ) {
00149                         cs.insert (cut);
00150                 }
00151 
00152                 delete cut;
00153                 delete [] coeff;
00154                 delete [] ind;
00155 }
00156 

Generated on Thu Sep 22 03:05:57 2011 by  doxygen 1.4.7