function [post nlZ dnlZ] = infGrid(hyp, mean, cov, lik, x, y, opt)

% Inference for a GP with Gaussian likelihood and covGrid covariance.
% The (Kronecker) covariance matrix used is given by:
%   K = kron( kron(...,K{2}), K{1} ) = K_p x .. x K_2 x K_1.
%
% Compute a parametrization of the posterior, the negative log marginal
% likelihood and its derivatives w.r.t. the hyperparameters.
% The result is exact for complete grids, otherwise results are approximate.
% See also "help infMethods".
%
% The function takes a specified covariance function (see covFunctions.m) and
% likelihood function (see likFunctions.m), and is designed to be used with
% gp.m and in conjunction with covFITC and likGauss.
%
% In case of equispaced data points, we use Toeplitz algebra. We offer to use a
% circulant embedding of the Toeplitz matrix and Fourier mechanics.
%
% Copyright (c) by Hannes Nickisch and Andrew Wilson, 2015-09-26.
%
% See also INFMETHODS.M, COVGRID.M.

if iscell(lik), likstr = lik{1}; else likstr = lik; end
if ~ischar(likstr), likstr = func2str(likstr); end
if ~strcmp(likstr,'likGauss')               % NOTE: no explicit call to likGauss
  error('Inference with infGrid only possible with Gaussian likelihood.');
end
if ~isfield(hyp,'mean'), hyp.mean = []; end  % set an empty mean hyper parameter
cov1 = cov{1}; if isa(cov1, 'function_handle'), cov1 = func2str(cov1); end
if ~strcmp(cov1,'covGrid'); error('Only covGrid supported.'), end    % check cov

