function y = tsne_vip(X, labels, no_dims, initial_dims, perplexity)
%TSNE Performs symmetric t-SNE on dataset X
%
%   mappedX = tsne(X, labels, no_dims, initial_dims, perplexity)
%   mappedX = tsne(X, labels, initial_solution, perplexity)
%
% The function performs symmetric t-SNE on the NxD dataset X to reduce its 
% dimensionality to no_dims dimensions (default = 2). The data is 
% preprocessed using PCA, reducing the dimensionality to initial_dims 
% dimensions (default = 30). Alternatively, an initial solution obtained 
% from an other dimensionality reduction technique may be specified in 
% initial_solution. The perplexity of the Gaussian kernel that is employed 
% can be specified through perplexity (default = 30). The labels of the
% data are not used by t-SNE itself, however, they are used to color
% intermediate plots. Please provide an empty labels matrix [] if you
% don't want to plot results during the optimization.
% The low-dimensional data representation is returned in mappedX.
%
%
% (C) Laurens van der Maaten, 2010
% University of California, San Diego


    if ~exist('labels', 'var')
        labels = [];
    end
    if ~exist('no_dims', 'var') || isempty(no_dims)
        no_dims = 2;
    end
     if ~exist('initial_dims', 'var') || isempty(initial_dims)
        initial_dims = min(50, size(X, 2));
    end
    if ~exist('perplexity', 'var') || isempty(perplexity)
        perplexity = 30;
    end
    
    % First check whether we already have an initial solution
    if numel(no_dims) > 1
        initial_solution = true;
        ydata = no_dims;
        no_dims = size(ydata, 2);
        perplexity = initial_dims;
    else
        initial_solution = false;
    end
   
normal_prep_X=1;    
if(normal_prep_X)
    % Normalize input data
    X = X - min(X(:));
    X = X / max(X(:));
    X = bsxfun(@minus, X, mean(X, 1));
    
    % Perform preprocessing using PCA
    if ~initial_solution
        disp('Preprocessing data using PCA...');
        if size(X, 2) < size(X, 1)
            C = X' * X;
        else
            C = (1 / size(X, 1)) * (X * X');
        end
        [M, lambda] = eig(C);
        [lambda, ind] = sort(diag(lambda), 'descend');
        M = M(:,ind(1:initial_dims));
        lambda = lambda(1:initial_dims);
        if ~(size(X, 2) < size(X, 1))
            M = bsxfun(@times, X' * M, (1 ./ sqrt(size(X, 1) .* lambda))');
        end
        X = bsxfun(@minus, X, mean(X, 1)) * M;
        clear M lambda ind
    end
    
    % Compute pairwise distance matrix
    sum_X = sum(X .^ 2, 2);
    D = bsxfun(@plus, sum_X, bsxfun(@plus, sum_X', -2 * (X * X')));
    sizeD = size(D)
else	
    % Compute pairwise distance matrix missing values
	
    [NX,DX] = size(X);
	D = zeros(NX);
	for i=1:NX
		i
	for j=1:NX
		v=0; c=0;
		for k=1:DX
			if (~ismissing(X(i,k)*X(j,k)))
				v = v + (X(i,k) - X(j,k))^2;
				c = c+1;
			end;
		end;
		D(i,j) = v / c;
	end;
	end;
end;    
    
   % ntrain=3681
   % ind1 = find(labels(1:ntrain)==1);
   % D(ind1, ind1) = D(ind1, ind1) / 10;
   %
   % ind0 = find(labels(1:ntrain)==0);
   % D(ind0, ind0) = D(ind0, ind0) / 10;
   %
   % D(ind0, ind1) = D(ind0, ind1) * 10;
   % D(ind1, ind0) = D(ind1, ind0) * 10;
    
        
    % Compute joint probabilities
    P = d2p(D, perplexity, 1e-5);                                           % compute affinities using fixed perplexity
    clear D
	
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
    
    % Run t-SNE
 


	    % if ~exist('labels', 'var')
	    %     labels = [];
	    % end
	    % if ~exist('no_dims', 'var') || isempty(no_dims)
	    %     no_dims = 2;
	    % end
    
	    % First check whether we already have an initial solution
	    % if numel(no_dims) > 1
	    %     initial_solution = true;
	    %     ydata = no_dims;
	    %     no_dims = size(ydata, 2);
	    % else
	    %     initial_solution = false;
	    % end
    
	    % Initialize some variables
	    n = size(P, 1);                                     % number of instances

	    stop_lying_iter = 100;                              % iteration at which lying about P-values is stopped
	    max_iter = 1000;                                    % maximum number of iterations
	    epsilon = 500;                                      % initial learning rate
	    min_gain = .01;                                     % minimum gain for delta-bar-delta
    
	    % Make sure P-vals are set properly
	    P(1:n + 1:end) = 0;                                 % set diagonal to zero
	    P = 0.5 * (P + P');                                 % symmetrize P-values
	    P = max(P ./ sum(P(:)), realmin);                   % make sure P-values sum to one
	    const = sum(P(:) .* log(P(:)));                     % constant in KL divergence
	    if ~initial_solution
	        P = P * 6;                                      % lie about the P-vals to find better local minima
	    end
    
	    % Initialize the solution
	        y = .0001 * randn(n, no_dims);
			y_incs  = zeros(size(y));
			gains = ones(size(y));
    
	    % Run the iterations
	    for iter=1:max_iter
        
	        % Compute joint probability that point i and j are neighbors
	        sum_y2 = sum(y .^ 2, 2);
	        Qnum = 1 ./ (1 + bsxfun(@plus, sum_y2, bsxfun(@plus, sum_y2', -2 * (y * y')))); 
			% Student-t distribution as 1 / (1 + ||y_i - y_j||^2 )
		
	        Qnum(1:n+1:end) = 0;                                                 % set diagonal to zero
	        Q = max(Qnum ./ sum(Qnum(:)), realmin);                               % normalize to get probabilities
        
	        % Compute the gradients (faster implementation)
	        L = (P - Q) .* Qnum;
	        y_grads = 4 * (diag(sum(L, 1)) - L) * y;
           
			%%%% Update the solution
			use_original_momentum_method=1;
			if(use_original_momentum_method)
				%%%% original momentum method
				%%%%% note that the y_grads are actually -y_grads
				gains = (gains + .2) .* (sign(y_grads) ~= sign(y_incs)) + (gains * .8) .* (sign(y_grads) == sign(y_incs));
				gains(gains < min_gain) = min_gain;
			 	y_incs =  - epsilon * (gains .* y_grads);
			else				   
			    %%%% simple gradient descent
				y_incs = -epsilon * y_grads;
			end;	
			
	        y = y + y_incs;
	        y = bsxfun(@minus, y, mean(y, 1));
        
	        % Update the momentum if necessary
	        if iter == stop_lying_iter && ~initial_solution
  	            P = P ./ 6;
  	        end
        
	        % Print out progress
	        if ~rem(iter, 10)
	            cost = const - sum(P(:) .* log(Q(:)));
	            disp(['Iteration ' num2str(iter) ': error is ' num2str(cost)]);
	        end
        
	        % Display scatter plot (maximally first three dimensions)
	        if ~rem(iter, 10) && ~isempty(labels)
				figure(3); clf;
	            if no_dims == 1
	                scatter(y, y, 19, labels, 'filled');
	            elseif no_dims == 2
	                scatter(y(:,1), y(:,2), 19, labels, 'filled');
	            else
	                scatter3(y(:,1), y(:,2), y(:,3), 40, labels, 'filled');
	            end
				grid on;
	            axis tight
	            axis on;
	            drawnow
	        end
	    end   
