function [ output_args ] = test_reconstruction_methods( distances_complete, savedir )
%Distances_complete may be a sparse incomplete distance matrix
rng(1)


if ~exist('savedir', 'var') | isempty(savedir)             
    savedir = pwd;
end

tempdir = [savedir filesep mfilename];
if ~exist(tempdir, 'dir')
    mkdir(tempdir)
end

methods = {struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 15), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 14), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 13), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 12), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 11), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 10), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 9), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 8), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 7), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 6), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 5), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 4), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 3), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 2), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', 1), ...
            struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', inf)}%, ...;
%             struct('method', 'shortest path',   'weight_factor', 0, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'shortest path',   'weight_factor', 1, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'shortest path',   'weight_factor', 2, 'desired_shape_space_dimensionality', inf)};%, ...
%             
%             struct('method', 'regress',         'weight_factor', 1, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'regress',         'weight_factor', 2, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'regress',         'weight_factor', 0, 'desired_shape_space_dimensionality', inf), ...
%             struct('method',  'shortest path',   'weight_factor', 0, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'shortest path',   'weight_factor', 1, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'shortest path',   'weight_factor', 2, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'regress',         'weight_factor', 0, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'regress',         'weight_factor', 1, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'regress',         'weight_factor', 2, 'desired_shape_space_dimensionality', 7), ...
%             };

% methods = {struct('method',  'shortest path',   'weight_factor', 0, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'shortest path',   'weight_factor', 1, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'shortest path',   'weight_factor', 2, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'regress',         'weight_factor', 0, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'regress',         'weight_factor', 1, 'desired_shape_space_dimensionality', inf), ...
%             struct('method', 'regress',         'weight_factor', 2, 'desired_shape_space_dimensionality', inf), ...
%             struct('method',  'shortest path',   'weight_factor', 0, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'shortest path',   'weight_factor', 1, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'shortest path',   'weight_factor', 2, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'regress',         'weight_factor', 0, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'regress',         'weight_factor', 1, 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'regress',         'weight_factor', 2, 'desired_shape_space_dimensionality', 7), ...
%             };
%                 struct('method', 'nystrom-euclidean', 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'landmark-kernel-svd', 'desired_shape_space_dimensionality', 7), ...
%             struct('method', 'landmark-distance-svd', 'desired_shape_space_dimensionality', 7)    

    sampling = {'column'};

    niter = 10;

    if ~exist('distances_complete', 'var')
        %simulate some gaussian distrubuteddata

        ndims = 10;
        npts = 250;
        pts = mvnrnd(zeros(1, ndims), eye(10), npts);

        distances_complete = squareform(pdist(pts));
    else
        distances_complete(eye(size(distances_complete))>0) = nan;

        keepinds = find(~all(isnan(distances_complete),2));

        distances_complete = distances_complete(keepinds,keepinds);
        distances_complete(eye(size(distances_complete))>0) = 0;

         nancounts = sum(isnan(distances_complete),2);

        while any(nancounts)
            [~, ind] = max(nancounts);

            distances_complete(ind,:) = [];
            distances_complete(:, ind) = [];
            keepinds(ind) = [];

            nancounts = sum(isnan(distances_complete),2);
        end

        npts = size(distances_complete,1);
    end

    eigval_best_temp = eig(distances_complete);
    eigval_best_temp = eigval_best_temp / max(eigval_best_temp);
    eigval_complete = zeros(1, npts);
    eigval_complete(1:length(eigval_best_temp)) = eigval_best_temp(end:-1:1);
    
    method = struct('method',  'regress',   'weight_factor', 0, 'desired_shape_space_dimensionality', inf);
    
    posfile = [tempdir filesep 'posfile.mat'];
    if ~exist(posfile, 'file')
        [positions_complete, ~, ~] = embed_partial_distance_matrix(distances_complete, method);
        save(posfile, 'positions_complete')
    else
        load(posfile)
    end
    
    distances_euclidean = squareform(pdist(positions_complete));

