Added ols and SD

miscmaths.cc

@@ -1407,6 +1407,27 @@ ReturnMatrix neq(const Matrix& mat1,const Matrix& mat2)
   return res;
 }
 
+ReturnMatrix SD(const Matrix& mat1,const Matrix& mat2) 
+{
+  if((mat1.Nrows() != mat2.Nrows()) ||
+     (mat1.Ncols() != mat2.Ncols()) ){
+    cerr <<"MISCMATHS::SD - matrices are of different dimensions"<<endl;
+    exit(-1);
+  }
+  Matrix ret(mat1.Nrows(),mat1.Ncols());
+  for (int r = 1; r <= mat1.Nrows(); r++) {
+    for (int c =1; c <= mat1.Ncols(); c++) {
+      if( mat2(r,c)==0)
+	ret(r,c)=0;
+      else
+	ret(r,c) = mat1(r,c)/mat2(r,c);
+    }
+  }
+
+  ret.Release();
+  return ret;
+}
+
 
 
 
@@ -1648,6 +1669,37 @@ ReturnMatrix read_vest(string p_fname)
   return p_mat;
 }
 
+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)
+  // tc is cons x ev (contrast matrix)
+  // cope and varcope will be cons x v
+  // but will be resized if they are wrong
+  // hence may be passed in uninitialised
+  // TB 2004
+  if(data.Nrows() != des.Nrows()){
+    cerr <<"MISCMATHS::ols - data and design have different number of time points"<<endl;
+    exit(-1);
+  }
+  if(des.Ncols() != tc.Ncols()){
+    cerr <<"MISCMATHS::ols - design and contrast matrix have different number of EVs"<<endl;
+    exit(-1);
+  }  
+  
+  Matrix pdes = pinv(des);
+  Matrix prevar=diag(tc*pdes*pdes.t()*tc.t());
+  Matrix R=Identity(des.Nrows())-des*pdes;
+  float tR=R.Trace();
+  Matrix pe=pdes*data;
+  cope=tc*pe;
+  Matrix res=data-des*pe;
+  Matrix sigsq=sum(SP(res,res))/tR;
+  varcope=prevar*sigsq;
+  
+  
+}
+
 int conjgrad(ColumnVector& x, const Matrix& A, const ColumnVector& b, int maxit,
 	     float reltol)
 {

miscmaths.h

@@ -186,7 +186,8 @@ namespace MISCMATHS {
   ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2); 
   ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2); 
   ReturnMatrix neq(const Matrix& mat1,const Matrix& mat2); 
-
+  ReturnMatrix SD(const Matrix& mat1,const Matrix& mat2); // Schur (element-wise) divide
+  
   void remmean(const Matrix& mat, Matrix& demeanedmat, Matrix& Mean,  const int dim = 1);
   ReturnMatrix remmean(const Matrix& mat, const int dim = 1);
   ReturnMatrix stdev(const Matrix& mat, const int dim = 1);
@@ -194,6 +195,15 @@ namespace MISCMATHS {
   ReturnMatrix corrcoef(const Matrix& mat, const int norm = 0);
   void symm_orth(Matrix &Mat);
   void powerspectrum(const Matrix &Mat1, Matrix &Result, bool useLog);
+  
+  // ols
+  // data is t x v
+  // des is t x ev (design matrix)
+  // 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);
+
 
   // Conjugate Gradient methods to solve for x in:   A * x = b 
   // A must be symmetric and positive definite