function [W Wn] = selftuning(inX, inK)
% self-tuning method to construct the KNN affinity graph without self-loop
% input: 
%       inX: n by p (each row is a data point)
%       inK: the parameter k for KNN graph
% output:
%       W: affinity matrix (without self-loop)
%       Wn: symmetric normalized affinity matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculate the pairwise square distance between all the data points
dist_square = L2_distance_1(inX',inX');
dist_square(find(dist_square<0)) = 0;
clear inX;
dist = sqrt(dist_square);
n = size(dist, 1);
[dumb idx] = sort(dist, 2); % sort each row

W = zeros(n, n);
for i = 1:n
    W(i, idx(i,2:inK+1)) = dist(i, idx(i,2:inK+1)) + eps;
end;
W = max(W, W');
%W = (W+W')/2;
clear dist_square dist idx dumb;

% Find the count of nonzero for each column
col_count = sum(W~=0, 1)';
col_sum = sum(W, 1)';
col_mean = col_sum ./ col_count;
[x y val] = find(W);
W = sparse(x, y, -val.*val./col_mean(x)./col_mean(y)./2);
clear col_count col_sum col_mean x y val;

% Do exp function sequentially because of memory limitation
num = 2000;
num_iter = ceil(n/num);
S = sparse([]);
for i = 1:num_iter
  start_index = 1 + (i-1)*num;
  end_index = min(i*num, n);
  S1 = spfun(@exp, W(:,start_index:end_index)); % sparse exponential func
  S = [S S1];
  clear S1;
end
W = real(S);
clear S;

W = W + 0.000000000001*eye(n);
W = max(W,W');
D = diag(sum(W,2));
Dd = diag(D);
Dn=diag(sqrt(1./Dd));
Wn = Dn*W*Dn; 
Wn = max(Wn,Wn');