1,120
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to searchBiostatistics 615/815 Fall 2011 (view source)
Revision as of 07:57, 31 March 2011
, 07:57, 31 March 2011→Problem Sets
** See below if you want to copy and paste example code
** See below if you want to copy and paste example code
* Problem Set 5 : Due on March 29th, 2011 -- [[Media:Biostat615-homework5.pdf | (Problem Set 5 - PDF)]]
* Problem Set 5 : Due on March 29th, 2011 -- [[Media:Biostat615-homework5.pdf | (Problem Set 5 - PDF)]]
* Problem Set 5 : Due on March 29th, 2011 -- [[Media:Biostat615-homework6.pdf | (Problem Set 6 - PDF)]]
=== Source code for Homework 6 - simplex615.h ===
#ifndef __SIMPLEX_615_H
#define __SIMPLEX_615_H
#include <vector>
#include <cmath>
#include <iostream>
#define ZEPS 1e-10
class optFunc {
public:
virtual double operator() (std::vector<double>& x) = 0;
};
// Simplex contains (dim+1)*dim points
class simplex615 {
protected:
std::vector<std::vector<double> > X;
std::vector<double> Y;
std::vector<double> midPoint;
std::vector<double> thruLine;
int dim, idxLo, idxHi, idxNextHi;
void evaluateFunction(optFunc& foo);
void evaluateExtremes();
void prepareUpdate();
bool updateSimplex(optFunc& foo, double scale);
void contractSimplex(optFunc& foo);
static int check_tol(double fmax, double fmin, double ftol);
public:
simplex615(double* p, int d);
void amoeba(optFunc& foo, double tol);
std::vector<double>& xmin();
double ymin();
};
simplex615::simplex615(double* p, int d) : dim(d) {
X.resize(dim+1);
Y.resize(dim+1);
midPoint.resize(dim);
thruLine.resize(dim);
for(int i=0; i < dim+1; ++i) {
X[i].resize(dim);
}
// set every point the same
for(int i=0; i < dim+1; ++i) {
for(int j=0; j < dim; ++j) {
X[i][j] = p[j];
}
}
// then increase each dimension by one except for the last point
for(int i=0; i < dim; ++i) {
X[i][i] += 1.;
}
}
void simplex615::evaluateFunction(optFunc& foo) {
for(int i=0; i < dim+1; ++i) {
Y[i] = foo(X[i]);
}
}
void simplex615::evaluateExtremes() {
if ( Y[0] > Y[1] ) {
idxHi = 0;
idxLo = idxNextHi = 1;
}
else {
idxHi = 1;
idxLo = idxNextHi = 0;
}
for(int i=2; i < dim+1; ++i) {
if ( Y[i] <= Y[idxLo] ) {
idxLo = i;
}
else if ( Y[i] > Y[idxHi] ) {
idxNextHi = idxHi;
idxHi = i;
}
else if ( Y[i] > Y[idxNextHi] ) {
idxNextHi = i;
}
}
}
void simplex615::prepareUpdate() {
for(int j=0; j < dim; ++j) {
midPoint[j] = 0;
}
for(int i=0; i < dim+1; ++i) {
if ( i != idxHi ) {
for(int j=0; j < dim; ++j) {
midPoint[j] += X[i][j];
}
}
}
for(int j=0; j < dim; ++j) {
midPoint[j] /= dim;
thruLine[j] = X[idxHi][j] - midPoint[j];
}
}
bool simplex615::updateSimplex(optFunc& foo, double scale) {
std::vector<double> nextPoint;
nextPoint.resize(dim);
for(int i=0; i < dim; ++i) {
nextPoint[i] = midPoint[i] + scale * thruLine[i];
}
double fNext = foo(nextPoint);
if ( fNext < Y[idxHi] ) { // exchange with maximum
for(int i=0; i < dim; ++i) {
X[idxHi][i] = nextPoint[i];
}
Y[idxHi] = fNext;
return true;
}
else {
return false;
}
}
void simplex615::contractSimplex(optFunc& foo) {
for(int i=0; i < dim+1; ++i) {
if ( i != idxLo ) {
for(int j=0; j < dim; ++j) {
X[i][j] = 0.5*( X[idxLo][j] + X[i][j] );
Y[i] = foo(X[i]);
}
}
}
}
void simplex615::amoeba(optFunc& foo, double tol) {
evaluateFunction(foo);
while(true) {
evaluateExtremes();
prepareUpdate();
if ( check_tol(Y[idxHi],Y[idxLo],tol) ) break;
updateSimplex(foo, -1.0); // reflection
if ( Y[idxHi] < Y[idxLo] ) {
updateSimplex(foo, -2.0); // expansion
}
else if ( Y[idxHi] >= Y[idxNextHi] ) {
if ( !updateSimplex(foo, 0.5) ) {
contractSimplex(foo);
}
}
}
}
std::vector<double>& simplex615::xmin() {
return X[idxLo];
}
double simplex615::ymin() {
return Y[idxLo];
}
int simplex615::check_tol(double fmax, double fmin, double ftol) {
double delta = fabs(fmax - fmin); double accuracy = (fabs(fmax) + fabs(fmin)) * ftol;
return (delta < (accuracy + ZEPS));
}
#endif // __SIMPLEX_615_H
=== Source code for Homework 4 - Matrix615.h ===
=== Source code for Homework 4 - Matrix615.h ===
