diff --git a/.gitignore b/.gitignore
index 3084eda..6d725d9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1 @@
-.m~
\ No newline at end of file
+*.m~
diff --git a/life/compute/feComputeEvidence.m b/life/compute/feComputeEvidence.m
index d81890b..c317fb2 100755
--- a/life/compute/feComputeEvidence.m
+++ b/life/compute/feComputeEvidence.m
@@ -1,4 +1,4 @@
-function se = feComputeEvidence(rmse1,rmse2)
+function se = feComputeEvidence(rmse1,rmse2,nbins)
 % Computes a series of distance metrics between two RMSE distribtions.
 %
 % Compute summary statistics on to characterize the lesion,:
@@ -20,7 +20,7 @@
 % 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.nbins     = 60;
+se.nbins     = nbins;
 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);
diff --git a/life/compute/feComputeEvidence_fixed.m b/life/compute/feComputeEvidence_fixed.m
new file mode 100644
index 0000000..ebf53b4
--- /dev/null
+++ b/life/compute/feComputeEvidence_fixed.m
@@ -0,0 +1,277 @@
+function se = feComputeEvidence_fixed(rmse1, rmse2, nbins, nmc, nboots)
+% 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 Earth Mover's distance
+%
+% Copyright (2020), Brent McPherson & Franco Pestilli, Indiana University
+%
+
+% fill in default parameters
+
+if(~exist('nbins', 'var') || isempty(nbins))
+    nbins = 128;
+end
+
+if(~exist('nmc', 'var') || isempty(nmc))
+    nmc = 5;
+end
+
+if(~exist('nboots', 'var') || isempty(nboots))
+    nboots = 250; % this almost certainly has to be higher, low for test
+end
+
+% pull the size of a distribution
+nobs1 = size(rmse1, 2); % nolesion
+nobs2 = size(rmse2, 2); % lesion
+
+% build the combined null distribution
+nrmse = [ rmse1, rmse2 ];
+
+% the min/max observed values for SOE / d-prime histogram computation
+minx = floor(mean(rmse1) - mean(rmse1)*.05);
+maxx = ceil(mean(rmse2) + mean(rmse2)*.05);
+
+% prepare the distribution of errors and the histograms describing the errors.
+se.nolesion.rmse.all = rmse1;
+se.nolesion.rmse.mean = mean(se.nolesion.rmse.all);
+
+se.lesion.rmse.all = rmse2;
+se.lesion.rmse.mean = mean(se.lesion.rmse.all);
+
+% compute histograms
+
+% define the bins based on the inputs / argument
+se.nbins = nbins;
+se.xrange = minmax([ se.lesion.rmse.all se.nolesion.rmse.all ]);
+se.bins = linspace(se.xrange(1), se.xrange(2), se.nbins);
+
+% store params for reference
+se.params.nmontecarlo = nmc;
+se.params.nboots = nboots;
+
+% compute the histograms on the data
+% the se.bins center have to be wrapped with inf to actually match hist
+% https://www.mathworks.com/matlabcentral/answers/377930-histcounts-error-in-place-of-histc
+[ se.lesion.hist, se.lesion.xhist ] = histcounts(se.lesion.rmse.all, [ se.bins inf ], 'Normalization', 'probability');
+[ se.nolesion.hist, se.nolesion.xhist ] = histcounts(se.nolesion.rmse.all, [ se.bins inf ], 'Normalization', 'probability');
+
+% drop inf from end to keep bin centers the same size, b/c this is "fixed"
+se.lesion.xhist = se.lesion.xhist(1:end-1);
+se.nolesion.xhist = se.nolesion.xhist(1:end-1);
+
+%% compute the true differences between the rmse / histograms
+
+% dissimilarity between raw RMSE values
+se.dis.name = 'Dissimilarity: https://nikokriegeskorte.org/category/representational-similarity-analysis/';
+se.dis.mean = pdist2(rmse1, rmse2, 'correlation');
+
+% cosine dissimilarity between raw RMSE values
+se.cos.name = 'Cosine Similarity: https://en.wikipedia.org/wiki/Cosine_similarity';
+se.cos.mean = pdist2(rmse1, rmse2, 'cosine');
+
+% Kullback-Leibler Divergence
+se.kl.name = sprintf('Kullback–Leibler divergence: http://en.wikipedia.org/wiki/Kullback-Leibler_divergence');
+tkl = se.nolesion.hist .* log2( (se.nolesion.hist) ./ (se.lesion.hist + eps) );
+se.kl.mean = nansum(tkl);
+clear tkl
+
+% Jeffrey's divergence
+se.j.name = sprintf('Jeffrey''s divergence: http://en.wikipedia.org/wiki/Divergence_(statistics)');
+tjd = 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(tjd); clear tjd
+
+% Earth Mover's Distance:
+% Note: This can be very slow and my require large amounts of memory for more than 1000 voxels
+se.em.name = sprintf('Earth Mover''s distance: http://en.wikipedia.org/wiki/Earth_mover''s_distance');
+
+try % Using Rubinov c-code fastest
+    
+    % estimate the distance between the bin centers
+    eDist = pdist2(se.nolesion.xhist', se.lesion.xhist', @(xi, xj) abs(xi-xj));
+    
+    % compute the emd
+    tmp_em = emd_mex(se.nolesion.hist, se.lesion.hist, eDist);
+
+catch % else
+
+    % compute EMD using slower MATLAB fxn
+    [ ~, tmp_em ] = emd(se.nolesion.xhist',se.lesion.xhist',se.nolesion.hist',se.lesion.hist',@gdf);
+    
+end
+se.em.mean = tmp_em;
+clear tmp_emp
+
+%% add the bootstrapped null / se measures for all metrics
+
+% define null histograms for SOE estimate
+mhist1 = nan(nbins, nmc);
+mhist2 = nan(nbins, nmc);
+
+% preallocate test vectors for every measure (SE)
+tcor = nan(nboots, nmc);
+tcos = nan(nboots, nmc);
+tkld = nan(nboots, nmc);
+tjfd = nan(nboots, nmc);
+temd = nan(nboots, nmc);
+tdp1 = nan(nboots, nmc);
+tdp2 = nan(nboots, nmc);
+
+% preallocate null vectors for every measure (pval)
+ncor = nan(nboots, nmc);
+ncos = nan(nboots, nmc);
+nkld = nan(nboots, nmc);
+njfd = nan(nboots, nmc);
+nemd = nan(nboots, nmc);
+ndpd = nan(nboots, nmc);
+
+for mc = 1:nmc
+    fprintf('.')
+    for boot = 1:nboots
+        
+        % pull random samples from each distribution for SE estimates
+        trmse1 = randsample(rmse1, nobs1, true); % nolesion
+        trmse2 = randsample(rmse2, nobs2, true); % lesion
+        
+        % pull random samples from combined distributions for p-value estimates
+        nrmse1 = randsample(nrmse, nobs1, true); % nolesion
+        nrmse2 = randsample(nrmse, nobs2, true); % lesion
+        
+        % build the boostrapped histograms for test and null
+
+        % compute the test histograms w/ the estimated bins
+        [ t1hist, t1xhist ] = histcounts(trmse1, [ se.bins inf ], 'Normalization', 'probability');
+        [ t2hist, t2xhist ] = histcounts(trmse2, [ se.bins inf ], 'Normalization', 'probability');
+        
+        % fix the "fixed" bins
+        t1xhist = t1xhist(1:end-1);
+        t2xhist = t2xhist(1:end-1);
+        
+        % build the the combined null histogram
+        
+        % compute the null histograms w/ the estimated bins
+        [ n1hist, n1xhist ] = histcounts(nrmse1, [ se.bins inf ], 'Normalization', 'probability');
+        [ n2hist, n2xhist ] = histcounts(nrmse2, [ se.bins inf ], 'Normalization', 'probability');
+
+        % fix the "fixed" bins
+        n1xhist = n1xhist(1:end-1);
+        n2xhist = n2xhist(1:end-1);
+
+        % estimate each SE / null throw on raw RMSE values
+
+        % SOE / d-prime test averages
+        tdp1(boot, mc) = mean(trmse1);
+        tdp2(boot, mc) = mean(trmse2);
+        
+        % SOE / d-prime null difference
+        ndp1 = mean(nrmse1);
+        ndp2 = mean(nrmse2);
+        ndpd(boot, mc) = ndp1 - ndp2;
+                
+        % estimate the test / null dissimilarity between rmse       
+        tcor(boot, mc) = pdist2(trmse1, trmse2, 'correlation');
+        ncor(boot, mc) = pdist2(nrmse1, nrmse2, 'correlation');
+        
+        % estimate the test / null cosine distance between rmse
+        tcos(boot, mc) = pdist2(trmse1, trmse2, 'cosine');
+        ncos(boot, mc) = pdist2(nrmse1, nrmse2, 'cosine');
+        
+        % build the SE / null estimates on the resampled / null histograms
+        
+        % kl divergence
+        tkld(boot, mc) = nansum(t1hist .* log2( (t1hist) ./ (t2hist + eps) ));
+        nkld(boot, mc) = nansum(n1hist .* log2( (n1hist) ./ (n2hist + eps) ));
+        
+        % jeffery's divergence
+        tjfd(boot, mc) = nansum(t1hist .* log2( (t1hist) ./ ((t2hist + t2hist + eps)./2)  ) + ...
+            t2hist .* log2( (t2hist) ./ ((t2hist + t2hist + eps)./2) ));
+        njfd(boot, mc) = nansum(n1hist .* log2( (n1hist) ./ ((n2hist + n2hist + eps)./2)  ) + ...
+            n2hist .* log2( (n2hist) ./ ((n2hist + n2hist + eps)./2) ));
+        
+        % Earth Mover's Distance; cased to failover to matlab fxb if MEX fxn fails
+        try
+            
+            % estimate EMD on resampled / null
+            temd(boot, mc) = emd_mex(t1hist, t2hist, eDist);
+            nemd(boot, mc) = emd_mex(n1hist, n2hist, eDist);
+            
+        % pass to slower full matlab fxn if MEX fails
+        catch
+            
+            % slower MATLAB fxn
+            [ ~, temd(boot, mc) ] = emd(t1xhist', t2xhist', t1hist', t2hist', @gdf);
+            [ ~, nemd(boot, mc) ] = emd(n1xhist', n2xhist', n1hist', n2hist', @gdf);
+            
+        end
+        
+    end
+    
+    % pull the estimated bins for SOE histogram
+    mbins = linspace(minx, maxx, nbins);
+    
+    % compute and normalize histogram of unlesioned rmse
+    [ mhist1(:, mc), mxhist1 ] = histcounts(tdp1(:, mc), [ mbins inf ], 'Normalization', 'probability');
+    [ mhist2(:, mc), mxhist2 ] = histcounts(tdp2(:, mc), [ mbins inf ], 'Normalization', 'probability');
+    
+    % fix the "fixed" bins
+    mxhist1 = mxhist1(1:end-1);
+    mxhist2 = mxhist2(1:end-1);
+    
+end
+disp(' done.')
+
+%% store the permutation resluts in output struct
+
+% permuted dissimilarity measure
+se.dis.se = mean(std(tcor, [], 1));
+se.dis.pval = min(sum(ncor > se.dis.mean) / nboots);
+
+% permuted cosine distance measure
+se.cos.se = mean(std(tcos, [], 1));
+se.cos.pval = min(sum(ncos > se.cos.mean) / nboots);
+
+% permuted kl values
+se.kl.se = mean(std(tkld, [], 1));
+se.kl.pval = min(sum(nkld > se.kl.mean) / nboots);
+
+% permuted jd values
+se.j.se = mean(std(tjfd, [], 1));
+se.j.pval = min(sum(njfd > se.j.mean) / nboots);
+
+% permuted emd values
+se.em.se = mean(std(temd, [], 1));
+se.em.pval = min(sum(nemd > se.em.mean) / nboots);
+
+% pull the tested SOE mean / std across MCs
+smn = diff([ mean(tdp1, 1); mean(tdp2, 1) ]);
+sse = sqrt(sum([ std(tdp1, [], 1); std(tdp2, [], 1) ].^2, 1));
+
+% collapse SOE across MC runs
+se.s.mean = mean(smn./sqrt(sse).^2);
+se.s.std = std(smn./sqrt(sse).^2);
+se.s.pval = min(sum(ndpd > smn) ./ nboots);
+
+% build mean histogram of SOE distributions
+mn_mhist1 = mean(mhist1, 2);
+mn_mhist2 = mean(mhist2, 2);
+
+% add 2 std around mean histogram estimates (confidence itervals)
+ci_mhist1 = [mn_mhist1, mn_mhist1] + 2*[-std(mhist1, [], 2), std(mhist1, [], 2)];
+ci_mhist2 = [mn_mhist2, mn_mhist2] + 2*[-std(mhist2, [], 2), std(mhist2, [], 2)];
+
+% store lesioned / unlesioned histograms
+se.s.min_x = minx;
+se.s.max_x = maxx;
+se.s.unlesioned_mn = mn_mhist1;
+se.s.unlesioned_ci = ci_mhist1;
+se.s.unlesioned_xbins = mxhist1';
+se.s.lesioned_mn = mn_mhist2;
+se.s.lesioned_ci = ci_mhist2;
+se.s.lesioned_xbins = mxhist2';
+
+end
diff --git a/life/compute/feComputeEvidence_norm.m b/life/compute/feComputeEvidence_norm.m
index 5398697..2f51eca 100644
--- a/life/compute/feComputeEvidence_norm.m
+++ b/life/compute/feComputeEvidence_norm.m
@@ -1,4 +1,4 @@
-function se = feComputeEvidence_norm(rmse1,rmse2)
+function se = feComputeEvidence_norm(rmse1, rmse2, nbins)
 % Computes a series of distance metrics between two RMSE distribtions.
 %
 % Compute summary statistics on to characterize the lesion,:
