Skip to content
Snippets Groups Projects
Commit 64152c77 authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

MNT: include cprob/libprob.h rather than libprob.h

parent db67b13c
No related branches found
No related tags found
1 merge request!9mnt/conda
Hide whitespace changes
Inline Side-by-side
Makefile
+ 6
− 7
View file @ 64152c77
......@@ -7,17 +7,17 @@ OPTFLAGS_alphaev6-dec-osf5.0-gcc2.95.2 = -O3 -mieee -mfp-trap-mode=sui
PROJNAME = melodic
USRINCFLAGS = -I${INC_NEWMAT} -I${INC_PROB} -I${INC_GD} -I${INC_GDC} -I${INC_PNG} -I${INC_ZLIB} -DCIFTILIB_USE_XMLPP -I${FSLEXTINC} -I${INC_XML2} -I${INC_XML++} -I${INC_XML++CONF} -I${INC_BOOST} -I${FSLDIR}/include/ciftiio
USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB} -L${LIB_GD} -L${LIB_GDC} -L${LIB_PNG} -L${LIB_ZLIB}
USRINCFLAGS = -I${INC_NEWMAT} -I${INC_CPROB} -I${INC_GD} -I${INC_GDC} -I${INC_PNG} -I${INC_ZLIB} -DCIFTILIB_USE_XMLPP -I${FSLEXTINC} -I${INC_XML2} -I${INC_XML++} -I${INC_XML++CONF} -I${INC_BOOST} -I${FSLDIR}/include/ciftiio
USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_CPROB} -L${LIB_GD} -L${LIB_GDC} -L${LIB_PNG} -L${LIB_ZLIB}
UNAME := $(shell uname)
ifeq (${UNAME},Darwin)
LIBS = -liconv
endif
LIBS += -lutils -lnewimage -lmiscplot -lmiscpic -lmiscmaths -lNewNifti -lcifti -lxml++-2.6 -lxml2 -lboost_filesystem -lboost_system -lznz -lnewmat -lprob -lm -lgdc -lgd -lpng -lz
LIBS += -lutils -lnewimage -lmiscplot -lmiscpic -lmiscmaths -lNewNifti -lcifti -lxml++-2.6 -lxml2 -lboost_filesystem -lboost_system -lznz -lcprob -lm -lgdc -lgd -lpng -lz
TEST_OBJS = test.o
TEST_OBJS = test.o
GGMIX_OBJS = ggmix.o
......@@ -33,7 +33,7 @@ FSL_REGFILT_OBJS = melhlprfns.o fsl_regfilt.o
FSL_SCHURPROD_OBJS = melhlprfns.o fsl_schurprod.o
MELODIC_OBJS = meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o melreport.o melodic.o
MELODIC_OBJS = meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o melreport.o melodic.o
TESTXFILES = test
......@@ -71,7 +71,7 @@ SCRIPTS = melodicreport dual_regression
all: ggmix fsl_regfilt fsl_glm melodic ${OTHERFILES}
ggmix: ${GGMIX_OBJS}
${AR} -r libggmix.a ${GGMIX_OBJS}
${AR} -r libggmix.a ${GGMIX_OBJS}
melodic: ${MELODIC_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${MELODIC_OBJS} ${LIBS}
......@@ -96,4 +96,3 @@ fsl_bwtf: ${FSL_BWTF_OBJS}
fsl_regfilt: ${FSL_REGFILT_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_REGFILT_OBJS} ${LIBS}
melhlprfns.cc
+ 89
− 89
View file @ 64152c77
/* MELODIC - Multivariate exploratory linear optimized decomposition into
/* MELODIC - Multivariate exploratory linear optimized decomposition into
independent components
melhlprfns.cc - misc functions
Christian F. Beckmann, FMRIB Analysis Group
Copyright (C) 1999-2013 University of Oxford */
/* CCOPYRIGHT */
#include "melhlprfns.h"
#include "libprob.h"
#include "cprob/libprob.h"
#include "miscmaths/miscprob.h"
#include "miscmaths/t2z.h"
#include "miscmaths/f2z.h"
......@@ -42,10 +42,10 @@ namespace Melodic{
if ( newmaskM(1,originalColumn)) {
for ( int64_t row=1;row<=nOriginalRows;row++)
newData(row,newColumn)=Data(row,originalColumn);
newColumn++;
newColumn++;
}
Data=newData;
mask = newmask[0];
mask = newmask[0];
}
}
......@@ -53,14 +53,14 @@ namespace Melodic{
{
for(int ctr=1; ctr<=howmany; ctr++){
in.deletevolume(ctr);
}
}
}
Matrix calc_FFT(const Matrix& Mat, const bool logpwr)
{
Matrix res;
for(int ctr=1; ctr <= Mat.Ncols(); ctr++){
ColumnVector tmpCol;
ColumnVector tmpCol;
tmpCol=Mat.Column(ctr);
ColumnVector FtmpCol_real;
ColumnVector FtmpCol_imag;
......@@ -84,19 +84,19 @@ namespace Melodic{
int ctr1 = temp.Nrows();
Matrix temp2(temp);
temp2=0;
temp = temp.Row(4) & temp.Row(3) & temp.Row(2) & temp & temp.Row(ctr1-1)
temp = temp.Row(4) & temp.Row(3) & temp.Row(2) & temp & temp.Row(ctr1-1)
& temp.Row(ctr1-2) &temp.Row(ctr1-3);
double kern[] ={0.0045 , 0.055, 0.25, 0.4, 0.25, 0.055, 0.0045};
double fac = 0.9090909;
for(int cc=1;cc<=temp2.Ncols();cc++){
for(int cr=1;cr<=temp2.Nrows();cr++){
temp2(cr,cc) = fac*( kern[0] * temp(cr,cc) + kern[1] * temp(cr+1,cc) +
kern[2] * temp(cr+2,cc) + kern[3] * temp(cr+3,cc) +
kern[4] * temp(cr+4,cc) + kern[5] * temp(cr+5,cc) +
temp2(cr,cc) = fac*( kern[0] * temp(cr,cc) + kern[1] * temp(cr+1,cc) +
kern[2] * temp(cr+2,cc) + kern[3] * temp(cr+3,cc) +
kern[4] * temp(cr+4,cc) + kern[5] * temp(cr+5,cc) +
kern[6] * temp(cr+6,cc));
}
}
......@@ -110,19 +110,19 @@ namespace Melodic{
float eps = 0.00001;
for(int ctr=1; ctr <= inp.Ncols(); ctr++){
if(meanimg(1,ctr) < eps)
if(meanimg(1,ctr) < eps)
meanimg(1,ctr) = eps;
}
for(int ctr=1; ctr <= inp.Nrows(); ctr++){
Matrix tmp;
tmp << inp.Row(ctr);
inp.Row(ctr) << 100 * SP((tmp - meanimg),pow(meanimg,-1));
inp.Row(ctr) << 100 * SP((tmp - meanimg),pow(meanimg,-1));
}
inp = remmean(inp);
return meanimg;
} //void convert_to_pbsc
} //void convert_to_pbsc
RowVector varnorm(Matrix& in, int dim, float level, int econ)
{
......@@ -137,16 +137,16 @@ namespace Melodic{
for(int ctr=1; ctr <=in.Nrows();ctr++)
in.Row(ctr) = SD(in.Row(ctr),vars);
}
RowVector varnorm(Matrix& in, SymmetricMatrix& Corr, int dim, float level, int econ)
{
{
Matrix tmpE, white, dewhite;
RowVector tmpD, tmpD2;
std_pca(remmean(in,2), Corr, tmpE, tmpD, econ);
calc_white(tmpE,tmpD, dim, white, dewhite);
Matrix ws = white * in;
for(int ctr1 = 1; ctr1<=ws.Ncols(); ctr1++)
for(int ctr2 = 1; ctr2<=ws.Nrows(); ctr2++)
......@@ -188,7 +188,7 @@ namespace Melodic{
{
for(int ctr=1; ctr <= in.Nrows(); ctr++){
in.Row(ctr) << SP(in.Row(ctr),weights.AsRow());
}
}
}
Matrix corrcoef(const Matrix& in1, const Matrix& in2){
......@@ -200,12 +200,12 @@ namespace Melodic{
}
Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part){
Matrix tmp1 = in1, tmp2 = in2, out;
Matrix tmp1 = in1, tmp2 = in2, out;
if(part.Storage()){
tmp1 = tmp1 - part * pinv(part) * tmp1;
tmp2 = tmp2 - part * pinv(part) * tmp2;
}
out = corrcoef(tmp1,tmp2);
return out;
}
......@@ -216,7 +216,7 @@ namespace Melodic{
// tmpD2= tmpD | tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols());
// cerr << tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols()) << endl;
//
//
// Matrix tmpPE;
// tmpPE = SP(param,ones(param.Nrows(),1)*pow(stdev(param,1)*std::sqrt((float)param.Nrows()),-1));
......@@ -243,9 +243,9 @@ namespace Melodic{
RD << abs(diag(tmpD2.t()));
// cerr << " " <<tmpD2.Ncols() << " " << N << " " << dim << endl;
RD << RD.SymSubMatrix(N-dim+1+param.Ncols(),N+param.Ncols());
RD << RD.SymSubMatrix(N-dim+1+param.Ncols(),N+param.Ncols());
float res = 1.0;
float res = 1.0;
white = sqrt(abs(RD.i()))*RE.t();
dewhite = RE*sqrt(RD);
......@@ -262,9 +262,9 @@ namespace Melodic{
dim = std::min(dim,N);
RE = tmpE.Columns(N-dim+1,N);
RD << abs(diag(tmpD.t()));
RD << RD.SymSubMatrix(N-dim+1,N);
RD << RD.SymSubMatrix(N-dim+1,N);
float res = 1.0;
float res = 1.0;
white = sqrt(abs(RD.i()))*RE.t();
dewhite = RE*sqrt(RD);
......@@ -277,14 +277,14 @@ namespace Melodic{
{
RowVector tmp(tmpE.Ncols());
float tmp2;
tmp2 = calc_white(tmpE,tmpD, tmp, dim, white, dewhite);
tmp2 = calc_white(tmpE,tmpD, tmp, dim, white, dewhite);
} //Matrix calc_white
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite)
{
RowVector tmp(tmpE.Ncols());
float tmp2;
tmp2 = calc_white(tmpE,tmpD, tmp, dim, param, paramS, white, dewhite);
tmp2 = calc_white(tmpE,tmpD, tmp, dim, param, paramS, white, dewhite);
} //Matrix calc_white
void calc_white(const SymmetricMatrix& Corr, int dim, Matrix& white, Matrix& dewhite)
......@@ -294,10 +294,10 @@ namespace Melodic{
RowVector tmp2;
EigenValues(Corr,RD,RE);
tmp2 = diag(RD).t();
calc_white(RE,tmp2, dim, white, dewhite);
calc_white(RE,tmp2, dim, white, dewhite);
} //Matrix calc_white
void std_pca(const Matrix& Mat, const Matrix& weights, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ)
{
if(weights.Storage()>0)
......@@ -346,7 +346,7 @@ namespace Melodic{
tmp2 = Mat*tmp.t()*tmp2;
C = tmp2;
}
symm_orth(C);
Matrix Evc;
SymmetricMatrix tmpC;
......@@ -358,7 +358,7 @@ namespace Melodic{
} //void em_pca
float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
{
{
SymmetricMatrix Corr;
Matrix Evecs, tmpWM, tmpDWM, tmp;
RowVector Evals;
......@@ -379,9 +379,9 @@ namespace Melodic{
for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++){
Matrix tmpVals(cols.Nrows(),rows.Nrows());
for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() *
tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() *
(ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1);
Matrix tmpcols, tmprows;
res(ctr1) =rankapprox(tmpVals, tmpcols, tmprows);
cols.Column(ctr1) = tmpcols;
......@@ -436,12 +436,12 @@ namespace Melodic{
for(int ctr=1;ctr<=CircleLaw.Ncols(); ctr++){
if(CircleLaw(ctr)<5*10e-10){CircleLaw(ctr) = 5*10e-10;}
}
}
/* float slope;
slope = CircleLaw.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
int(AdjEV.Ncols()/4)).Sum() /
AdjEV.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
slope = CircleLaw.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
int(AdjEV.Ncols()/4)).Sum() /
AdjEV.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
int(AdjEV.Ncols()/4)).Sum();*/
RowVector PercEV(AdjEV);
......@@ -452,7 +452,7 @@ namespace Melodic{
SortDescending(AdjEV);
int maxEV = 1;
float threshold = 0.98;
for(int ctr_i = 1; ctr_i < PercEV.Ncols(); ctr_i++){
for(int ctr_i = 1; ctr_i < PercEV.Ncols(); ctr_i++){
if((PercEV(ctr_i)<threshold)&&(PercEV(ctr_i+1)>=threshold)){maxEV=ctr_i;}
}
......@@ -475,7 +475,7 @@ namespace Melodic{
RowVector AdjEV, PercEV;
AdjEV = in.Reverse();
SortDescending(AdjEV);
PercEV = Melodic::cumsum(AdjEV / sum(AdjEV,2).AsScalar());
AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
out1 = AdjEV;
......@@ -484,7 +484,7 @@ namespace Melodic{
RowVector Feta(int n1, int n2)
{
float nu = (float) n1/n2;
float nu = (float) n1/n2;
float bm = pow((1-sqrt(nu)),2.0);
float bp = pow((1+sqrt(nu)),2.0);
......@@ -492,21 +492,21 @@ namespace Melodic{
float urange = 1.1*bp;
RowVector eta(30*n1);
float rangestepsize = (urange - lrange) / eta.Ncols();
for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
float rangestepsize = (urange - lrange) / eta.Ncols();
for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
eta(ctr_i) = lrange + rangestepsize * (ctr_i);
}
RowVector teta(10*n1);
teta = 0;
float stepsize = (bp - bm) / teta.Ncols();
for(int ctr_i = 1; ctr_i <= teta.Ncols(); ctr_i++){
for(int ctr_i = 1; ctr_i <= teta.Ncols(); ctr_i++){
teta(ctr_i) = stepsize*(ctr_i);
}
}
RowVector feta(teta);
feta = SP(pow(2*M_PI*nu*(teta + bm),-1), pow(SP(teta, bp-bm-teta),0.5));
teta = teta + bm;
RowVector claw(eta);
......@@ -518,7 +518,7 @@ namespace Melodic{
tmpval += feta(ctr_j);
claw(ctr_i) = max(n1*(1-stepsize*tmpval),0.0);
}
RowVector Res(n1); //invert the CDF
Res = 0;
for(int ctr_i = 1; ctr_i < eta.Ncols(); ctr_i++){
......@@ -526,7 +526,7 @@ namespace Melodic{
Res(int(floor(claw(ctr_i)))) = eta(ctr_i);
}
}
return Res;
} //RowVector Feta
......@@ -540,23 +540,23 @@ namespace Melodic{
} //RowVector cumsum
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);
int maxEV = 1;
RowVector NewEV;
adj_eigspec(tmpD.AsRow(),AdjEV,PercEV,NewEV,maxEV,in.Ncols(),resels);
int res;
PPCA = ppca_est(NewEV, in.Ncols(),resels);
ColumnVector tmp = ppca_select(PPCA, res, maxEV, which);
PPCA = tmp | PPCA;
return res;
} //int ppca_dim
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, float resels, string which)
{
{
RowVector tmpD;
Matrix tmpE;
SymmetricMatrix Corr;
......@@ -577,14 +577,14 @@ namespace Melodic{
{
RowVector estimators(5);
estimators = 1.0;
for(int ctr=1; ctr<=PPCAest.Ncols(); ctr++){
PPCAest.Column(ctr) = (PPCAest.Column(ctr) -
min(PPCAest.Column(ctr)).AsScalar()) /
( max(PPCAest.Column(ctr)).AsScalar() -
PPCAest.Column(ctr) = (PPCAest.Column(ctr) -
min(PPCAest.Column(ctr)).AsScalar()) /
( max(PPCAest.Column(ctr)).AsScalar() -
min(PPCAest.Column(ctr)).AsScalar());
}
int ctr_i = 1;
while((ctr_i< PPCAest.Nrows()-1)&&
(PPCAest(ctr_i,2) < PPCAest(ctr_i+1,2))&&(ctr_i<maxEV))
......@@ -612,7 +612,7 @@ namespace Melodic{
PercEV = Melodic::cumsum(PercEV / sum(PercEV,2).AsScalar());
if(which == string("aut")) {
if(int(estimators(2)) < int(estimators(1)) &&
if(int(estimators(2)) < int(estimators(1)) &&
float(PercEV(int(estimators(2))))>0.8){
res=int(estimators(2));
PPCA << PPCAest.Column(3);
......@@ -620,7 +620,7 @@ namespace Melodic{
res = int(estimators(1));
PPCA << PPCAest.Column(2);
}
}
}
if(which == string("lap")){
res = int(estimators(1));
PPCA << PPCAest.Column(2);
......@@ -642,11 +642,11 @@ namespace Melodic{
PPCA << PPCAest.Column(6);
}
if(which == string("median")){
RowVector unsorted(estimators);
RowVector unsorted(estimators);
SortAscending(unsorted);
ctr_i=1;
res=int(unsorted(3));
while(res != int(estimators(ctr_i)))
while(res != int(estimators(ctr_i)))
ctr_i++;
PPCA << PPCAest.Column(ctr_i);
}
......@@ -660,7 +660,7 @@ namespace Melodic{
} //RowVector ppca_select
Matrix ppca_est(const RowVector& eigenvalues, const int N1, const float N2)
{
{
Matrix Res;
Res = ppca_est(eigenvalues, (int) floor(N1/(2.5*N2)));
return Res;
......@@ -677,16 +677,16 @@ namespace Melodic{
for(int ctr_i = 1; ctr_i <=d; ctr_i++){
k(ctr_i)=ctr_i;
}
RowVector m(d);
m=d*k-0.5*SP(k,k+1);
m=d*k-0.5*SP(k,k+1);
RowVector loggam(d);
loggam = 0.5*k.Reverse();
for(int ctr_i = 1; ctr_i <=d; ctr_i++){
loggam(ctr_i)=lgam(loggam(ctr_i));
}
loggam = Melodic::cumsum(loggam);
loggam = Melodic::cumsum(loggam);
RowVector l_probU(d);
l_probU = -log(2)*k + loggam - Melodic::cumsum(0.5*log(M_PI)*k.Reverse());
......@@ -706,7 +706,7 @@ namespace Melodic{
tmp3(d) = 1.0;
RowVector tmp4;
tmp4 = SP(tmp1,pow(tmp3,-1));
tmp4 = SP(tmp1,pow(tmp3,-1));
for(int ctr_i = 1; ctr_i <=d; ctr_i++){
if(tmp4(ctr_i)<0.01){tmp4(ctr_i)=0.01;}
if(tmp3(ctr_i)<0.01){tmp3(ctr_i)=0.01;}
......@@ -737,13 +737,13 @@ namespace Melodic{
for(int ctr_i = 1; ctr_i <= t1.Ncols(); ctr_i++){
for(int ctr_j = 1; ctr_j <= t1.Nrows(); ctr_j++){
if(t1(ctr_j,ctr_i)<=0){t1(ctr_j,ctr_i)=1;}
}
}
}
for(int ctr_i = 1; ctr_i <= t2.Ncols(); ctr_i++){
for(int ctr_j = 1; ctr_j <= t2.Nrows(); ctr_j++){
if(t2(ctr_j,ctr_i)<=0){t2(ctr_j,ctr_i)=1;}
}
}
}
t1 = Melodic::cumsum(sum(log(t1),2).AsRow());
t2 = Melodic::cumsum(sum(log(t2),2).AsRow());
......@@ -752,7 +752,7 @@ namespace Melodic{
RowVector l_lap;
l_lap = l_probU + l_nu +l_Az + l_lam + 0.5*log(2*M_PI)*(m+k)-0.5*log(N)*k;
RowVector l_BIC;
l_BIC = l_lam + l_nu - 0.5*log(N)*(m+k);
......@@ -775,8 +775,8 @@ namespace Melodic{
Res |= l_MDL.t();
Res |= l_RRN.t();
Res |= l_AIC.t();
return Res;
} //Matrix ppca_est
......@@ -785,7 +785,7 @@ namespace Melodic{
double meanval;
meanval = mean(in).AsScalar();
int tpoints = in.Nrows();
ColumnVector y, res;
Matrix X, tmp;
......@@ -806,11 +806,11 @@ namespace Melodic{
ColumnVector res;
res = acf(in, maxorder);
for(int ctr1 = 1; ctr1 <= maxorder; ctr1++)
if ( res.Column(ctr1).AsScalar() < (1/tpoint) + 2/(float)std::pow(tpoint,0.5))
if ( res.Column(ctr1).AsScalar() < (1/tpoint) + 2/(float)std::pow(tpoint,0.5))
res.Column(ctr1) = 0;
return res;
} //ColumnVector pacf
Matrix est_ar(const Matrix& Mat, int maxorder)
{
Matrix res;
......@@ -845,7 +845,7 @@ namespace Melodic{
for(int ctr=1; ctr <= in.Ncols(); ctr++){
tmp = in.Column(ctr);
res.Column(ctr) = gen_ar(tmp, maxorder);
}
}
return res;
} //Matrix gen_ar
......@@ -862,13 +862,13 @@ namespace Melodic{
return res;
} //Matrix gen_arCorr
void basicGLM::olsfit(const Matrix& data, const Matrix& design,
void basicGLM::olsfit(const Matrix& data, const Matrix& design,
const Matrix& contrasts, int requestedDOF) {
beta = zeros(design.Ncols(),1);
residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1);
beta = zeros(design.Ncols(),1);
residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1);
t = zeros(1); z = zeros(1); p=-1.0*ones(1);
dof = (int)-1; cbeta = -1.0*ones(1);
if(data.Nrows()==design.Nrows()){
dof = (int)-1; cbeta = -1.0*ones(1);
if(data.Nrows()==design.Nrows()){
if ( requestedDOF>0)
dof = requestedDOF;
else
......@@ -878,9 +878,9 @@ namespace Melodic{
Matrix fullFit(design*beta);
residu= data-fullFit;
sigsq = sum(SP(residu,residu))/(float)dof; //residuals now hold sigmasquared
fullFit = sum(SP(fullFit,fullFit))/(float)design.Ncols();
fullFit = sum(SP(fullFit,fullFit))/(float)design.Ncols();
f_fmf= SD(fullFit,sigsq);
pf_fmf = f_fmf.Row(1);
pf_fmf = f_fmf.Row(1);
for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),dof,f_fmf.Column(ctr1).AsScalar());
if(contrasts.Storage()>0 && contrasts.Ncols()==beta.Nrows()){
......@@ -888,7 +888,7 @@ namespace Melodic{
Matrix tmp = contrasts*pinvModel*pinvModel.t()*contrasts.t();
varcb = diag(tmp)*sigsq;
t = SP(cbeta,pow(varcb,-0.5));
z = p = t;
z = p = t;
T2z& t2z = T2z::getInstance();
double logp;
for(int64_t col=1;col<=t.Ncols();col++)
......
melpca.cc
+ 24
− 26
View file @ 64152c77
/* MELODIC - Multivariate exploratory linear optimized decomposition into
/* MELODIC - Multivariate exploratory linear optimized decomposition into
independent components
melpca.cc - PCA and whitening
melpca.cc - PCA and whitening
Christian F. Beckmann, FMRIB Analysis Group
Copyright (C) 1999-2013 University of Oxford */
/* CCOPYRIGHT */
......@@ -17,19 +17,19 @@
#include "newmatio.h"
#include "melpca.h"
#include "melhlprfns.h"
#include "libprob.h"
#include "cprob/libprob.h"
using namespace Utilities;
namespace Melodic{
void MelodicPCA::perf_pca(Matrix& in, Matrix& weights){
void MelodicPCA::perf_pca(Matrix& in, Matrix& weights){
message("Starting PCA ... ");
SymmetricMatrix Corr;
Matrix tmpE;
RowVector tmpD, AdjEV, PercEV;
if(opts.paradigmfname.value().length()>0)
{
basicGLM tmpglm;
......@@ -37,28 +37,28 @@ namespace Melodic{
std_pca(tmpglm.get_residu(),weights,Corr,tmpE,tmpD);
}
else{
std_pca(in,weights,Corr,tmpE,tmpD);
}
std_pca(in,weights,Corr,tmpE,tmpD);
}
if(opts.tsmooth.value()){
message(endl << " temporal smoothing of Eigenvectors " << endl);
tmpE=smoothColumns(tmpE);
}
adj_eigspec(tmpD,AdjEV,PercEV);
melodat.set_pcaE(tmpE);
melodat.set_pcaD(tmpD);
melodat.set_EVP(PercEV);
melodat.set_EVP(PercEV);
melodat.set_EV((AdjEV));
write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV);
message("done" << endl);
}
void MelodicPCA::perf_white(){
int N = melodat.get_pcaE().Ncols();
int N = melodat.get_pcaE().Ncols();
if(opts.pca_dim.value() > N){
message("dimensionality too large - using -dim " << N <<
message("dimensionality too large - using -dim " << N <<
" instead " << endl);
opts.pca_dim.set_T(N);
}
......@@ -89,10 +89,10 @@ namespace Melodic{
float varp = 1.0;
if (opts.paradigmfname.value().length()>0)
varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
melodat.get_EVP(),opts.pca_dim.value(),melodat.get_param(),melodat.get_paramS(),tmpWhite,tmpDeWhite);
else
varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
melodat.get_EVP(),opts.pca_dim.value(),tmpWhite,tmpDeWhite);
melodat.set_white(tmpWhite);
......@@ -102,24 +102,22 @@ namespace Melodic{
}
int MelodicPCA::pcadim()
{
{
message("Estimating the number of sources from the data (PPCA) ..." << endl);
ColumnVector PPCA; RowVector AdjEV, PercEV;
int res = ppca_dim(melodat.get_Data(),melodat.get_RXweight(), PPCA,
ColumnVector PPCA; RowVector AdjEV, PercEV;
int res = ppca_dim(melodat.get_Data(),melodat.get_RXweight(), PPCA,
AdjEV, PercEV, melodat.get_resels(), opts.pca_est.value());
write_ascii_matrix(logger.appendDir("PPCA"),PPCA);
write_ascii_matrix(logger.appendDir("PPCA"),PPCA);
write_ascii_matrix(logger.appendDir("eigenvalues_adjusted"),AdjEV.t());
write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV.t());
melodat.set_EVP(PercEV);
melodat.set_EV(AdjEV);
melodat.set_PPCA(PPCA);
return res;
}
}
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