xg = cov{3}; p = numel(xg);                 % extract underlying grid parameters
if nargin<=6, opt = []; end                        % make opt variable available
if isfield(opt,'proj_ortho'), proj_ortho = opt.proj_ortho;    % enforce P*P' = I
else proj_ortho = false; end, hP = 1;                     % no projection at all
if isfield(hyp,'P')                  % apply transformation matrix P if provided
  hP = hyp.P; if proj_ortho, hP = sqrtm(hP*hP')\hP; end  % orthonormal projector
end

% START
additive_mode = numel(cov)>3;
if additive_mode
  Q = cov{4}; if numel(cov)>4, w=cov{5}; else w={@meanZero}; end,cov = cov(1:3);
  fprintf('%d additive scaling terms\n',Q)
end
% END

[K,M] = feval(cov{:}, hyp.cov, x*hP');    % evaluate covariance mat constituents
[xx,ng,Dg] = covGrid('idx2dat',xg,x);             % turn into data vector if idx
kronmvm = K.kronmvm;
if isfield(opt,'stat'), stat = opt.stat; else stat = false; end  % report status
if stat      % show some information about the nature of the p Kronecker factors
  s = 'K = ';
  for i=1:p
    if isnumeric(K.kron(i).factor)
      si = sprintf('mat(%d)',size(K.kron(i).factor,1));
    else
      sz = num2str(size(xg{i}{1},1));
      for j=2:numel(xg{i}), sz = [sz,'x',num2str(size(xg{i}{j},1))]; end
      si = sprintf('%s(%s)',K.kron(i).factor.descr(1:4),sz);
    end
    if i<p, si = [si,' x ']; end
    s = sprintf('%s%s',s,si);
  end
  fprintf('%s\n',s)
end
n = size(x,1); D = sum(Dg); N = prod(ng);                                 % dims
mid = n==N; if mid, q = randn(N,1); mid = max(abs(q-M*q))<1e-10; end % M==eye(N)
m = feval(mean{:}, hyp.mean, xx*hP');                     % evaluate mean vector

if isfield(opt,'cg_tol'), cgtol = opt.cg_tol;         % stop conjugate gradients
else cgtol = 1e-5; end
if isfield(opt,'cg_maxit'), cgmit = opt.cg_maxit;      % number of cg iterations
else cgmit = min(n,2000); end

sn2 = exp(2*hyp.lik);                               % noise variance of likGauss
circ_emb = -5;                                % Whittle embedding overlap factor
[V,e,ir] = eigkron(K,xg,ng,circ_emb);                 % eigenvalue decomposition
vnum = true;              % flag indicating whether all entries in V are numeric
for j=1:p, vnum = vnum && isnumeric(V{j}); end

% START
vnum = vnum && ~additive_mode;
if additive_mode   % set up K.mvm() for inference
  [Kq,Sq,dKq,dWq] = additive_scaled_mvm_prec(N,Q,D,w,cov,hyp.cov,x,xg,hP);
  Vq = cell(Q,1); Eq = zeros(N,Q);
  for j=1:Q, [Vq{j},Eq(:,j)] = eigkron(Kq{j},xg,ng,circ_emb); end
  mvm = @(t) additive_scaled_mvm(t,Kq,Sq);
  K = struct('mvm',mvm);
  nw = eval(feval(w{:})); nc = eval(feval(cov{:}));
  dKi = @(t,i) additive_scaled_mvm(t,Kq,Sq,dKq,dWq,nc,nw,i);         % deriv mvm

  % See also Prob.III.6.14 in Matrix Analysis, Bhatia 1997.
  % lam(A*B) < lam(A)*lam(B), where lam(.) returns sorted eigvals
  [junk,Si] = sort(Sq,1,'descend'); Si = Si + N*ones(N,1)*(0:Q-1);
  [junk,Ei] = sort(Eq,1,'descend'); Ei = Ei + N*ones(N,1)*(0:Q-1);
  SES = Eq(Ei).*Sq(Si).^2;
  [ldu,I] = logdet_weyl(SES); Ei = Ei(I); Si = Si(I);   % propagate Weyl indices
  e = sum(SES(I),2);                              % e = sum(Eq(Ei).*Sq(Si).^2,2)
end
% END

if mid
  s = 1./(e+sn2); ord = 1:N;                  % V*diag(s)*V' = inv(K+sn2*eye(N))
else
  [eord,ord] = sort(e,'descend');              % sort long vector of eigenvalues
  if n>N, eord = [eord;zeros(n-N,1)]; ord = [ord;(N+1:n)']; end   % special case
  s = 1./((n/N)*eord(1:n) + sn2);               % approx using top n eigenvalues
end

if mid && vnum                    % decide between Toeplitz or Kronecker algebra
  L=@(k) real(kronmvm(V,repmat(1/sn2-s(ir),1,size(k,2)).*kronmvm(V,k,1)))-k/sn2;
else                                                % go for conjugate gradients
  mvm = @(t) mvmK(t,K,sn2,M);                             % single parameter mvm
  L = @(k) -solveMVM(k,mvm,cgtol,cgmit);
end
alpha = -L(y-m);

post.alpha = alpha;                            % return the posterior parameters
post.sW = ones(n,1)/sqrt(sn2);                         % sqrt of noise precision
post.L = L;                           % function to compute inv(K+sn2*eye(N))*Ks
Mtal = M'*alpha;                              % blow up alpha vector from n to N
post.fmu = @(xs) pred_mu(xs,K.mvm(Mtal),xg,hyp,mean);  % interpolate latent mean
post.ymu = post.fmu;       % predictive mean is the same for Gaussian likelihood
if isfield(opt,'pred_var') && opt.pred_var
  ns = 20;                          % number of samples for perturb-and-MAP ...
                     % ... for details see George Papandreou and Alan L. Yuille:
  % "Efficient Variational Inference in Large-Scale Bayesian Compressed Sensing"
  % explained variance on the grid vg=diag(Ku*M'* inv(M*Ku*M'+sn2*eye(n)) *M*Ku)
  % relative accuracy r = std(vg_est)/vg_exact = sqrt(2/ns)
  z = kronmvm(V,repmat(sqrt(e(ir)),1,ns).*kronmvm(V,randn(N,ns),1)); % z~N(0,Ku)
  z = M*z + sqrt(sn2)*randn(n,ns);              % z~N(0,A), A=M*Ku*M'+sn2*eye(n)
  z = K.mvm(M'*L(z)); vg = sum(z.*z,2)/ns;            % z~N(0,Ku*M'*inv(A)*M*Ku)
else
  vg = zeros(N,1);                                       % no variance explained
end
post.fs2 = @(xs) pred_s2(xs,vg,xg,hyp,cov);                    % latent variance
post.ys2 = @(xs) post.fs2(xs) + sn2;      % predictive = latent + noise variance

if nargout>1                               % do we want the marginal likelihood?
  lda = -sum(log(s));
  % exact: lda = 2*sum(log(diag(chol(M*kronmvm(K,eye(N))*M'+sn2*eye(n)))));
  nlZ = (y-m)'*alpha/2 + n*log(2*pi)/2 + lda/2;
  if nargout>2                                         % do we want derivatives?

    %size(V{1})
    %size(e(e>0))
    %size(e(e<0))
    %size(s)
    %inv(K.kron.factor);
    %exit

    dnlZ = hyp;     % allocate space for derivatives, define Q=inv(M*K*M'+sn2*I)

    % START
    if additive_mode
      Ei = Ei(1:min(n,N),:); Si = Si(1:min(n,N),:);  % restrict to relevant part
      Eqi = Eq(Ei); Sqi = Sq(Si);
      ee = sum(Eqi.*Sqi.^2,2);
      assert( norm(ee-e(1:min(n,N)))<1e-12 )
      t = ee + N/n*sn2;
      o = sum(log(t))/2;               % complexity penalty up to constant terms
      assert( abs(o+min(n,N)*log(n/N)/2 - lda/2)<1e-12 )
      for i = 1:numel(hyp.cov)
        ii = mod(i-1,nc+nw)+1;       % hyperparameter index within additive term
        ij = 1+(i-ii)/(nc+nw);                      % index of the additive term
        if ii<=nw              % derivative w.r.t. squashing function parameters
          sj_dWij = Sq(:,ij).*dWq{ij}(ii);
          dSq = zeros(size(Sq)); dSq(:,ij) = 1;
          dSq = (dSq - Sq).*(sj_dWij*ones(1,Q));
          u = 2*sum(Eqi.*Sqi.*dSq(Si),2);
        else            % derivatives w.r.t. ordinary covariance hyperparameters
          % d lam(K) = diag(V'*dK*V), for psd matrix K = V*diag(lam)*V'.
          dKjf = dKq{ij}(ii-nw); dKj = struct2cell(dKjf.kron);      % mvm for dK
          Vj = Vq{ij};                                               % mvm for V
          dEqi = 1;                                    % dEqi = diag(Vj'*dKj*Vj)
          for j=1:p    % first compute orthogonal complement, then grow diagonal
            if isnumeric(Vj{j})
              dgj = diag(Vj{j}'*dKj{j}*Vj{j}); nj = size(Vj{j},1);
            else
              [xj,nj] = covGrid('expand',xg{j});          % extract subgrid size
              dkj = circ(dKjf.kron(j).factor.kii,nj,circ_emb);
              dgj = real(Vj{j}.mvmt(dkj))*sqrt(prod(nj));
            end
            if numel(dgj)<prod(nj), dgj(prod(nj),1) = 0; end
            dEqi = kron(dgj,dEqi);
          end
          dEqi = dEqi(Ei(:,ij)-N*(ij-1));         % restrict to relevant indices
          u = dEqi.*Sqi(:,ij).^2;              % only one single term ij is left
        end
        dnlZ.cov(i) = (sum(u./t) - Mtal'*dKi(Mtal,i))/2;
      end
    end
    % END

    for i = 1:numel(hyp.cov)

      % START
      if additive_mode, break, end
      % END

      dK = feval(cov{:}, hyp.cov, x*hP', [], i); dk = cell(size(dK));
      P = cell(size(dK));               % auxiliary Kronecker matrix P = V'*dK*V
      for j=1:p                    % perform the circulant embedding as for K{j}
        if ~isnumeric(dK.kron(j).factor)
          [xj,nj] = covGrid('expand',xg{j});              % extract subgrid size
          dk{j} = circ(dK.kron(j).factor.kii,nj,circ_emb);
        end
        if isstruct(V{j})                   % V{j} is a partial Fourier operator
          P{j} = real(V{j}.mvmt(dk{j}))*sqrt(prod(nj));
        else                                            % V{j} is a dense matrix
          P{j} = real(sum((V{j}'*dK.kron(j).factor).*V{j}.',2));
        end
      end
      dp = 1; for j=1:p, dp = kron(zfill(P{j},ng(j)),dp); end     % dp = diag(P)
      if n>N, dp = [dp; zeros(n-N,1)]; end                        % special case
      dp = (n/N)*dp(ord(1:n));      % approximate if incomplete grid observation
      dnlZ.cov(i) = (dp'*s - Mtal'*dK.mvm(Mtal))/2; % complexity + data fit term
    end

    dnlZ.lik = sn2*(sum(s) - alpha'*alpha);                  % sum(s) = trace(Q)
    for i = 1:numel(hyp.mean)
      dnlZ.mean(i) = -feval(mean{:}, hyp.mean, xx*hP', i)'*alpha;
    end
    if isfield(hyp,'P')
      dnlZ.P = deriv_P(alpha,hP,K,covGrid('flatten',xg),m,mean,hyp,x);
      if proj_ortho, dnlZ.P = chain_ortho(hyp.P,dnlZ.P); end
    end
  end
end

% START
% precompute terms for mvm with additive scaled covariance
function [Kq,Sq,dKq,dWq] = additive_scaled_mvm_prec(N,Q,D,w,cov,hyp,x,xg,hP)
  nw = eval(feval(w{:})); nc = eval(feval(cov{:}));           % number of hypers
  Kq = cell(Q,1); dKq = cell(Q,1); Wq = zeros(N,Q); dWq = cell(Q,1);  % allocate
  for j=1:Q
    hypwq = hyp((j-1)*(nw+nc)+   (1:nw));      % scaling weight hyper parameters
    hypcq = hyp((j-1)*(nw+nc)+nw+(1:nc)); % covariance function hyper parameters
    Wq(:,j) = feval(w{:}, hypwq, covGrid('expand',xg));        % scaling weights
    dWq{j} = @(i) feval(w{:}, hypwq, covGrid('expand',xg), i); % derivative sc-w
    Kq{j} = feval(cov{:}, hypcq, x*hP');                     % covariance matrix
    dKq{j} = @(i) feval(cov{:}, hypcq, x*hP', [], i);    % derivative cov matrix
  end
  Sq = softmax(Wq);

% mvm with additive scaled covariance
function r = additive_scaled_mvm(t,Kq,Sq,dKq,dWq,nc,nw,i)
  r = zeros(size(t)); o = ones(1,size(r,2)); Q = numel(Kq);
  if nargin>7                                                      % derivatives
    ii = mod(i-1,nc+nw)+1;           % hyperparameter index within additive term
    ij = 1+(i-ii)/(nc+nw);                          % index of the additive term
    sj_dWij = Sq(:,ij).*dWq{ij}(ii);                            % factor for dsj
    if ii<=nw                  % derivative w.r.t. squashing function parameters
      for j=1:Q
        dsj = ((double(ij==j)-Sq(:,j)).*sj_dWij)*o; sj = Sq(:,j)*o;  % dsj/dh_ii
        r = r + dsj.*Kq{j}.mvm(sj.*t) + sj.*Kq{j}.mvm(dsj.*t);
      end
    else                % derivatives w.r.t. ordinary covariance hyperparameters
      dKij = dKq{ij}(ii-nw);                   % only term ij is left in the sum
      sj = Sq(:,ij)*o; r = r + sj.*dKij.mvm(sj.*t);
    end
  else                                                        % plain evaluation
    for j=1:Q, sj = Sq(:,j)*o; r = r + sj.*Kq{j}.mvm(sj.*t); end
  end

% softmax along the second dimension
function s = softmax(w)
  Q = size(w,2); oq = ones(1,Q);
  s = exp(w-max(w,[],2)*oq);
  s = s./(sum(s,2)*oq);
  assert(norm(sum(s,2)-1)<1e-12)

% Upper bound the log determinant of a sum of p symmetric positive semi-definite
% matrices of size nxn represented in terms of their eigenvalues E using Weyl's
% inequalities:
% Let a(nx1) and b(nx1) be the orderered eigenvalues of the Hermitian matrices
% A(nxn) and B(nxn). Then, the eigenvalues c(nx1) of the matrix C = A+B are
% upper bounded by c(i+j-1) <= a(i)+b(j).
%
% Each of the p columns of the matrix E of size nxp contains the
% eigenvalues of the respective matrix.
% The string mod contains the method of upper bounding the log determinant:
%  'halve'  - select eigenvalues of roughly equal indices i,j
%  'greedy' - select the locally smallest bound; the integer s is the number of
%             search steps ahead
%
% The index matrix I can be used to reconstruct the contributions of the
% individual eigenvalues to the overall log determinant upper bound given by:
%   ldu = sum(log( sum(E(I),2) ));
function [ldu,I] = logdet_weyl(E,mod,s)
if nargin<2, mod = 'halve'; end
[n,p] = size(E); I = zeros(n,p);                         % get dims, init result
[e,I(:,1)] = sort(E(:,1),'descend');                    % get upper bound so far
i = ceil((1:n)'/2); j = (1:n)'-i+1;        % default: similar index size (halve)
for k=2:p
  [fs,js] = sort(E(:,k),'descend');                     % sort both k-th summand
  if strcmp(mod,'greedy')     % set up i and j adaptively using greedy heuristic
    if nargin<3, s = 1; end                                  % set default value
    i(1) = 1; j(1) = 1;                                                   % init
    for m=2:n
      ss = s;                           % reduce steps if it violates boundaries
      if ss>  i(m-1) || ss>  j(m-1), ss = min(  i(m-1),   j(m-1)); end
      if ss>n-i(m-1) || ss>n-j(m-1), ss = min(n-i(m-1), n-j(m-1)); end
      i_ss = i(m-1)+(1-ss:ss); j_ss = m-i_ss+1;                % left/right step
      [junk,idx] = min(e(i_ss) + fs(j_ss));   % keep min of steps left and right
      i(m) = i_ss(idx); j(m) = j_ss(idx);
    end
  end
  I(:,1:k-1) = I(i,1:k-1);               % sort previous summands according to i
  I(:,k) = js(j) + (k-1)*n;              % bound contributions from k-th summand
  e = e(i) + fs(j);                              % compute Weyl eigenvalue bound
end
ldu = sum(log(e));
assert( norm(e-sum(E(I),2))<1e-10 )
if ~isinf(ldu), assert( abs(ldu-sum(log(sum(E(I),2))))<1e-10 ), end
% END

% compute derivative of neg log marginal likelihood w.r.t. projection matrix P
function dP = deriv_P(alpha,P,K,xg,m,mean,hyp,x)
  xP = x*P'; [M,dM] = covGrid('interp',xg,xP); % grid interp derivative matrices
  beta = K.mvm(M'*alpha);                          % dP(i,j) = -alpha'*dMij*beta
  dP = zeros(size(P)); h = 1e-4;               % allocate result, num deriv step
  for i=1:size(P,1)
    if equi(xg,i), wi = max(xg{i})-min(xg{i}); else wi = 1; end % scaling factor
    xP(:,i) = xP(:,i)+h;
    dmi = (feval(mean{:},hyp.mean,xP)-m)/h;         % numerically estimate dm/di
    xP(:,i) = xP(:,i)-h;
    betai = dmi + dM{i}*beta/wi;
    for j=1:size(P,2), dP(i,j) = -alpha'*(x(:,j).*betai); end
  end

function y = zfill(z,n)                                      % fill z with zeros
  y = accumarray((1:numel(z))',z,[n,1]);

function eq = equi(xg,i)                        % grid along dim i is equispaced
  ni = size(xg{i},1);
  if ni>1                              % diagnose if data is linearly increasing
    dev = abs(diff(xg{i})-ones(ni-1,1)*(xg{i}(2,:)-xg{i}(1,:)));
    eq = max(dev(:))<1e-9;
  else
    eq = true;
  end

% chain rule for the function Q = sqrtm(P*P')\P;  for d sqrtm(X) see the website
function dQ = chain_ortho(P,dP) % http://math.stackexchange.com/questions/540361
  [V,F] = eig(P*P'); sf = sqrt(diag(F)); S = V*diag(sf)*V';         % eig-decomp
  H = dP'/S; G = H'*(P'/S); o = ones(size(dP,1),1);                 % chain rule
  dQ = (H - P'*V*((V'*(G+G')*V)./(sf*o'+o*sf'))*V')';

function q = solveMVM(p,mvm,varargin) % solve q = mvm(p) via conjugate gradients
  [q,flag,relres,iter] = conjgrad(mvm,p,varargin{:});                 % like pcg
  if ~flag,error('Not converged after %d iterations, r=%1.2e\n',iter,relres),end

function q = mvmK(p,K,sn2,M) % mvm q = M*K*M'*p + sn2*p using the Kronecker repr
  q = M*K.mvm(M'*p) + sn2*p;

function fmu = pred_mu(xs,Kalpha,xg,hyp,mean) % predict interpolated latent mean
  Ms = covGrid('interp',xg,xs);                    % obtain interpolation matrix
  xs = covGrid('idx2dat',xg,xs);                        % deal with index vector
  ms = feval(mean{:},hyp.mean,xs);                         % evaluate prior mean
  fmu = ms + Ms*Kalpha;                 % combine and perform grid interpolation

function fs2 = pred_s2(xs,vg,xg,hyp,cov)  % predict interpolated latent variance
  if norm(vg,1)>1e10
    ve = covGrid('interp',xg,xs)*vg;       % interpolate grid variance explained
  else
    ve = 0;
  end
  xs = covGrid('idx2dat',xg,xs);                        % deal with index vector
  ks = feval(cov{:},hyp.cov,xs,'diag');                % evaluate prior variance
  fs2 = max(ks-ve,0);                % combine, perform grid interpolation, clip

% Solve x=A*b with symmetric A(n,n), b(n,m), x(n,m) using conjugate gradients.
% The method is along the lines of PCG but suited for matrix inputs b.
function [x,flag,relres,iter,r] = conjgrad(A,b,tol,maxit)
if nargin<3, tol = 1e-10; end
if nargin<4, maxit = min(size(b,1),20); end
x0 = zeros(size(b)); x = x0;
if isnumeric(A), r = b-A*x; else r = b-A(x); end, r2 = sum(r.*r,1); r2new = r2;
nb = sqrt(sum(b.*b,1)); flag = 0; iter = 1;
relres = sqrt(r2)./nb; todo = relres>=tol; if ~any(todo), flag = 1; return, end
on = ones(size(b,1),1); r = r(:,todo); d = r;
for iter = 2:maxit
  if isnumeric(A), z = A*d; else z = A(d); end
  a = r2(todo)./sum(d.*z,1);
  a = on*a;
  x(:,todo) = x(:,todo) + a.*d;
  r = r - a.*z;
  r2new(todo) = sum(r.*r,1);
  relres = sqrt(r2new)./nb; cnv = relres(todo)<tol; todo = relres>=tol;
  d = d(:,~cnv); r = r(:,~cnv);                           % get rid of converged
  if ~any(todo), flag = 1; return, end
  b = r2new./r2;                                               % Fletcher-Reeves
  d = r + (on*b(todo)).*d;
  r2 = r2new;
end

% Eigendecomposition of a Kronecker matrix with dense, Toeplitz or BTTB factors.
% K.mvm(z) == kronmvm(V,e.*kronmvm(V,z,1)), approximate for Toeplitz/BTTB
function [V,e,ir] = eigkron(K,xg,ng,circ_emb)
isbttb = @(Ki) isstruct(Ki) && (strcmp (Ki.descr,'toep') ...
                             || strncmp(Ki.descr,'bttb',4));   % BTTB covariance
p = numel(K.kron); V = cell(p,1); E = cell(p,1);     % sizes and allocate memory
e = 1; ir = 1;      % compute eigenvalue diagonal matrix full/reduced rank index
for j=1:p
  if isbttb(K.kron(j).factor)
    [xj,nj] = covGrid('expand',xg{j}); Nj = prod(nj);     % extract subgrid size
    ej = max(real(fftn( circ(K.kron(j).factor.kii,nj,circ_emb) )),0);
    sNj = sqrt(Nj);                            % eigenvalues of circulant matrix
    ep = ej>0; Ej = ej(ep); sub = @(s,ep) s(ep);    % non-negative eigenvalues
    V{j}.mvmt = @(v) sub(fftn(v)/sNj,ep);       % V{j} is partial Fourier matrix
    V{j}.mvm = @(v) ifftn(reshape(accumarray(find(ep),v,[Nj,1]),[nj,1]))*sNj;
    V{j}.size = [numel(ep),sum(ep)];
  else [V{j},Ej] = eigr(K.kron(j).factor,0); % eigenvalues of non-Toeplitz mat
  end
  E{j} = Ej;
end
for j=1:p
  de = zfill(E{j},ng(j)); e = kron(de,e); ir = kron(de>0,ir);
end
ir = ir~=0;

% Real eigenvalues and eigenvectors up to the rank of a real symmetric matrix.
% Decompose A into V*D*V' with orthonormal matrix V and diagonal matrix
% D = diag(d)represented by a column vector.
% Entries of d with index larger tha the rank r of the matrix A as returned by
% the call rank(A,tol) are zero.
function [V,d] = eigr(A,econ,tol)
if nargin<2, econ = false; end                          % assign a default value
if isnumeric(econ), econ = econ==0; end                      % turn into boolean
if ischar(econ), econ = strcmp(econ,'econ'); end             % turn into boolean
if nargout==0, return, end                     % avoid computation in limit case
if nargout==1
  d = sort(eig((A+A')/2),'descend');       % eigenvalues of strictly symmetric A
else
  [V,D] = eig((A+A')/2);                 % decomposition of strictly symmetric A
  d = max(real(diag(D)),0); [d,ord] = sort(d,'descend');      % tidy up and sort
end
n = size(A,1);                                                  % dimensionality
if nargin<3, tol = n*eps(max(d)); end, r = sum(d>tol);             % get rank(A)
d(r+1:n) = 0;                             % set junk eigenvalues to strict zeros
if econ, d = d(1:r); end                                   % truncate if desired
if nargout==1                                        % only eigenvalues required
  V = d;
else                                             % entire decomposition required
  V(:,1:r) = real(V(:,ord(1:r)));                    % eigenvectors up to rank r
  if econ
    V = V(:,1:r);
  else                                            % ortho completion if required
    V(:,r+1:n) = null(V(:,1:r)');
  end
end

% Construct a circular embedding c(nx1) from a covariance function k.
%  - k is a function and the call k(1:n) returns a vector of length n
%  - s is the setting for the embedding with values s={..,-2,-1,0,1,2,3,4,5}.
%
% s<0  Whittle embedding [2] as described by Guinness & Fuentes [3] in
%      equation (5) with N = |s|.
% s=0  No embedding c = k(n).
% s=1  Generic embedding c = [k(1:n/2), k(n/2:-1:2)] without smoothing
%      as in the Strang preconditioner.
% s=2  Variant of smoothing inspired by [1] respecting the conditions
%      c(1) = k(1), c(i)=c(n-i+2) for i=2..n and idempotence. We use the generic
%      scheme c(2:n) = w.*k(2:end) + flipud(w.*k(2:end)), where w are sigmoidal
%      interpolation weights.
% s=3  Linear smoothing as in T. Chan's preconditioner.
% s=4  Tyrtyshnikov preconditioner.
% s=5  R. Chan's preconditioner.
%
% [1] Helgason, Pipiras & Abry, Smoothing windows for the synthesis of
%     Gaussian stationary random fields using circulant matrix embedding,
%     Journal of Computational and Graphical Statistics, 2014, 23(3).
% [2] Whittle, On stationary processes in the plane, Biometrika, 1954, 41(3/4).
% [3] Guinness & Fuentes, Circulant embedding of approximate covariances for
%     inference from Gaussian data on large lattices, 2014, preprint,
%     http://www4.stat.ncsu.edu/~guinness/circembed.html.
function c = circ(k,n,s)
p = numel(n); n = n(:)';                             % dimensions and row vector
if nargin<3, s = -2; end                                         % default value
if s==0                                                    % no embedding at all
  xg = cell(p,1); for i=1:p, xg{i} = (1:n(i))'; end              % standard grid
  xc = covGrid('expand',xg);
  c = reshape(k(xc),[n,1]);
elseif s<0                                           % Whittle/Guinness aliasing
  N = abs(s);
  xg = cell(p,1); for i=1:p, xg{i} = (1-N*n(i):N*n(i))'; end
  sz = [n; 2*N*ones(1,p)]; sz = sz(:)';
  c = reshape(k(covGrid('expand',xg)),sz);
  for i=1:p, c = sum(c,2*i); end, c = squeeze(c);
else
  if     s==1                                                 % Strang embedding
    xg = cell(p,1);
    for i=1:p, n2 = floor(n(i)/2)+1; xg{i} = [1:n2, n2-n(i)+1:0]'; end
    xc = covGrid('expand',xg);
    c = reshape(k(xc),[n,1]);
  elseif s==2                                               % Helgason smoothing
    if numel(n)>1, error('Only 1d allowed'), end
    k = k(1:n); k = k(:);
    k0 = k(2:n); % smooth by sigmoid interpolation between the two flipped parts
    wk = k0./(1+exp( 30*(((1:n-1)-n/2)'/n) ));           % sigmoid-weighted part
    c = [k(1); wk + flipud(wk)];
  elseif s==3                                       % T. Chan's linear smoothing
    if numel(n)>1, error('Only 1d allowed'), end
    k = k(1:n); k = k(:);
    w = (1:n-1)'/n;
    c = k; c(2:n) = (1-w).*c(2:n) + w.*flipud(c(2:n));
  elseif s==4                                                     % Tyrtyshnikov
    if numel(n)>1, error('Only 1d allowed'), end
    k = k(1:n); k = k(:);
    w = (1:n-1)'/n;
    c = k; c2 = c(2:n).*c(2:n);
    c(2:n) = ((1-w).*c2 + w.*flipud(c2)) ./ ((1-w).*c(2:n) + w.*flipud(c(2:n)));
  elseif s==5                                                          % R. Chan
    if numel(n)>1, error('Only 1d allowed'), end
    k = k(1:n); k = k(:);
    c = k; c(2:n) = c(2:n) + flipud(c(2:n));                   % mirrored add-up
  end
end
assert(numel(c)==prod(n))                                               % length
c2 = c;
for i=1:p
  c2 = reshape(c2,[prod(n(1:i-1)-1),n(i),prod(n(i+1:end))]); c2 = c2(:,2:end,:);
end
c2f = c2; for i=1:p, c2f = flipdim(c2f,i); end
if s~=0, assert(max(reshape(abs(c2-c2f),[],1))<1e-10), end         % circularity