@@ -9,24 +9,33 @@
 % - The Earth Mover's distance
 %
 % Copyright (2016), Franco Pestilli, Indiana University, frakkopesto@gmail.com.
+%
 
-% Prepare the distribution of errors and the histograms describing the
-% erros.
+% Prepare the distribution of errors and the histograms describing the errors.
 se.nolesion.rmse.all  = rmse1;
-se.lesion.rmse.all    = rmse2;
 se.nolesion.rmse.mean = mean(se.nolesion.rmse.all);
+
+se.lesion.rmse.all    = rmse2;
 se.lesion.rmse.mean   = mean(se.lesion.rmse.all);
 
 % Histograms
+
+% define the bins based on the inputs / argument
+se.nbins     = nbins;
 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.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);
+
+% compute the histograms on the data
+[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);
+
+% normalize the histograms
+se.lesion.hist     = se.lesion.hist ./ sum(se.lesion.hist);
 se.nolesion.hist   = se.nolesion.hist ./ sum(se.nolesion.hist);
 
+%% compute the true differences between the histograms
+
 % 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) );
@@ -62,7 +71,7 @@
 clear tmp_emp
 
 % Strenght of evidence (effect size)
-fprintf('[%s] Computing the Strength of Evidence... ',mfilename)
+fprintf('[%s] Computing the Strength of Evidence... ', 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;
@@ -76,10 +85,33 @@
 for inm = 1:se.s.nmontecarlo
     fprintf('.')
     for 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));
