MOCHA 0.9
|
00001 // Copyright 2009 by Jesus De Loera, David Haws, Jon Lee, Allison O'Hair, University 00002 // of California at Davis, IBM, Inc. All rights reserved. 00003 // 00004 // Distributed under the Eclipse Public License v 1.0. See ../LICENSE 00005 00006 // $Rev$ $Date$ 00007 00009 #ifndef MATHPROG_H 00010 #define MATHPROG_H 00011 00012 #include <iostream> 00013 #include <set> 00014 #include <list> 00015 #include <ctime> 00016 #include "matrix.h" 00017 #include "matroid.h" 00018 00019 using namespace::std; 00020 00021 class Functional { 00022 public: 00023 Functional() {}; 00024 ~Functional() {}; 00025 virtual double Evaluate (Matrix M) = 0; 00026 protected: 00027 }; 00028 00029 class FunctionalGeneral : public Functional { 00030 public: 00031 FunctionalGeneral (double (*tempFct) (Matrix M)); 00032 ~FunctionalGeneral (); 00033 double Evaluate (Matrix M); 00034 protected: 00035 double (*BalFct) (Matrix M) ; 00036 }; 00037 00038 class FunctionalMinVariance : public Functional { 00039 public: 00040 FunctionalMinVariance (Matrix M); 00041 ~FunctionalMinVariance (); 00042 double Evaluate (Matrix M); 00043 protected: 00044 Matrix WeightSum; 00045 }; 00046 00047 class FunctionalQuadratic : public Functional { 00048 public: 00049 FunctionalQuadratic (Matrix M); 00050 ~FunctionalQuadratic (); 00051 double Evaluate (Matrix M); 00052 protected: 00053 Matrix Apex; 00054 }; 00055 00056 class FunctionalLinear : public Functional { 00057 public: 00058 FunctionalLinear (Matrix M); 00059 ~FunctionalLinear (); 00060 protected: 00061 double Evaluate (Matrix M); 00062 Matrix C; 00063 }; 00064 00065 class MathProg { 00066 public: 00067 MathProg (); 00068 MathProg (std::istream &in); 00069 virtual ~MathProg (); 00070 friend std::ostream& operator<< (std::ostream& o, MathProg &someMathProg); 00071 friend std::istream& operator>> (std::istream& in, MathProg &someMathProg); 00072 00073 protected: 00074 virtual void printMathProg (std::ostream &o) { o << "MathProg called\n";}; 00075 virtual void getMathProg (std::istream &in) { }; 00076 }; 00077 00078 class MatroidOpt : public MathProg { 00079 public: 00080 MatroidOpt () {M = 0;}; 00081 protected: 00082 unsigned MatroidType; 00083 Matroid *M; 00084 }; 00085 00086 class ProjMatroidOpt : public MatroidOpt{ 00087 public: 00088 // This will calculate all projected spanning trees 00089 // only if the underlying matroid is a graph. Uses 00090 // Matsui algorithm to find all spanning trees and the 00091 // weighting to project 00092 set <Matrix, ltcolvec> calcAllProjBases (); 00093 protected: 00094 Matrix Weight; 00095 }; 00096 00097 class ProjBalMatroidOpt : public ProjMatroidOpt { 00098 public: 00099 ProjBalMatroidOpt (); 00100 ProjBalMatroidOpt (std::istream &in); 00101 ProjBalMatroidOpt (std::istream &in, double (*tempFct) (Matrix)); 00102 void setBalFct (double (*tempFct) (Matrix)); 00103 00104 // This assumes BalFct is a convex function and we wish to minimize 00105 // Actually, BalFct need not be convex, but if not the algorithms behavior 00106 // is unexpected. 00107 // This implements the steepest decent pivoting algorithm 00108 // Also known as local search 00109 list <set <unsigned> > LocalSearch (set <unsigned> firstBasis); 00110 list <set <unsigned> > LocalSearchRandomStart (); 00111 00112 //list <set <unsigned> > LocalSearch (set <unsigned> firstBasis, double (*tempFct) (Matrix)); 00113 list <set <unsigned> > LocalSearch (set <unsigned> firstBasis, Functional *F); 00114 //list <set <unsigned> > LocalSearchRandomStart (double (*tempFct) (Matrix)); 00115 list <set <unsigned> > LocalSearchRandomStart (Functional *F); 00116 00117 00118 // This assumes BalFct is a convex function and we wish to minimize 00119 // Actually, BalFct need not be convex, but if not the algorithms behavior 00120 // is unexpected. 00121 // This implements the steepest decent pivoting algorithm 00122 // Also known as local search 00123 list <set <unsigned> > FirstComeFirstServe (set <unsigned> firstBasis); 00124 list <set <unsigned> > FirstComeFirstServeRandomStart (); 00125 00126 //list <set <unsigned> > FirstComeFirstServe (set <unsigned> firstBasis, double (*tempFct) (Matrix)); 00127 list <set <unsigned> > FirstComeFirstServe (set <unsigned> firstBasis, Functional *F); 00128 //list <set <unsigned> > FirstComeFirstServeRandomStart (double (*tempFct) (Matrix)); 00129 list <set <unsigned> > FirstComeFirstServeRandomStart (Functional *F); 00130 00131 00132 // This implements the Tabu search method. It always pivots to the minimum 00133 // neighbor under the BalFct. It will also never return to a previously seen basis. 00134 // pivotLimit species the threshold of pivots allowed without improvement 00135 // Returns the bases pivoted to. 00136 list <set <unsigned> > TabuSearchHeuristic (set <unsigned> firstBasis, unsigned pivotLimit, list <set <unsigned> > &tabuBasesList); 00137 list <set <unsigned> > TabuSearchHeuristicRandomStart (unsigned pivotLimit, list <set <unsigned> > &tabuBasesList); 00138 00139 //list <set <unsigned> > TabuSearchHeuristic (set <unsigned> firstBasis, double (*tempFct) (Matrix), unsigned pivotLimit, list <set <unsigned> > &tabuBasesList); 00140 list <set <unsigned> > TabuSearchHeuristic (set <unsigned> firstBasis, Functional *F, unsigned pivotLimit, list <set <unsigned> > &tabuBasesList); 00141 //list <set <unsigned> > TabuSearchHeuristicRandomStart (double (*tempFct) (Matrix), unsigned pivotLimit, list <set <unsigned> > &tabuBasesList); 00142 list <set <unsigned> > TabuSearchHeuristicRandomStart (Functional *F, unsigned pivotLimit, list <set <unsigned> > &tabuBasesList); 00143 00144 // This function takes in a list of subsets and outputs the 00145 // min bases under the functional and weighting for this object 00146 set <unsigned> FindMin (list <set <unsigned> > &); 00147 00148 // This function takes in a list of subsets and outputs the 00149 // min bases under the given functional and weighting 00150 set <unsigned> FindMin (list <set <unsigned> > &, Functional *); 00151 00152 // This function takes in a list of subsets and outputs the 00153 // min bases under the functional for this object 00154 Matrix FindMin (set <Matrix, ltcolvec> &); 00155 00156 // This function takes in a list of subsets and outputs the 00157 // min bases under the given functional 00158 Matrix FindMin (set <Matrix, ltcolvec> &, Functional *); 00159 00160 //list <set <unsigned> > ProjBalMatroidOpt::SimulatedAnnealing(set <unsigned> firstBasis, list < double > temperatures, list < unsigned > times, list <set <unsigned> > &minBases); 00161 list <set <unsigned> > SimulatedAnnealing(set <unsigned> firstBasis, list < double > temperatures, list < unsigned > times, list <set <unsigned> > &minBases); 00162 list <set <unsigned> > SimulatedAnnealing(set <unsigned> firstBasis, list < double >, list < unsigned >, list <set <unsigned> > &, Functional *F); 00163 00164 // This takes in a basis, a temperature T, and a function tempFct. It outputs an adjacent basis 00165 // according to the Boltzmann distribution using the Metropolis Chain method. 00166 //set <unsigned> MetropolisBoltzmannUpdateFunction(set <unsigned>, double T, double (*tempFct) (Matrix)); 00167 set <unsigned> MetropolisBoltzmannUpdateFunction(set <unsigned>, double T, Functional *F); 00168 // This function assumes tempFct is BalFct 00169 set <unsigned> MetropolisBoltzmannUpdateFunction(set <unsigned>, double T); 00170 00171 // This will enumerate a superset of the boundary of the convex hull starting with firstBasis 00172 // Currently only works when the projected dimension is 2 00173 list <set <unsigned> > Boundary(set <unsigned> firstBasis, set <Matrix, ltcolvec> &CH); 00174 // Does the same, except will calculate the firstBasis with a linear program 00175 list <set <unsigned> > Boundary(set <Matrix, ltcolvec> &CH); 00176 00177 00178 // This function takes in pareto Boundary points (in 2d) and will calculate triangular regions 00179 // where other pareto optimum may lie. 00180 set <Matrix, ltcolvec> BoundaryTrianglesTwoDim ( set <Matrix, ltcolvec> &CH); 00181 00182 set <Matrix, ltcolvec> BFSDifferentFiber(set <unsigned> firstBasis); 00183 set <Matrix, ltcolvec> BFSDifferentFiberRandomStart(); 00184 00185 // This will attempt up to numTests using LocalSearch to minimize 00186 // a convex function that zeros only at the weightedSet point. 00187 // Returns 1 if weighted set is a projected point, 0 if can not determine 00188 int PivotTestLocalSearch(Matrix weightedSet, int numTests); 00189 set <Matrix, ltcolvec> PivotTestLocalSearch(set <Matrix, ltcolvec> &, int numTests); 00190 00191 // This will attempt up to numTests using TabuSearch to minimize 00192 // a convex function that zeros only at the weightedSet point. 00193 // Returns 1 if weighted set is a projected point, 0 if can not determine 00194 int PivotTestTabuSearch(Matrix weightedSet, int numTests, int pivotLimit); 00195 set <Matrix, ltcolvec> PivotTestTabuSearch(set <Matrix, ltcolvec> &, int numTests, int pivotLimit); 00196 00197 00198 // This will tests all points in the box defined by lowerCorner and upperCorner 00199 // using PivotTestLocalSearch numTests times. It also will not test points in projPoints. 00200 void BoxPivotTestLocalSearch(const Matrix &lowerCorner, const Matrix &upperCorner, set <Matrix, ltcolvec> &projPoints, int numTests, set <Matrix, ltcolvec> &newPoints); 00201 00202 // This will tests all points in the box defined by lowerCorner and upperCorner 00203 // using PivotTestTabuSearch numTests times, with pivotLimit. It also will not test points in projPoints. 00204 void BoxPivotTestTabuSearch(const Matrix &lowerCorner, const Matrix &upperCorner, set <Matrix, ltcolvec> &projPoints, int numTests, set <Matrix, ltcolvec> &newPoints, int pivotLimit); 00205 00206 // This will use local search to find upper and lower bounds for each dimension. 00207 // Then it will call BoxPivotTestLocalSearch 00208 // using PivotTestLocalSearch numTests times. It also will not test points in projPoints. 00209 void AutoBoundsPivotTestLocalSearch( set <Matrix, ltcolvec> &projPoints, int numTests, set <Matrix, ltcolvec> &newPoints); 00210 00211 // This will use local search to find upper and lower bounds for each dimension. 00212 // Then it will call BoxPivotTestTabuSearch 00213 // using PivotTestTabuSearch numTests times, with pivotLimit. It also will not test points in projPoints. 00214 void AutoBoundsPivotTestTabuSearch( set <Matrix, ltcolvec> &projPoints, int numTests, set <Matrix, ltcolvec> &newPoints, int pivotLimit); 00215 00216 // This will perform numSearches many BFS random starts 00217 // It will always try to find a new random basis to start mod projection 00218 // from. If it can not find a new random basis before newRandTolerance attempts 00219 // it returns what it currently has calculated. -1 =: infinity 00220 // findAllBoundary == 1 means do not find random boundary elements but 00221 // compute the entire boundary and run bfs on each point 00222 void MultiBFSRandomStarts(int numSearches, int BFSSearchDepth, int newRandToleranceBoundary, int findAllBoundary, int newRandTolerance, set <Matrix, ltcolvec> &); 00223 Matrix projectSet(set <unsigned>); 00224 00225 // Allisons work. 00226 // This takes in points in R^2 and outputs the Pareto Optimum. 00227 set <Matrix, ltcolvec> ParetoOptimum( set <Matrix, ltcolvec> &); 00228 00229 // This takes in points in R^2 and outputs the minmax Optimum. 00230 set <Matrix, ltcolvec> MinMax(set <Matrix, ltcolvec> &); 00231 00232 // This will create a random vector and pivot to a terminal basis and return it. 00233 // By basis we mean a projected basis. 00234 set <unsigned> randomLinearBasis(); 00235 00236 static void printPivots(const list <set <unsigned> > &); 00237 static void printBFSList (const set <Matrix, ltcolvec> &); 00238 static void printBFSListMatlab (const set <Matrix, ltcolvec> &,string); 00239 void printPivotsMin(const list <set <unsigned> > &); 00240 //void printPivotsMin(const list <set <unsigned> > &,double (*tempFct) (Matrix)); 00241 void printPivotsMin(const list <set <unsigned> > &, Functional *F); 00242 void printPivotsMatlab(const list <set <unsigned> > &,string); 00243 00244 static void writePivots(const list <set <unsigned> > &, std::ostream &); 00245 static void writeBFSList (const set <Matrix, ltcolvec> &, std::ostream &); 00246 static void writeBFSListMatlab (const set <Matrix, ltcolvec> &, std::ostream &,string); 00247 void writePivotsMin(const list <set <unsigned> > &, std::ostream &); 00248 //void writePivotsMin(const list <set <unsigned> > &, std::ostream &,double (*tempFct) (Matrix)); 00249 void writePivotsMin(const list <set <unsigned> > &, std::ostream &, Functional *F); 00250 void writePivotsMatlab(const list <set <unsigned> > &, std::ostream &,string); 00251 00252 double evalFct(set <unsigned>); 00253 static unsigned BFSLevel; 00254 static int BFSTerminateLevel; // -1 signifies infinity 00255 static unsigned BFSPrinted; 00256 static unsigned pivotsPrinted; 00257 unsigned projDim (); //Returns the number of rows of the weighting matrix 00258 set <unsigned> randomBasis (); 00259 protected: 00260 unsigned totalBoxTests; 00261 time_t BoxTestsModTime; 00262 unsigned BoxTestsModValue; 00263 unsigned BoxTests; 00264 // This is the function that optimizes. 00265 // It takes in a list since we do not know the proper dimensions are run time 00266 //void BFSDifferentFiberInternal(set <unsigned>,set <Matrix, ltcolvec> &,set <Matrix, ltcolvec> &); 00267 void BFSDifferentFiberInternal(set <unsigned>,set <Matrix, ltcolvec> &); 00268 00269 // Internal recursive function to go through lowerCorner - upperCorner points and test 00270 void BoxPivotTestLocalSearchRec(int colIndex, Matrix ¤tMatrix, const Matrix &lowerCorner, const Matrix &upperCorner, set <Matrix, ltcolvec> &projPoints, set <Matrix, ltcolvec> &newPoints); 00271 // Internal recursive function to go through lowerCorner - upperCorner points and test 00272 void BoxPivotTestTabuSearchRec(int colIndex, Matrix ¤tMatrix, const Matrix &lowerCorner, const Matrix &upperCorner, set <Matrix, ltcolvec> &projPoints, set <Matrix, ltcolvec> &newPoints, int pivotLimit); 00273 double (*BalFct) (Matrix M) ; 00274 Functional *BalanceFunction; 00275 virtual void printMathProg (std::ostream &o); 00276 void getMathProg (std::istream &in); 00277 }; 00278 00279 class MinVarianceBalClustering : public ProjBalMatroidOpt { 00280 public: 00281 MinVarianceBalClustering (); 00282 // This constructor reads in the data points. There is an implied uniform matroid structure. 00283 // If n points are read in, there is a rank n/2 uniform matroid on n elements. 00284 MinVarianceBalClustering (std::istream &in); 00285 00286 // This takes in a subset of columns of Weights (the data points) 00287 // and returns two matrices; one that is a submatrix of weights using S (M1) and 00288 // the remaining (M2) 00289 //void getClusters(set <unsigned> S, Matrix &M1, Matrix &M2); 00290 00291 // S is the basis that represents half of the points. 00292 // label is the matlab label to use. It will use label1 and label2 00293 void writeClustersMatlab (set <unsigned> S, std::ostream &,string label); 00294 protected: 00295 virtual void printMathProg(std::ostream &o); 00296 }; 00297 00298 00299 class TwoDSetAngle 00300 { 00301 public: 00302 set <unsigned> S; //This holds a set 00303 double angle; // This is the angle this set with some fixed set makes with the origin 00304 unsigned dely; 00305 unsigned delx; 00306 }; 00307 00308 struct ltTwoDSetAngle { 00309 bool operator()(const TwoDSetAngle &S1, const TwoDSetAngle &S2) const 00310 { 00311 if (S1.angle < S2.angle) 00312 { 00313 return 1; 00314 } 00315 return 0; 00316 } 00317 00318 }; 00319 00320 class TwoDAngleOrderedList 00321 { 00322 public: 00323 //TwoDAngleOrderedList(); 00324 TwoDAngleOrderedList(set <unsigned> initialSet, ProjBalMatroidOpt *PBMO); 00325 ~TwoDAngleOrderedList(); 00326 // Inserts S into the list using refPBMO for angle calculation 00327 void insert(set <unsigned> S); 00328 // Returns the two sets with the largest angle between them and the angle 00329 // Returns 0 if there are less than 2 sets in the list. 00330 double largestAngle (set <unsigned> &S1, set <unsigned> &S2); 00331 00332 void print (); 00333 protected: 00334 ProjBalMatroidOpt *refPBMO; 00335 double x,y; //x and y values for the initial set 00336 set < TwoDSetAngle, ltTwoDSetAngle> orderedRays; 00337 }; 00338 00339 #endif