Skip to content

Commit

Permalink
ENH: Adding files to the initial empty repository these fiels work on…
Browse files Browse the repository at this point in the history 
… mac64 and LINUX64
  • Loading branch information
@francopestilli
francopestilli committed Jul 24, 2015
0 parents commit 90a46ad
Show file tree
Hide file tree
Showing 334 changed files with 34,344 additions and 0 deletions.
113 changes: 113 additions & 0 deletions compute/feComputeEvidence.m
@@ -0,0 +1,113 @@
function se = feComputeEvidence(rmse1,rmse2)
% Computes a series of distance metrics between two RMSE distribtions.
%
% Compute summary statistics on to characterize the lesion,:
% We compute:
% - The strength of evidence, namely the effect size of the lesion
% - The Kullback-Leibler Divergence
% - Jeffrey's Divergence
% - The Eath Mover's distance
%
% Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.

% Prepare the distribution of errors and the histograms describing the
% erros.
se.nolesion.rmse.all = rmse1;
se.lesion.rmse.all = rmse2;
se.nolesion.rmse.mean = mean(se.nolesion.rmse.all);
se.lesion.rmse.mean = mean(se.lesion.rmse.all);

% Histograms
se.xrange(1) = ceil(min([se.lesion.rmse.all se.nolesion.rmse.all]));
se.xrange(2) = floor(max([se.lesion.rmse.all se.nolesion.rmse.all]));
se.nbins = 60;
se.bins = linspace(se.xrange(1),se.xrange(2),se.nbins);
[se.lesion.hist, se.lesion.xhist] = hist(se.lesion.rmse.all, se.bins);
[se.nolesion.hist, se.nolesion.xhist] = hist(se.nolesion.rmse.all,se.bins);
se.lesion.hist = se.lesion.hist ./ sum(se.lesion.hist);
se.nolesion.hist = se.nolesion.hist ./ sum(se.nolesion.hist);

% Kullback-Leibler Divergence
se.kl.name = sprintf('Kullback–Leibler divergence: http://en.wikipedia.org/wiki/Kullback-Leibler_divergence');
tmp = se.nolesion.hist .* log2( (se.nolesion.hist) ./ (se.lesion.hist + eps) );
se.kl.mean = nansum(tmp);clear tmp
se.kl.std = nan;

% Jeffrey's divergence
se.j.name = sprintf('Jeffrey''s divergence: http://en.wikipedia.org/wiki/Divergence_(statistics)');
tmp = se.nolesion.hist .* log2( (se.nolesion.hist) ./ ((se.lesion.hist + se.lesion.hist + eps)./2) ) + ...
se.lesion.hist .* log2( (se.lesion.hist) ./ ((se.lesion.hist + se.lesion.hist + eps)./2) );
se.j.mean = nansum(tmp); clear tmp
se.j.std = nan;

% Earth Mover's Distance:
% Note: This can be very slow and my require large amounts of memory for more than 1000 voxels
fprintf('[%s] Computing the Earth Mover''s distance... \n',mfilename)
se.em.name = sprintf('Earth Mover''s distance: http://en.wikipedia.org/wiki/Earth_mover''s_distance');

try
%if (exist('emd_mex.m','file') == 2) % Using Rubinov c-code fastest
pairwiseDist = zeros(size(se.lesion.xhist,2));
for i=1:size(se.nolesion.xhist,2)
for j=1:size(se.lesion.xhist,2)
pairwiseDist(i,j) = abs(se.lesion.xhist(i)-se.nolesion.xhist(j));
end
end
tmp_em = emd_mex(se.nolesion.hist,se.lesion.hist,pairwiseDist);
catch ME %else
fprintf('[%s] Cannot find compiled c-code for Earth Movers Distance.\nUsing the slower and less reliable MatLab implementation.',mfilename)
[~,tmp_em] = emd(se.nolesion.xhist',se.lesion.xhist',se.nolesion.hist',se.lesion.hist',@gdf);
end
se.em.mean = tmp_em;
se.em.std = nan;
clear tmp_emp

% Strenght of evidence (effect size)
fprintf('[%s] Computing the Strength of Evidence... \n',mfilename)
se.s.name = sprintf('strength of evidence, d-prime: http://en.wikipedia.org/wiki/Effect_size');
se.s.nboots = 5000;
se.s.nmontecarlo = 5;
se.s.nbins = 200;
sizeunlesioned = length(se.nolesion.rmse.all);
nullDistributionW = nan(se.s.nboots,se.s.nmontecarlo);
nullDistributionWO = nan(se.s.nboots,se.s.nmontecarlo);
min_x = floor(mean([se.nolesion.rmse.all]) - mean([se.nolesion.rmse.all])*.05);
max_x = ceil( mean([se.lesion.rmse.all]) + mean([se.lesion.rmse.all])*.05);

for inm = 1:se.s.nmontecarlo
fprintf('.')
parfor ibt = 1:se.s.nboots
nullDistributionW(ibt,inm) = mean(randsample(se.nolesion.rmse.all, sizeunlesioned,true));
nullDistributionWO(ibt,inm) = mean(randsample(se.lesion.rmse.all,sizeunlesioned,true));
end

% Distribution unlesioned
[y(:,inm),xhis] = hist(nullDistributionW(:,inm),linspace(min_x,max_x,se.s.nbins));
y(:,inm) = y(:,inm)./sum(y(:,inm));

