function [F,I,N] = dpg_cal( varargin );

% [F,I,N ] = dpg_cal( Wk, Sk, Y, [opt args]);
% 
% opt args:
%  'v', {0,1}              - verbosity off/on
%  'kernel', string        - grappa kernel size, x by y. e.g. '2x5', '4x1'
%  'N', (xy)-by-num_patterns-L matrix   - grappa reconstruction parameters
%  'acs', array of indices - index locations to use for ACS lines
% 
% This implementation is derived from 
% "Generalized autocalibrating partially parallel acquisitions"
%   by M. A. Griswold, P. M. Jakob, et. al..
%  Mag. Reson. Med, 47(6):1202-1210, Jun 2002
%

%
% Copyright 2004-2016  W. Scott Hoge  (wsh032580 at proton dot me)
%
% Licensed under the terms of the MIT License 
% (https://opensource.org/licenses/MIT)

%
% GRAPPA: data from all coils used to patch data in one coil
%
%  coil 1     o * o * a * o *
%  coil 2     o * o * a * o *
%  coil 3     o * o * a * o *
%  coil 4     o * o * x * o *
%
%  *: acquired line, o: unacquired line, a: auto-callibration line
%  x: reconstructed line
%
%  S_j (k_y - m\Delta k_y ) = \sum_{l=1}^L \sum_{b=0}^{N_b-1}
%                                  n(j,b,l,m) S_l(k_y - b A \Delta k_y)
%

Wk  = varargin{1};
kin = varargin{2};
Y   = varargin{3};

acs = []; N = []; verbose=0; kernel = '2x5'; dks = []; prnt = 0;
eta = 1.0;
chi = 0.0001;

vararg = varargin(4:end);
for cnt=1:2:length(vararg),
  if strcmp( vararg{cnt}, 'v' ),
    verbose = cell2mat(vararg(cnt+1));
  elseif strcmp( vararg{cnt}, 'N' ),
    N = vararg{cnt+1};
  elseif strcmp( vararg{cnt}, 'dks' ),
    dks = vararg{cnt+1};
  elseif strcmp( vararg{cnt}, 'kernel' ),
    kernel = vararg{cnt+1};
  elseif strcmp( vararg{cnt}, 'acs' ),
    acs = vararg{cnt+1};
  elseif (strcmp( vararg{cnt}, 'report' ) || strcmp( vararg{cnt}, 'r' )),
    prnt = vararg{cnt+1};
  end;
end;

if length(size(Wk)) == 3, 
  [Wm,Wn,ncoils] = size(Wk);
else,
  Wm = Wk(1); Wn = Wk(2); ncoils = Wk(3);
end;


%% zero-pad the k-space data
Sk = zeros( Wm, Wn, ncoils );
if iscell(kin)
  dualks = 1;
  s = kin{3};
else
  dualks = 0;
  s = kin; clear kin;
end;
Sk(Y,:,:) = s;

rptupd(prnt,['* ' mfilename ': kernel size = ' kernel '\n']);

d = diff(Y(:));
if ( isempty(acs) ),
  prfl = max( sum(abs( Sk(:,:,:)),3).' );

  if length( d == 1) == 1,
    acs = find( prfl == max(prfl(Y(find(d==1)))) );
  else
    acs = Y([ min(find(d==1)); find( d == 1)+1 ]);
  end;
  rptupd(prnt,['* ' mfilename ': acs =']);
  rptupd(prnt,sprintf(' %d', acs) );
  rptupd(prnt,'\n');
end;

if ( isempty(dks) ),
  %% determine which delta-k's to repair.
  dks = [];
  for cnt=max(d):-1:2,
    if (sum( cnt == d ) > 0),
      dks = [ cnt; dks ];
    end;
  end;
end;

% set the kernel spacings based on 'kernel' string
nky = str2double(kernel(1));           
nkx = str2double(kernel(3:end)); % should be odd

indx3 = ceil((0-nkx)/2) -1 + (1:nkx); % [-2:2];

pattern = {};

if ~isempty(N),
  for cnt=1:length(N)
    pattern{cnt} = N(cnt).pattern;
  end;
end;

if isempty(pattern)
%% foreach dk spacing
for cnt=1:length(dks),
  
  if ( dks(cnt) == 1 )
    % for no acceleration, just specify the pattern
    tmppat = repmat('*',[1 nky]);
    for patcnt = 1:nky,
      pattern{length(pattern)+1} = tmppat;
      pattern{length(pattern)}(patcnt) = '<';
    end;

    for patcnt = 1:nky,
      pattern{length(pattern)+1} = tmppat;
      pattern{length(pattern)}(patcnt) = '>';
    end;

    continue;
  end;

  tmp = [ repmat([ '*' repmat('o',1,dks(cnt)-1) ],1,nky-1) '*' ];
  z2 = strfind(tmp,'o');
  if isempty(z2),                       
    % pad with zeros, if the kernel is 1-by-N
    tmp = [ 'o' tmp 'o' ];   
    z2 = strfind(tmp,'o');
  end;
  [tmp2,indx] = sort( abs(z2 - z2( floor(median(1:length(z2))) )) );

  olength = length(pattern);
  if (dualks == 1)
    % for dual-kernels, we want to genereate coefficients for the sampled
    % lines as well.
    npat = 2*length(tmp);
  else
    npat = length(z2);
  end
  for cnt2 = 1:npat, % (dks(cnt)-1), % 

    if (dualks == 0 )
      cnt3 = indx(cnt2);
      pattern{length(pattern)+1} = tmp;
      pattern{length(pattern)}(z2(cnt3)) = 'x';
    else
      cnt3 = mod( cnt2-1, length(tmp) ) + 1;
      pattern{length(pattern)+1} = tmp;
      if ( cnt2 <= npat/2 )
        pattern{length(pattern)}( cnt3 ) = '<';
      else
        pattern{length(pattern)}( cnt3 ) = '>';
      end;
    end;
  
  end;

end;
end;


% find the recon patterns we can use 
plist = cell(length(pattern),1);

for cnt = 1:length(pattern),
  v = extrapolate_pattern(pattern{cnt},nky);     

  % if cnt == 14, keyboard; end;

  for cnt2 = 1:length(acs);
    %% check if the GRAPPA pattern...
    indx = acs(cnt2) + v;

    %% can be accomodated by the sampling pattern 
    tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
    if ( length(tmp) == length(indx) ),
      plist{cnt} =  [ plist{cnt} acs(cnt2) ];
    end;
  end;
end;

rptupd(prnt,['* ' mfilename ': sampling patterns in use = \n' ]);
for cnt1=1:length(pattern);
  rptupd(prnt,sprintf('  %s: %d acs lines =', pattern{cnt1}, length(plist{cnt1})) );
  rptupd(prnt,sprintf(' %3d', plist{cnt1}) );
  % for cnt2 = 1:length( plist{cnt1} ),
  %   rptupd(prnt,('  %s', pattern{plist{cnt1}} );
  %   % v = extrapolate_pattern(pattern{plist{cnt1}(cnt2)},nky);
  %   % norm(vec(Sk( acs(cnt1)+v, :, : )));
  % end;
  rptupd(prnt,'\n');
end;
if (length(plist) == 0), 
  error(['sorry, acquisition pattern is incompatible with available GRAPPA recon patterns']);
end;

st = clock; % tic;

if verbose, rptupd(prnt,'\n'); end;

mss = zeros( size(Sk,1), 1 );         % set the missing lines index vector 
mss(Y) = 1;

if isempty(N),
  % N = zeros( length(indx3)*nky*size(s,3),length(pattern),size(s,3));

  rptupd(prnt,'calc recon param\n');
  for cnt = 1:length(pattern), 
    v = extrapolate_pattern(pattern{cnt},nky);

    %% and for each pattern ... 
    %% A2 = [];     b2 = [];
    
    %%  A is data from the expanded input data matrix.
    %%     size(A) = [ kxnum * nacs, ncoils * nky * nkx ]
    %%
    %% in the matlab (and c-code) structure is:
    %%    for each element of indx3
    %%                   coil 1   |    coil 2   |    coil 3   | ...
    %%                -3 -1  1  3 | -3 -1  1  3 | -3 -1  1  3 | ...
    %%                 o  o  o  o |  o  o  o  o |  o  o  o  o | ...
    %%             ^
    %%  ACS(1) xres|
    %%             v
    %%             ^
    %%  ACS(2) xres|
    %%             v
    %%             ^
    %%  ACS(3) xres|
    %%             v
    %% 
    %%  the older version of matlab code sorted the columns by
    %% 
    %%     | c1 c2 c3 ... | c1 c2 c3 ... | c1 c2 c3 ...  | c1 c2 c3 ...
    %%     | -3 -3 -3 ... | -1 -1 -1 ... |  1  1  1 ...  |  3  3  3 ...
    %%     |  o  o  o ... |  o  o  o ... |  o  o  o ...  |  o  o  o ...
    %% 
    %%  the order is immaterial, as long at it is consistent...
    %%  UPDATE (12 Jan 2011): the order does matter for GROG, so we have
    %%                        now changed it here so that the c-code and matlab code match.
    %%   
    %%  The old way to switch between the two was:
    %%    Accode(:, permute( reshape(1:size(Accode,2),YKS,Ncoils,XKS), [3 2 1]) ) == Amatlab
    %%    Accode = Amatlab(:, permute( reshape(1:size(Accode,2),XKS,Ncoils,YKS), [3 2 1]))
    A = zeros( length(plist{cnt})*size(s,2), length(indx3)*nky*size(s,3) );
    
    rptupd(prnt,sprintf('%3d  %s: ', cnt, pattern{cnt}));

    % v = extrapolate_pattern(pattern{cnt});
    N(length(N)+1).pattern = pattern{cnt};

    n = zeros( nky*size(s,3)*length(indx3),size(s,3));

    for cnt2 = 1:length(plist{cnt}),
      %% and ACS line ...

      %% if the sampling pattern accomodates the GRAPPA pattern...
      indx  = plist{cnt}(cnt2) + v;         % feedforward indices
      % indx2 = plist{cnt}(cnt2) + f;         % feedback indices

      %% determine the k-space filling coefficients 
      for cnt3=1:length(indx3),
        y0 = mod( indx3(cnt3) + [1:size(Sk,2)] , size(Sk,2) );
        y0( y0 == 0 ) = size(Sk,2);
        
        % generate a hybrid source matrix
        if (dualks)
          tmpSksub = zeros( size( Sk(indx,y0,:) ) );
          if ( sum(pattern{cnt}=='<') ), 
            tmpSksub( 1:2:length(indx), :, : ) = kin{1}( indx(1:2:end),y0,: );
            tmpSksub( 2:2:length(indx), :, : ) = kin{2}( indx(2:2:end),y0,: );
%            tmpSksub( (v < 0), :, : )  = kin{2}( indx( (v < 0) ),y0,: );
%            tmpSksub( (v >= 0), :, : ) = kin{3}( indx( (v >= 0) ),y0,: );
          else
            tmpSksub( 1:2:length(indx), :, : ) = kin{2}( indx(1:2:end),y0,: );
            tmpSksub( 2:2:length(indx), :, : ) = kin{1}( indx(2:2:end),y0,: );
%            tmpSksub( (v >= 0), :, : ) = kin{2}( indx( (v >= 0) ),y0,: );
%            tmpSksub( (v < 0), :, : )  = kin{3}( indx( (v < 0) ),y0,: );
          end
        else
          tmpSksub = Sk(indx,y0,:);
        end;
        % A( size(s,2)*(cnt2-1) + (1:size(s,2)), cnt3:length(indx3):size(A,2) ) = ...
        A( size(s,2)*(cnt2-1) + (1:size(s,2)), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
            reshape( permute( tmpSksub, [ 2 1 3 ]), size(Sk,2), size(Sk,3)*length(indx) );

      end;
      %% for an N-by-1 kernel, these next lines are enough...
      % A2 = [ A2; ...
      %         reshape( permute(squeeze( Sk(indx,:,:) ),[ 2 3 1 ]), ...
      %                  size(s,2), size(s,3)*length(indx) )
      %       ];
      % 
      % b2 = [b2; squeeze( Sk( plist{cnt}(cnt2), :, l ) ).' ];
    end;

    % normalize, following to Algo 4.2 of M. Schneider, "GPGPU for Accelerated GRAPPA Autocalibration in
    %    Magnetic Resonance Imaging," Masters Thesis, U of Erlangen, Apr 2008
    % A_bak = A;
    p_src = sum(abs(A),2).^(-eta/2);    % compute the power in each line of the linsys
    A = tmult(A,diag(p_src),1);         % normalize by the power in each linsys line

    % Apinv = pinv(A);
    % AA = A'*A;
    AA = [A ]'*[A ];
    At = [A ]';
    clear b2;
    for l=1:size(Sk,3),
      rptupd(prnt,'.');
      % A is the same for each coil ...
      b = vec( squeeze(Sk( plist{cnt}, :, l )).' ); % for dual-kernel, this defaults
                                                    % to the target data

      b = b.*p_src(:);                  % normalize by the power in each linsys line

      b2(:,l) = b;
      n(:,l) = A \ b;

      if (verbose>2), pltcmplx( At'*n(:,l), b ); keyboard; end;
    end;
    b1 = b2;
    N(length(N)).n = n;
    
    rptupd(prnt,'\n');
  end;  
end;

rptupd(prnt,'recon missing lines\n');
%% now go recon the missing lines of k-space:

if (verbose); keyboard; end;

kypts = 1:length(mss);
kypts( 1:floor(length(mss)/2) ) = kypts( (floor(length(mss)/2)):-1:1 );

linerpt = repmat(' ',1,length(mss));

for kycnt=1:length(mss),
  cnt2 = kypts(kycnt);
  
  if verbose, rptupd(prnt,sprintf(' line %3d, from lines ',cnt2)); end;

  if ( find( Y == cnt2) ); 
    % rptupd(prnt,'*'); 
    linerpt(cnt2) = '*';
    if verbose, rptupd(prnt,'\n'); end;
    continue; 
  end;

  tmp = [];
  %% determine which recon pattern we can use:
  for cnt4 = 1:length(pattern),
    v = extrapolate_pattern(pattern{cnt4},nky);

    if ( (kycnt < length(mss)/2) && sum(find( pattern{cnt4} == '-')) ), continue, end;
    if ( (kycnt > length(mss)/2) && sum(find( pattern{cnt4} == '+')) ), continue, end;

    indx  = cnt2 + v;
    % indx2 = cnt2 + f;
    %     tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
    %     if  ( ... % ( (min(indx) > 1) & (max(indx) < size(Sk,1)) ) & ...
    %         ( length(tmp) == length(indx) ) ),
    %       break;                          % jump out of the loop
    %     end;
    if ( indx(1) < 1) | (indx(end) > size(Sk,1)); continue; end;
    
    tmp = (mss( indx ) == ones(length(indx),1));
    if ( sum(tmp) == length(tmp)  ),
      break; 
    end;

  end;  
  
  if ( ( (min(indx) < 1) | (max(indx) > size(Sk,1)) ) | ...
       ( sum(tmp) ~= length(indx) ) ),
    % 
    linerpt(cnt2) = 'o';
    if verbose, rptupd(prnt,'o'); rptupd(prnt,'\n'); end;
    continue; 
  end;
  
  %% recon k-space line based on that pattern:
  rsz = size(Sk,2);
  A = zeros( rsz, length(indx3)*length(v)*size(s,3) );
  B = []; % zeros( rsz, length(indx3)*length(f)*size(s,3) );
  
  for cnt3=1:length(indx3),
    y0 = mod( indx3(cnt3) + [(size(Sk,2)-rsz)/2+(1:rsz)] , size(Sk,2) );
    y0( y0 == 0 ) = size(Sk,2);
    
    A( (1:rsz), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
        reshape( permute( Sk(indx,y0,:), [ 2 1 3 ]), ...
                 rsz, size(Sk,3)*length(indx) );
  end;

  if (verbose), % ( sum(Sk(Y(mss(cnt2))+1,:,l)) ~= 0 ); 
    rptupd(prnt,sprintf(' %3d',[ indx(:) ]));
    rptupd(prnt,sprintf(' (pattern: %s)', pattern{cnt4}) );
    rptupd(prnt,'\n');
    if (verbose>1),
      figure(1); imagesc( log10(abs(Sk(:,:,1)) + 1e-1) ); pause(0.5);
      keyboard;
    end;
  end;
  linerpt(cnt2) = 'x';

  for cnt3=1:length(N),
    if strcmp( N(cnt3).pattern, pattern{cnt4} ),
      break;                            % jump out of the loop
    end;
  end;
  n = N(cnt3).n;
  
  for l=1:size(Sk,3),
    % for each coil ...

    % if cnt2==109, keyboard; end;
    % rptupd(prnt,(' line %3d, pattern %3d\n', cnt2, plist(cnt3) );
    Sk( cnt2, :, l ) = [A B] * n(:,l);
  end;

  if (verbose>1), keyboard; end;
end;
rptupd(prnt,linerpt);
rptupd(prnt,'  \n');

% Sk(1:2:size(Sk,1),:,:) = -Sk(1:2:size(Sk,1),:,:);
% Sk(:,1:2:size(Sk,2),:) = -Sk(:,1:2:size(Sk,2),:);

F = Sk;
I = sqrt( sum(abs(ifft2(Sk)).^2,3) ); 
I = I./(prod(size(I)));
mask = 0;

rpttime = etime(clock,st);
rptupd(prnt,sprintf(' GRAPPA recon time: %g (avg per kernel: %g)\n',  ... 
        rpttime, rpttime/length(pattern)) ); 

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rptupd( prnt, str )

if (prnt)
  fprintf( str );
end;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = extrapolate_pattern(p,YKS)
 
r = find( (p == 'x') | (p == '<') | (p == '>') );
v = find( (p == '*') | (p == '<') | (p == '>') ) - r;

if length(v)>YKS,
  v = v( find( v~= 0 ) );
end;
return;


