function [nn, kk] = idxSIDmap(sid, i, j, d, verbose)
% idxSIDmap       Index map between SID dictionaries of different size
%                 i.e. Di = Dj(nn,kk), where i+d <= j, and  0 <= d
% Di and Dj are dictionaries of different size depending on i and j,
% each consists of a "basic"-block, D0, and i or j repetitions of a small
% block, Ds. The block Di may be shifted d small blocks in Dj.
% The values in the SID are not relevant here, only the structure.
%
% use:   [nn, kk] = idxSIDmap(sid, i, j, d, v);
% -------------------------------------------------------------------------
%   nn    Row indexes for Di in Dj
%   kk    Column indexes for Di in Dj
%   sid   The Shift Invariant Dictionary, as in initSID
%   i     Gives size of Di, i = 1,2,...
%   j     Gives size of Dj, j = i,i+1,i+2,...   (j >= i)
%   d     How many small blocks Di is "shifted" in Dj, d = 0,1,..
%   v     verbose
% -------------------------------------------------------------------------
% ex:
%    Qm = [35,25,20,10,10];  Pm = [ 5, 2, 2, 1, 1]; 
%    sid = initSID(Qm, Pm, randn(sum(Qm),1));
%    [nn, kk] = idxSIDmap(sid, 17, 1997, 0, 1); 
%    D2 = full(makeSIDmatrix(sid,50));
%    D5 = full(makeSIDmatrix(sid,80)); 
%    [nn, kk] = idxSIDmap(sid, 2, 5, 1); disp(norm(D2-D5(nn,kk)));

% -------------------------------------------------------------------------
% Made by Karl Skretting, University of Stavanger (UiS), 
% Department of Electrical Engineering and Computer Science (IDE)
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  20.11.2015  KS: made file
% -------------------------------------------------------------------------

nn = [];
kk = [];
if ~isstruct(sid) || ~isfield(sid,'Type') || ~strcmp(sid.Type, 'SID')
    disp([mfilename,': First argument not a SID structure.']);
    return
end
if (nargin < 3)
    disp([mfilename,': There should be at least three input arguments.']);
    return
end
if (i > j)
    disp([mfilename,': should have j >= i.']);
    return
end
if (nargin < 4)
    d = 0;
end
if nargin < 5
    verbose = false;
end
if (d < 0)
    disp([mfilename,': negative d is set to 0']);
    d = 0; 
end
if ((i+d) > j) 
    disp([mfilename,': too large d, ',int2str(d),', is set to max value ',int2str(j-i),'.']);
    d = j-i; 
end

M = sid.M;    %  number of atoms in dictionary
Q = sid.Q;    %  total number of free variables, i.e. sum of Qm
Qm = sid.Qm;  %  length of each atom, size Mx1
Pm = sid.Pm;  %  shift for each atom, size Mx1
Qmax = max(Qm);
Ns = thisLCM(Pm);

Ksm = Ns./Pm;
Ks = sum(Ksm);
N0 = ceil(max(Qm)/Ns)*Ns - Ns;    % = N1-Ns, N0 is a multiple of Ns
K0 = sum(1+floor((N0+Ns-Qm)./Pm)) - Ks;

% ******  OLD and probably wrong code
% Ksm = Ns./Pm;
% Ks = sum(Ksm);
% N1m = (Ksm-1).*Pm + Qm;
% N0 = max(N1m) - Ns;
% K1m = floor( (N0 + Ns + Pm - Qm)./Pm );
% K0m = K1m - Ksm;
% K0 = sum(K1m) - Ks;
%    fprintf('    Ksm : %7i',Ksm(1)); for m=2:M; fprintf(', %7i',Ksm(m)); end; fprintf('\n');
%    fprintf('    K0m : %7i',K0m(1)); for m=2:M; fprintf(', %7i',K0m(m)); end; fprintf(',  sum = %7i.\n',sum(K0m));
% ****** 

Ni = N0 + i*Ns;
Ki = K0 + i*Ks;
Nj = N0 + j*Ns;
Kj = K0 + j*Ks;
Kim = floor( (Ni+Pm-Qm)./Pm );
Kjm = floor( (Nj+Pm-Qm)./Pm );
k2i = cumsum(Kim);
k1i = [1; k2i(1:(M-1))+1];
k2j = cumsum(Kjm);
k1j = [1; k2j(1:(M-1))+1];
nn = (1+d*Ns):(Ni+d*Ns);
kk = zeros(1,Ki);
for m = 1:M
    a = (k1j(m) + d*Ksm(m));
    kk(k1i(m):k2i(m)) = a:(a+Kim(m)-1);
end

if verbose
    fprintf('%s: SID has M = %i atoms, Q = %i, and Qmax = %i.\n', mfilename, M, Q, Qmax);
    fprintf('    Input parameters:  i = %i,  j = %i,  d = %i.\n', i, j, d);
    fprintf('    Small block Ds is %i x %i,  "basic block" D0 is %i x %i.\n', Ns, Ks, N0, K0);
    fprintf('    Di is %i x %i,  Dj is %i x %i.\n', Ni, Ki, Nj, Kj);
    fprintf('    Qm  : %7i',Qm(1)); for m=2:M; fprintf(', %7i',Qm(m)); end; fprintf('\n');
    fprintf('    Pm  : %7i',Pm(1)); for m=2:M; fprintf(', %7i',Pm(m)); end; fprintf('\n');
    fprintf('    Kim : %7i',Kim(1)); for m=2:M; fprintf(', %7i',Kim(m)); end; fprintf(',  sum = %7i.\n',sum(Kim));
    fprintf('    Kjm : %7i',Kjm(1)); for m=2:M; fprintf(', %7i',Kjm(m)); end; fprintf(',  sum = %7i.\n',sum(Kjm));
    fprintf('    k1i : %7i',k1i(1)); for m=2:M; fprintf(', %7i',k1i(m)); end; fprintf('\n');
    fprintf('    k2i : %7i',k2i(1)); for m=2:M; fprintf(', %7i',k2i(m)); end; fprintf('\n');
    fprintf('    k1j : %7i',k1j(1)); for m=2:M; fprintf(', %7i',k1j(m)); end; fprintf('\n');
    fprintf('    k2j : %7i',k2j(1)); for m=2:M; fprintf(', %7i',k2j(m)); end; fprintf('\n');
    % fprintf('    N1m : %7i',N1m(1)); for m=2:M; fprintf(', %7i',N1m(m)); end; fprintf('\n');
end

end

%% subfunctions
function y = thisLCM(x)
% x may be a vector of several elements (positive integers)
% matlab lcm function allows only two input arguments, i.e it works
% elementwise on vectors and matrix inputs!
M = numel(x);
if (M == 1) || (min(x) == max(x))
    y = x(1);
else
    y = lcm(x(1), x(2));
    for m = 3:M
        y = lcm(y, x(m));
    end
end
end

