global DISTM;
global X;
global L;
global N;
global D;

%% Mnist dataset
% W = load ('../../../data/mnist-original.mat');
% X = W.data';
% L = W.label';
% X = single(X);
% [N,D] = size(X)


% 20NG dataset
W = load ('../../../data/20news-bydate-matlab/train.data');
N = max(W(:,1))
XS = spconvert(W);
X = full(XS);
L = load ('../../../data/20news-bydate-matlab/train.label');
 [N,D] = size(X)

% W = load ('../../../data/20news-bydate-matlab/test.data');
% TS = spconvert(W);
% TS(:,53976:end) = [];
% T = full(TS);function ydata = 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
			%
			%     % 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')));
			%
			%     % 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 = 200;                                    % 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 * 4;                                      % lie about the P-vals to find better local minima
			% 	    end
			%
			% 	    % Initialize the solution
			% 	        ydata = .0001 * randn(n, no_dims);
			% y_incs  = zeros(size(ydata));
			% gains = ones(size(ydata));
			%
			% 	    % Run the iterations
			% 	    for iter=1:max_iter
			%
			% 	        % Compute joint probability that point i and j are neighbors
			% 	        sum_ydata = sum(ydata .^ 2, 2);
			% 	        num = 1 ./ (1 + bsxfun(@plus, sum_ydata, bsxfun(@plus, sum_ydata', -2 * (ydata * ydata')))); % Student-t distribution as 1 / (1 + ||x_i - x_j||^2 )
			%
			% 	        num(1:n+1:end) = 0;                                                 % set diagonal to zero
			% 	        Q = max(num ./ sum(num(:)), realmin);                               % normalize to get probabilities
			%
			% 	        % Compute the gradients (faster implementation)
			% 	        L = (P - Q) .* num;
			% 	        y_grads = 4 * (diag(sum(L, 1)) - L) * ydata;
			%
			% 	        % Update the solution
			% 	        gains = (gains + .2) .* (sign(y_grads) ~= sign(y_incs)) ...         % note that the y_grads are actually -y_grads
			% 	              + (gains * .8) .* (sign(y_grads) == sign(y_incs));
			% 	        gains(gains < min_gain) = min_gain;
			% 	        y_incs =  - epsilon * (gains .* y_grads);
			% 	        ydata = ydata + y_incs;
			% 	        ydata = bsxfun(@minus, ydata, mean(ydata, 1));
			%
			% 	        % Update the momentum if necessary
			% 	        if iter == stop_lying_iter && ~initial_solution
			% 	            P = P ./ 4;
			% 	        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(2); clf;
			% 	            if no_dims == 1
			% 	                scatter(ydata, ydata, 19, labels, 'filled');
			% 	            elseif no_dims == 2
			% 	                scatter(ydata(:,1), ydata(:,2), 19, labels, 'filled');
			% 	            else
			% 	                scatter3(ydata(:,1), ydata(:,2), ydata(:,3), 40, labels, 'filled');
			% 	            end
			% 	grid on
			% 	            axis tight
			% 	            axis off
			% 	            drawnow
			% 	        end
			% 	    end 
% LT = load ('../../../data/20news-bydate-matlab/test.label');
% %%%%% weight features by IDF
% DF = sum (X~=0,1);
% IDF = log2(11269 ./ DF);
% IDF(find(IDF<3)) =0;
% X = bsxfun(@times, X, IDF) ;
%T = bsxfun(@times, T, IDF)

%% sample dataset
Nnew= 10000;
r = randperm(N);
ind = r(1:Nnew)';
X = X(ind, :);
L = L(ind);
N=Nnew

whos; 



numGroups = length(unique(L));
clr = hsv(numGroups)



 no_dims=2;
 PCA_dims=800;
 perplexity=20; % # of effective neighbors

% Xtnse = tsne(X,'Algorithm','exact','NumPCAComponents',PCA_dims, 'NumDimensions',2,'Perplexity',30, 'Verbose',1 );
  
 Xtnse = tsne_vip (X, L, no_dims, PCA_dims, perplexity);


 
 figure(1);clf;
 gscatter(Xtnse(:,1),Xtnse(:,2),L,clr);
 
 %%%%%%% Euclidian distance fast
  X2 = dot(Xtnse,Xtnse,2);
 %%%% euclidian
 DISTM = bsxfun(@plus, X2, X2') - 2 * Xtnse * Xtnse';

 avd_dist = mean(mean(DISTM))
 
 esp = avd_dist/500
 minpts = N/500
%idx = dbscan(Xtnse, epsilon ,minpts)
idx = dbscan(DISTM, esp, minpts,'Distance','precomputed');
NUM_COLORS_MATL = length(unique(idx))
[hc_MATL, edges_MATL] = histcounts(idx,NUM_COLORS_MATL+1);
hc_MATL

[counts, groupnames] = groupcounts(idx);[groupnames counts]

CONFMAT = confusionmat (idx, L);CONFMAT (:,1:10)

pause;

%
%
% labels=L
 
 
 
 
 % X = Xtnse_maat;
 % X = single(X);
 % [N,D]= size(X)
 %
 % X2 = dot(X,X,2);
 % DISTM = bsxfun(@plus, X2, X2') - 2 * X * X';
 %
 % avd_dist = mean(mean(DISTM))
 % pause;
 %
 %  dbcolor = dbscan_alg(6, 40, 4);
 
 
 
 % figure(2);clf;
 % scatter3(Xtnse_maat(:,1),Xtnse_maat(:,2),Xtnse_maat(:,3),L);


% % Y3 = tsne(X,'Algorithm','barneshut','NumPCAComponents',50,'numdimensions',3, 'perplexity', 20);
% Y3 = tsne(X,'Algorithm','barneshut','NumPCAComponents',300,'numdimensions',3,'distance','cosine','perplexity',50);
% figure(4);clf;
% %scatter3 ( Y3(:,1),Y3(:,2),Y3(:,3), 10 ,L)
% scatter3(Y3(:,1),Y3(:,2),Y3(:,3),15,clr(L,:),'filled');
