Skip to content

Instantly share code, notes, and snippets.

/units.cpp

Created Apr 21, 2014
Embed
What would you like to do?
template<int cols>
class Matrix {
typedef double Row[cols];
Row*const mat;
const int r;
Matrix(const Matrix&); // need copy constructor, but I'm cheating
public:
int rows() const { return r; }
Matrix(int rows):mat(new Row[rows]), r(rows) {}
double& operator()(int row, int col) {
if(row < 0 || row >= r || col <0 || col >= cols) throw "bug!";
return mat[row][col];
}
};
template<int rows>
struct TransposedMatrix {
int cols() const { return mat.rows(); }
TransposedMatrix(Matrix<rows>& m) : mat(m) {}
double operator()(int row, int col) { return mat(col,row); }
private:
Matrix<rows>& mat;
};
template<int n>
TransposedMatrix<n> transpose(Matrix<n>& m) { return m; }
template<int rows, int cols>
struct SmallMatrix {
double mat[rows][cols];
double& operator()(int row, int col) {
if(row < 0 || row >= rows || col <0 || col >= cols) throw "bug!";
return mat[row][col];
}
SmallMatrix(){
for(int i=0; i<rows; i++)
for(int j=0; j<cols; j++)
mat[i][j] = 0;
}
};
template<int rows, int cols>
SmallMatrix<rows,cols> operator*(TransposedMatrix<rows>& a, Matrix<cols>& b) {
if(a.cols() != b.rows()) throw "bug!";
return matmul<rows,cols>(a, b, a.cols());
}
template<int rows, int cols, int n>
SmallMatrix<rows,cols> operator*(SmallMatrix<rows, n>& a, SmallMatrix<n, cols>& b) {
return matmul<rows,cols>(a, b, n);
}
template<int rows, int cols, typename M1, typename M2>
SmallMatrix<rows,cols> matmul(M1& a, M2& b, int n){
SmallMatrix<rows, cols> c;
for(int i=0; i<rows; i++)
for(int j=0; j<cols; j++)
for(int k=0; k<n; k++)
c(i,j) += a(i,k) * b(k,j);
return c;
}
template<int rows, int cols> // Gaussian elemination alg simplified from Wikipedia, with different variable names
void row_reduce(SmallMatrix<rows, cols>& mat) {
int n = rows < cols ? rows : cols;
for(int i=0; i<n; i++) {
int p = i;
while(p<n && mat(p,i)==0) p++;
for(int j=0; j<cols; j++){
double temp = mat(i,j);
mat(i,j) = mat(p,j);
mat(p,j) = temp;
}
double scale = mat(i,i);
for(int j=i; j<cols; j++)
mat(i,j) /= scale;
for(int j=0; j<rows; j++) {
if(j==i) continue;
double factor = mat(j,i);
for(int k=i; k<cols; k++)
mat(j,k) -= mat(i,k) * factor;
}
}
}
template<int n>
SmallMatrix<n, n> inverse(SmallMatrix<n, n> mat) {
SmallMatrix<n, 2*n> temp;
for(int i=0; i<n; i++) // glue on the identity matrix
for(int j=0; j<n; j++)
temp(i,j) = mat(i,j),
temp(i,j+n) = i == j;
row_reduce(temp);
SmallMatrix<n,n> ans;
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
ans(i,j) = temp(i,j+n);
return ans;
}
struct Line { double slope, intercept; };
Line least_squares_fit(const double x[], const double y[], int n) {
Matrix<1> Y(n);
for(int i=0; i<n; i++) Y(i,0) = y[i];
Matrix<2> X(n);
for(int i=0; i<n; i++)
X(i,0) = x[i],
X(i,1) = 1;
TransposedMatrix<2> Xt = transpose(X);
SmallMatrix<2,1> ans = inverse(Xt*X)*(Xt*Y);
Line ret;
ret.slope = ans(0,0);
ret.intercept = ans(1,0);
return ret;
}
#include<cstdio>
int main() {
double x[] = {1,2,3}, y[] = {7,8,9}; // input
Line line = least_squares_fit(x, y, 3);
printf("y=%gx+%g\n", line.slope, line.intercept);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment