Merge pull request #59 from brain-life/revert-58-toshare_sub01

Revert "Toshare sub01"

README.md

@@ -1,96 +1,7 @@
-# Encode: Multidimensional encoding of brain connectomes
+# life
 
-![alt tag](https://cloud.githubusercontent.com/assets/11638664/18485100/66313a68-79b9-11e6-8b04-80bec2f8530e.png)
+This version is the one used for the arxiv.org paper:
 
-# About
-This software implements a framework to encode structural brain connectomes into multidimensional arrays. These arrays are commonly referred to as [tensors](https://arxiv.org/abs/1403.4462). Encoding Connectomes provides an agile framework for computing over connectome edges and nodes efficiently. We provide several examples of operations that can be performed using the framework.
-
-One major application of the tensor encoding is the implementaion of the [Linear Fascicle Evaluation method](http://francopestilli.github.io/life/), in short [LiFE](http://www.nature.com/nmeth/journal/v11/n10/abs/nmeth.3098.html). The tensor encoding method allows implementing LiFE with dramatic reduction in storage requirements, up to 40x compression factors. Furtheremore, connectome encoding allows performing multiple computational neuroanatomy operations such as tract-dissections, virtual lesions, and connectivity estimates very efficiently using the machine-friendly array operators. 
-
-We provide demos to expain how to:
- (1) Load and encode diffusion-weighted data and tractography models of white matter fascicles, as well as perform multidimensional arrays operations. 
- (2) Build and optimize a Linear Fascicle Evaluation model. 
- (4) Perform neuronatomical segmentations, computational neuroanatomy operations and virtual lesions using the connectome encoding framework.
- (4) Reproduce some fo the figures of article describing the method implemented in thsi toolbox: Caiafa and Pestilli, forthcoming.
-
-## Application.
-* Encoding of brain conenctome and associated phenotypes into multidimensional arrays.
-* Evaluate the evidence supporting white-matter connectomes generated using [magnetic resonance diffusion-weighted imaging](http://en.wikipedia.org/wiki/Diffusion_MRI) and [computational tractography ](http://en.wikipedia.org/wiki/Tractography).
-* Perform statistical inference on white-matter connectomes: Compare white-matter connectomes, show the evidence for white-matter tracts and connections between brain areas.
-
-## License.
-#### Copyright (2016), [Franco Pestilli](http://francopestilli.com/), frakkopesto@gmail.com, [Cesar Caiafa](http://web.fi.uba.ar/~ccaiafa), ccaiafa@gmail.com
-
-## [Documentation](TBA).
-
-## [Stable code release](TBA).
-
-## How to cite the software.
-[Caiafa, C. and Pestilli, F.](Multidimensional encoding of brain connectomes) Multidimensional encoding of brain connectomes (forthcoming.)
-
-## Funding.
-This work was supported by grants by the Indiana Clinical and Translational Institute (CTSI, NIH ULTTR001108).
-
-## Installation.
-1. Download (https://github.com/brain-life/encode).
-2. [Start MatLab](http://www.mathworks.com/help/matlab/startup-and-shutdown.html).
-3. Add repository to the [matlab search path](http://www.mathworks.com/help/matlab/ref/addpath.html).
-
-## Dependencies.
-* [MatLab](http://www.mathworks.com/products/matlab/).
-* [vistasoft](https://github.com/vistalab/vistasoft).
-* [Matlab Brain Anatomy (MBA)](https://github.com/francopestilli/mba).
-
-## Getting started.
-
-### 1. [Download the repository](https://github.com/brain-life/encode).
-* Download the Encode repository from the TAR/ZIP files linked [here](https://github.com/brain-life/encode/archive/master.zip).
-* UNZIP/UNTAR the file.
-* Add the encode folder to your matlab search path. To do so in the MatLab prompt type: 
-```
-   >> addpath(genpath('/my/path/to/the/encode/folder/'))
-```
-
-### 2. [Download the vistasoft repository](https://github.com/vistalab/vistasoft).
-* Download the VISTASOFT repository from the TAR/ZIP files linked [here](https://github.com/vistalab/vistasoft/archive/master.zip).
-* UNZIP/UNTAR the file.
-* Add the VISTASOFT folder to your matlab search path. To do so in the MatLab prompt type: 
-```
-   >> addpath(genpath('/my/path/to/the/VISTASOFT/folder/'))
-```
-### 3. [Download the MBA repository](https://github.com/francopestilli/mba).
-* Download the MBA repository from the TAR/ZIP files linked [here](https://github.com/francopestilli/mba/archive/master.zip).
-* UNZIP/UNTAR the file.
-* Add the MBA folder to your matlab search path. To do so in the MatLab prompt type: 
-```
-   >> addpath(genpath('/my/path/to/the/MBA/folder/'))
-```
-
-### 4. [Download the Demo Datasets](http://XXXXXXXXXXX).
-* Download the demo datasets from the repository [here](https://XXXXXX.tar.gz).
-* UNZIP/UNTAR the file.
-* Add the unzipped/untarred Data folder to your matlab search path. To do so in the MatLab prompt type:
-```
-   >> addpath(genpath('/my/path/to/the/data_demo/folder/'))
-```
-### 5. [Run the demo_connectome_encoding code](/scripts/demos/demo_connectome_encoding.m).
-Here you will learn about creating the tensor representation of a connectoms and perform basic operations such as identifying fascicles having a particular spatial orientation in a small voxel area. 
-```
-  >>  demo_connectome_encoding.m
-```
-### 6. [Run the demo_connectome_data_comparison code](/scripts/demos/demo_connectome_data_comparison.m).
-This code reproduce Fig. 3 of the paper "Multidimensional encoding of brain connectomes", by C. Caiafa and F. Pestilli. 
-```
-  >>  demo_connectome_data_comparison.m
-```
-### 7. [Run the demo_virtual_lesion code](/scripts/demos/demo_virtual_lesion.m).
-This code allows you to compute virtual lesions on a particular brain dataset and visualize particular major tracts together with their path-neighborhood, i.e. fascicles sharing same voxels. 
-```
-  >>  demo_virtual_lesion.m
-```
-### 8. [Run the demo_LiFE code](/scripts/demos/demo_LiFE.m).
-This code allows you to compute compute the fascicles weights for two different tractography methods, probabilistic and deterministic tractographies, on a same brain. This is similar to the original LiFE demo in  https://github.com/francopestilli/life but here a full brain dataset is used. The optimization (fitting fascicles weights) runs in about 3 hours on a modern Intel processor with 8GB of RAM. This code has been tested with MatLab 2015b on Ubuntu 15+ and Mac OSX 10.11.
-```
-  >>  demo_LiFE.m
-```
+C. Caiafa and F. Pestilli Sparse multiway decomposition for analysis and modelling of diffusion imagign and tractogrpahy.
 
+http://arxiv.org/abs/1505.07170

life/compute/feComputeEvidence.m → compute/feComputeEvidence.m

@@ -6,9 +6,9 @@
 % - The strength of evidence, namely the effect size of the lesion
 % - The Kullback-Leibler Divergence
 % - Jeffrey's Divergence
-% - The Earth Mover's distance
+% - The Eath Mover's distance
 %
-% Copyright (2016), Franco Pestilli, Indiana University, frakkopesto@gmail.com.
+% Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
 
 % Prepare the distribution of errors and the histograms describing the
 % erros.
@@ -18,8 +18,8 @@
 se.lesion.rmse.mean   = mean(se.lesion.rmse.all);
 
 % Histograms
-se.xrange(1) = min([se.lesion.rmse.all se.nolesion.rmse.all]);
-se.xrange(2) = max([se.lesion.rmse.all se.nolesion.rmse.all]);
+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);
@@ -45,15 +45,16 @@
 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 % Using Rubinov c-code fastest
+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 
+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
@@ -62,7 +63,7 @@
 clear tmp_emp
 
 % Strenght of evidence (effect size)
-fprintf('[%s] Computing the Strength of Evidence... ',mfilename)
+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;
@@ -75,7 +76,7 @@
 
 for inm = 1:se.s.nmontecarlo
     fprintf('.')
-    for ibt = 1:se.s.nboots
+    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

life/compute/feFitModel.m → compute/feFitModel.m

@@ -1,11 +1,5 @@
-function [fit, w, R2] = feFitModel(varargin)
-M = varargin{1};
-dSig = varargin{2};
-fitMethod = varargin{3};
-Niter = varargin{4};
-preconditioner = varargin{5};
-
-
+function [fit, w, R2] = feFitModel(M,dSig,fitMethod)
+% 
 % feFitModel() function in LiFE but restricted to the
 % 
 % BBNNLS algorithm and using the Factorization model.
@@ -54,30 +48,10 @@
 % 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.
 
-% Check for mexfiles and generate them if necessary
-% Mtransp_times_b function
-checkMexCompiled('-largeArrayDims', '-output', 'Mtransp_times_b', '-DNDEBUG','Mtransp_times_b.c', 'Mtransp_times_b_sub.c')
-% M_times_w function
-checkMexCompiled('-largeArrayDims', '-output', 'M_times_w', '-DNDEBUG', 'M_times_w.c', 'M_times_w_sub.c')
-% compute_diag function
-checkMexCompiled('-largeArrayDims', '-output', 'compute_diag', '-DNDEBUG', 'compute_diag.c', 'compute_diag_sub.c')
-
-if nargin <6 % no initial w0 is provided
-    [nFibers] = size(M.Phi,3);
-    w0 = zeros(nFibers,1);
-else
-    w0 = varargin{6};
-end
-
-
-if strcmp(preconditioner,'preconditioner') 
-    [nFibers] = size(M.Phi,3); %feGet(fe,'nfibers');
-    h = compute_diag(M.Phi.subs(:,1), M.Phi.subs(:,3), M.Phi.vals, M.DictSig,nFibers);
-    vals = M.Phi.vals./h(M.Phi.subs(:,3));
-    M.Phi = sptensor(M.Phi.subs,vals,size(M.Phi));
-end
-
+mycomputer = computer();
+release = version('-release');
 
 switch fitMethod
    case {'bbnnls'}
@@ -89,16 +63,19 @@
     tic
     fprintf('\nLiFE: Computing least-square minimization with BBNNLS...\n')
     opt = solopt;
-    opt.maxit = Niter;
+    opt.maxit = 5000;
     opt.use_tolo = 1;
-    opt.tolg = 1e-5;
-
-    out_data = bbnnls(M,dSig,w0,opt);
 
-    if strcmp(preconditioner,'preconditioner')
-        out_data.x = out_data.x./h;
+    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)