+        nullDistributionW(ibt, inm) = mean(randsample(se.nolesion.rmse.all, sizeunlesioned, true));      
+        nullDistributionWO(ibt, inm) = mean(randsample(se.lesion.rmse.all, sizeunlesioned, true));
+    end
+    
+    % build the combined null distribution
+    ndist = [ se.lesion.rmse.all, se.nolesion.rmse.all ];
+    
+    % build the true difference to compare after null
+    empDiff = se.lesion.rmse.mean - se.nolesion.rmse.mean;
+    
+    % for every boot strap
+    for ibt = 1:se.s.nboots
+        
+        % pull a random sample from combined distributions
+        tnullDistributionW  = mean(randsample(ndist, sizeunlesioned, true));      
+        tnullDistributionWO = mean(randsample(ndist, sizeunlesioned, true));
+        
+        % pull the null difference for comparison
+        diffDist(ibt) = tnullDistributionW - tnullDistributionWO;
+        
     end
     
+    % get the count of nulls that are greater than observed
+    se.s.pval(inm) = sum(diffDist > empDiff) / se.s.nboots; 
+    % find the code in CCA and copy it here to estimate this.
+    %pval(ii) = sum(grotRr(1, ii, :) > cca.cca.hocorrs(ii)) / Nperm;    
+    
     % Distribution unlesioned
     [y(:,inm),xhis] = hist(nullDistributionW(:,inm),linspace(min_x,max_x,se.s.nbins));
     y(:,inm) = y(:,inm)./sum(y(:,inm));
@@ -87,8 +119,9 @@
     % Distribution lesioned
     [woy(:,inm),woxhis] = hist(nullDistributionWO(:,inm),linspace(min_x,max_x,se.s.nbins));
     woy(:,inm) = woy(:,inm)./sum(woy(:,inm));
-end
     
+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); ...
@@ -109,4 +142,66 @@
 se.s.min_x   = min_x;
 se.s.max_x = max_x;
 
