From e4848afb9e6f5bf28579cc325a2f39b2d9594c01 Mon Sep 17 00:00:00 2001
From: Mark Woolrich <woolrich@fmrib.ox.ac.uk>
Date: Wed, 14 Apr 2004 14:54:54 +0000
Subject: [PATCH] *** empty log message ***
---
Makefile | 4 ++--
minimize.cc | 4 +---
minimize.h | 43 +++++++++++++++++++++++++++----------------
miscmaths.cc | 2 +-
miscmaths.h | 2 +-
quick.cc | 26 ++++++++++++++++++++++++--
sparse_matrix.cc | 47 +++++++++++++++++++++++++++++++++++++++++++++++
sparse_matrix.h | 2 ++
8 files changed, 105 insertions(+), 25 deletions(-)
diff --git a/Makefile b/Makefile
index b2fb9e0..73aa54b 100644
--- a/Makefile
+++ b/Makefile
@@ -18,13 +18,13 @@ LIBS = -lutils -lnewmat -lprob -lm
all: libmiscmaths.a
quick:${OBJS} quick.o
- ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} quick.o ${LIBS}
+ ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} quick.o ${LIBS} -lavwio
libmiscmaths.a: ${OBJS}
${AR} -r libmiscmaths.a ${OBJS}
testprog: testprog.o
- ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} testprog.o ${LIBS}
+ ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} testprog.o ${LIBS} -lavwio
diff --git a/minimize.cc b/minimize.cc
index eb199a8..7380983 100644
--- a/minimize.cc
+++ b/minimize.cc
@@ -24,8 +24,6 @@ using namespace NEWMAT;
using namespace std;
///////////////////////////////////////////////////////
-//fminsearch.m
-
namespace MISCMATHS {
float diff1(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord)
@@ -193,7 +191,7 @@ ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func, float h,in
}
-void fminsearch(ColumnVector& x, const EvalFunction& func){
+void minsearch(ColumnVector& x, const EvalFunction& func){
//perform generic function minimization without gradient info
int n=x.Nrows(), maxiter=200*n,iter=0;
int func_evals=0;
diff --git a/minimize.h b/minimize.h
index 39de572..8716aef 100644
--- a/minimize.h
+++ b/minimize.h
@@ -41,6 +41,23 @@ public:
}
};
+ class EvalFunction;
+ class gEvalFunction;
+
+float diff1(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff derivative
+
+float diff2(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff 2nd derivative
+
+float diff2(const ColumnVector& x, const EvalFunction& func, int i,int j,float h,int errorord=4);// finite diff cross derivative
+
+ReturnMatrix gradient(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff derivative vector
+
+ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff hessian
+
+void minsearch(ColumnVector& x, const EvalFunction& func);
+
+void scg(ColumnVector& x, const gEvalFunction& func, float tol = 0.0000001, float eps=1e-16, int niters=100);
+
class EvalFunction
{//Function where gradient is not analytic (or you are too lazy to work it out) (required for fminsearch)
public:
@@ -48,6 +65,11 @@ public:
virtual float evaluate(const ColumnVector& x) const = 0; //evaluate the function
virtual ~EvalFunction(){};
+ virtual void minimize(ColumnVector& x)
+ {
+ minsearch(x,*this);
+ }
+
private:
const EvalFunction& operator=(EvalFunction& par);
EvalFunction(const EvalFunction&);
@@ -62,31 +84,20 @@ public:
virtual ReturnMatrix g_evaluate(const ColumnVector& x) const = 0; //evaluate the gradient
virtual ~gEvalFunction(){};
+ virtual void minimize(ColumnVector& x)
+ {
+ scg(x,*this);
+ }
+
private:
const gEvalFunction& operator=(gEvalFunction& par);
gEvalFunction(const gEvalFunction&);
};
-
-float diff1(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff derivative
-
-float diff2(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff 2nd derivative
-
-float diff2(const ColumnVector& x, const EvalFunction& func, int i,int j,float h,int errorord=4);// finite diff cross derivative
-
-ReturnMatrix gradient(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff derivative vector
-
-ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff hessian
-
-void fminsearch(ColumnVector& x, const EvalFunction& func);
-
-void scg(ColumnVector& x, const gEvalFunction& func, float tol = 0.0001, float eps=1e-16, int niters=100);
}
-
-
#endif
diff --git a/miscmaths.cc b/miscmaths.cc
index e6f0237..53ee0b5 100644
--- a/miscmaths.cc
+++ b/miscmaths.cc
@@ -1669,7 +1669,7 @@ ReturnMatrix read_vest(string p_fname)
return p_mat;
}
-void ols(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Matrix& varcope){
+void OLS(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Matrix& varcope){
// ols
// data is t x v
// des is t x ev (design matrix)
diff --git a/miscmaths.h b/miscmaths.h
index 123d4ba..a75c198 100644
--- a/miscmaths.h
+++ b/miscmaths.h
@@ -202,7 +202,7 @@ namespace MISCMATHS {
// tc is cons x ev (contrast matrix)
// cope and varcope will be cons x v
// but will be resized if they are wrong
- void ols(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Matrix& varcope);
+ void OLS(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Matrix& varcope);
int ols_dof(const Matrix& des);
diff --git a/quick.cc b/quick.cc
index 1217739..976093c 100644
--- a/quick.cc
+++ b/quick.cc
@@ -6,6 +6,7 @@
#define WANT_MATH
#include "miscmaths.h"
+#include "sparse_matrix.h"
#include "libprob.h"
using namespace MISCMATHS;
@@ -13,9 +14,30 @@ using namespace MISCMATHS;
int main(int argc, char *argv[])
{
try{
- double p = MISCMATHS::fdtr(1.4449e+12, 1, 1);
+
+ SparseMatrix x(3,4);
+ SparseMatrix y(3,4);
+ x.insert(1,1,1);
+ y.insert(1,1,2);
+ y.insert(1,2,3);
+
+ x.insert(2,3,1.5);
+
+ y.insert(2,4,3);
+ x.insert(2,4,-1);
+
+ y.insert(3,4,6);
+
+ SparseMatrix ret;
+
+
+ OUT(x.AsMatrix());
+ OUT(y.AsMatrix());
+
+ add(x,y,ret);
+
+ OUT(ret.AsMatrix());
- OUT(p);
}
catch(Exception p_excp)
{
diff --git a/sparse_matrix.cc b/sparse_matrix.cc
index 5980f6b..4f6c6c7 100644
--- a/sparse_matrix.cc
+++ b/sparse_matrix.cc
@@ -253,6 +253,53 @@ namespace MISCMATHS {
}
}
+ void add(const SparseMatrix& lm, const SparseMatrix& rm, SparseMatrix& ret)
+ {
+ Tracer_Plus tr("SparseMatrix::add");
+
+ int nrows = lm.Nrows();
+ int ncols = lm.Ncols();
+
+ if(lm.Ncols() != rm.Ncols() || lm.Nrows() != rm.Nrows()) throw Exception("Rows and cols don't match in SparseMatrix::add");
+
+ ret.ReSize(nrows,ncols);
+
+ for(int j = 1; j<=nrows; j++)
+ {
+ const SparseMatrix::Row& lmrow = lm.row(j);
+ const SparseMatrix::Row& rmrow = rm.row(j);
+
+ SparseMatrix::Row::const_iterator lmit = lmrow.begin();
+ SparseMatrix::Row::const_iterator rmit = rmrow.begin();
+ int lmc = (*lmit).first+1;
+ int rmc = (*rmit).first+1;
+
+ while(lmit!=lmrow.end() || rmit!=rmrow.end())
+ {
+ if((lmc<rmc && lmit!=lmrow.end()) || rmit==rmrow.end())
+ {
+ ret.insert(j,lmc,(*lmit).second+rm(j,lmc));
+ lmit++;
+ lmc = (*lmit).first+1;
+ }
+ else if((rmc<lmc && rmit!=rmrow.end()) || lmit==lmrow.end())
+ {
+ ret.insert(j,rmc,(*rmit).second+lm(j,rmc));
+ rmit++;
+ rmc = (*rmit).first+1;
+ }
+ else
+ {
+ //lmc==rmc
+ ret.insert(j,rmc,(*lmit).second+(*rmit).second);
+ lmit++;
+ lmc = (*lmit).first+1;
+ rmit++;
+ rmc = (*rmit).first+1;
+ }
+ }
+ }
+ }
}
diff --git a/sparse_matrix.h b/sparse_matrix.h
index c6e720f..c793bbe 100644
--- a/sparse_matrix.h
+++ b/sparse_matrix.h
@@ -133,6 +133,8 @@ namespace MISCMATHS {
void multiply(const SparseMatrix& lm, const SparseMatrix::Row& rm, ColumnVector& ret);
+ void add(const SparseMatrix& lm, const SparseMatrix& rm, SparseMatrix& ret);
+
void colvectosparserow(const ColumnVector& col, SparseMatrix::Row& row);
}
--
GitLab