for k = 1:niter
    rand_columns = randperm(npts);

    for i = 1:length(methods)

        tempfile = [tempdir filesep 'temp' num2str(k) '_' num2str(i) ];

        [can_start, final_name, final_exists, ~] = chunk_start(tempfile, '.mat');

        if can_start & ~final_exists
            disp(methods{i}.method);


            eignorm_vs_prev = ones(1,npts-1) * -1;
            eignorm_vs_best = ones(1,npts-1) * -1;
            matnorm_vs_prev = ones(1,npts-1) * -1;
            matnorm_vs_best = ones(1,npts-1) * -1;
            
            pos_vs_prev = ones(1,npts-1) * -1;
            pos_vs_best = ones(1,npts-1) * -1;

            corr_vs_best = ones(1,npts-1) * -1;
            
            this_data = false(npts, npts);
            this_data(logical(eye(npts))) = true;
            %add a column-worth of data each iteration

            this_data(rand_columns(1),:) = true;
            this_data(:,rand_columns(1)) = true;       

            positions = cell(1, npts-1);
            
            try

                for j = 1:(npts-1)
                    disp([num2str(j) filesep num2str(npts)]);

                    inds = find(~all(this_data,1));

                    if j == 1
                        endind = 1;
                    else
                        endind = 0;
                    end


                    switch sampling{1}
                        case 'random' %this is currently not implemented correctly
                            npts = sum(~this_data(inds(1:endind),:));

                            unsampled_points = find(~this_data);

                            inds = randperm(length(unsampled_points));

                            this_data(inds(npts)) = true;
                        case 'column'
                            rand_column = rand_columns(j+1);

                            this_data(rand_column,:) = true;
                            this_data(:,rand_column) = true;                
                        otherwise
                    end

                    fract_data(j) = sum(this_data(:))/ (npts^2);

                    distances_temp = distances_complete;
                    distances_temp(~this_data) = nan;

    %                 methods{i} = ml_initparam(methods{i}, struct

                    [positions{j}, ~, ~] = embed_partial_distance_matrix(distances_temp, methods{i});
                    dists = squareform(pdist(real(positions{j})));

                    eigval_tmp = eig(dists);
                    eigval_tmp = eigval_tmp / max(eigval_tmp);

                    eigval = zeros(1, npts);
                    eigval(1:length(eigval_tmp)) = eigval_tmp(end:-1:1);

                    if ~exist('dists_prev', 'var')
                        dists_prev = dists;
                        eigval_prev = eigval;
                        positions_prev = positions{j};
                    end

                    eignorm_vs_prev(j) = sum((eigval(:) - eigval_prev(:)).^2);
                    eignorm_vs_best(j) = sum((eigval(:) - eigval_complete(:)).^2);

                    matnorm_vs_prev(j) = sum((dists(:) - dists_prev(:)).^2);
                    matnorm_vs_best(j) = sum((dists(:) - distances_complete(:)).^2);

                    pos_vs_prev(j) = procrustes_l2(positions{j}, positions_prev);
                    pos_vs_best(j) = procrustes_l2(positions{j}, positions_complete);

                    corr = corrcoef(dists(:), distances_complete(:));
                    corr_vs_best(j) = corr(1,2).^2;
                    
                    corr_euc = corrcoef(dists(:), distances_euclidean(:));
                    corr_euc_vs_best = corr_euc(1,2).^2;
                    

                    dists_prev = dists;
                    eigval_prev = eigval;
                    positions_prev = positions{j};
                end


                save(tempfile, 'eignorm_vs_prev', 'eignorm_vs_best', 'matnorm_vs_prev', 'matnorm_vs_best', 'pos_vs_prev', 'pos_vs_best', 'positions', 'corr_vs_best', 'corr_euc_vs_best');
            catch
                disp(['Method ' methods{i}.method ' failed.']);
            end 
        end 
    end
end


eignorm_vs_prev_all = zeros((npts-1), length(methods), niter);
eignorm_vs_best_all = zeros((npts-1), length(methods), niter);

matnorm_vs_prev_all = zeros((npts-1), length(methods), niter);
matnorm_vs_best_all = zeros((npts-1), length(methods), niter);

pos_vs_prev_all = zeros((npts-1), length(methods), niter);
pos_vs_best_all = zeros((npts-1), length(methods), niter);

corecoeff_vs_best_all = zeros((npts-1), length(methods), niter);
corecoeff_euc_vs_best_all = zeros((npts-1), length(methods), niter);

for k = 1:niter   
    for i = 1:length(methods)

        tempfile = [tempdir filesep 'temp' num2str(k) '_' num2str(i) ];

        if exist([tempfile '.mat'], 'file')
            load(tempfile)
            
            eignorm_vs_prev_all(:,i,k) = eignorm_vs_prev;
            eignorm_vs_best_all(:,i,k) = eignorm_vs_best;

            matnorm_vs_prev_all(:,i,k) = matnorm_vs_prev;
            matnorm_vs_best_all(:,i,k) = matnorm_vs_best;

            pos_vs_prev_all(:,i,k) = pos_vs_prev;
            pos_vs_best_all(:,i,k) = pos_vs_best;
            
            corecoeff_vs_best_all(:,i,k) = 1- corr_vs_best;
            corecoeff_euc_vs_best_all(:,i,k) = 1 - corr_euc_vs_best;
        end
    end
end


keepinds = ~all(all(matnorm_vs_best_all == -1,3),1);
% keepinds([2,13,14]) = 0;

c = 1;
for i = 1:length(methods)
    if keepinds(i)
        
        functioning_methods{c} = ['D = ' num2str(methods{i}.desired_shape_space_dimensionality)];
        
%         if isfield(methods{i}, 'weight_factor')
% %             functioning_methods{c} = [functioning_methods{c} ' w=' num2str(methods{i}.weight_factor)];
%               
%         end
        c = c + 1;
    end
end


matnorm_vs_best = matnorm_vs_best_all(:,keepinds,:);
eignorm_vs_best = eignorm_vs_best_all(:,keepinds,:);
eignorm_vs_prev = eignorm_vs_prev_all(:,keepinds,:);
matnorm_vs_prev = matnorm_vs_prev_all(:,keepinds,:);
pos_vs_best = pos_vs_best_all(:,keepinds,:);
pos_vs_prev = pos_vs_prev_all(:,keepinds,:);

corecoeff_vs_best = corecoeff_vs_best_all(:,keepinds,:);
corecoeff_euc_vs_best = corecoeff_euc_vs_best_all(:,keepinds,:);


xdat = (2:npts)./npts;

colormap(jet)

figure('color', 'w')
h = plot_mean_and_stdev(xdat, matnorm_vs_best);
title('L2 distance matrix vs complete data')
xlabel('Fraction of completed columns')
ylabel('L2 vs complete data');
legend(h, functioning_methods)

if exist(savedir, 'dir')
    saveas(gcf, [savedir filesep mfilename '_dmat_vs_complete.tif'], 'tif')
end

figure('color', 'w')
h = plot_mean_and_stdev(xdat, corecoeff_vs_best);
title('correlation coefficient vs complete data')
xlabel('Fraction of completed columns')
ylabel('1 - R^2');
legend(h, functioning_methods, 'Location', 'SouthWest')

axis([0,1,0,1])

if exist(savedir, 'dir')
    saveas(gcf, [savedir filesep mfilename '_corrcoef_vs_complete.tif'], 'tif')
end





% 
% [~, ind] = sort(sum(sum(matnorm_vs_best,3),1));
% functioning_methods(ind)'
% 
% 
% figure('color', 'w')
% h = plot_mean_and_stdev(xdat, pos_vs_best);
% title('L2 positions vs complete data')
% xlabel('Fraction of completed columns')
% ylabel('L2 vs complete data');
% legend(h, functioning_methods)
% 
% if exist(savedir, 'dir')
%     saveas(gcf, [savedir filesep mfilename '_positions_vs_complete.tif'], 'tif')
% end
% 
% 
% figure('color', 'w')
% h = plot_mean_and_stdev(xdat, eignorm_vs_best);
% title('L2 eigenvalues vs complete data')
% xlabel('Fraction of completed columns')
% ylabel('L2 vs complete data');
% legend(h, functioning_methods)
% 
% if exist(savedir, 'dir')
%     saveas(gcf, [savedir filesep mfilename '_eig_vs_complete.tif'], 'tif')
% end
% 
% figure('color', 'w')
% h = plot_mean_and_stdev(xdat, matnorm_vs_prev);
% title('L2 distance matrix vs previous iteration')
% xlabel('Fraction of completed columns')
% ylabel('L2 vs complete data');
% legend(h, functioning_methods)
% 
% if exist(savedir, 'dir')
%     saveas(gcf, [savedir filesep mfilename '_dmat_vs_prev.tif'], 'tif')
% end
% 
% figure('color', 'w')
% h = plot_mean_and_stdev(xdat, eignorm_vs_prev);
% hold on
% title('L2 eigenvaluse vs previous iteration')
% xlabel('Fraction of completed columns')
% ylabel('L2 vs complete data');
% legend(h, functioning_methods)
% 
% if exist(savedir, 'dir')
%     saveas(gcf, [savedir filesep mfilename '_eig_vs_prev.tif'], 'tif')
% end
% 
% 
% figure('color', 'w')
% h = plot_mean_and_stdev(xdat, pos_vs_prev);
% title('L2 positions vs previous iteration')
% xlabel('Fraction of completed columns')
% ylabel('L2 vs complete data');
% legend(h, functioning_methods)
% 
% if exist(savedir, 'dir')
%     saveas(gcf, [savedir filesep mfilename '_pos_vs_prev.tif'], 'tif')
% end
% 

end

function h = plot_mean_and_stdev(xdat, ydat)
    m = mean(ydat,3);
    st = std(ydat,[], 3);
    
    colors = jet(size(m,2))*0.8;
    
    h = zeros(1, size(m,2));
    
    for i = 1:size(m,2)

        h(i) = plot(xdat, (m(:,i)), 'color', colors(i,:));
        hold on
        plot(xdat, (m(:,i)+st(:,i)), '--', 'color', colors(i,:))
        plot(xdat, (m(:,i)-st(:,i)), '--', 'color', colors(i,:))

        
    end
    
%     axis([0, 1, prctile(real(log10(m(m~=0))), 0.5), prctile(real(log10(m(m~=0))), 99.5)])
    
end