+%% add the bootstrapped null / std measures for all metrics
+
+% pull the size of a distribution
+nobs1 = size(rmse1, 2);
+nobs2 = size(rmse2, 2);
+
+% build the combined null distribution
+nrmse = [ rmse1, rmse2 ];
+
+% preallocate null vector for every measure
+ncor = nan(se.s.nboots, 1);
+ncos = nan(se.s.nboots, 1);
+nkld = nan(se.s.nboots, 1);
+
+for ibt = 1:se.s.nboots
+    
+    % pull random samples from each distribution
+    trmse1 = randsample(rmse1, nobs1, true); % nolesion
+    trmse2 = randsample(rmse2, nobs2, true); % lesion
+    
+    % pull random samples from combined distributions
+    nrmse1 = randsample(nrmse, nobs1, true);
+    nrmse2 = randsample(nrmse, nobs2, true);
+    
+    % estimate each null throw on raw RMSE values
+    ncor(ibt) = pdist2(nrmse1, nrmse2, 'correlation');
+    ncos(ibt) = pdist2(nrmse1, nrmse2, 'cosine');
+    
+    % build the null histogram
+    
+    % define the bins based on the inputs / argument
+    tse.xrange = minmax([ trmse1 trmse2 ]);
+    tse.bins = linspace(tse.xrange(1), tse.xrange(2), nbins);
+    
+    % compute the histograms on the data
+    [ tse.lesion.hist, tse.lesion.xhist ] = hist(trmse2, nbins);
+    [ tse.nolesion.hist, tse.nolesion.xhist ] = hist(trmse1, nbins);
+    
+    % normalize the histograms
+    tse.lesion.hist = tse.lesion.hist ./ sum(tse.lesion.hist);
+    tse.nolesion.hist = tse.nolesion.hist ./ sum(tse.nolesion.hist);
+    
+    % kl divergence estimated on null
+    tkld = tse.nolesion.hist .* log2( (tse.nolesion.hist) ./ (tse.lesion.hist + eps) );
+    nkld(ibt) = nansum(tkld);
+    
+end
+
+% dissimilarity between raw RMSE values
+se.dis.name = 'Dissimilarity';
+se.dis.mean = pdist2(rmse1, rmse2, 'correlation');
+se.dis.pval = 1 - (sum(ncor > se.dis.mean) / se.s.nboots);
+
+% cosine dissimilarity between raw RMSE values
+se.cos.name = 'Cosine Distance';
+se.cos.mean = pdist2(rmse1, rmse2, 'cosine');
+se.cos.pval = 1 - (sum(ncos > se.cos.mean) / se.s.nboots);
+
+se.nkl.name = "New KL Divergence";
+se.nkl.mean = se.kl.mean;
+se.nkl.pval = 1 - (sum(nkld > se.kl.mean) / se.s.nboots);
+
 end
diff --git a/life/compute/feFitModel.m b/life/compute/feFitModel.m
index 54f7edb..a9b4539 100755
--- a/life/compute/feFitModel.m
+++ b/life/compute/feFitModel.m
@@ -56,12 +56,14 @@
 
 
 % 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 ~isdeployed
+    % 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')
+end
 
 if nargin <6 % no initial w0 is provided
     [nFibers] = size(M.Phi,3);
diff --git a/life/fe/feBuildDictionaries.m b/life/fe/feBuildDictionaries.m
index 3f84152..2db5c9a 100755
--- a/life/fe/feBuildDictionaries.m
+++ b/life/fe/feBuildDictionaries.m
@@ -14,7 +14,7 @@
 %  Cesar F. Caiafa ccaiafa@gmail.com
 
 tic
-fprintf(['\n[%s] Computing demeaned and non-demeaned difussivities dictionaries in a (',num2str(Nphi),'x',num2str(Ntheta),')-grid', ' ...'],mfilename); 
+fprintf(['\n[%s] Computing demeaned and non-demeaned diffusivities dictionaries in a (',num2str(Nphi),'x',num2str(Ntheta),')-grid', ' ...'],mfilename); 
 fprintf('took: %2.3fs.\n',toc)
 
 % Compute orientation vectors
@@ -43,29 +43,125 @@
 DictSig = zeros(nBvecs,Norient); % Initialize Signal Dictionary matrix
 %DictTensors = zeros(9,Norient); % Initialize Tensors Dictionary matrix
 
-D = diag(fe.life.modelTensor); % diagonal matix with diffusivities
-
-%[ shells, ~, shellindex ] = unique(round(bvals)); % round to get around noise in shells - FRAGILE FIX
+% catch Kurtosis estimates for debugging
+%KDict = zeros(nBvecs,Norient); 
+%KDictSig = zeros(nBvecs,Norient); 
 
+% pull the shell information
 ubv = feGet(fe, 'nshells');
 ubi = feGet(fe, 'shellindex');
+ubl = unique(ubi);
+
+% pull diffusion tensor
+dt = feGet(fe, 'model tensor');
+
+% pull kurtosis fit
+kt = feGet(fe, 'model Kurtosis');
+
+% if data is single shell
+if ubv == 1 || all(isnan(kt))
+    
+    % akc is invalid, set to nan
+    akc = nan(nBvecs,1); 
+
+% if data is multishell
+else
+    
+    % preallocate akc parameters to 0
+    akc = zeros(nBvecs,1);
+            
+    % build apparent kurtosis coefficient - akc
+    for i=1:ubv
+        
+        % pull tensor parameters for hard indexing of equations
+        sdt = dt(i,:);
+        
+        % find the indices for the shell
+        si = ubi == ubl(i);
+        
+        % mean diffusivity of tensor from forward model
+        md = mean(sdt);
+        
+        % apparent diffusion coefficient - matches dipy
+        adc = bvecs(si,1) .* bvecs(si,1) * sdt(1) + ...
+            2 * bvecs(si,1) .* bvecs(si,2) * sdt(2) + ...
+            bvecs(si,2) .* bvecs(si,2) .* sdt(3) + ...
+            2 * bvecs(si,1) .* bvecs(si,3) * 0 + ... % only store primary eigenvalues in
+            2 * bvecs(si,2) .* bvecs(si,3) * 0 + ... % LiFE prediction. These are 0.
+            bvecs(si,3) .* bvecs(si,3) * 0;
+        
+        % apparent diffusion variance - matches dipy
+        adv = ...
+            bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* kt(1) + ...       % xxxx
+            bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(2) + ...       % yyyy
+            bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(3) + ...       % zzzz
+            4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* kt(4) + ...   % xxxy
+            4 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* kt(5) + ...   % xxxz
+            4 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* kt(6) + ...   % xyyy
+            4 * bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(7) + ...   % xzzz
+            4 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(8) + ...   % yyyz
+            4 * bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* bvecs(si,3) .* kt(9) + ...   % yzzz
+            6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* kt(10) + ...  % xxyy
+            6 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,3) .* bvecs(si,3) .* kt(11) + ...  % xxzz
+            6 * bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(12) + ...  % yyzz
+            12 * bvecs(si,1) .* bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* kt(13) + ... % xxyz
+            12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,2) .* bvecs(si,3) .* kt(14) + ... % xyyz
+            12 * bvecs(si,1) .* bvecs(si,2) .* bvecs(si,3) .* bvecs(si,3) .* kt(15);      % xyzz
+        
+        % zero out noisy (bad) adc estimates
+        adc(adc < 0) = 0;
+        
+        % estimate apparent kutosis coefficient - matches dipy
+        akc(si) = adv .* ((md ./ adc).^2);
+        
+        % zero out noisy (bad) akc estimates
+        akc(akc < -3/7) = -3/7;
+        
+    end
+    
+end
 
 % Compute each dictionary column for a different kernel orientation
 for j=1:Norient
