Skip to content
Snippets Groups Projects
Commit 49644732 authored by Matthew Webster's avatar Matthew Webster
Browse files

Changed Matrix to SymmetricMatrix for Covariance

parent 7c7f613e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
meldata.cc
+ 5
− 3
View file @ 49644732
......@@ -367,7 +367,8 @@ namespace Melodic{
//estimate model order
Matrix tmpPPCA;
RowVector AdjEV, PercEV;
Matrix Corr, tmpE;
Matrix tmpE;
SymmetricMatrix Corr;
int order;
order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
......@@ -506,7 +507,8 @@ namespace Melodic{
//reduce dim down to manageable level
if(Data.Nrows() > opts.migp_factor.value()*opts.migpN.value() || ctr==numfiles-1){
message(" Reducing data matrix to a " << opt.migpN.value() << " dimensional subspace " << endl);
Matrix pcaE, Corr;
Matrix pcaE;
SymmetricMatrix Corr;
RowVector pcaD;
std_pca(Data, RXweight, Corr, pcaE, pcaD, opts.econ.value());
pcaE = pcaE.Columns(pcaE.Ncols()-opts.migpN.value()+1,pcaE.Ncols());
......@@ -518,7 +520,7 @@ namespace Melodic{
//update mask
if(opts.update_mask.value()){
message(endl<< "Excluding voxels with constant value ...");
message(endl<< "Excluding voxels with constant value ...");
update_mask(Mask, Data);
message(" done" << endl);
}
......
melhlprfns.cc
+ 45
− 86
View file @ 49644732
......@@ -127,21 +127,19 @@ namespace Melodic{
RowVector varnorm(Matrix& in, int dim, float level, int econ)
{
Matrix Corr;
Corr = calc_corr(in,econ);
SymmetricMatrix Corr(cov_r(in,false,econ));
RowVector out;
out = varnorm(in,Corr,dim,level, econ);
return out;
} //RowVector varnorm
void varnorm(Matrix& in, const RowVector& vars)
{
for(int ctr=1; ctr <=in.Nrows();ctr++){
in.Row(ctr) = SD(in.Row(ctr),vars);
}
}
void varnorm(Matrix& in, const RowVector& vars)
{
for(int ctr=1; ctr <=in.Nrows();ctr++)
in.Row(ctr) = SD(in.Row(ctr),vars);
}
RowVector varnorm(Matrix& in, Matrix& Corr, int dim, float level, int econ)
RowVector varnorm(Matrix& in, SymmetricMatrix& Corr, int dim, float level, int econ)
{
Matrix tmpE, white, dewhite;
......@@ -211,43 +209,6 @@ namespace Melodic{
return out;
}
Matrix calc_corr(const Matrix& in, int econ)
{
Matrix Res;
int nrows=in.Nrows();
int ncols=in.Ncols();
Res = zeros(nrows,nrows);
if(econ>0){
ColumnVector rowmeans(nrows);
for (int n=1; n<=nrows; n++) {
rowmeans(n)=0;
for (int m=1; m<=ncols; m++) {
rowmeans(n)+=in(n,m);
}
rowmeans(n)/=ncols;
}
int dcol = econ;
Matrix suba;
for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
int scolmax = suba.Ncols();
for (int n=1; n<=nrows; n++) {
for (int m=1; m<=scolmax; m++) {
suba(n,m)-=rowmeans(n);
}
}
Res += suba*suba.t() / ncols;
}
}
else
Res = cov(in.t());
return Res;
} //Matrix calc_corr
Matrix calc_corr(const Matrix& in, const Matrix& weights, int econ)
{
Matrix Res;
......@@ -259,36 +220,36 @@ namespace Melodic{
localweights = ones(1,ncols);
else
localweights = weights;
if(econ>0){
ColumnVector rowmeans(nrows);
for (int n=1; n<=nrows; n++) {
rowmeans(n)=0;
for (int m=1; m<=ncols; m++) {
rowmeans(n)+=in(n,m);
}
rowmeans(n)/=ncols;
}
int dcol = econ;
Matrix suba;
for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
int scolmax = suba.Ncols();
for (int n=1; n<=nrows; n++) {
double locw =localweights(ctr + n - 1);
for (int m=1; m<=scolmax; m++) {
suba(n,m)-=rowmeans(n);
suba(n,m)*=locw;
if(econ>0){
RowVector colmeans(ncols);
for (int n=1; n<=ncols; n++) {
colmeans(n)=0;
for (int m=1; m<=nrows; m++) {
colmeans(n)+=in(m,n);
}
}
Res += suba*suba.t() / ncols;
}
colmeans(n)/=nrows;
}
int dcol = econ;
Matrix suba;
for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
int scolmax = suba.Ncols();
for (int n=1; n<=scolmax; n++) {
double cmean=colmeans(ctr + n - 1);
double locw =localweights(ctr + n - 1);
for (int m=1; m<=nrows; m++) {
suba(m,n)-=cmean;
suba(m,n)*= locw;
}
}
Res += suba*suba.t() / ncols;
}
}
else{
Res = SP2(in,localweights);
Res = calc_corr(Res, 0);
Res = cov_r(Res,false,econ);
}
return Res;
} //Matrix calc_corr
......@@ -372,33 +333,29 @@ namespace Melodic{
tmp2 = calc_white(tmpE,tmpD, tmp, dim, param, paramS, white, dewhite);
} //Matrix calc_white
void calc_white(const Matrix& Corr, int dim, Matrix& white, Matrix& dewhite)
void calc_white(const SymmetricMatrix& Corr, int dim, Matrix& white, Matrix& dewhite)
{
Matrix RE;
DiagonalMatrix RD;
SymmetricMatrix tmp;
RowVector tmp2;
tmp << Corr;
EigenValues(tmp,RD,RE);
EigenValues(Corr,RD,RE);
tmp2 = diag(RD).t();
calc_white(RE,tmp2, dim, white, dewhite);
} //Matrix calc_white
void std_pca(const Matrix& Mat, const Matrix& weights, Matrix& Corr, Matrix& evecs, RowVector& evals, int econ)
void std_pca(const Matrix& Mat, const Matrix& weights, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ)
{
if(weights.Storage()>0)
Corr << calc_corr(Mat, weights, econ);
else
Corr << calc_corr(Mat, econ);
SymmetricMatrix tmp;
tmp << Corr;
Corr = cov_r(Mat,false,econ);
DiagonalMatrix tmpD;
EigenValues(tmp,tmpD,evecs);
EigenValues(Corr,tmpD,evecs);
evals = tmpD.AsRow();
} //void std_pca
void std_pca(const Matrix& Mat, Matrix& Corr, Matrix& evecs, RowVector& evals, int econ)
void std_pca(const Matrix& Mat, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ)
{
Matrix weights;
std_pca(Mat,weights,Corr,evecs,evals, econ);
......@@ -436,7 +393,8 @@ namespace Melodic{
}
symm_orth(C);
Matrix Evc, tmpC;
Matrix Evc;
SymmetricMatrix tmpC;
RowVector Evl;
tmp = C.t() * Mat;
std_pca(tmp,tmpC,Evc,Evl);
......@@ -446,7 +404,8 @@ namespace Melodic{
float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
{
Matrix Corr, Evecs, tmpWM, tmpDWM, tmp;
SymmetricMatrix Corr;
Matrix Evecs, tmpWM, tmpDWM, tmp;
RowVector Evals;
std_pca(Mat.t(), Corr, Evecs, Evals);
calc_white(Corr, dim, tmpWM, tmpDWM);
......@@ -626,7 +585,7 @@ namespace Melodic{
return Res;
} //RowVector cumsum
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, Matrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which)
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, SymmetricMatrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which)
{
std_pca(in,weights,Corr,tmpE,tmpD);
......@@ -646,7 +605,7 @@ namespace Melodic{
{
RowVector tmpD;
Matrix tmpE;
Matrix Corr;
SymmetricMatrix Corr;
int res = ppca_dim(in, weights, PPCA, AdjEV, PercEV, Corr, tmpE, tmpD, resels, which);
return res;
......
melhlprfns.h
+ 6
− 7
View file @ 49644732
......@@ -44,9 +44,9 @@ namespace Melodic{
Matrix convert_to_pbsc(Matrix& Mat);
RowVector varnorm(Matrix& in, int dim = 30, float level = 1.6, int econ = 0);
RowVector varnorm(Matrix& in, int dim = 30, float level = 1.6, int econ = 20000);
void varnorm(Matrix& in, const RowVector& vars);
RowVector varnorm(Matrix& in, Matrix& Corr, int dim = 30, float level = 1.6, int econ = 0);
RowVector varnorm(Matrix& in, SymmetricMatrix& Corr, int dim = 30, float level = 1.6, int econ = 20000);
Matrix SP2(const Matrix& in, const Matrix& weights, int econ = 20000);
void SP3(Matrix& in, const Matrix& weights);
......@@ -56,17 +56,16 @@ namespace Melodic{
Matrix corrcoef(const Matrix& in1, const Matrix& in2);
Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part);
Matrix calc_corr(const Matrix& in, int econ = 20000);
Matrix calc_corr(const Matrix& in, const Matrix& weights, int econ = 20000);
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite);
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& white, Matrix& dewhite);
void calc_white(const Matrix& Corr, int dim, Matrix& white, Matrix& dewhite);
void calc_white(const SymmetricMatrix& Corr, int dim, Matrix& white, Matrix& dewhite);
void std_pca(const Matrix& Mat, Matrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
void std_pca(const Matrix& Mat, const Matrix& weights, Matrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
void std_pca(const Matrix& Mat, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
void std_pca(const Matrix& Mat, const Matrix& weights, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
......@@ -79,7 +78,7 @@ namespace Melodic{
void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2, RowVector& out3, int& out4, int num_vox, float resels);
void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2);
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, Matrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which);
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, SymmetricMatrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which);
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, float resels, string which);
int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which);
ColumnVector ppca_select(Matrix& PPCAest, int& dim, int maxEV, string which);
......
melica.cc
+ 3
− 5
View file @ 49644732
......@@ -34,7 +34,6 @@ namespace Melodic {
//initialize matrices
Matrix redUMM_old, rank1_old;
Matrix tmpU;
//srand((unsigned int)timer(NULL));
redUMM = melodat.get_white()*
unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
......@@ -57,10 +56,9 @@ namespace Melodic {
if(opts.approach.value() == string("tica"))
opts.maxNumItt.set_T(opts.rank1interval.value());
rank1_old = melodat.get_dewhite()*redUMM;
rank1_old = melodat.expand_dimred(rank1_old);
rank1_old = krapprox(rank1_old,int(rank1_old.Nrows()/melodat.get_numfiles()));
rank1_old = melodat.get_dewhite()*redUMM;
rank1_old = melodat.expand_dimred(rank1_old);
rank1_old = krapprox(rank1_old,int(rank1_old.Nrows()/melodat.get_numfiles()));
do{// TICA loop
itt_ctr = 1;
do{ // da loop!!!
......
melpca.cc
+ 1
− 2
View file @ 49644732
......@@ -28,7 +28,7 @@ namespace Melodic{
void MelodicPCA::perf_pca(Matrix& in, Matrix& weights){
message("Starting PCA ... ");
Matrix Corr;
SymmetricMatrix Corr;
Matrix tmpE;
RowVector tmpD, AdjEV, PercEV;
......@@ -37,7 +37,6 @@ namespace Melodic{
basicGLM tmpglm;
tmpglm.olsfit(in,melodat.get_param(),IdentityMatrix(melodat.get_param().Ncols()));
std_pca(tmpglm.get_residu(),weights,Corr,tmpE,tmpD);
// std_pca(in,weights,Corr,tmpE,tmpD,melodat.get_param());
}
else{
std_pca(in,weights,Corr,tmpE,tmpD);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment