vb_sparse_inv_distance

 Compute sparse inverse distance matrix with (distance) < Rmax

 Dinv = vb_sparse_distance(X, Y, Rmax, Nblk);
 Dinv = vb_sparse_distance(X, Y, Rmax);

 X : (M x D) Data matrix
 Y : (N x D) Data matrix
 X(m,:) : D-dim vector for m-th point 
 Y(n,:) : D-dim vector for n-th point 

 M, N : total number of points for X, Y
 D : data dimension

 Rmax : Max distance
 Nblk : block length for block-wise distance matrix calculation

 Sparse inverse distance matrix with (distance) < Rmax
   Dinv = sparse(I, J, 1./dd, M, N);
   
 dd(n) = ||X(I(n),:) - Y(J(n),:)||  < Rmax

   2017-3-6 Masa-aki Sato

0001 function    Dinv = vb_sparse_inv_distance(X, Y, Rmax, Nblk);
0002 % Compute sparse inverse distance matrix with (distance) < Rmax
0003 %
0004 % Dinv = vb_sparse_distance(X, Y, Rmax, Nblk);
0005 % Dinv = vb_sparse_distance(X, Y, Rmax);
0006 %
0007 % X : (M x D) Data matrix
0008 % Y : (N x D) Data matrix
0009 % X(m,:) : D-dim vector for m-th point
0010 % Y(n,:) : D-dim vector for n-th point
0011 %
0012 % M, N : total number of points for X, Y
0013 % D : data dimension
0014 %
0015 % Rmax : Max distance
0016 % Nblk : block length for block-wise distance matrix calculation
0017 %
0018 % Sparse inverse distance matrix with (distance) < Rmax
0019 %   Dinv = sparse(I, J, 1./dd, M, N);
0020 %
0021 % dd(n) = ||X(I(n),:) - Y(J(n),:)||  < Rmax
0022 %
0023 %   2017-3-6 Masa-aki Sato
0024 
0025 if nargin < 4 || isempty(Nblk), Nblk = 100; end;
0026 if nargin < 3,  help vb_sparse_inv_distance; return; end
0027 if (size(X,2) ~= size(Y,2)), error('size(X,2) ~= size(Y,2)'); end
0028 
0029 % Rmax : Max distance
0030 Rmax = Rmax^2;
0031 
0032 % total number of points
0033 M = size(X,1);
0034 N = size(Y,1);
0035 
0036 % Index for Sparse distance matrix
0037 II = zeros(M,1);
0038 JJ = zeros(M,1);
0039 DD = zeros(M,1);
0040 
0041 ns = 0;
0042 
0043 if M < N
0044     % Loop for X
0045     for i1=1:Nblk:M
0046         i2 = i1+Nblk-1;
0047         if (i2> M), i2=M; end;
0048         
0049         % block-wise distance matrix
0050         XX = X(i1:i2,:);  
0051         dd = vb_sq_distance(XX,Y,1);
0052         
0053         % find neighbor points with distance < Rmax
0054         sw = (dd < Rmax);
0055         [I,J] = find( sw );
0056         D = dd(sw);
0057         
0058         n = length(I);
0059         ix = (1:n) + ns;
0060         
0061         II(ix) = I(:) + i1 - 1;
0062         JJ(ix) = J(:);
0063         DD(ix) = D(:);
0064         
0065         ns = ns + n;
0066         
0067 %       II = I + i1 - 1;
0068 %       ds = ds + sparse(II(:), J(:), DD(:), M, N);
0069     end;
0070 else
0071     % Loop for Y
0072     for i1=1:Nblk:N
0073         i2 = i1+Nblk-1;
0074         if (i2> N), i2=N; end;
0075         
0076         % block-wise distance matrix
0077         YY = Y(i1:i2,:);  
0078         dd = vb_sq_distance(X,YY,1);
0079         
0080         % find neighbor points with distance < Rmax
0081         sw = (dd < Rmax);
0082         [I,J] = find( sw );
0083         D = dd(sw);
0084         
0085         n = length(I);
0086         ix = (1:n) + ns;
0087         
0088         II(ix) = I(:);
0089         JJ(ix) = J(:) + i1 - 1;
0090         DD(ix) = D(:);
0091         
0092         ns = ns + n;
0093         
0094 %       JJ = J + i1 - 1;
0095 %       ds = ds + sparse(I(:), JJ(:), DD(:), M, N);
0096     end;
0097 end;
0098 
0099 
0100 % matrix index with (distance) < Rmax
0101 II = II(1:ns);
0102 JJ = JJ(1:ns);
0103 
0104 % distance with (distance) < Rmax
0105 DD = DD(1:ns);
0106 
0107 %fprintf('rate of nonzero/(MxN) = %g in sparce inverse distance\n', ...
0108 %   100*ns/(M*N))
0109 
0110 % set min distance for inverse distance
0111 MinVal = Rmax * 1e-10;
0112 DD = max(DD, MinVal);
0113 
0114 % inverse distance
0115 DD = sqrt(1./DD);
0116 
0117 % Sparse inverse distance matrix with (distance) < Rmax
0118 Dinv = sparse(II, JJ, DD, M, N);
0119