-    [Rot,~, ~] = svd(orient(:,j)); % Compute the eigen vectors of the kernel orientation
     
-    Q = Rot*D*Rot';
-    %DictTensors(:,j) = Q(:);
-    Dict(:,j) = exp(- bvals .* diag(bvecs*Q*bvecs')); % Compute the signal contribution of a fiber in the kernel orientation divided S0
+    % Compute the eigen vectors of the kernel orientation
+    [Rot,~, ~] = svd(orient(:,j));
     
     % for every shell
-    for k=1:size(ubv, 1) 
-        si = ubi == ubv(k); % find the indices for the shell
-        DictSig(si,j) = Dict(si,j) - median(Dict(si,j)); % demedianed signal by shell (used to demean)
+    for k=1:ubv
+        
+        % find the indices for the shell
+        si = ubi == ubl(k); 
+        
+        % create diagonal matix with diffusivities for current shell
+        % this assumes tensor fits for shell are entered in the order this will parse them in
+        D = diag(dt(k,:));
+        
+        % estimate Q for tensor values in shell
+        Q = Rot*D*Rot';
+        
+        % STORE SEPARATE? THIS WILL GREATLY INCREASE THE SIZE
+        
+        % Compute the signal contribution of a fiber in the kernel orientation divided S0
+        if all(kt == 0) || all(isnan(kt)) % if kurtosis is zeroed b/c it's invalid / requested
+            Dict(si,j)  = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)')); 
+        else % if kurtosis parameters are passed
+            Dict(si,j) = exp(-bvals(si) .* diag(bvecs(si,:)*Q*bvecs(si,:)') + (-bvals(si).^2 .* diag(bvecs(si,:)*Q*bvecs(si,:)').^2 .* akc(si))/6);
+        end
+        
+        % demeaned signal by shell
+        DictSig(si,j) = Dict(si,j) - mean(Dict(si,j)); 
+        %KDictSig(si,j) = KDict(si,j) - mean(KDict(si,j));
+        
     end
     
 end
 
 fe = feSet(fe,'dictionary parameters',{Nphi,Ntheta,orient,Dict,DictSig});
+fe = feSet(fe,'akc',akc);
+
+% hard add kurtosis dictionaries / akc for debugging
+%fe.life.M.KDict = KDict;
+%fe.life.M.KDictSig = KDictSig;
 
 end
\ No newline at end of file
diff --git a/life/fe/feConnectomeEncoding.m b/life/fe/feConnectomeEncoding.m
index f4e3f38..e29bb66 100755
--- a/life/fe/feConnectomeEncoding.m
+++ b/life/fe/feConnectomeEncoding.m
@@ -62,7 +62,9 @@
 nBatch       = max([5,ceil(nTotalNodes/nNodesMax)]); % Number of batches
 nFib_Batch   = ceil(nFibers/nBatch); % Number of fibers per batch
 
-disp(['Max number of nodes per batch',num2str(nNodesMax),', Total number of nodes = ',num2str(nTotalNodes),'.  Building the model in ',num2str(nBatch),' batches'])
+disp(['Max number of nodes per batch: ',num2str(nNodesMax)]);
+disp(['Total number of nodes = ',num2str(nTotalNodes)]);
+disp(['Building the model in ',num2str(nBatch),' batches.']);
 
 fprintf('\n');
 for n=1:nBatch
diff --git a/life/fe/feConnectomeInit.m b/life/fe/feConnectomeInit.m
index 6a36de7..9d2b738 100755
--- a/life/fe/feConnectomeInit.m
+++ b/life/fe/feConnectomeInit.m
@@ -1,4 +1,4 @@
-function fe = feConnectomeInit(dwiFile,fgFileName,feFileName,savedir,dwiFileRepeated,anatomyFile,varargin)
+function fe = feConnectomeInit(dwiFile,fgFileName,dwiFileRepeated,anatomyFile,N,tensor,kurtosis,reset)
 % Initialize a new connectome (fe) structure. 
 %
 %    fe = feConnectomeInit(dwiFile,fgFileName,feFileName,savedir,dwiFileRepeated,anatomyFile,varargin);
@@ -13,34 +13,54 @@
 
 %feOpenLocalCluster
 
+if(~exist('N', 'var') || isempty(N))
+    N = 360;
+end
+
+if(~exist('tensor', 'var') || isempty(tensor))
+    axialDiffusion = 1; radialDiffusion = 0;
+    tensor = [ axialDiffusion radialDiffusion radialDiffusion ];
+end
+
+if(~exist('kurtosis', 'var') || isempty(kurtosis))
+    kurtosis = zeros(15,1);
+end
+
+if(~exist('reset', 'var') || isempty(reset))
+    reset = true;
+end
+
 % Intialize the fe structure.
 fe = feCreate;
 
-% Set the based dir for fe, this dire will be used 
-if notDefined('savedir')
-    if isstruct(fgFileName)
-        savedir = fullfile(fileparts(fgFileName.name),'life');
-    else
-        savedir = fullfile(fileparts(fgFileName),'life');
-    end
-end
-fe = feSet(fe,'savedir',savedir);
+% Set the based dir for fe, this dir will be used 
+% if notDefined('savedir')
+%     if isstruct(fgFileName)
+%         savedir = fullfile(fileparts(fgFileName.name),'life');
+%     else
+%         savedir = fullfile(fileparts(fgFileName),'life');
+%     end
+% end
+% fe = feSet(fe,'savedir',savedir);
 
 % Set the xforms (transformations from diffusion data to acpc)
 tempNi = niftiRead(dwiFile);
 fe = feSet(fe, 'img2acpc xform', tempNi.qto_xyz);
 fe = feSet(fe, 'acpc2img xform', inv(tempNi.qto_xyz));