% Distribution lesioned
[woy(:,inm),woxhis] = hist(nullDistributionWO(:,inm),linspace(min_x,max_x,se.s.nbins));
woy(:,inm) = woy(:,inm)./sum(woy(:,inm));
end

se.s.mean = mean(diff([mean(nullDistributionW,1); ...
mean(nullDistributionWO,1)])./sqrt(sum([std(nullDistributionW,[],1);std(nullDistributionWO,[],1)].^2,1)));
se.s.std = std(diff([mean(nullDistributionW,1); ...
mean(nullDistributionWO,1)])./sqrt(sum([std(nullDistributionW,[],1);std(nullDistributionWO,[],1)].^2,1)));
disp(' done.')

% Extract the mean and error of the botstrapped disributions of mean errors
y_m = mean(y,2);
ywo_m = mean(woy,2);
y_e = [y_m, y_m] + 2*[-std(y,[],2),std(y,[],2)];
ywo_e = [ywo_m, ywo_m] + 2*[-std(woy,[],2),std(woy,[],2)];
se.s.lesioned_e = ywo_e;
se.s.lesioned_m = ywo_m;
se.s.unlesioned_e = y_e;
se.s.unlesioned_m = y_m;
se.s.lesioned.xbins = woxhis;
se.s.unlesioned.xbins = xhis;
se.s.min_x = min_x;
se.s.max_x = max_x;

end
100 changes: 100 additions & 0 deletions compute/feFitModel.m
@@ -0,0 +1,100 @@
function [fit, w, R2] = feFitModel(M,dSig,fitMethod)
%
% feFitModel() function in LiFE but restricted to the
%
% BBNNLS algorithm and using the Factorization model.
% M is the factorization model composed by:
%
% M.DictSig Dictionary
%
% Copyright (2015), Franco Pestilli (Indiana Univ.) - Cesar F. Caiafa (CONICET)
% email: pestillifranco@gmail.com and ccaiafa@gmail.com


% Below are the old comments which are not longer correct in general
% Fit the LiFE model.
%
% Finds the weights for each fiber to best predict the directional
% diffusion signal (dSig)
%
% fit = feFitModel(M,dSig,fitMethod)
%
% dSig: The diffusion weighted signal measured at each
% voxel in each direction. These are extracted from
% the dwi data at some white-matter coordinates.
% M: The LiFE difusion model matrix, constructed
% by feConnectomeBuildModel.m
%
% fitMethod:
% - 'bbnnls' - DEFAULT and best, faster large-scale solver.
%
% See also: feCreate.m, feConnectomeBuildModel.m, feGet.m, feSet.m
%
% Example:
%
% Copyright (2015), Franco Pestilli (Indiana Univ.) - Cesar F. Caiafa (CONICET)
% email: pestillifranco@gmail.com and ccaiafa@gmail.com
%
% Notes about the LiFE model:
%
% The rows of the M matrix are nVoxels*nBvecs. We are going to predict the
% diffusion signal in each voxel for each direction.
%
% The columns of the M matrix are nFibers + nVoxels. The diffusion signal
% for each voxel is predicted as the weighted sum of predictions from each
% fibers that passes through a voxel plus an isotropic (CSF) term.
%
% In addition to M, we typically return dSig, which is the signal measured
% at each voxel in each direction. These are extracted from the dwi data
% and knowledge of the roiCoords.

% fit the model, by selecting the proper toolbox.

mycomputer = computer();
release = version('-release');

switch fitMethod
case {'bbnnls'}
[nFibers] = size(M.Phi,3); %feGet(fe,'nfibers');
%[nTheta] = size(M.DictSig,1);
[nAtoms] = size(M.DictSig,2); %feGet(fe,'natoms');
[Nvoxels] = size(M.Phi,2); %feGet(fe,'nvoxels');

tic
fprintf('\nLiFE: Computing least-square minimization with BBNNLS...\n')
opt = solopt;
opt.maxit = 5000;
opt.use_tolo = 1;

switch strcat(mycomputer,'_',release)
case {'GLNXA64_2015a'}
out_data = bbnnls_GLNXA64(M,dSig,zeros(nFibers,1),opt);
case {'MACI64_2014b'}
out_data = bbnnls_MACI64(M,dSig,zeros(nFibers,1),opt);
otherwise
sprintf('WARNING: currently LiFE is optimized for an efficient usage of memory \n using the Sparse Tucker Decomposition aproach (Caiafa&Pestilli, 2015) \n ONLY for Linux (MatlabR2015a) and MacOS (MatlabR2014b). \n If you have a different system or version you can still \n use the old version of LiFE (memory intensive). \n\n')
sprintf('\n Starting using old version of LiFE...\n')
out_data = bbnnls_OLD(M.MmatrixM,dSig,zeros(nFibers,1),opt);
end
fprintf('BBNNLS status: %s\nReason: %s\n',out_data.status,out_data.termReason);
w = out_data.x;
fprintf(' ...fit process completed in %2.3fminutes\n',toc/60)
% Save the state of the random generator so that the stochasit cfit can be recomputed.
defaultStream = RandStream.getGlobalStream; %RandStream.getDefaultStream;
fit.randState = defaultStream.State;

% Save out some results
fit.results.R2 = [];
fit.results.nParams = size(M,2);
fit.results.nMeasures = size(M,1);
R2=[];

otherwise
error('Cannot fit LiFE model using method: %s.\n',fitMethod);
end

% Save output structure.
fit.weights = w;
fit.params.fitMethod = fitMethod;

end

0 comments on commit 90a46ad

Please sign in to comment.