function img = demo3D20
%demo3D20
%
% Trains a generative model of the framework using one diffeomorphic model

% Created: Gregory Johnson
%
% Copyright (C) 2013 Murphy Lab
% Lane Center for Computational Biology
% School of Computer Science
% Carnegie Mellon University
%
% Aug 26, 2013 Updated to display shape space and error associated with 
%               incomplete distance matricies
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published
% by the Free Software Foundation; either version 2 of the License,
% or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA.
%
% For additional information visit http://murphylab.web.cmu.edu or
% send email to murphy@cmu.edu

curr_path = which('demo3D20.m');
curr_path = curr_path(1:end-10);
cd(curr_path);

 frameworkflag = 0;

%icaoberg 27/7/2012
if nargin == 0 | nargin == 1
 if strcmpi( lower(frameworkflag), 'yes' )
    frameworkflag = true;
 else
    frameworkflag = false;
 end
end

if ~isa( frameworkflag, 'logical' )
   warning( 'Input variable framework flag must be boolean. Exiting demo.' );
   return
end



pattern = 'LAM';
directory = [ pwd filesep 'LAM' ];
param.protein.class = 'lysosome';
param.protein.name = 'lamp2';
param.nucleus.name = 'lamp2';
param.cell.model = 'lamp2';


tarball = [ '3DHela_' pattern '.tgz' ];
url = 'http://murphylab.web.cmu.edu/data/Hela/3D/multitiff';
   
if ~exist( [pwd filesep tarball] )
 disp('Downloading dataset');
 urlwrite([url filesep tarball], tarball);
 disp('Extracting files');
 untar( tarball );
else
    disp('Pattern found on disk. Extracting files from existing tarball');
    if ~exist( directory )
        disp('Extracting files');
        system(['tar -xvf ' tarball] );
    else
        disp('Directory found. Skipping decompression');
    end
end

dimensionality = '3D';
answer = false;


for i = 16:50
    delete([directory filesep 'cell' num2str(i) '_*']);
end

dna = [directory filesep 'cell*ch0*.tif'] ;
cell = [directory filesep 'cell*ch1*.tif'];
protein = [directory filesep 'cell*ch2*.tif'];

param.masks = [directory filesep 'cell*mask*.tif'];

% param.masks = [directory filesep 'cell*mask*.tif'];

% generic model options
% ---------------------
% model.name                (optional) Holds the name of the model. Default is empty.
param.model.name = 'all';

% model.id                  (optional) Holds the id of the model. Default is empty.
try
   [status,result] = system( 'uuidgen' ); 
   param.model.id = result(1:end-1);
catch
   param.model.id = num2str(now);
end

% model.filename            Holds the output filename.
param.model.filename = [ lower(pattern) '.xml' ];

% nucleus.id                (optional) Holds the id of the nuclear model. Default is empty.
param.nucleus.type = 'diffeomorphic';
try
   [status,result] = system( 'uuidgen' ); 
   param.nucleus.id = result(1:end-1);
catch
   param.nucleus.id = num2str(now);
end

% cell shape model options
% ------------------------
% cell.type                 Holds the cell model type. Default is "ratio".
param.cell.type = 'diffeomorphic';

% cell.id                   (optional) Holds the id the cell model. Default is empty.
try
   [status,result] = system( 'uuidgen' ); 
   param.cell.id = result(1:end-1);
catch
   param.cell.id = num2str(now);
end

% protein shape model options
% ---------------------------
% protein.type              (optional) Holds the protein model type. The default is "vesicle".
param.protein.type = 'vesicle';

% protein.id                (optional) Holds the id of the protein model. The default is empty.
try
   [status,result] = system( 'uuidgen' ); 
   param.protein.id = result(1:end-1);
catch
   param.protein.id = num2str(now);
end

% protein.cytonuclearflag   (optional) Determines whether the protein pattern will be generated in
%                           the cytosolic space ('cyto'), nuclear space ('nuc') or everywhere ('all').
%                           Default is cyto.
if strcmpi( param.protein.class, 'nuc' )
   param.protein.cytonuclearflag = 'nuc';
else
   param.protein.cytonuclearflag = 'cyto';
end

% other options
% -------------
% verbose                   (optional) Displays messages to screen. The default is true.
param.verbose = true;

% debug                     (optional) Reports errors and warnings. Default is false.
param.debug = true;

% train.flag                (optional) Selects what model is going to be trained ('nuclear',
%                           'framework', or 'all'). Default is 'all'
param.train.flag = 'framework';

%documentation
param.documentation.description = 'This model has been trained using demo3D20 from CellOrganizer';

%model resolution
param.model.downsampling = [10,10,1];
param.model.proteinUpsampleZ = 5;

%this is the resolution of the datasets
param.model.resolution = [0.049, 0.049, 0.2000];

success = img2slml( dimensionality, dna, cell, protein, param );

load('lam.mat')

% plotReconstructionError(model)
img = showShapeSpace(model);

imwrite(img, ['demo3D20ShapeSpace.png'], 'png')

param = [];
param.targetDirectory = curr_path;
param.prefix = 'synthesizedImages';
param.numberOfSynthesizedImages = 1;
param.compression = 'lzw';
param.microscope = 'none';

param.verbose = true;
param.debug = true;
param.synthesis = 'framework';
param.synthesis.diffeomorphic.maximum_iterations = 40;

img = model2img({model});

imwrite(im2projection_RGB(img), ['demo3D20SampledCell.png'], 'png')

end


function plotReconstructionError(model)
    desired_shape_space_dimensionality = 3;
    distmat = model.nuclearShapeModel.distances;
    nimgs = model.nuclearShapeModel.numimgs;
    
    ninds = (nimgs^2-nimgs) / 2;
    
    niter = 25;
    
    template = triu(ones(size(distmat))) .* ~eye(size(distmat));
    inds = find(template);
    
    
    nrmv = round(ninds/2);
    
    sse = zeros(nrmv+1, niter);
    
    for i = nrmv:-1:0
        for j = 1:niter
            newinds = randsample(inds,i);
            
            dist_temp = triu(distmat);
            dist_temp(newinds) = NaN;
            dist_temp = dist_temp+dist_temp';
            
            
            
            distances = reconstruct_distance_matrix2(dist_temp, desired_shape_space_dimensionality);
            
            se = (distmat.*template - distances.*template).^2;
                
            sse(i+1,j) = sum(se(:));

            
        end
    end
    
    
    figure, 
    
    plot(fliplr((nrmv:-1:0) ./ ninds), mean(sse,2))
    hold on
    plot(fliplr((nrmv:-1:0) ./ ninds), mean(sse,2)+ std(sse,[],2), '--')
    plot(fliplr((nrmv:-1:0) ./ ninds), mean(sse,2)- std(sse,[],2), '--')
    
    xlabel('Percent of Distances Missing')
    ylabel('SSE')
    axis tight 
end