-clear tempNi
 
-% Set up the fe name
-if isstruct(fgFileName),  n  = fgFileName.name;
-else                   [~,n] = fileparts(fgFileName);
-end
+% set the offset for files as well, then clear it
+%fe = feSet(fe, 'dwi offset', tempNi.qto_xyz(1:3,4)');
+clear tempNi
 
-if notDefined('feFileName')
-  feFileName = sprintf('%s-%s', datestr(now,30),n);
-end
-fe = feSet(fe, 'name',feFileName);
+% % Set up the fe name
+% if isstruct(fgFileName),  n  = fgFileName.name;
+% else                   [~,n] = fileparts(fgFileName);
+% end
+%
+% if notDefined('feFileName')
+%   feFileName = sprintf('%s-%s', datestr(now,30),n);
+% end
+% fe = feSet(fe, 'name',feFileName);
 
 %% Load a connectome from disk
 if isstruct(fgFileName), fg = fgFileName; clear fgFileName
@@ -60,9 +80,78 @@
   fe = feConnectomeSetDwi(fe,dwiFileRepeated,1);
 end
 
-%% Anatomy Install the path tot he anatomical high-resolution file.
+%% Anatomy Install the path to the anatomical high-resolution file.
 if ~notDefined('anatomyFile')
-  fe = feSet(fe,'anatomy file',anatomyFile);
+    
+    fe = feSet(fe,'anatomy file',anatomyFile);
+    
+    try
+        tempNi = niftiRead(anatomyFile);
+        fe = feSet(fe, 'anat2acpcxform', tempNi.qto_xyz);
+        fe = feSet(fe, 'acpc2anatxform', inv(tempNi.qto_xyz));
+        clear tempNi
+    catch
+        error('Anatomy file unable to be opened with niftiRead()');
+    end
+       
+    % add check to file headers for consistent orientation between dwi / anat
+    
+    % load minimum data for check - ignore repeat dwi for now
+    doff = feGet(fe, 'dwi offset');
+    aoff = feGet(fe, 'anat offset');
+    
+    % check x dimension is consistently signed
+    if ~(all([ doff(1) aoff(1) ] < 0) || all([ doff(1) aoff(1) ] > 0))
+        if reset
+            warning('POTENTIAL X-FLIP - bvecs flipped on x-axis to correct.');
+            bvecs = feGet(fe, 'bvecs');
+            bvecs(:,1) = -bvecs(:,1);
+            fe = feSet(fe, 'bvecs', bvecs);
+        else
+            warning('POTENTIAL X-FLIP - qoffset_x is different between the dwi and anatomy files.');
+        end
+    end
+    
+    % check y dimension is consistently signed
+    if ~(all([ doff(2) aoff(2) ] < 0) || all([ doff(2) aoff(2) ] > 0))
+        if reset
+            warning('POTENTIAL Y-FLIP - bvecs flipped on y-axis to correct.');
+            bvecs = feGet(fe, 'bvecs');
+            bvecs(:,2) = -bvecs(:,2);
+            fe = feSet(fe, 'bvecs', bvecs);
+        else
+            warning('POTENTIAL Y-FLIP - qoffset_y is different between the dwi and anatomy files.');
+        end
+    end
+    
+    % check z dimension is consistently signed
+    if ~(all([ doff(3) aoff(3) ] < 0) || all([ doff(3) aoff(3) ] > 0))
+        if reset
+            warning('POTENTIAL Z-FLIP - bvecs flipped on z-axis to correct.');
+            bvecs = feGet(fe, 'bvecs');
+            bvecs(:,3) = -bvecs(:,3);
+            fe = feSet(fe, 'bvecs', bvecs);
+        else
+            warning('POTENTIAL Z-FLIP - qoffset_z is different between the dwi and anatomy files.');
+        end
+    end
+    
+end
+
+%% check / fix the tensor model for multishell input
+
+nshells = feGet(fe, 'nshells');
+if (size(tensor, 1) == 1) && nshells > 1
+    tensorParams = repmat(tensor, [ nshells 1 ]);
+else
+    tensorParams = tensor;
+end
+
+% if data is single shell, zero out kurtosis model regardless of user input
+if nshells == 1
+    kurtosisParams = nan(15,1);
+else
+    kurtosisParams = kurtosis;
 end
 
 %% Precompute Dictionary of orientations and canonical demeaned signals
@@ -70,23 +159,13 @@
 % These numbers represents the steps in which the range [0,pi] is divided
 % for spherical coordinates phi: azimuth and theta: pi/2-elevetion
 
-% Set model for canonical tensor
-%tic
-if ~isempty(varargin)
-  N = varargin{1}(1);
-  axialDiffusion  = varargin{2}(1);
-  radialDiffusion = varargin{2}(2);  
-else % Default to stick and ball
-  N = 360;   
-  axialDiffusion  = 1;
-  radialDiffusion = 0;
-end
-dParms(1) =  axialDiffusion; 
-dParms(2) = radialDiffusion; 
-dParms(3) = radialDiffusion;
 Nphi = N;
 Ntheta = N;
-fe = feSet(fe,'model tensor',dParms);
+fe = feSet(fe,'model tensor',tensorParams);
+fe = feSet(fe,'model kurtosis',kurtosisParams);
+% if single shell data is passed (impossible to fit kurtosis), 
+% feBuildDictionaries will zero out any submitted kurtosis parameters
+
 fe = feBuildDictionaries(fe,Nphi,Ntheta);
 
 %% Encode Connectome in multidimensional tensor
diff --git a/life/fe/feConnectomeSetDwi.m b/life/fe/feConnectomeSetDwi.m
index c876175..a4d50ce 100755
--- a/life/fe/feConnectomeSetDwi.m
+++ b/life/fe/feConnectomeSetDwi.m
@@ -48,18 +48,18 @@
             dwiGet(dwi, 'diffusionimagenums'));
             
 % Handling multishell data. 
-% The original model was only able to handl single-shell. 
+% The original model was only able to handle single-shell. 
 % This new version will handle also multishell data.
 % 
 % To handle multishell data we will store the number of shells encoded
 % for the data (1,2,3, etc) and the indices to the bvecs for each shell 
 % (this will NOT assume that the number of BVECS is equal for each shell
-% (say 30 diffusion directiosn for each shell).
+% (say 30 diffusion directions for each shell).
 
 % We extract the bvalues and find the unique of them and assign them to each shell.
 bvals = feGet(fe,'bvals');
 [nshells, ~, shellindex] = unique(round(bvals));
-fe  = feSet(fe, 'nshells', nshells);
+fe  = feSet(fe, 'nshells', length(nshells));
 fe  = feSet(fe, 'shellindex', shellindex);
 
 dim = dwi.nifti.dim;
diff --git a/life/fe/feGet.m b/life/fe/feGet.m
index 2194e17..f14b5fd 100755
--- a/life/fe/feGet.m
+++ b/life/fe/feGet.m
@@ -506,13 +506,27 @@
     % This is the number of unique shells for multishell data
     %
     % val = feGet(fe,'nshells');
-    val = fe.life.shells_n;
+    
+    % if nshells isn't present, it's an old single shell structure
+    if ~isfield(fe.life, 'nshells')
+        val = 1; % fill in structure and set the value to 1
+        feSet(fe, 'nshells', val); % add it to the structure
+    else % otherwise pull the value set during encoding
+        val = fe.life.nshells;
+    end
 
   case {'shellindex','indextoeachshellbvecs'}
   % This is the index to the bvecs associated to each unique shell  
     %
     % val = feGet(fe,'nshshellindexells');
-    val = fe.life.shells_index;
+    
+    % if shellsindex isn't present, it's an old single shell structure
+    if ~isfield(fe.life, 'shells_index') 
+        val = ones(feGet(fe, 'nbvals'), 1); % fill in all ones for the index
+        feSet(fe, 'shellindex', val); % add it to the structure
+    else % otherwise pull the value set during encoding
+        val = fe.life.shells_index;
+    end
     
   case {'bvecs'}
     % Diffusion directions.
@@ -566,7 +580,8 @@
     voxelIndices = feGet(fe,'voxelsindices',varargin);
     val = fe.life.diffusion_signal_img(voxelIndices,:) - repmat(mean(fe.life.diffusion_signal_img(voxelIndices,:), 2),1,nBvecs);
     keyboard
-    % THis seems to be wrong
+    % THIS IS WRONG - DOES NOT DO MULTISHELL CORRECTLY
+    % BUT IS IT USED ANYWHERE?
     
   case {'b0signalimage','b0vox'}
     % Get the diffusion signal at 0 diffusion weighting (B0) for this voxel
@@ -929,11 +944,26 @@
     % dSig = feGet(fe,'dsigdemeaned',coords);
     nVoxels = feGet(fe,'nVoxels');
     nBvecs  = feGet(fe,'nBvecs');
+    shells  = feGet(fe,'shellindex');
+    nshell  = feGet(fe,'nshells');
+    ushell  = unique(shells);
+    
+    % this works now - data needs to be transposed from how it's stored to
+    % match the unwound data used during prediction
+    dSig = fe.life.diffusion_signal_img';
+    
+    % loop over shells for demeaing
+    for shell = 1:nshell
+        s = ushell(shell); % index into listed shell for good logic
+        dSigMean = mean(dSig(shells == s, :));
+        dSig(shells == s, :) = dSig(shells == s, :) - repmat(dSigMean, [ sum(shells == s) 1 ]);
+    end
+    val = reshape(dSig, [nVoxels*nBvecs,1]);
+    
+    % % original demeaning operation
+    %dSig = reshape(fe.life.diffusion_signal_img',[1,nVoxels*nBvecs]);
+    %val = (dSig - reshape(repmat(mean(reshape(dSig, nBvecs, nVoxels),1),nBvecs,1), size(dSig)))';
     
-    dSig = reshape(fe.life.diffusion_signal_img',[1,nVoxels*nBvecs]);
-    val     = (dSig - reshape(repmat( ...
-      mean(reshape(dSig, nBvecs, nVoxels),1),...
-      nBvecs,1), size(dSig)))';
     % Return a subset of voxels
     if ~isempty(varargin)
       % voxelIndices     = feGet(fe,'voxelsindices',varargin);
@@ -1140,7 +1170,7 @@
     % res = feGet(fe,'res sig full');
     % res = feGet(fe, 'res sig full',coords);
     % res = feGet(fe, 'res sig full',voxelIndex);
-    val = (feGet(fe,'dsig full')    - feGet(fe,'psig full')');
+    val = (feGet(fe,'dsig full') - feGet(fe,'psig full')');
     if ~isempty(varargin)
       % voxelIndices     = feGet(fe,'voxelsindices',varargin);
       % voxelRowsToKeep  = feGet(fe,'voxel rows',voxelIndices);
@@ -1310,11 +1340,11 @@
     val       = sqrt(mean((measured - predicted).^2,1));
     val       = val(feGet(fe,'voxelsindices',varargin));
    
-  case {'voxrmses0norm'}
+  case {'voxrmses0norm'} % CAN THIS BE OPTIMIZED FOR A MULTISHELL?
       % A volume of RMSE normalized by the S0 value in each voxel.
       rmse = feGet(fe,'vox rmse');
       s0   = feGet(fe,'b0signalimage')';
-      % Some voxels can have a S0=0. We replace the S0 in these voxles with
+      % Some voxels can have a S0=0. We replace the S0 in these voxels with
       % a NaN. 
       idx = (s0 == 0);
       rmse(idx) = nan(size(find(idx)));
@@ -1651,11 +1681,29 @@
         
         s0 = fe.life.s0;
         val = ones(nTheta,1)*s0' + D*B; % with iso          
-        
+  
+  case {'modelkurtosis'}
+      val = fe.life.modelKurtosis;
+  case {'akc'}
+      val = fe.life.M.akc;
+  case {'dwi2acpc'}
+      val = fe.life.xform.img2acpc;
+  case {'dwi2img'}
+      val = fe.life.xform.acpc2img;      
+  case {'anat2acpc'}
+      val = fe.life.xform.anat2acpc;
+  case {'anat2img'}
+      val = fe.life.xform.anat2img;      
+  case 'dwioffset'
+      tmp = feGet(fe, 'dwi2acpc');
+      val = tmp(1:3,4)';
+  case 'anatoffset'
+      tmp = feGet(fe, 'anat2acpc');
+      val = tmp(1:3,4)';
+      
   otherwise
     help('feGet')
     fprintf('[feGet] Unknown parameter << %s >>...\n',param);
-    keyboard
 end
 
 end % END MAIN FUNCTION
diff --git a/life/fe/feSet.m b/life/fe/feSet.m
index 9d5b352..ad6e2f5 100755
--- a/life/fe/feSet.m
+++ b/life/fe/feSet.m
@@ -90,7 +90,7 @@
     fe.life.bvals = val;
   case {'nshells','numberofshells'}
   % This is the number of unique shells for multishell data
-    fe.life.shells_n = val;
+    fe.life.nshells = val;
   case {'shellindex','indextoeachshellbvecs'}
   % This is the index to the bvecs associated to each unique shell  
     fe.life.shells_index = val;
@@ -170,7 +170,7 @@
     fe.rep.bvals = val;
   case {'nshellsrepeat','numberofshellsrepeat'}
   % This is the number of unique shells for multishell data
-    fe.rep.shells_n = val;
+    fe.rep.nshells = val;
   case {'shellindexrepeat','indextoeachshellbvecsrepeat'}
   % This is the index to the bvecs associated to each unique shell  
     fe.rep.shells_index = val;
@@ -217,6 +217,19 @@
     fe.life.M.tracts{Ntracts+1}.ind = find(val.index==0);
     fe.life.M.tracts{Ntracts+1}.name = 'not a tract'; 
             
+  case {'xformanat2acpc','anat2acpc','anat2acpcxform'}
+    fe.life.xform.anat2acpc = val;
+  case {'xformacpc2anat','acpc2anat','acpc2anatxform'}
+    fe.life.xform.acpc2anat = val;
+%   case 'dwioffset'
+%     fe.life.xform.dwi_offset = val;
+%   case 'anatoffset'
+%     fe.life.xform.anat_offset = val;
+  case {'modelkurtosis'}
+    fe.life.modelKurtosis = val;
+  case {'akc'}
+    fe.life.M.akc = val;
+
   otherwise
     error('Unknown parameter %s\n',param);
 end
diff --git a/life/utility/feSparseMatrixCoords.m b/life/utility/feSparseMatrixCoords.m
new file mode 100644
index 0000000..adf1ecf
--- /dev/null
+++ b/life/utility/feSparseMatrixCoords.m
@@ -0,0 +1,93 @@
+function [ out ] = feSparseMatrixCoords(M)
+%[ mat ] = feSparseMatrixCoords(Phi)
+%   Convert the ENCODE model object (Phi) and the dictionary signal (D) to 
+%       a set of coordinates to be used in the 2d sparse matrix version of 
+%       LiFE that incorporates streamline groupings during optimization.
+%
+% INPUTS:
+%   M   - The model object from an initialized ENCODE object
+%         example: M = feGet(fe, 'model');
+%
+% OUTPUT:
+%   out - The conversion of the 3d sparse tensor and dictionary to the 
+%         coordinates / values required for the 2d sparse matrix LiFE model
+%
+% Brent McPherson put this function together, but
+% Andy Womack and Dan McDonald did most of the driving.
+%
+% Brent McPherson, Indiana University (c) 2021
+%
+
+% pull the number of streamlines
+nfib = size(M.Phi, 3);
+
+% Each slice (streamline) of Phi needs to be multiplied by the dictionary
+% signal to create the 2d matrix for of each streamlines prediction. The
+% streamlines can be stacked into a single 2d matrix by adding a column
+% indicating which streamline the entries (subscripts, values) they
+% correspond to.
+
+fprintf('Converting %d streamlines in Phi to sparse matrix coordinates...\n', nfib);
+
+% initialize the output
+%out = [];
+out = cell(nfib, 1);
+
+% for every streamline
+for fib = 1:nfib
+    
+    % print out every 1000th streamline as an update
+    if mod(fib, 1000) == 0
+        fprintf('Running streamline %d\n', fib);
+    end
+    
+    % pull the slice of the model tensor for the streamline
+    slc = M.Phi(:,:,fib);
+    
+    % use the sptensor toolbox to multiply the slice by the dictionary for
+    % making each streamlines prediction
+    mm = ttm(slc, M.DictSig, 1);
+    
+    % create the index to identify the streamline in the new matrix
+    idx = ones(size(mm.vals)) * fib;
+    
+    % for this streamline, store the values from the sparse product:
+    % - x (dictionary) / y (voxel) subscripts into the sparse model matrix
+    % - the index of the streamline for these values
+    % - the value stored at each entry
+    iter = [ mm.subs, idx, mm.vals ];
+    
+    % append this streamline to the output matrix
+    %out = [ out; iter ];
+    out{fib} = iter;
+    
+end
+
+fprintf('Stacking the streamline data into sparse matrix coordinates...\n');
+
+% merge the output
+out = cat(1, out{:});
+
+%
+% concerns / things to note
+%
+% mm.subs is type double from ttm - why?
+% - this is dumb, these should be int - ought to fix
+% - indices are integers, why store them as doubles?
+% - BE CAREFUL WITH THE PRECISION WHEN SAVING THIS OUTPUT
+% -- because they're saved as doubles, the integer indices will round if
+%    the precision isn't high enough, resulting in errors.
+%
+% In theory, the sptensor toolbox should be able to transform the tensor
+% to a matrix with a function distributed in the tool. However, that 
+% doesn't appear to be the case (their documentation / notation is bad).
+%
+% Further, the tensor-matrix multiplication doesn't work for these types 
+% (handling sparse objects sparsely). So attempts to work on the 3d tensor 
+% directly resulted in comically large memory overflows. Which is weird, 
+% because they have a @sptensor/ttm function that ought to handle sparse 
+% tensors correctly. Maybe it worked in an older version of matlab?
+%
+
+end
+