text:
function [g, gdata, gprior] = gpla_g(w, gp, x, y, varargin)
%GPLA_G   Evaluate gradient of Laplace approximation's marginal
%         log posterior estimate (GPLA_E)
%
%  Description
%    G = GPLA_G(W, GP, X, Y, OPTIONS) takes a full GP parameter
%    vector W, structure GP a matrix X of input vectors and a
%    matrix Y of target vectors, and evaluates the gradient G of
%    EP's marginal log posterior estimate. Each row of X
%    corresponds to one input vector and each row of Y corresponds
%    to one target vector.
%
%    [G, GDATA, GPRIOR] = GPLA_G(W, GP, X, Y, OPTIONS) also returns
%    the data and prior contributions to the gradient.
%
%    OPTIONS is optional parameter-value pair
%      z - optional observed quantity in triplet (x_i,y_i,z_i)
%          Some likelihoods may use this. For example, in case of
%          Poisson likelihood we have z_i=E_i, that is, expected
%          value for ith case.
%
%  See also
%    GP_SET, GP_G, GPLA_E, GPLA_PRED

% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GPLA_G';
  ip.addRequired('w', @(x) isvector(x) && isreal(x) && all(isfinite(x)));
  ip.addRequired('gp',@isstruct);
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.parse(w, gp, x, y, varargin{:});
  z=ip.Results.z;

  gp = gp_unpak(gp, w);       % unpak the parameters
  ncf = length(gp.cf);
  n=size(x,1);

  g = [];
  gdata = [];
  gprior = [];

  if isfield(gp, 'savememory') && gp.savememory
    savememory=1;
  else
    savememory=0;
  end

  % First Evaluate the data contribution to the error
  switch gp.type
    case 'FULL'
    % ============================================================
    % FULL
    % ============================================================

    if ~isfield(gp.lik, 'nondiagW')
      % Likelihoos with diagonal Hessian

      % Calculate covariance matrix and the site parameters
      K = gp_trcov(gp,x);
      if isfield(gp,'meanf')
        [H,b_m,B_m]=mean_prep(gp,x,[]);
        K=K+H'*B_m*H;
      end

      [e, edata, eprior, f, L, a, W, p] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
      if isnan(e)
        g=NaN; gdata=NaN; gprior=NaN;
        return;
      end

      if W >= 0              % This is the usual case where likelihood is log concave
        % for example, Poisson and probit
        if issparse(K)                               % use sparse matrix routines

          % permute
          y = y(p);
          x = x(p,:);
          K = K(p,p);
          if ~isempty(z)
            z = z(p,:);
          end

          sqrtW = sqrt(W);

          R = sqrtW*spinv(L,1)*sqrtW;
          sqrtWK = sqrtW*K;
          C = ldlsolve(L,sqrtWK);
          C2 = diag(K) - sum(sqrtWK.*C,1)';
          s2 = 0.5*C2.*gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
        else                                         % evaluate with full matrices
          sqrtW = diag(sqrt(W));
          c = L\sqrtW;
          R = c'*c;
          C2 = diag(K) - sum((c*K).^2,1)' ;
          s2 = 0.5*C2.*gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
        end
      else                         % We might end up here if the likelihood is not log-concave
        % For example Student-t likelihood.
        C = L;
        V = L*diag(W);
        R = diag(W) - V'*V;
        C2 = sum(C.^2,1)';
        s2 = 0.5*C2.*gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
      end

      % =================================================================
      % Gradient with respect to covariance function parameters
      if ~isempty(strfind(gp.infer_params, 'covariance'))
        % Evaluate the gradients from covariance functions
        for i=1:ncf
          i1=0;
          if ~isempty(gprior)
            i1 = length(gprior);
          end

          gpcf = gp.cf{i};
          if savememory
            % If savememory option is used, just get the number of
            % hyperparameters and calculate gradients later
            np=gpcf.fh.cfg(gpcf,[],[],[],0);
          else
            DKffc = gpcf.fh.cfg(gpcf, x);
            np=length(DKffc);
          end
          gprior_cf = -gpcf.fh.lpg(gpcf);

          g1 = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);
          for i2 = 1:np
            if savememory
              DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
            else
              DKff=DKffc{i2};
            end
            i1 = i1+1;
            if ~isfield(gp,'meanf')
              s1 = 0.5 * a'*DKff*a - 0.5*sum(sum(R.*DKff));
            else
              s1 = 0.5 * (a-K\(H'*b_m))'*DKff*(a-K\(H'*b_m)) - 0.5*sum(sum(R.*DKff));
            end
            b = DKff * g1;
            if issparse(K)
              s3 = b - K*(sqrtW*ldlsolve(L,sqrtW*b));
            else
              s3 = b - K*(R*b);
              %s3 = (1./W).*(R*b);
            end
            gdata(i1) = -(s1 + s2'*s3);
            gprior(i1) = gprior_cf(i2);
          end

          % Set the gradients of hyperparameter
          if length(gprior_cf) > np
            for i2=np+1:length(gprior_cf)
              i1 = i1+1;
              gdata(i1) = 0;
              gprior(i1) = gprior_cf(i2);
            end
          end
        end

      end

      % =================================================================
      % Gradient with respect to likelihood function parameters
      if ~isempty(strfind(gp.infer_params, 'likelihood')) ...
          && ~isempty(gp.lik.fh.pak(gp.lik))

        gdata_lik = 0;
        lik = gp.lik;

        g_logPrior = -lik.fh.lpg(lik);
        if ~isempty(g_logPrior)

          DW_sigma = lik.fh.llg3(lik, y, f, 'latent2+param', z);
          DL_sigma = lik.fh.llg(lik, y, f, 'param', z);
          b = K * lik.fh.llg2(lik, y, f, 'latent+param', z);
          s3 = b - K*(R*b);
          nl= size(DW_sigma,2);

          gdata_lik = - DL_sigma - 0.5.*sum(repmat(C2,1,nl).*DW_sigma) - s2'*s3;

          % set the gradients into vectors that will be returned
          gdata = [gdata gdata_lik];
          gprior = [gprior g_logPrior];
          i1 = length(g_logPrior);
          i2 = length(gdata_lik);
          if i1  > i2
            gdata = [gdata zeros(1,i1-i2)];
          end
        end
      end

      g = gdata + gprior;

    else
      % Likelihoods with non-diagonal Hessian

      [n,nout]=size(y);
      if isfield(gp, 'comp_cf')  % own covariance for each ouput component
        multicf = true;
        if length(gp.comp_cf) ~= nout && nout > 1
          error('GPLA_ND_G: the number of component vectors in gp.comp_cf must be the same as number of outputs.')
        end
      else
        multicf = false;
      end

      % Get help parameters
      [e, edata, eprior, f, L, a, E, M] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
      if isnan(e)
        return
      end

      switch gp.lik.type

        case {'LGP', 'LGPC'}

          nl=n;
          nlp=length(nl); % number of latent processes

          if isfield(gp.latent_opt, 'kron') && gp.latent_opt.kron==1
            gptmp=gp; gptmp.jitterSigma2=0;
            Ka = gp_trcov(gptmp, unique(x(:,1)));
            wtmp=gp_pak(gptmp); wtmp(1)=0; gptmp=gp_unpak(gptmp,wtmp);
            Kb = gp_trcov(gptmp, unique(x(:,2)));
            clear gptmp
            n1=size(Ka,1);
            n2=size(Kb,1);
          else
            K = gp_trcov(gp,x);
          end

          if isfield(gp,'meanf')
            [H,b_m,B_m]=mean_prep(gp,x,[]);
            Hb_m=H'*b_m;
            if isfield(gp.latent_opt, 'kron') && gp.latent_opt.kron==1
              % only zero mean function implemented for Kronecker
              % approximation
              iKHb_m=zeros(n,1);
            else
              K=K+H'*B_m*H;
              iKHb_m=K\Hb_m;
            end
          end

          g2 = gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
          ny=sum(y);

          g3=gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);

          if isfield(gp.latent_opt, 'kron') && gp.latent_opt.kron==1

            [Va,Da]=eig(Ka); [Vb,Db]=eig(Kb);
            % eigenvalues of K matrix
            Dtmp=kron(diag(Da),diag(Db));
            [sDtmp,istmp]=sort(Dtmp,'descend');

            % Form the low-rank approximation. Exclude eigenvalues
            % smaller than gp.latent_opt.eig_tol or take
            % gp.latent_opt.eig_prct*n eigenvalues at most.
            nlr=min([sum(sDtmp>gp.latent_opt.eig_tol) round(gp.latent_opt.eig_prct*n)]);
            sDtmp=sDtmp+gp.jitterSigma2;

            itmp1=meshgrid(1:n1,1:n2);
            itmp2=meshgrid(1:n2,1:n1)';
            ind=[itmp1(:) itmp2(:)];

            % included eigenvalues
            Dlr=sDtmp(1:nlr);
            % included eigenvectors
            Vlr=zeros(n,nlr);
            for i1=1:nlr
              Vlr(:,i1)=kron(Va(:,ind(istmp(i1),1)),Vb(:,ind(istmp(i1),2)));
            end

            Lb=gp_trvar(gp,x)-sum(bsxfun(@times,Vlr.*Vlr,Dlr'),2);

            if isfield(gp,'meanf')
              Dt=[Dlr; diag(B_m)];
              Vt=[Vlr H'];
            else
              Dt=Dlr;
              Vt=Vlr;
            end

            Lbt=ny*(g2)+1./Lb;

            St=[diag(1./Dt)+Vt'*bsxfun(@times,1./Lb,Vt) zeros(size(Dt,1),1); ...
              zeros(1,size(Dt,1)) 1];
            Pt=[bsxfun(@times,1./Lb,Vt) sqrt(ny)*g2];
            Ptt=bsxfun(@times,1./sqrt(Lbt),Pt);

            StL=chol(St-Ptt'*Ptt,'lower');
            iStL=StL\(bsxfun(@times,Pt',1./Lbt'));

            dA=(1./Lbt+sum(iStL.*iStL)');
            iStLg3=iStL*g3;
            const1=( 0.5*ny*(sum( dA.*g3))-ny*(g3'*(g3./Lbt)+iStLg3'*iStLg3) );
            const2=(g3./Lbt)+iStL'*iStLg3;

            s2=const1.*g3 - 0.5*ny*dA.*g3 + ny*const2.*g3;

          else
            if strcmpi(gp.lik.type,'LGPC')
              n1=gp.lik.gridn(1); n2=gp.lik.gridn(2);
              ny2=sum(reshape(y,fliplr(gp.lik.gridn)));
              g2sq=sqrt(g2);

              R=zeros(n);
              RR=zeros(n,n2);
              for k1=1:n1
                R((1:n2)+(k1-1)*n2,(1:n2)+(k1-1)*n2)=sqrt(ny2(k1))*(diag(g2sq((1:n2)+(k1-1)*n2))-g2((1:n2)+(k1-1)*n2)*g2sq((1:n2)+(k1-1)*n2)');
                RR((1:n2)+(k1-1)*n2,:)=R((1:n2)+(k1-1)*n2,(1:n2)+(k1-1)*n2)*R((1:n2)+(k1-1)*n2,(1:n2)+(k1-1)*n2)';
              end
              KW=K;
              for k1=1:n1
                KW(:,(1:n2)+(k1-1)*n2)=KW(:,(1:n2)+(k1-1)*n2)*RR((1:n2)+(k1-1)*n2,:);
              end
              %KW=K*(R*R');

              KW(1:(n+1):end)=KW(1:(n+1):end)+1;
              iKW=KW\eye(n);
              A=iKW*K;

              s2=zeros(n,1);
              for k1=1:n1
                if ny2(k1)~=0
                  g3tmp=g3((1:n2)+(k1-1)*n2);
                  Atmp=A((1:n2)+(k1-1)*n2,(1:n2)+(k1-1)*n2);
                  for ind2=1:n2
                    g3dtmp=-g3tmp*g3tmp(ind2);
                    g3dtmp(ind2)=g3dtmp(ind2)+g3tmp(ind2);
                    s2( ind2+(k1-1)*n2 ) = -ny2(k1)*0.5*sum(diag(Atmp).*g3dtmp) ...
                      + ny2(k1)*sum(sum(Atmp.*(bsxfun(@times,g3tmp,g3dtmp'))));
                  end
                end
              end

            else
              KW=-(K*(sqrt(ny)*g2))*(sqrt(ny)*g2)'- bsxfun(@times, K, (-ny*g2)');

              KW(1:(n+1):end)=KW(1:(n+1):end)+1;
              iKW=KW\eye(n);
              A=iKW*K;

              const1=( 0.5*ny*(sum(A(1:(n+1):end).*g3'))-ny*sum(sum(A.*(g3*g3'))) );
              const2=sum(bsxfun(@times,A,g3));
              s2=const1.*g3 - 0.5*ny*diag(A).*g3 + ny*const2'.*g3;
            end
          end

          % =================================================================
          % Gradient with respect to covariance function parameters
          if ~isempty(strfind(gp.infer_params, 'covariance'))
            % Evaluate the gradients from covariance functions
            for i=1:ncf

              i1=0;
              if ~isempty(gprior)
                i1 = length(gprior);
              end

              gpcf = gp.cf{i};
              if isfield(gp.latent_opt, 'kron') && gp.latent_opt.kron==1

                gptmp=gp; gptmp.jitterSigma2=0;
                DKa = gpcf.fh.cfg(gptmp.cf{1}, unique(x(:,1)));
                wtmp=gp_pak(gptmp); wtmp(1)=0; gptmp=gp_unpak(gptmp,wtmp);
                DKb = gpcf.fh.cfg(gptmp.cf{1}, unique(x(:,2)));
                clear gptmp

                for j1=1:length(DKa)
                  [DVa{j1},DDa{j1}]=eig(DKa{j1});
                end
                for j1=1:length(DKb)
                  [DVb{j1},DDb{j1}]=eig(DKb{j1});
                end

                % low-rank approximation of the derivative matrix w.r.t.
                % magnitude hyperparameter, low-rank + diagonal correction
                Datmp=kron(diag(DDa{1}),diag(Db));
                [sDatmp,isatmp]=sort(Datmp,'descend');
                nlr=min([sum(sDtmp>gp.latent_opt.eig_tol) round(gp.latent_opt.eig_prct*n)]);
                Dmslr=sDatmp(1:nlr);
                Vmslr=zeros(n,nlr);
                for j1=1:nlr
                  Vmslr(:,j1)=kron(DVa{1}(:,ind(isatmp(j1),1)),Vb(:,ind(isatmp(j1),2)));
                end
                % diagonal correction
                dc=gpcf.fh.cfg(gpcf, x(1,:));
                Lms=dc{1}*ones(n,1)-sum(bsxfun(@times,Vmslr.^2,Dmslr'),2);

                % low-rank approximation of the derivative matrix w.r.t.
                % lengthscale hyperparameter, low-rank1 + lowr-rank2 + diagonal correction
                %
                % first low-rank part
                Datmp=kron(diag(DDa{2}),diag(Db));
                [sDatmp,isatmp]=sort(abs(Datmp),'descend');
                nlr=min([sum(sDtmp>gp.latent_opt.eig_tol) round(gp.latent_opt.eig_prct*n)]);
                sDatmp=Datmp(isatmp);
                Dlslr1=sDatmp(1:nlr);
                Vlslr1=zeros(n,nlr);
                for j1=1:nlr
                  Vlslr1(:,j1)=kron(DVa{2}(:,ind(isatmp(j1),1)),Vb(:,ind(isatmp(j1),2)));
                end
                % second low-rank part
                Dbtmp=kron(diag(Da),diag(DDb{2}));
                [sDbtmp,isbtmp]=sort(abs(Dbtmp),'descend');
                nlr=min([sum(sDtmp>gp.latent_opt.eig_tol) round(gp.latent_opt.eig_prct*n)]);
                sDbtmp=Dbtmp(isbtmp);
                Dlslr2=sDbtmp(1:nlr);
                Vlslr2=zeros(n,nlr);
                for j1=1:nlr
                  Vlslr2(:,j1)=kron(Va(:,ind(isbtmp(j1),1)),DVb{2}(:,ind(isbtmp(j1),2)));
                end
                % diagonal correction
                Lls=dc{2}*ones(n,1)-sum(bsxfun(@times,Vlslr1.^2,Dlslr1'),2)-sum(bsxfun(@times,Vlslr2.^2,Dlslr2'),2);

              else
                if savememory
                  % If savememory option is used, just get the number of
                  % hyperparameters and calculate gradients later
                  np=gpcf.fh.cfg(gpcf,[],[],[],0);
                else
                  DKffc = gpcf.fh.cfg(gpcf, x);
                  np=length(DKffc);
                end
              end

              gprior_cf = -gpcf.fh.lpg(gpcf);
              g1 = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);

              if isfield(gp.latent_opt, 'kron') && gp.latent_opt.kron==1

                for i2 = 1:length(DDa)
                  i1 = i1+1;

                  if ~isfield(gp,'meanf')
                    if i2==1
                      % derivative wrt magnitude hyperparameter
                      Vmsa = Vmslr'*a;
                      s1a = 0.5*(Vmsa'*(bsxfun(@times,Dmslr,Vmsa)) + a'*(Lms.*a));
                    elseif i2==2
                      % derivative wrt lengthscale hyperparameter
                      Vls1 = Vlslr1'*a;
                      Vls2 = Vlslr2'*a;
                      s1a = 0.5*( Vls1'*bsxfun(@times,Dlslr1,Vls1) + Vls2'*bsxfun(@times,Dlslr2,Vls2) + a'*(Lls.*a));
                    end
                  else
                    if i2==1
                      % derivative wrt magnitude hyperparameter
                      Vmsa = Vmslr'*(a-iKHb_m);
                      s1a = 0.5*(Vmsa'*(bsxfun(@times,Dmslr,Vmsa)) + (a-iKHb_m)'*(Lms.*(a-iKHb_m)));
                    elseif i2==2
                      % derivative wrt lengthscale hyperparameter
                      Vls1 = Vlslr1'*(a-iKHb_m);
                      Vls2 = Vlslr2'*(a-iKHb_m);
                      s1a = 0.5*( Vls1'*bsxfun(@times,Dlslr1,Vls1) + Vls2'*bsxfun(@times,Dlslr2,Vls2) + (a-iKHb_m)'*(Lls.*(a-iKHb_m)));
                    end
                  end

                  % DKg2=DKff{i2}*g2;
                  if i2==1
                    DKg2 = Vmslr*bsxfun(@times,Dmslr,Vmslr'*g2) + Lms.*g2;
                  elseif i2==2
                    DKg2 =  Vlslr1*bsxfun(@times,Dlslr1,Vlslr1'*g2) + Vlslr2*bsxfun(@times,Dlslr2,Vlslr2'*g2) + Lls.*g2;
                  end

                  WDKg2=ny*(g2.*DKg2-(g2*(g2'*DKg2)));
                  s1b = -0.5*(ny)*( ( - (DKg2-((WDKg2./Lbt)+(iStL'*(iStL*WDKg2)))))'*(g2) );

                  if i2==1 % magnitude hyperparameter

                    % low-rank
                    WDVa=ny*( bsxfun(@times,g2,Vmslr)-g2*(g2'*Vmslr) );
                    stmp=Vmslr-(bsxfun(@times,(1./Lbt),WDVa)+(iStL'*(iStL*WDVa)));
                    s1clr = 0.5*sum( (sum(bsxfun(@times,stmp,Dmslr').*Vmslr,2))' .*(-ny*g2)' );

                    % diagonal correction
                    s1cdtmp = Lms - ( ny*( (g2.*Lms)./Lbt  - (g2./Lbt).*(g2'.*Lms')' ) + ...
                      sum(iStL.* (ny*( bsxfun(@times,iStL,(g2.*Lms)') - (iStL*g2)*(g2'.*Lms') )) )' );
                    s1cd=0.5*sum( s1cdtmp' .*(-ny*g2)' );
                    s1c=s1clr+s1cd;
                    DKg = Vmslr*bsxfun(@times,Dmslr,Vmslr'*g1) + Lms.*g1;

                  elseif i2==2 % lengthscale hyperparameter

                    % low-rank 1
                    WDVa=ny*( bsxfun(@times,g2,Vlslr1)-g2*(g2'*Vlslr1) );
                    stmp=Vlslr1-(bsxfun(@times,(1./Lbt),WDVa)+(iStL'*(iStL*WDVa)));
                    s1clr1 = 0.5*sum( (sum(bsxfun(@times,stmp,Dlslr1').*Vlslr1,2))' .*(-ny*g2)' );

                    % low-rank 2
                    WDVb=ny*( bsxfun(@times,g2,Vlslr2)-g2*(g2'*Vlslr2) );
                    stmp=Vlslr2-(bsxfun(@times,(1./Lbt),WDVb)+(iStL'*(iStL*WDVb)));
                    s1clr2 = 0.5*sum( (sum(bsxfun(@times,stmp,Dlslr2').*Vlslr2,2))' .*(-ny*g2)' );

                    % diagonal correction
                    s1cdtmp = Lls - ( ny*( (g2.*Lls)./Lbt  - (g2./Lbt).*(g2'.*Lls')' ) + ...
                      sum(iStL.* (ny*( bsxfun(@times,iStL,(g2.*Lls)') - (iStL*g2)*(g2'.*Lls') )) )' );
                    s1cd=0.5*sum( s1cdtmp' .*(-ny*g2)' );

                    s1c=s1clr1+s1clr2+s1cd;

                    DKg = Vlslr1*bsxfun(@times,Dlslr1,Vlslr1'*g1) + Vlslr2*bsxfun(@times,Dlslr2,Vlslr2'*g1) + Lls.*g1;

                  end

                  s1=s1a+s1b+s1c;

                  %DKg=DKff{i2}*g1;
                  WDKg=ny*(g2.*DKg-(g2*(g2'*DKg)));
                  s3=DKg-((WDKg./Lbt)+(iStL'*(iStL*WDKg)));

                  gdata(i1) = -(s1 + s2'*s3);
                  gprior(i1) = gprior_cf(i2);
                end

              else

                for i2 = 1:np
                  i1 = i1+1;
                  if savememory
                    DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
                  else
                    DKff=DKffc{i2};
                  end
                  if ~isfield(gp,'meanf')
                    if strcmpi(gp.lik.type,'LGPC')
                      %s1 = 0.5 * a'*DKff{i2}*a - 0.5*((-iKW*(DKff{i2}*(sqrt(ny)*g2)))'*(sqrt(ny)*g2)) + 0.5*sum(sum(iKW'.*DKff{i2}).*(-ny*g2)');
                      s1 = 0.5 * a'*DKff*a - 0.5*sum(diag( R*(R'*(iKW*DKff))));
                    else
                      s1 = 0.5 * a'*DKff*a - 0.5*((-iKW*(DKff*(sqrt(ny)*g2)))'*(sqrt(ny)*g2)) + 0.5*sum(sum(iKW'.*DKff).*(-ny*g2)');
                    end
                  else
                    if strcmpi(gp.lik.type,'LGPC')
                      s1 = 0.5 * (a-iKHb_m)'*DKff*(a-iKHb_m) - 0.5*sum(diag( R*(R'*(iKW*DKff))));
                    else
                      s1 = 0.5 * (a-iKHb_m)'*DKff*(a-iKHb_m) - 0.5*((-iKW*(DKff*(sqrt(ny)*g2)))'*(sqrt(ny)*g2)) + 0.5*sum(sum(iKW'.*DKff).*(-ny*g2)');
                    end
                  end
                  %b = DKff{i2} * g1;
                  if issparse(K)
                    s3 = b - K*(sqrtW*ldlsolve(L,sqrtW*b));
                  else
                    s3=iKW*(DKff*g1);
                  end

                  gdata(i1) = -(s1 + s2'*s3);
                  gprior(i1) = gprior_cf(i2);
                end
              end

              if isfield(gp.latent_opt, 'kron') && gp.latent_opt.kron==1
                % Set the gradients of hyperparameter
                if length(gprior_cf) > length(DKa)
                  for i2=length(DKff)+1:length(gprior_cf)
                    i1 = i1+1;
                    gdata(i1) = 0;
                    gprior(i1) = gprior_cf(i2);
                  end
                end
              else
                % Set the gradients of hyperparameter
                if length(gprior_cf) > np
                  for i2=np+1:length(gprior_cf)
                    i1 = i1+1;
                    gdata(i1) = 0;
                    gprior(i1) = gprior_cf(i2);
                  end
                end
              end
            end
          end

          % =================================================================
          % Gradient with respect to likelihood function parameters
          if ~isempty(strfind(gp.infer_params, 'likelihood')) ...
              && ~isempty(gp.lik.fh.pak(gp.lik))

            gdata_lik = 0;
            lik = gp.lik;

            g_logPrior = -lik.fh.lpg(lik);
            if ~isempty(g_logPrior)

              DW_sigma = lik.fh.llg3(lik, y, f, 'latent2+param', z);
              DL_sigma = lik.fh.llg(lik, y, f, 'param', z);
              b = K * lik.fh.llg2(lik, y, f, 'latent+param', z);
              s3 = iKW*b;

              gdata_lik = - DL_sigma - 0.5.*sum(sum((A.*DW_sigma))) - s2'*s3;

              % set the gradients into vectors that will be returned
              gdata = [gdata gdata_lik];
              gprior = [gprior g_logPrior];
              i1 = length(g_logPrior);
              i2 = length(gdata_lik);
              if i1  > i2
                gdata = [gdata zeros(1,i1-i2)];
              end
            end
          end

          % g = gdata + gprior;

        case {'Softmax', 'Multinom'}

          K = zeros(n,n,nout);
          if multicf
            for i1=1:nout
              K(:,:,i1) = gp_trcov(gp, x, gp.comp_cf{i1});
            end
          else
            Ktmp=gp_trcov(gp, x);
            for i1=1:nout
              K(:,:,i1) = Ktmp;
            end
          end

          % softmax
          f2=reshape(f,n,nout);

          llg = gp.lik.fh.llg(gp.lik, y, f2, 'latent', z);
          [pi2_vec, pi2_mat] = gp.lik.fh.llg2(gp.lik, y, f2, 'latent', z);
          % W = -diag(pi2_vec) + pi2_mat*pi2_mat', where
          % W_ij = -d^2(log(p(y|f)))/(df_i)(df_j)
          R = repmat(1./pi2_vec,1,n).*pi2_mat;
          RE = zeros(n,n*nout);
          for i1=1:nout
            RE(:,(1:n)+(i1-1)*n) = R((1:n)+(i1-1)*n,:)'*E(:,:,i1);
          end

          inv_iWK=zeros(n,n,nout);

          % Matrices for computing the derivative of determinant term w.r.t. f
          A=zeros(nout, nout, n);
          Minv=M\(M'\eye(n));
          Minv=(Minv+Minv')./2;
          for cc1=1:nout
            EMinv=RE(:,(1:n)+(cc1-1)*n)'*Minv;
            KEMinv=K(:,:,cc1)*EMinv;
            for cc2=1:nout
              if cc2>=cc1
                if cc1==cc2
                  EMtmp = - EMinv*RE(:,(1:n)+(cc2-1)*n);
                  EMtmp = EMtmp + E(:,:,cc1);
                  inv_iWK(:,:,cc1) = EMtmp;
                  A(cc1,cc1,:) = diag(K(:,:,cc1))-sum((K(:,:,cc1)*EMtmp).*K(:,:,cc1),2);
                else
                  EMtmp = - KEMinv*RE(:,(1:n)+(cc2-1)*n);
                  A(cc1,cc2,:) = -sum(EMtmp.*K(:,:,cc2),2);
                  A(cc2,cc1,:) = -sum(EMtmp.*K(:,:,cc2),2);
                end
              end
            end
          end

          % Derivative of determinant term w.r.t. f
          s2=zeros(n*nout,1);
          dw_mat = gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
          for cc3=1:nout
            for ii1=1:n
              % s2(i)=-0.5*trace(inv(inv(K)+W)*dW/df_i)
              s2(ii1+(cc3-1)*n) = -0.5*trace(A(:,:,ii1)*dw_mat(:,:,cc3,ii1));
            end
          end

          % Loop over the covariance functions
          for i=1:ncf
            DKllg=zeros(size(a));
            EDKllg=zeros(size(a));
            DKffba=zeros(n*nout,1);

            % check in which components the covariance function is present
            do = false(nout,1);
            if multicf
              for z1=1:nout
                if any(gp.comp_cf{z1}==i)
                  do(z1) = true;
                end
              end
            else
              do = true(nout,1);
            end

            i1=0;
            if ~isempty(gprior)
              i1 = length(gprior);
            end

            % Gradients from covariance functions
            gpcf = gp.cf{i};
            % DKff{j} = dK(x,x)/dtheta_j
            if savememory
              % If savememory option is used, just return number of
              % hyperparameters and calculate gradients later
              np=gpcf.fh.cfg(gpcf,[],[],[],0);
            else
              DKff = gpcf.fh.cfg(gpcf, x);
              np=length(DKff);
            end
            gprior_cf = -gpcf.fh.lpg(gpcf);

            for i2 = 1:np
              i1 = i1+1;
              if savememory
                DKffb=gpcf.fh.cfg(gpcf,x,[],[],i2);
              else
                DKffb=DKff{i2};
              end

              % Derivative of explicit terms
              trace_sum_tmp=0;
              for z1=1:nout
                if do(z1)
                  DKffba((1:n)+(z1-1)*n)=DKffb*a((1:n)+(z1-1)*n);
                  trace_sum_tmp = trace_sum_tmp + sum(sum( inv_iWK(:,:,z1) .* DKffb ));
                end
              end
              % s1=0.5*f'*inv(K)*dK/dtheta_j*inv(K)*f - 0.5*trace(inv(inv(W)+K)*dK/dtheta_j)
              s1 = 0.5 * a'*DKffba - 0.5.*trace_sum_tmp;

              % Derivative of f w.r.t. theta
              for z1=1:nout
                if do(z1)
                  DKllg((1:n)+(z1-1)*n)=DKffb*llg((1:n)+(z1-1)*n);
                  EDKllg((1:n)+(z1-1)*n)=E(:,:,z1)*DKllg((1:n)+(z1-1)*n);
                end
              end
              s3 = EDKllg - RE'*(M\(M'\(RE*DKllg)));
              for z1=1:nout
                s3((1:n)+(z1-1)*n)=K(:,:,z1)*s3((1:n)+(z1-1)*n);
              end
              % s3=inv(I+KW)*dK/dtheta_j*d(log(p(y|f)))/df
              s3 = DKllg - s3;

              gdata(i1) = -(s1 + s2'*s3);
              gprior(i1) = gprior_cf(i2);

            end

            % Set the gradients of hyper-hyperparameter
            if length(gprior_cf) > np
              for i2=np+1:length(gprior_cf)
                i1 = i1+1;
                gdata(i1) = 0;
                gprior(i1) = gprior_cf(i2);
              end
            end
          end

          %         % =================================================================
          %         % Gradient with respect to likelihood function parameters
          %         if ~isempty(strfind(gp.infer_params, 'likelihood')) ...
          %             && ~isempty(gp.lik.fh.pak(gp.lik))
          %
          %             gdata_likelih = 0;
          %             lik = gp.lik;
          %
          %             g_logPrior = feval(lik.fh.gprior, lik);
          %             if ~isempty(g_logPrior)
          %
          %                 DW_sigma = feval(lik.fh.llg3, lik, y, f, 'latent2+hyper', z);
          %                 DL_sigma = feval(lik.fh.llg, lik, y, f, 'hyper', z);
          %                 b = K * feval(lik.fh.llg2, lik, y, f, 'latent+hyper', z);
          %                 s3 = b - K*(R*b);
          %                 nl= size(DW_sigma,2);
          %
          %                 gdata_lik = - DL_sigma - 0.5.*sum(repmat(C2,1,nl).*DW_sigma) - s2'*s3;
          %
          %                 % set the gradients into vectors that will be returned
          %                 gdata = [gdata gdata_lik];
          %                 gprior = [gprior g_logPrior];
          %                 i1 = length(g_logPrior);
          %                 i2 = length(gdata_lik);
          %                 if i1  > i2
          %                     gdata = [gdata zeros(1,i1-i2)];
          %                 end
          %             end
          %         end

        otherwise

          if isfield(gp.lik,'xtime')
            xtime=gp.lik.xtime;
            if isfield(gp.lik, 'stratificationVariables')
              ebc_ind=gp.lik.stratificationVariables;
              ux = unique(x(:,ebc_ind), 'rows');
              gp.lik.n_u = size(ux,1);
              for i1=1:size(ux,1)
                gp.lik.stratind{i1}=(x(:,ebc_ind)==ux(i1));
              end
              [xtime1, xtime2] = meshgrid(ux, xtime);
              xtime = [xtime2(:) xtime1(:)];
              if isfield(gp.lik, 'removeStratificationVariables') && gp.lik.removeStratificationVariables
                x(:,ebc_ind)=[];
              end
            end
            ntime = size(xtime,1);
            nl=[ntime n];

            % Second derivatives of log-likelihood
            [llg2vec, llg2mat] = gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
            % W = [diag(Wdiag(1:ntime)) Wmat; Wmat' diag(Wdiag(ntime+1:end))]
            Wdiag=-llg2vec; Wmat=-llg2mat;
          else
            nl=repmat(n,1,length(gp.comp_cf));

            % Second derivatives of log-likelihood
            Wvec=-gp.lik.fh.llg2(gp.lik, y, f, 'latent',z);
            % W = [diag(Wvec(1:n,1)) diag(Wvec(1:n,2)); diag(Wvec(n+1:end,1)) diag(Wvec(n+1:end,2))]
            Wdiag=[Wvec(1:nl(1),1); Wvec(nl(1)+(1:nl(2)),2)];
          end
          nlp=length(nl); % Number of latent processes

          % K is block-diagonal covariance matrix where blocks correspond to
          % latent processes
          K = zeros(sum(nl));
          if isfield(gp.lik,'xtime')
            K(1:ntime,1:ntime)=gp_trcov(gp, xtime, gp.comp_cf{1});
            K((1:n)+ntime,(1:n)+ntime) = gp_trcov(gp, x, gp.comp_cf{2});
          else
            for i1=1:nlp
              K((1:n)+(i1-1)*n,(1:n)+(i1-1)*n) = gp_trcov(gp, x, gp.comp_cf{i1});
            end
          end

          if isfield(gp,'meanf')
            [H,b_m,B_m]=mean_prep(gp,x,[]);
            Hb_m=H'*b_m;
            K=K+H'*B_m*H;
            iKHb_m=K\Hb_m;
          end

          KW=zeros(sum(nl));
          KW(1:nl(1),1:nl(1))=bsxfun(@times, K(1:nl(1),1:nl(1)), Wdiag(1:nl(1))');
          KW(nl(1)+(1:nl(2)),nl(1)+(1:nl(2)))=bsxfun(@times, K(nl(1)+(1:nl(2)),nl(1)+(1:nl(2))), Wdiag(nl(1)+(1:nl(2)))');
          if isfield(gp.lik,'xtime')
            KW(1:nl(1),nl(1)+(1:nl(2)))=K(1:nl(1),1:nl(1))*Wmat;
            KW(nl(1)+(1:nl(2)),1:nl(1))=K(nl(1)+(1:nl(2)),nl(1)+(1:nl(2)))*Wmat';
          else
            KW(1:nl(1),nl(1)+(1:nl(2)))=bsxfun(@times, K((1:nl(1)),(1:nl(1))), Wvec((nl(1)+1):2*n,1)');
            KW(nl(1)+(1:nl(2)),1:nl(1))=bsxfun(@times, K(nl(1)+(1:nl(2)),nl(1)+(1:nl(2))), Wvec(1:n,2)');
          end

          % B = (I + K*W)
          B=KW; B(1:(nl(1)+nl(2)+1):end)=B(1:(nl(1)+nl(2)+1):end)+1;

          iB=B\eye(sum(nl));

          % A = inv(I+K*W)*K
          A=iB*K;

          s2=zeros(sum(nl),1);

          if isfield(gp.lik,'xtime')
            A_diag=diag(A);
            A_mat=A(1:ntime,ntime+(1:n));
            for i1=1:sum(nl)
              % Third derivatives
              [dw_diag,dw_mat]=gp.lik.fh.llg3(gp.lik, y, f, 'latent', z, i1);
              % s2(i) = -0.5*trace(inv(inv(K)+W)*dW/df_i)
              s2(i1) = 0.5*(sum(A_diag.*dw_diag)+2*sum(sum(A_mat.*dw_mat)));
            end
          else
            % Third derivatives
            dw_mat = gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
            for i1=1:n
              % s2(i) = -0.5*trace(inv(inv(K)+W)*dW/df_i)
              s2(i1) = 0.5*trace(A(i1:n:(i1+n),i1:n:(i1+n))*dw_mat(:,:,1,i1));
              s2(i1+n) = 0.5*trace(A(i1:n:(i1+n),i1:n:(i1+n))*dw_mat(:,:,2,i1));
            end
          end

          % =================================================================
          % Gradient with respect to covariance function parameters
          if ~isempty(strfind(gp.infer_params, 'covariance'))
            % Evaluate the gradients from covariance functions
            for i=1:ncf

              i1=0;
              if ~isempty(gprior)
                i1 = length(gprior);
              end

              gpcf = gp.cf{i};

              % check in which components the covariance function is present
              do = false(nlp,1);
              for z1=1:nlp
                if any(gp.comp_cf{z1}==i)
                  do(z1) = true;
                end
              end

              if isfield(gp.lik,'xtime')
                if ~isempty(intersect(gp.comp_cf{1},i))
                  if savememory
                    % If savememory option is used, just get the number of
                    % hyperparametrs and calculate gradients later
                    np=gpcf.fh.cfg(gpcf,[],[],[],0);
                  else
                    DKffc = gpcf.fh.cfg(gpcf, xtime);
                    np=length(DKffc);
                  end
                else
                  if savememory
                    % If savememory option is used, just get the number of
                    % hyperparametrs and calculate gradients later
                    np=gpcf.fh.cfg(gpcf,[],[],[],0);
                  else
                    DKffc = gpcf.fh.cfg(gpcf, x);
                    np=length(DKffc);
                  end
                end
              else
                if savememory
                  % If savememory option is used, just get the number of
                  % hyperparametrs and calculate gradients later
                  np=gpcf.fh.cfg(gpcf,[],[],[],0);
                else
                  DKffc = gpcf.fh.cfg(gpcf, x);
                  np=length(DKffc);
                end
              end
              gprior_cf = -gpcf.fh.lpg(gpcf);
              g1 = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);

              WiB11=bsxfun(@times, Wdiag(1:nl(1)),iB(1:nl(1),1:nl(1)));
              WiB12=bsxfun(@times, Wdiag(1:nl(1)),iB(1:nl(1),nl(1)+(1:nl(2))));
              WiB22=bsxfun(@times, Wdiag(nl(1)+(1:nl(2))),iB(nl(1)+(1:nl(2)),nl(1)+(1:nl(2))));
              if isfield(gp.lik,'xtime')
                WiB11=WiB11 + Wmat*iB(nl(1)+(1:nl(2)),1:nl(1));
                WiB12=WiB12 + Wmat*iB(nl(1)+(1:nl(2)),nl(1)+(1:nl(2)));
                WiB22=WiB22 + Wmat'*iB(1:nl(1),nl(1)+(1:nl(2)));
              else
                WiB11=WiB11 + bsxfun(@times,Wvec(1:n,2),iB(nl(1)+(1:nl(2)),1:nl(1)));
                WiB12=WiB12 + bsxfun(@times,Wvec(1:n,2),iB(nl(1)+(1:nl(2)),nl(1)+(1:nl(2))));
                WiB22=WiB22 + bsxfun(@times,Wvec(nl(1)+(1:nl(2)),1),iB(1:nl(1),nl(1)+(1:nl(2))));
              end
              WiB=[WiB11 WiB12; WiB12' WiB22];
              % WiB=W*inv(I+KW)

              for i2 = 1:np
                i1 = i1+1;
                if ~isfield(gp,'meanf')
                  dKnl = zeros(sum(nl));
                  if isfield(gp.lik,'xtime')
                    if ~isempty(intersect(gp.comp_cf{1},i)) %do(indnl)
                      if savememory
                        DKff=gpcf.fh.cfg(gpcf,xtime,[],[],i2);
                      else
                        DKff=DKffc{i2};
                      end
                      dKnl(1:ntime,1:ntime) = DKff;
                      %end
                    else
                      if savememory
                        DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
                      else
                        DKff=DKffc{i2};
                      end
                      %if do(indnl)
                      dKnl(ntime+(1:n),ntime+(1:n)) = DKff;
                      %end
                    end
                  else
                    if savememory
                      DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
                    else
                      DKff=DKffc{i2};
                    end
                    for indnl=1:nlp
                      if do(indnl)
                        dKnl((1:n)+(indnl-1)*n,(1:n)+(indnl-1)*n) = DKff;
                      end
                    end
                  end
                  % s1 = 0.5*f'*inv(K)*dK/dtheta_j*inv(K)*f -
                  % 0.5*trace(inv(inv(W)+K)*dK/dtheta_j)
                  s1 = 0.5 * a'*dKnl*a - 0.5*sum(sum((dKnl.*WiB)));
                else
                  %s1 = 0.5 * (a-K\(H'*b_m))'*DKff{i2}*(a-K\(H'*b_m)) - 0.5*sum(sum(R.*DKff{i2}));
                end
                %b = DKff{i2} * g1;
                b = dKnl*g1;
                % s3 = inv(I+KW)*dK/dtheta_j*d(log(p(y|f)))/df
                s3=iB*b;

                gdata(i1) = -(s1 + s2'*s3);
                gprior(i1) = gprior_cf(i2);
              end
            end

            % Set the gradients of hyperparameter
            if length(gprior_cf) > np
              for i2=np+1:length(gprior_cf)
                i1 = i1+1;
                gdata(i1) = 0;
                gprior(i1) = gprior_cf(i2);
              end
            end
          end

          % =================================================================
          % Gradient with respect to likelihood function parameters
          if ~isempty(strfind(gp.infer_params, 'likelihood')) ...
              && ~isempty(gp.lik.fh.pak(gp.lik))

            gdata_lik = 0;
            lik = gp.lik;

            g_logPrior = -lik.fh.lpg(lik);
            if ~isempty(g_logPrior)

              DW_sigma = lik.fh.llg3(lik, y, f, 'latent2+param', z);
              DL_sigma = lik.fh.llg(lik, y, f, 'param', z);
              b = K * lik.fh.llg2(lik, y, f, 'latent+param', z);
              s3 = iB*b;

              gdata_lik = - DL_sigma - 0.5.*sum(sum((A.*DW_sigma))) - s2'*s3;

              % set the gradients into vectors that will be returned
              gdata = [gdata gdata_lik];
              gprior = [gprior g_logPrior];
              i1 = length(g_logPrior);
              i2 = length(gdata_lik);
              if i1  > i2
                gdata = [gdata zeros(1,i1-i2)];
              end
            end
          end

      end

      g = gdata + gprior;
    end

    case 'FIC'
      % ============================================================
      % FIC
      % ============================================================
      g_ind = zeros(1,numel(gp.X_u));
      gdata_ind = zeros(1,numel(gp.X_u));
      gprior_ind = zeros(1,numel(gp.X_u));

      u = gp.X_u;
      m = size(u,1);

      [e, edata, eprior, f, L, a, La1] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);

      K_fu = gp_cov(gp, x, u);         % f x u
      K_uu = gp_trcov(gp, u);          % u x u, noiseles covariance K_uu
      Luu = chol(K_uu);
      B=Luu'\(K_fu');       % u x f
      iKuuKuf = Luu\B;

      W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
      sqrtW = sqrt(W);

      % Components for trace( inv(inv(W) + K) * dK) )
      Lah = 1 + sqrtW.*La1.*sqrtW;
      sWKfu = (repmat(sqrtW,1,m).*K_fu);
      B3 = repmat(Lah,1,m).\sWKfu;
      A = K_uu + sWKfu'*B3; A=(A+A')/2;
      L2 = repmat(sqrtW,1,m).*(B3/chol(A));
      iLa2W = sqrtW.*(Lah.\sqrtW);

      LL = sum(L2.*L2,2);
      BB = sum(B.^2)';

      % Evaluate s2
      C1 = L2'*B'*B;
      C2 = L2'.*repmat(La1',m,1);

      s2t = La1 + BB;
      s2t = s2t - (La1.*iLa2W.*La1 - sum(C2.^2)' + sum(B'.*((B*(repmat(iLa2W,1,m).*B'))*B)',2)...
                   - sum(C1.^2)' + 2*La1.*iLa2W.*BB - 2*La1.*sum(L2.*C1',2));

      s2 = 0.5*s2t.*gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
      b3 = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);

      % =================================================================
      % Gradient with respect to covariance function parameters
      if ~isempty(strfind(gp.infer_params, 'covariance'))
        for i=1:ncf
          i1=0;
          if ~isempty(gprior)
            i1 = length(gprior);
          end

          % Get the gradients of the covariance matrices
          % and gprior from gpcf_* structures
          gpcf = gp.cf{i};
          if savememory
            % If savememory option is used, just get the number of
            % hyperparameters and calculate gradients later
            np=gpcf.fh.cfg(gpcf,[],[],[],0);
          else
            DKffc = gpcf.fh.cfg(gpcf, x, [], 1);
            DKuuc = gpcf.fh.cfg(gpcf, u);
            DKufc = gpcf.fh.cfg(gpcf, u, x);
            np=length(DKuuc);
          end
          gprior_cf = -gpcf.fh.lpg(gpcf);

          for i2 = 1:np
            i1 = i1+1;
            if savememory
              DKff=gpcf.fh.cfg(gpcf,x,[],1,i2);
              DKuu=gpcf.fh.cfg(gpcf,u,[],[],i2);
              DKuf=gpcf.fh.cfg(gpcf,u,x,[],i2);
            else
              DKff=DKffc{i2};
              DKuu=DKuuc{i2};
              DKuf=DKufc{i2};
            end
            % 0.5* a'*dK*a, where a = K\f
            KfuiKuuKuu = iKuuKuf'*DKuu;
            gdata(i1) = -0.5.*((2.*a'*DKuf'-(a'*KfuiKuuKuu))*(iKuuKuf*a) + (a'.*DKff')*a...
                               - (2.*a'.*sum(DKuf'.*iKuuKuf',2)'*a-a'.*sum(KfuiKuuKuu.*iKuuKuf',2)'*a) );

            % trace( inv(inv(W) + K) * dQ) )
            gdata(i1) = gdata(i1) - 0.5.*sum(sum(L2'.*(2.*L2'*DKuf'*iKuuKuf - L2'*KfuiKuuKuu*iKuuKuf)));
            gdata(i1) = gdata(i1) + 0.5.*sum(DKff.*iLa2W - LL.*DKff);
            gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*sum(DKuf'.*iKuuKuf',2)) - sum(LL.*sum(KfuiKuuKuu.*iKuuKuf',2)));

            b = (2*DKuf' - KfuiKuuKuu)*(iKuuKuf*b3) + DKff.*b3 - sum((2.*DKuf'- KfuiKuuKuu).*iKuuKuf',2).*b3;
            bb = sqrtW.*(Lah.\(sqrtW.*b)) - L2*(L2'*b);
            s3 = b - (La1.*bb + B'*(B*bb));
            gdata(i1) = gdata(i1) - s2'*s3;

            gprior(i1) = gprior_cf(i2);
          end

          % Set the gradients of hyperparameter
          if length(gprior_cf) > np
            for i2=np+1:length(gprior_cf)
              i1 = i1+1;
              gdata(i1) = 0;
              gprior(i1) = gprior_cf(i2);
            end
          end
        end

      end

      % =================================================================
      % Gradient with respect to likelihood function parameters

      if ~isempty(strfind(gp.infer_params, 'likelihood')) && ~isempty(gp.lik.fh.pak(gp.lik))
        gdata_lik = 0;
        lik = gp.lik;

        DW_sigma = lik.fh.llg3(lik, y, f, 'latent2+param', z);
        DL_sigma = lik.fh.llg(lik, y, f, 'param', z);
        DL_f_sigma = lik.fh.llg2(lik, y, f, 'latent+param', z);
%         b = La1.*DL_f_sigma + B'*(B*DL_f_sigma);
%         bb = (iLa2W.*b - L2*(L2'*b));
%         s3 = b - (La1.*bb + B'*(B*bb));
        b = repmat(La1,1,size(DL_f_sigma,2)).*DL_f_sigma + B'*(B*DL_f_sigma);
        bb = (repmat(iLa2W,1,size(DL_f_sigma,2)).*b - L2*(L2'*b));
        s3 = b - (repmat(La1,1,size(DL_f_sigma,2)).*bb + B'*(B*bb));

%         gdata_lik = - DL_sigma - 0.5.*sum(s2t.*DW_sigma) - s2'*s3;
        gdata_lik = - DL_sigma - 0.5.*sum(repmat(s2t,1,size(DL_f_sigma,2)).*DW_sigma) - s2'*s3;

        % evaluate prior contribution for the gradient
        if isfield(gp.lik, 'p')
          g_logPrior = -lik.fh.lpg(lik);
        else
          g_logPrior = zeros(size(gdata_lik));
        end
        % set the gradients into vectors that will be returned
        gdata = [gdata gdata_lik];
        gprior = [gprior g_logPrior];
        i1 = length(gdata);
      end

      % =================================================================
      % Gradient with respect to inducing inputs

      if ~isempty(strfind(gp.infer_params, 'inducing'))
        if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
          m = size(gp.X_u,2);
          st=0;
          if ~isempty(gprior)
            st = length(gprior);
          end

          gdata(st+1:st+length(gp.X_u(:))) = 0;
          i1 = st+1;
          for i = 1:size(gp.X_u,1)
            if iscell(gp.p.X_u) % Own prior for each inducing input
              pr = gp.p.X_u{i};
              gprior(i1:i1+m) = pr.fh.lpg(gp.X_u(i,:), pr);
            else % One prior for all inducing inputs
              gprior(i1:i1+m-1) = gp.p.X_u.fh.lpg(gp.X_u(i,:), gp.p.X_u);
            end
            i1 = i1 + m;
          end

          for i=1:ncf
            i1=st;

            gpcf = gp.cf{i};
            if savememory
              % If savememory option is used, just get the number of
              % covariates in X and calculate gradients later
              np=gpcf.fh.ginput(gpcf,u,[],0);
            else
              np=1;
              DKuu = gpcf.fh.ginput(gpcf, u);
              DKuf = gpcf.fh.ginput(gpcf, u, x);
            end

            for i3 = 1:np
              if savememory
                DKuu=gpcf.fh.ginput(gpcf,u,[],i3);
                DKuf=gpcf.fh.ginput(gpcf,u,x,i3);
              end
              for i2 = 1:length(DKuu)
                i1 = i1+1;

                % 0.5* a'*dK*a, where a = K\f
                KfuiKuuKuu = iKuuKuf'*DKuu{i2};
                gdata(i1) = gdata(i1) -0.5.*((2.*a'*DKuf{i2}'-(a'*KfuiKuuKuu))*(iKuuKuf*a) + ...
                                           - (2.*a'.*sum(DKuf{i2}'.*iKuuKuf',2)'*a-a'.*sum(KfuiKuuKuu.*iKuuKuf',2)'*a) );

                % trace( inv(inv(W) + K) * dQ) )
                gdata(i1) = gdata(i1) - 0.5.*(sum(sum(L2'.*(2.*L2'*DKuf{i2}'*iKuuKuf - L2'*KfuiKuuKuu*iKuuKuf))));
                gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*sum(DKuf{i2}'.*iKuuKuf',2)) - sum(LL.*sum(KfuiKuuKuu.*iKuuKuf',2)));

                % b2*dK*b3
                b = (2*DKuf{i2}'-KfuiKuuKuu)*(iKuuKuf*b3)  - sum((2.*DKuf{i2}'- KfuiKuuKuu).*iKuuKuf',2).*b3;
                bb = (iLa2W.*b - L2*(L2'*b));
                s3 = b - (La1.*bb + B'*(B*bb));
                gdata(i1) = gdata(i1) - s2'*s3;
              end
            end
          end
        end
      end

      g = gdata + gprior;

    case {'PIC' 'PIC_BLOCK'}
      % ============================================================
      % PIC
      % ============================================================
      g_ind = zeros(1,numel(gp.X_u));
      gdata_ind = zeros(1,numel(gp.X_u));
      gprior_ind = zeros(1,numel(gp.X_u));

      u = gp.X_u;
      m = size(u,1);
      ind = gp.tr_index;

      [e, edata, eprior, f, L, a, La1] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);

      K_fu = gp_cov(gp, x, u);         % f x u
      K_uu = gp_trcov(gp, u);          % u x u, noiseles covariance K_uu
      K_uu = (K_uu+K_uu')./2;          % ensure the symmetry of K_uu
      Luu = chol(K_uu);
      iKuuKuf = Luu\(Luu'\K_fu');
      B=Luu'\(K_fu');       % u x f

      W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
      sqrtW = sqrt(W);

      % Components for trace( inv(inv(W) + K) * dK) )
      B2 = (repmat(sqrtW,1,m).*K_fu);
      for i=1:length(ind)
        La2{i} = eye(size(La1{i})) + diag(sqrtW(ind{i}))*La1{i}*diag(sqrtW(ind{i}));
        LLa2{i} = chol(La2{i});
        B3(ind{i},:) = LLa2{i}\(LLa2{i}'\B2(ind{i},:));
      end
      A2 = K_uu + B2'*B3; A2=(A2+A2')/2;
      L2 = repmat(sqrtW,1,m).*B3/chol(A2);
      for i=1:length(ind)
        iLa2W{i} = diag(sqrtW(ind{i}))*(LLa2{i}\(LLa2{i}'\diag(sqrtW(ind{i}))));
      end

      LL = sum(L2.*L2,2);
      BB = sum(B.^2)';

      % Evaluate s2
      C1 = L2'*B'*B;
      s2t = BB;
      for i=1:length(ind)
        C2(:,ind{i}) = L2(ind{i},:)'*La1{i};
        s2t1(ind{i},:) = diag(La1{i}*iLa2W{i}*La1{i});
        s2t2(ind{i},:) = La1{i}*iLa2W{i}*B(:,ind{i})';
        s2t3(ind{i},:) = La1{i}*L2(ind{i},:);
        Bt(ind{i},:) = iLa2W{i}*B(:,ind{i})';
        s2t(ind{i}) = s2t(ind{i}) + diag(La1{i});
      end

      s2t = s2t - (s2t1 - sum(C2.^2)' + sum(B'.*((B*Bt)*B)',2)...
                   - sum(C1.^2)' + 2*sum(s2t2.*B',2) - 2*sum(s2t3.*C1',2));

      s2 = 0.5*s2t.*gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
      b3 = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);

      % =================================================================
      % Gradient with respect to covariance function parameters
      if ~isempty(strfind(gp.infer_params, 'covariance'))
        for i=1:ncf
          i1=0;
          if ~isempty(gprior)
            i1 = length(gprior);
          end

          % Get the gradients of the covariance matrices
          % and gprior from gpcf_* structures
          gpcf = gp.cf{i};

          if savememory
            % If savememory option is used, just get the number of
            % hyperparameters and calculate gradients later
            np=gpcf.fh.cfg(gpcf,[],[],[],0);
          else
            DKuuc = gpcf.fh.cfg(gpcf, u);
            DKufc = gpcf.fh.cfg(gpcf, u, x);
            for kk = 1:length(ind)
              DKffc{kk} = gpcf.fh.cfg(gpcf, x(ind{kk},:));
            end
            np=length(DKuuc);
          end
          gprior_cf = -gpcf.fh.lpg(gpcf);

          for i2 = 1:np
            i1 = i1+1;
            if savememory
              DKuu=gpcf.fh.cfg(gpcf,u,[],[],i2);
              DKuf=gpcf.fh.cfg(gpcf,u,x,[],i2);
            else
              DKuu=DKuuc{i2};
              DKuf=DKufc{i2};
            end
            KfuiKuuKuu = iKuuKuf'*DKuu;
            gdata(i1) = -0.5.*((2.*a'*DKuf'-(a'*KfuiKuuKuu))*(iKuuKuf*a) );
            gdata(i1) = gdata(i1) - 0.5.*(sum(sum(L2'.*(2.*L2'*DKuf'*iKuuKuf - L2'*KfuiKuuKuu*iKuuKuf))));

            b = (2*DKuf'-KfuiKuuKuu)*(iKuuKuf*b3);
            for kk=1:length(ind)
              if savememory
                DKff=gpcf.fh.cfg(gpcf, x(ind{kk},:),[],[],i2);
              else
                DKff=DKffc{kk}{i2};
              end
              gdata(i1) = gdata(i1) -0.5.*(a(ind{kk})'*DKff*a(ind{kk})...
                                           - (2.*a(ind{kk})'*DKuf(:,ind{kk})'*iKuuKuf(:,ind{kk})*a(ind{kk})...
                                              -a(ind{kk})'*KfuiKuuKuu(ind{kk},:)*iKuuKuf(:,ind{kk})*a(ind{kk})) );

              % trace( inv(inv(W) + K) * dQ) )
              gdata(i1) = gdata(i1) + 0.5.*(sum(sum(iLa2W{kk}.*DKff)) - sum(sum(L2(ind{kk},:)'.*(L2(ind{kk},:)'*DKff))));
              gdata(i1) = gdata(i1) + 0.5.*(2.*sum(sum(L2(ind{kk},:)'.*(L2(ind{kk},:)'*DKuf(:,ind{kk})'*iKuuKuf(:,ind{kk})))) - ...
                                            sum(sum(L2(ind{kk},:)'.*((L2(ind{kk},:)'*KfuiKuuKuu(ind{kk},:))*iKuuKuf(:,ind{kk})))));

              b(ind{kk}) = b(ind{kk}) + DKff*b3(ind{kk})...
                  - (2.*DKuf(:,ind{kk})'- KfuiKuuKuu(ind{kk},:))*iKuuKuf(:,ind{kk})*b3(ind{kk});
              bbt(ind{kk},:) = iLa2W{kk}*b(ind{kk});
            end

            % b2*dK*b3
            bb = (bbt - L2*(L2'*b));
            for kk=1:length(ind)
              s3t(ind{kk},:) = La1{kk}*bb(ind{kk});
            end
            s3 = b - (s3t + B'*(B*bb));
            gdata(i1) = gdata(i1) - s2'*s3;

            gprior(i1) = gprior_cf(i2);
          end

          % Set the gradients of hyperparameter
          if length(gprior_cf) > np
            for i2=np+1:length(gprior_cf)
              i1 = i1+1;
              gdata(i1) = 0;
              gprior(i1) = gprior_cf(i2);
            end
          end
        end

      end

      % =================================================================
      % Gradient with respect to likelihood function parameters

      if ~isempty(strfind(gp.infer_params, 'likelihood')) && ~isempty(gp.lik.fh.pak(gp.lik))
        gdata_lik = 0;
        lik = gp.lik;

        DW_sigma = lik.fh.llg3(lik, y, f, 'latent2+param', z);
        DL_sigma = lik.fh.llg(lik, y, f, 'param', z);
        DL_f_sigma = lik.fh.llg2(lik, y, f, 'latent+param', z);
        b = B'*(B*DL_f_sigma);
        for kk=1:length(ind)
          b(ind{kk}) = b(ind{kk}) + La1{kk}*DL_f_sigma(ind{kk});
          bbt(ind{kk},:) = iLa2W{kk}*b(ind{kk});
        end
        bb = (bbt - L2*(L2'*b));
        for kk=1:length(ind)
          s3t(ind{kk},:) = La1{kk}*bb(ind{kk});
        end
        s3 = b - (s3t + B'*(B*bb));

        gdata_lik = - DL_sigma - 0.5.*sum(s2t.*DW_sigma) - s2'*s3;

        % evaluate prior contribution for the gradient
        if isfield(gp.lik, 'p')
          g_logPrior = -lik.fh.lpg(lik);
        else
          g_logPrior = zeros(size(gdata_lik));
        end
        % set the gradients into vectors that will be returned
        gdata = [gdata gdata_lik];
        gprior = [gprior g_logPrior];
        i1 = length(gdata);
      end

      % =================================================================
      % Gradient with respect to inducing inputs

      if ~isempty(strfind(gp.infer_params, 'inducing'))
        if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
          m = size(gp.X_u,2);

          st=0;
          if ~isempty(gprior)
            st = length(gprior);
          end
          gdata(st+1:st+length(gp.X_u(:))) = 0;

          i1 = st+1;
          for i = 1:size(gp.X_u,1)
            if iscell(gp.p.X_u) % Own prior for each inducing input
              pr = gp.p.X_u{i};
              gprior(i1:i1+m) = pr.fh.lpg(gp.X_u(i,:), pr);
            else % One prior for all inducing inputs
              gprior(i1:i1+m-1) = gp.p.X_u.fh.lpg(gp.X_u(i,:), gp.p.X_u);
            end
            i1 = i1 + m;
          end

          % Loop over the  covariance functions
          for i=1:ncf
            i1=st;
            gpcf = gp.cf{i};
            if savememory
              % If savememory option is used, just get the number of
              % covariates in X and calculate gradients later
              np=gpcf.fh.ginput(gpcf,u,[],0);
            else
              DKuu = gpcf.fh.ginput(gpcf, u);
              DKuf = gpcf.fh.ginput(gpcf, u, x);
              np=1;
            end

            for i3 = 1:np
              if savememory
                DKuu=gpcf.fh.ginput(gpcf,u,[],i3);
                DKuf=gpcf.fh.ginput(gpcf,u,x,i3);
              end
              for i2 = 1:length(DKuu)
                i1 = i1+1;

                KfuiKuuKuu = iKuuKuf'*DKuu{i2};
                gdata(i1) = -0.5.*((2.*a'*DKuf{i2}'-(a'*KfuiKuuKuu))*(iKuuKuf*a) );
                gdata(i1) = gdata(i1) - 0.5.*(sum(sum(L2'.*(2.*L2'*DKuf{i2}'*iKuuKuf - L2'*KfuiKuuKuu*iKuuKuf))));

                b = (2*DKuf{i2}'-KfuiKuuKuu)*(iKuuKuf*b3);
                for kk=1:length(ind)
                  gdata(i1) = gdata(i1) -0.5.*(- (2.*a(ind{kk})'*DKuf{i2}(:,ind{kk})'*iKuuKuf(:,ind{kk})*a(ind{kk})...
                                                -a(ind{kk})'*KfuiKuuKuu(ind{kk},:)*iKuuKuf(:,ind{kk})*a(ind{kk})) );

                  % trace( inv(inv(W) + K) * dQ) )
                  gdata(i1) = gdata(i1) + 0.5.*(2.*sum(sum(L2(ind{kk},:)'.*(L2(ind{kk},:)'*DKuf{i2}(:,ind{kk})'*iKuuKuf(:,ind{kk})))) - ...
                                              sum(sum(L2(ind{kk},:)'.*((L2(ind{kk},:)'*KfuiKuuKuu(ind{kk},:))*iKuuKuf(:,ind{kk})))));

                  b(ind{kk}) = b(ind{kk}) + ...
                    - (2.*DKuf{i2}(:,ind{kk})'- KfuiKuuKuu(ind{kk},:))*iKuuKuf(:,ind{kk})*b3(ind{kk});
                  bbt(ind{kk},:) = iLa2W{kk}*b(ind{kk});
                end

                % b2*dK*b3
                bb = (bbt - L2*(L2'*b));
                for kk=1:length(ind)
                  s3t(ind{kk},:) = La1{kk}*bb(ind{kk});
                end
                s3 = b - (s3t + B'*(B*bb));
                gdata(i1) = gdata(i1) - s2'*s3;
              end
            end
          end
        end
      end

      g = gdata + gprior;

    case 'CS+FIC'
      % ============================================================
      % CS+FIC
      % ============================================================
      u = gp.X_u;
      m = size(u,1);

      [e, edata, eprior, f, L, a, La1] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);

      cf_orig = gp.cf;

      cf1 = {};
      cf2 = {};
      j = 1;
      k = 1;
      for i = 1:ncf
        if ~isfield(gp.cf{i},'cs')
          cf1{j} = gp.cf{i};
          j = j + 1;
        else
          cf2{k} = gp.cf{i};
          k = k + 1;
        end
      end
      gp.cf = cf1;

      % First evaluate needed covariance matrices
      % v defines that parameter is a vector
      [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % f x 1  vector
      K_fu = gp_cov(gp, x, u);         % f x u
      K_uu = gp_trcov(gp, u);    % u x u, noiseles covariance K_uu
      K_uu = (K_uu+K_uu')./2;     % ensure the symmetry of K_uu
      gp.cf = cf_orig;

      W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);

      % Find fill reducing permutation and permute all the
      % matrices
      p = analyze(La1);
      if ~isempty(z)
        z = z(p,:);
      end
      f = f(p);
      y = y(p);
      La1 = La1(p,p);
      K_fu = K_fu(p,:);
      L = L(p,:);
      x = x(p,:);
      W = W(p);
      a = a(p);

      Luu = chol(K_uu)';
      B=Luu\(K_fu');       % u x f
      iKuuKuf = Luu'\B;
      sW = sqrt(W);
      sqrtW = sparse(1:n,1:n,sW,n,n);
      Inn = sparse(1:n,1:n,1,n,n);

      % Components for trace( inv(inv(W) + K) * dK) )
      Lah = Inn + sqrtW*La1*sqrtW;
      LD2 = ldlchol(Lah);
      B2 = (repmat(sW,1,m).*K_fu);
      %B3 = Lah\B2;
      B3 = ldlsolve(LD2,B2);
      A2 = K_uu + B2'*B3; A2=(A2+A2')/2;
      L2 = repmat(sW,1,m).*B3/chol(A2);

      siLa2 = spinv(LD2,1);
      dsiLa2 = diag(siLa2);

      LL = sum(L2.*L2,2);
      BB = sum(B.^2)';

      % Evaluate s2
      C1 = L2'*B'*B;
      C2 = L2'*La1;
      C3 = repmat(sW,1,m).*ldlsolve(LD2,repmat(sW,1,m).*B');

      s2t = diag(La1) + BB;
      %diag(La1*sqrtW*ldlsolve(LD2,sqrtW*La1))
      s2t = s2t - (diag(La1) - sum(La1*sqrtW.*siLa2',2)./sW - sum(C2.^2)' + sum(B'.*(B*C3*B)',2)...
                   - sum(C1.^2)' + 2*sum((La1*C3).*B',2) - 2*sum(C2'.*C1',2));

      s2 = 0.5*s2t.*gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);

      b3 = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);

      % =================================================================
      % Gradient with respect to covariance function parameters
      if ~isempty(strfind(gp.infer_params, 'covariance'))
        for i=1:ncf
          i1=0;
          if ~isempty(gprior)
            i1 = length(gprior);
          end

          gpcf = gp.cf{i};

          % Evaluate the gradient for FIC covariance functions
          if ~isfield(gpcf,'cs')
            % Get the gradients of the covariance matrices
            % and gprior from gpcf_* structures
            if savememory
              % If savememory option is used, just get the number of
              % hyperparameters and calculate gradients later
              np=gpcf.fh.cfg(gpcf,[],[],[],0);
            else
              DKffc = gpcf.fh.cfg(gpcf, x, [], 1);
              DKuuc = gpcf.fh.cfg(gpcf, u);
              DKufc = gpcf.fh.cfg(gpcf, u, x);
              np=length(DKuuc);
            end
            gprior_cf = -gpcf.fh.lpg(gpcf);

            for i2 = 1:np
              i1 = i1+1;
              if savememory
                DKff=gpcf.fh.cfg(gpcf,x,[],1,i2);
                DKuu=gpcf.fh.cfg(gpcf,u,[],[],i2);
                DKuf=gpcf.fh.cfg(gpcf,u,x,[],i2);
              else
                DKff=DKffc{i2};
                DKuu=DKuuc{i2};
                DKuf=DKufc{i2};
              end
              % 0.5* a'*dK*a, where a = K\f
              KfuiKuuKuu = iKuuKuf'*DKuu;
              gdata(i1) = -0.5.*((2.*a'*DKuf'-(a'*KfuiKuuKuu))*(iKuuKuf*a) + (a'.*DKff')*a...
                                 - (2.*a'.*sum(DKuf'.*iKuuKuf',2)'*a-a'.*sum(KfuiKuuKuu.*iKuuKuf',2)'*a) );

              % trace( inv(inv(W) + K) * dQ) )
              gdata(i1) = gdata(i1) - 0.5.*(sum(sum(L2'.*(2.*L2'*DKuf'*iKuuKuf - L2'*KfuiKuuKuu*iKuuKuf))));
              gdata(i1) = gdata(i1) + 0.5.*(sum(DKff.*dsiLa2.*W - LL.*DKff));
              gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*sum(DKuf'.*iKuuKuf',2)) - sum(LL.*sum(KfuiKuuKuu.*iKuuKuf',2)));

              gdata(i1) = gdata(i1) + 0.5.*sum(sum(sqrtW*ldlsolve(LD2,repmat(sW,1,m).*(2.*DKuf' - KfuiKuuKuu)).*iKuuKuf',2));
              gdata(i1) = gdata(i1) - 0.5.*sum(sW.*dsiLa2.*sW.*sum((2.*DKuf' - KfuiKuuKuu).*iKuuKuf',2) );

              % b2*dK*b3
              b = (2*DKuf'-KfuiKuuKuu)*(iKuuKuf*b3) + DKff.*b3 - sum((2.*DKuf'- KfuiKuuKuu).*iKuuKuf',2).*b3;
              bb = (sW.*ldlsolve(LD2,sW.*b) - L2*(L2'*b));
              s3 = b - (La1*bb + B'*(B*bb));
              gdata(i1) = gdata(i1) - s2'*s3;

              gprior(i1) = gprior_cf(i2);
            end

            % Evaluate the gradient for compact support covariance functions
          else
            % Get the gradients of the covariance matrices
            % and gprior from gpcf_* structures
            if savememory
              % If savememory option is used, just get the number of
              % hyperparameters and calculate gradients later
              np=gpcf.fh.cfg(gpcf,[],[],[],0);
            else
              DKffc = gpcf.fh.cfg(gpcf, x);
              np=length(DKffc);
            end
            gprior_cf = -gpcf.fh.lpg(gpcf);

            for i2 = 1:np
              i1 = i1+1;
              if savememory
                DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
              else
                DKff=DKffc{i2};
              end

              % Evaluate the gradient with respect to magnSigma
              gdata(i1) = 0.5*(sum(sW.*sum(siLa2.*(sqrtW*DKff)',2)) - sum(sum(L2.*(L2'*DKff')')) - a'*DKff*a);
              b = DKff*b3;
              bb = (sW.*ldlsolve(LD2,sW.*b) - L2*(L2'*b));
              s3 = b - (La1*bb + B'*(B*bb));
              gdata(i1) = gdata(i1) - s2'*s3;
              gprior(i1) = gprior_cf(i2);

            end
          end

          % Set the gradients of hyperparameter
          if length(gprior_cf) > np
            for i2=np+1:length(gprior_cf)
              i1 = i1+1;
              gdata(i1) = 0;
              gprior(i1) = gprior_cf(i2);
            end
          end
        end

      end

      % =================================================================
      % Gradient with respect to likelihood function parameters

      if ~isempty(strfind(gp.infer_params, 'likelihood')) && ~isempty(gp.lik.fh.pak(gp.lik))
        gdata_lik = 0;
        lik = gp.lik;

        DW_sigma = lik.fh.llg3(lik, y, f, 'latent2+param', z);
        DL_sigma = lik.fh.llg(lik, y, f, 'param', z);
        DL_f_sigma = lik.fh.llg2(lik, y, f, 'latent+param', z);
        b = La1*DL_f_sigma + B'*(B*DL_f_sigma);
        bb = (sW.*ldlsolve(LD2,sW.*b) - L2*(L2'*b));
        s3 = b - (La1*bb + B'*(B*bb));

        gdata_lik = - DL_sigma - 0.5.*sum(s2t.*DW_sigma) - s2'*s3;

        % evaluate prior contribution for the gradient
        if isfield(gp.lik, 'p')
          g_logPrior = -lik.fh.lpg(lik);
        else
          g_logPrior = zeros(size(gdata_lik));
        end
        % set the gradients into vectors that will be returned
        gdata = [gdata gdata_lik];
        gprior = [gprior g_logPrior];
        i1 = length(gdata);
      end

      % =================================================================
      % Gradient with respect to inducing inputs

      if ~isempty(strfind(gp.infer_params, 'inducing'))
        if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
          m = size(gp.X_u,2);
          st=0;
          if ~isempty(gprior)
            st = length(gprior);
          end

          gdata(st+1:st+length(gp.X_u(:))) = 0;
          i1 = st+1;
          for i = 1:size(gp.X_u,1)
            if iscell(gp.p.X_u) % Own prior for each inducing input
              pr = gp.p.X_u{i};
              gprior(i1:i1+m) = pr.fh.lpg(gp.X_u(i,:), pr);
            else % One prior for all inducing inputs
              gprior(i1:i1+m-1) = gp.p.X_u.fh.lpg(gp.X_u(i,:), gp.p.X_u);
            end
            i1 = i1 + m;
          end

          for i=1:ncf
            i1=st;
            gpcf = gp.cf{i};
            if ~isfield(gpcf,'cs')
              if savememory
                % If savememory option is used, just get the number of
                % covariates in X and calculate gradients later
                np=gpcf.fh.ginput(gpcf,u,[],0);
              else
                DKuu = gpcf.fh.ginput(gpcf, u);
                DKuf = gpcf.fh.ginput(gpcf, u, x);
                np=1;
              end

              for i3 = 1:np
                if savememory
                  DKuu=gpcf.fh.ginput(gpcf,u,[],i3);
                  DKuf=gpcf.fh.ginput(gpcf,u,x,i3);
                end
                for i2 = 1:length(DKuu)
                  i1=i1+1;

                  % 0.5* a'*dK*a, where a = K\f
                  KfuiKuuKuu = iKuuKuf'*DKuu{i2};
                  gdata(i1) = -0.5.*((2.*a'*DKuf{i2}'-(a'*KfuiKuuKuu))*(iKuuKuf*a) ...
                                   - (2.*a'.*sum(DKuf{i2}'.*iKuuKuf',2)'*a-a'.*sum(KfuiKuuKuu.*iKuuKuf',2)'*a) );

                  % trace( inv(inv(W) + K) * dQ) )
                  gdata(i1) = gdata(i1) - 0.5.*(sum(sum(L2'.*(2.*L2'*DKuf{i2}'*iKuuKuf - L2'*KfuiKuuKuu*iKuuKuf))));
                  gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*sum(DKuf{i2}'.*iKuuKuf',2)) - sum(LL.*sum(KfuiKuuKuu.*iKuuKuf',2)));

                  gdata(i1) = gdata(i1) + 0.5.*sum(sum(sqrtW*ldlsolve(LD2,repmat(sW,1,size(u,1)).*(2.*DKuf{i2}' - KfuiKuuKuu)).*iKuuKuf',2));
                  gdata(i1) = gdata(i1) - 0.5.*( sum(sW.*dsiLa2.*sW.*sum((2.*DKuf{i2}' - KfuiKuuKuu).*iKuuKuf',2)) );

                  % b2*dK*b3
                  b = (2*DKuf{i2}'-KfuiKuuKuu)*(iKuuKuf*b3) - sum((2.*DKuf{i2}'- KfuiKuuKuu).*iKuuKuf',2).*b3;
                  bb = (sW.*ldlsolve(LD2,sW.*b) - L2*(L2'*b));
                  s3 = b - (La1*bb + B'*(B*bb));
                  gdata(i1) = gdata(i1) - s2'*s3;
                end
              end
            end
          end
        end
      end

      g = gdata + gprior;

    case {'DTC', 'VAR', 'SOR'}
      % ============================================================
      % DTC, VAR, SOR
      % ============================================================
      g_ind = zeros(1,numel(gp.X_u));
      gdata_ind = zeros(1,numel(gp.X_u));
      gprior_ind = zeros(1,numel(gp.X_u));

      u = gp.X_u;
      m = size(u,1);

      [e, edata, eprior, f, L, a, La2] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);

      K_fu = gp_cov(gp, x, u);         % f x u
      K_uu = gp_trcov(gp, u);          % u x u, noiseles covariance K_uu
      Luu = chol(K_uu);
      B=Luu'\(K_fu');       % u x f
      iKuuKuf = Luu\B;

      W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
      sqrtW = sqrt(W);

      % Components for trace( inv(inv(W) + K) * dK) )
      sWKfu = bsxfun(@times, sqrtW, K_fu);
      % L = chol(K_uu + sWKfu'*sWKfu)
      L2 = repmat(sqrtW,1,m).*(sWKfu/L);

      LL = sum(L2.*L2,2);
      BB = sum(B.^2)';

      % Evaluate s2
      C1 = L2'*B'*B;
      C2 = zeros(size(L2'));

      s2t = BB;
      s2t = s2t - (sum(C2.^2)' + sum(B'.*((B*(repmat(W,1,m).*B'))*B)',2)...
        - sum(C1.^2)');

      g3 = gp.lik.fh.llg3(gp.lik, y, f, 'latent', z);
      s2 = 0.5*s2t.*g3;
      b3 = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);

      % =================================================================
      % Gradient with respect to covariance function parameters
      if ~isempty(strfind(gp.infer_params, 'covariance'))
        for i=1:ncf
          i1=0;
          if ~isempty(gprior)
            i1 = length(gprior);
          end

          % Get the gradients of the covariance matrices
          % and gprior from gpcf_* structures
          gpcf = gp.cf{i};
          if savememory
            % If savememory option is used, just get the number of
            % hyperparameters and calculate gradients later
            np=gpcf.fh.cfg(gpcf,[],[],[],0);
          else
            %             DKffc = gpcf.fh.cfg(gpcf, x, [], 1);
            DKuuc = gpcf.fh.cfg(gpcf, u);
            DKufc = gpcf.fh.cfg(gpcf, u, x);
            if isequal(gp.type, 'VAR')
              DKffc=gpcf.fh.cfg(gpcf,x,[],1);
            end
            np=length(DKuuc);
          end
          gprior_cf = -gpcf.fh.lpg(gpcf);

          for i2 = 1:np
            i1 = i1+1;
            if savememory
              DKuu=gpcf.fh.cfg(gpcf,u,[],[],i2);
              DKuf=gpcf.fh.cfg(gpcf,u,x,[],i2);
            else
              DKuu=DKuuc{i2};
              DKuf=DKufc{i2};
            end
            % 0.5* a'*dK*a, where a = K\f
            KfuiKuuDKuu = iKuuKuf'*DKuu;
            gdata(i1) = -0.5.*((2.*a'*DKuf'-(a'*KfuiKuuDKuu))*(iKuuKuf*a));

            % trace( inv(inv(W) + K) * dQ) )
            gdata(i1) = gdata(i1) - 0.5.*sum(sum(L2'.*(2.*L2'*DKuf'*iKuuKuf - L2'*KfuiKuuDKuu*iKuuKuf)));
            gdata(i1) = gdata(i1) + 0.5.*(2.*sum(W.*sum(DKuf'.*iKuuKuf',2)) - sum(W.*sum(KfuiKuuDKuu.*iKuuKuf',2)));

            b = (2*DKuf' - KfuiKuuDKuu)*(iKuuKuf*b3);
            bb = sqrtW.*sqrtW.*b - L2*(L2'*b);
            s3 = b - (B'*(B*bb));
            gdata(i1) = gdata(i1) - s2'*s3;

            if isequal(gp.type, 'VAR')
              % Derivative of tr(diag(K-Q).*W) (Titsias, 2009) wrt th
              % Here we have to also split in explicit and implicit
              % derivatives
              if savememory
                DKff=gpcf.fh.cfg(gpcf,x,[],1,i2);
              else
                DKff=DKffc{i2};
              end
              gdata(i1) = gdata(i1) + 0.5*sum(DKff.*W) - 0.5.*(2.*sum(W.*sum(DKuf'.*iKuuKuf',2)) - sum(W.*sum(KfuiKuuDKuu.*iKuuKuf',2)));
              % La2 = diag(K - Q);
              gdata(i1) = gdata(i1) - 0.5*(La2.*g3)'*s3;
            end

            gprior(i1) = gprior_cf(i2);
          end

          % Set the gradients of hyperparameter
          if length(gprior_cf) > np
            for i2=np+1:length(gprior_cf)
              i1 = i1+1;
              gdata(i1) = 0;
              gprior(i1) = gprior_cf(i2);
            end
          end
        end

      end

      % =================================================================
      % Gradient with respect to likelihood function parameters

      if ~isempty(strfind(gp.infer_params, 'likelihood')) && ~isempty(gp.lik.fh.pak(gp.lik))
        gdata_lik = 0;
        lik = gp.lik;

        DW_sigma = lik.fh.llg3(lik, y, f, 'latent2+param', z);
        DL_sigma = lik.fh.llg(lik, y, f, 'param', z);
        DL_f_sigma = lik.fh.llg2(lik, y, f, 'latent+param', z);
        %         b = La1.*DL_f_sigma + B'*(B*DL_f_sigma);
        %         bb = (iLa2W.*b - L2*(L2'*b));
        %         s3 = b - (La1.*bb + B'*(B*bb));

        % b = Kfu*inv(Kuu)*Kfu'*dW/dth
        b = B'*(B*DL_f_sigma);
        % bb = Kfu*inv(Kuu+Kfu'*W*Kfu)*Kfu'*W
        bb = repmat(1./W,1,size(DL_f_sigma,2)).*(L2*(L2'*b));
        % s3 = df/dth
        s3 = b - bb;
        %bb = (repmat(W,1,size(DL_f_sigma,2)).*b - L2*(L2'*b));
        %s3 = b + B'*(B*bb);

        %         gdata_lik = - DL_sigma - 0.5.*sum(s2t.*DW_sigma) - s2'*s3;
        gdata_lik = - DL_sigma - 0.5.*sum(repmat(s2t,1,size(DL_f_sigma,2)).*DW_sigma) - s2'*s3;

        if isequal(gp.type, 'VAR')
          % Derivative of the trace term tr(diag(K-Q).*W)
          % Explicit dW/dth
          % La2 = diag(K-Q)
          gdata_lik = gdata_lik - 0.5*sum(La2.*DW_sigma);
          % Implicit dW/df*df/dth
          gdata_lik = gdata_lik + 0.5.*(La2.*g3)'*s3;
        end

        % evaluate prior contribution for the gradient
        if isfield(gp.lik, 'p')
          g_logPrior = -lik.fh.lpg(lik);
        else
          g_logPrior = zeros(size(gdata_lik));
        end
        % set the gradients into vectors that will be returned
        gdata = [gdata gdata_lik];
        gprior = [gprior g_logPrior];
        i1 = length(gdata);
      end

      % =================================================================
      % Gradient with respect to inducing inputs

      if ~isempty(strfind(gp.infer_params, 'inducing'))
        if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
          m = size(gp.X_u,2);
          st=0;
          if ~isempty(gprior)
            st = length(gprior);
          end

          gdata(st+1:st+length(gp.X_u(:))) = 0;
          i1 = st+1;
          for i = 1:size(gp.X_u,1)
            if iscell(gp.p.X_u) % Own prior for each inducing input
              pr = gp.p.X_u{i};
              gprior(i1:i1+m) = pr.fh.lpg(gp.X_u(i,:), pr);
            else % One prior for all inducing inputs
              gprior(i1:i1+m-1) = gp.p.X_u.fh.lpg(gp.X_u(i,:), gp.p.X_u);
            end
            i1 = i1 + m;
          end

          for i=1:ncf
            i1=st;

            gpcf = gp.cf{i};
            if savememory
              % If savememory option is used, just get the number of
              % covariates in X and calculate gradients later
              np=gpcf.fh.ginput(gpcf,u,[],0);
            else
              np=1;
              DKuu = gpcf.fh.ginput(gpcf, u);
              DKuf = gpcf.fh.ginput(gpcf, u, x);
            end

            for i3 = 1:np
              if savememory
                DKuu=gpcf.fh.ginput(gpcf,u,[],i3);
                DKuf=gpcf.fh.ginput(gpcf,u,x,i3);
              end

              for i2 = 1:length(DKuu)
                i1 = i1+1;

                % 0.5* a'*dK*a, where a = K\f
                KfuiKuuDKuu = iKuuKuf'*DKuu{i2};
                gdata(i1) = -0.5.*((2.*a'*DKuf{i2}'-(a'*KfuiKuuDKuu))*(iKuuKuf*a));

                % trace( inv(inv(W) + K) * dQ) )
                gdata(i1) = gdata(i1) - 0.5.*sum(sum(L2'.*(2.*L2'*DKuf{i2}'*iKuuKuf - L2'*KfuiKuuDKuu*iKuuKuf)));
                tt=0.5.*(2.*sum(W.*sum(DKuf{i2}'.*iKuuKuf',2)) - sum(W.*sum(KfuiKuuDKuu.*iKuuKuf',2)));
                gdata(i1) = gdata(i1) + tt;

                b = (2*DKuf{i2}' - KfuiKuuDKuu)*(iKuuKuf*b3);
                bb = sqrtW.*sqrtW.*b - L2*(L2'*b);
                s3 = b - (B'*(B*bb));
                gdata(i1) = gdata(i1) - s2'*s3;

                if isequal(gp.type, 'VAR')
                  % Derivative of 0.5.*tr(diag(K-Q).*W) (Titsias, 2009) wrt X_u
                  % Here we have to also split in explicit and implicit
                  % derivatives
                  gdata(i1) = gdata(i1) + 0 - tt;
                  % La2 = diag(K - Q);
                  gdata(i1) = gdata(i1) - 0.5*(La2.*g3)'*s3;
                end
              end
            end
          end
        end
      end

      g = gdata + gprior;

  end

  assert(isreal(gdata))
  assert(isreal(gprior))
end


text:
function [row, col, mu] = isolated_peaks_new(S1, ops)
% takes a matrix of timepoints by channels S1
% outputs threshold crossings that are relatively isolated from other peaks
% outputs row, column and magnitude of the threshold crossing
loc_range = getOr(ops, 'loc_range', [5 4]);
long_range = getOr(ops, 'long_range', [30 6]);
Th = ops.spkTh;
nt0 = ops.nt0;

% loc_range = [3  1];
% long_range = [30  6];

% finding the local minimum in a sliding window within plus/minus loc_range extent
% across time and across channels
smin = my_min(S1, loc_range, [1 2]);
peaks = single(S1<smin+1e-3 & S1<Th); % the peaks are samples that achieve this local minimum, AND have negativities less than a preset threshold

% only take local peaks that are isolated from other local peaks
sum_peaks = my_sum(peaks, long_range, [1 2]); % if there is another local peak close by, this sum will be at least 2
peaks = peaks .* (sum_peaks<1.2) .* S1; % set to 0 peaks that are not isolated, and multiply with the voltage values

% exclude temporal buffers
peaks([1:nt0 end-nt0:end], :) = 0;

[row, col, mu] = find(peaks); % find the non-zero peaks, and take their amplitudes

mu = - mu; % invert the sign of the amplitudes

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [g1, g2, g3] = sdlfmXsdrbfKernGradientBlockIGJ(lfmKern, rbfKern, ...
    t1, t2, i, j, generalConst, generalConstGrad, covGrad, typeMeanGrad)

% SDLFMXSDRBFKERNGRADIENTBLOCKIGJ
% FORMAT
% DESC computes the gradients of the parameters for output system and
% latent force system when i is greater than j.
% ARG lfmKern : structure containing parameters for the output system
% ARG rbfKern : structure containing parameters for the latent force system
% ARG t1 : times at which the output system is evaluated
% ARG t2 : times at which the LF system is evaluated
% ARG i : interval in which the output is evaluated
% ARG j : interval in which the latent force is being evaluated
% ARG i : interval to be evaluated for system 1
% ARG j : interval to be evaluated for system 2
% ARG generalConstant : constants evaluated with sdlfmKernComputeConstant.m
% ARG generalConstGrad : derivatives of the constants computed with
% sdlfmKernGradientConstant.m
% ARG covGrad : partial derivatives of the objective function wrt portion
% of the corresponding kernel matrix
% ARG typeMeanGrad : specify the mean functions used to compute this part
% of the kernel
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin < 10
    typeMeanGrad = {'lfm', 'lfmv'};
end

fhandleMeanParam = str2func(['sd' typeMeanGrad{1} 'MeanCompute']);
fhandleMeanGradParam = str2func(['sd' typeMeanGrad{1} 'MeanGradient']);
fhandleMeanGradSwitching = str2func(['sd' typeMeanGrad{2} 'MeanCompute']);

c1 = fhandleMeanParam(lfmKern, t1, 'Pos');
e1 = fhandleMeanParam(lfmKern, t1, 'Vel');

if isempty(generalConst{i,j})
    coeffPosRbf = c1;
    coeffVelRbf = e1;
else
    coeffPosRbf = generalConst{i,j}(1,1)*c1 + generalConst{i,j}(2,1)*e1;
    coeffVelRbf = generalConst{i,j}(1,2)*c1 + generalConst{i,j}(2,2)*e1;
end

PosRbf = lfmXrbfKernCompute(lfmKern, rbfKern, rbfKern.limit, t2);
[g1Kern, g2Kern] = lfmXrbfKernGradient(lfmKern, rbfKern, ...
        rbfKern.limit, t2, covGrad.', coeffPosRbf);
VelRbf = lfmvXrbfKernCompute(lfmKern, rbfKern, rbfKern.limit, t2);
[g1KLocal, g2KLocal] = lfmvXrbfKernGradient(lfmKern, rbfKern, ...
        rbfKern.limit, t2, covGrad.', coeffVelRbf);
g1Kern = g1Kern + g1KLocal;
g2Kern = g2Kern + g2KLocal;

[gcAlpha, gcOmega] = fhandleMeanGradParam(lfmKern, t1, 'Pos');
[geAlpha, geOmega] = fhandleMeanGradParam(lfmKern, t1, 'Vel');
[gradAlpha, gradOmega] = getLocalGradAlphaOmega(lfmKern);

g1Pos = fhandleMeanGradSwitching(lfmKern, t1, 'Pos');
h1Vel = fhandleMeanGradSwitching(lfmKern, t1, 'Vel');

% This means that are only two switching points involved, t1
% and t0 appear in k, like k_{ff}(t1-t0,t'-t0) and
% k_{mf}(t1-t0,t'-t0). The other points appear in c1(t - t1)
% and e1(t - t1)
gsp1Kern = 0;
gsp2Kern = 0;

temp = lfmvXrbfKernCompute(lfmKern, rbfKern, rbfKern.limit, t2);
covPartial = sum(sum((coeffPosRbf*temp).*covGrad));
gsp2Kern = gsp2Kern - covPartial;
gsp1Kern = gsp1Kern + covPartial;
temp = lfmXrbfvKernCompute(lfmKern, rbfKern, rbfKern.limit, t2);
covPartial = sum(sum((coeffPosRbf*temp).*covGrad));
gsp2Kern = gsp2Kern - covPartial;
temp = lfmaXrbfKernCompute(lfmKern, rbfKern, rbfKern.limit, t2);
covPartial = sum(sum((coeffVelRbf*temp).*covGrad));
gsp2Kern = gsp2Kern - covPartial;
gsp1Kern = gsp1Kern + covPartial;
temp = lfmvXrbfvKernCompute(lfmKern, rbfKern, rbfKern.limit, t2);
gsp2Kern = gsp2Kern - sum(sum((coeffVelRbf*temp).*covGrad));

% Organize gradients accordingly
if isempty(generalConst{i,j})
    matGradAlpha = gcAlpha*PosRbf + geAlpha*VelRbf;
    matGradOmega = gcOmega*PosRbf + geOmega*VelRbf;
    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    matGradSp1 = g1Pos*PosRbf + h1Vel*VelRbf;
    gsp1 = - sum(sum(matGradSp1.*covGrad));
    g3(i) = sum(gsp1Kern) + gsp1; % switching point 1
    g3(j) = sum(gsp2Kern);        % switching point 2
else
    constGradAlpha = generalConstGrad{1}{i,j};
    constGradOmega = generalConstGrad{2}{i,j};
    constVal = generalConst{i,j};
    % Derivative wrt parameters
    matGradAlpha = (constVal(1,1)*gcAlpha + constGradAlpha(1,1)*c1 ...
        + constVal(2,1)*geAlpha + constGradAlpha(2,1)*e1)*PosRbf ...
        + (constVal(1,2)*gcAlpha + constGradAlpha(1,2)*c1 ...
        + constVal(2,2)*geAlpha + constGradAlpha(2,2)*e1)*VelRbf;
    matGradOmega = (constVal(1,1)*gcOmega + constGradOmega(1,1)*c1 ...
        + constVal(2,1)*geOmega + constGradOmega(2,1)*e1)*PosRbf ...
        + (constVal(1,2)*gcOmega + constGradOmega(1,2)*c1 ...
        + constVal(2,2)*geOmega + constGradOmega(2,2)*e1)*VelRbf;
    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    % Derivative wrt switching points
    % Firts, notice that gsp1Kern doesn't correspond to the switching point
    % of the interval i (or say 1). It's just the derivative of the inner
    % switching point that leads to the actual switching point for the
    % interval 1. So we rename it first.
    gspInt = sum(gsp1Kern);
    % Compute the derivative of the swicthing point i, not appearing in the
    % constant
    matGradSp1 = (constVal(1,1)*g1Pos + constVal(2,1)*h1Vel)*PosRbf + ...
        (constVal(1,2)*g1Pos + constVal(2,2)*h1Vel)*VelRbf;
    gsp1 = - sum(sum(matGradSp1.*covGrad));
    % Compute the derivatives for all the switching points appearing in the
    % constant. The order is descending, i.e., t_3_, t_2, t_1
    constGradSPoint = generalConstGrad{3}{i,j};
    numberSP = size(constGradSPoint,2);
    gspInBetween = zeros(1, numberSP);
    for k=1:numberSP
        temp = (constGradSPoint(1,k)*c1 + constGradSPoint(3,k)*e1)*PosRbf + ...
            (constGradSPoint(2,k)*c1 + constGradSPoint(4,k)*e1)*VelRbf;
        gspInBetween(k) = sum(sum(temp.*covGrad));
    end
    % Assign derivatives wrt all other switching points, with a correction
    % for the derivative of the innermost switching point in the constant
    gspInBetween(end) =  gspInBetween(end) + gspInt;
    g3(i) = gsp1 + gspInBetween(1);
    g3(j) = sum(gsp2Kern);
    g3(j+1:i-1) = fliplr(gspInBetween(2:end));
end
g1 = zeros(1,5);
g2 = zeros(1,2);
% Assign derivatives wrt first system
g1(1:3) = g1Kern(1:3) + gCoeff;          % mass, spring, damper
g1(4) = g1Kern(4);                            % inverse width
g1(5) = g1Kern(5);                            % sensitivity
% Assign derivatives wrt second system
g2(1) = g2Kern(1);                      % variance
g2(2) = g2Kern(2);                   % inverse widths

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ y ] = tapas_rdcm_subsample(Y, r_dt)
% [ y ] = tapas_rdcm_subsample(Y, r_dt)
%
% Subsamples signal Y (in frequency domain) with a rate r_dt
%
%   Input:
%       Y       - original signal
%       r_dt    - time step (delta_t)
%
%   Output:
%       y       - subsampled signal
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2022 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% dimensionality of data
[N, dim] = size(Y);
p = factor(r_dt);
y = zeros(N/r_dt,dim);

for i = 1:dim
    y_tmp = Y(:,i);
    for j = 1:length(p)
        y_tmp = decimate(y_tmp, p(j),12); % decimating (subsampling) u
    end
    y(:,i) = y_tmp;
end

% subsample the data
y = Y(1:r_dt:end,:)/r_dt;

end

author: translationalneuromodeling

repo: tapas


text:
function [model,labels] = fgplvmReadFromFile(fileName)

% FGPLVMREADFROMFILE Load a file produced by the C++ implementation.
% FORMAT
% DESC loads a GP-LVM model from a file produced by the C++
% implementation of the GP-LVM.
% ARG fileName : the file name written by the C++ software.
% RETURN model : a MATLAB GP-LVM model structure containing the
% model from the file.
% RETURN labels : any available labels from the GP-LVM file.
%
% SEEALSO : fgplvmReadFromFID, fgplvmCreate
%
% COPYRIGHT : Neil D. Lawrence, 2005

% FGPLVM

FID = fopen(fileName);
if FID==-1
  error(['Cannot find file ' fileName])
end
[model, labels] = fgplvmReadFromFID(FID);
fclose(FID);


text:
function testBranin(varargin)

run_high_dim = 0;
embedding = 1;
embedding_consecutive = 0;
test_percentage = 0;

total_iter = 500;
high_dim = 25;
rotate = 0;

%% Test the percentage of success
if test_percentage
    dim = 2;                                         % Intrinsic Dimensionality.
    used_dim = 2;                                              % Dimension used.

    maximizers = trueMaximizer();
    scale = max(1.5*log(used_dim));
    test_bounds = stardardBounds(dim)*scale;                  % Standard Bounds.

    prct = success_prctg(high_dim, 1000, dim, used_dim, test_bounds, ...
        maximizers);                         % Success Percentage by simulation.

    fprintf('Success Percentage by using %d dimensions is approximately: \n%f.\n', ...
        used_dim, prct);
end

%% Random embedding with true intrisic dimensions
%  True maximum is ensured to fall in to the bounds
if embedding
    dim = 2;
    model = rembo(total_iter, dim, high_dim, rotate, 1);
    ditance_log_plot(total_iter, model.f);
end

%% Random embedding with true intrisic dimensions
%  True maximum is NOT ensured to fall in to the bounds
%  Runs are repeated to ensure covering of the true maximum
if embedding_consecutive
    dim = 2;
    rotate = 0;
    run_iter = 125; num_trial = 4;
    f_values = zeros(run_iter*num_trial,1);

    for i = 1:num_trial
        model = rembo(run_iter, dim, high_dim, rotate);
        f_values((i-1)*run_iter+1:i*run_iter) = model.f;
    end

    ditance_log_plot(run_iter*num_trial, f_values);
end

%% High Dim
if run_high_dim
    rotate = 0;
    model = rembo(total_iter, high_dim, high_dim, rotate, 0, 0);

    ditance_log_plot(total_iter, model.f);
end

%% Run rembo.
    function model = rembo(total_iter, dim, high_dim, rotate, force_in_bounds, ...
        embed)
        % total_iter: total number of iterations.
        % dim: embedding dimension.
        % high_dim: ambient dimension.
        % roate: whether to randomly rotate the objective function.
        % force_in_bounds: to force an optimizer in bound by repeatly drawing
        %                  random embedding matrices.
        % embed: Whether to use a randome embedding matrix. (If not
        %        then we effectively use regular BO.)

        if nargin < 5
            force_in_bounds = 0;       % Whether to force an optimizer in bound.
        end

        if nargin < 6
            embed = 1;                        % Whether to use embedding or not.
        end

        if rotate
            % Rotate the objective function.
            [rm, ~] = qr(randn(high_dim, high_dim), 0);
        else
            % Do not rotate the objective function.
            rm = eye(high_dim);
        end

        if embed
            % Generate random projection matrix A.
            [in_bounds, A] = test_fall_in_bound(dim, rm);
            while ~in_bounds && force_in_bounds
                % Ensure that at least one maximizer fall in bound by
                % generating as many random projection matrix A as needed.
                [in_bounds, A] = test_fall_in_bound(dim, rm);
            end
        else
            % By setting A to be identity we do not use embedding here.
            A = eye(high_dim, dim);
        end

        scale = max(1.5*log(dim), 1);
        bounds = stardardBounds(dim)*scale;                 % Initialize bounds.
        obj_fct = @(x)-branin((A*x')');     % Initialize the objective function.

        init_pt = zeros(1, dim);                                % Initial point.
        init_f = obj_fct(init_pt);                     % Evaluate initial point.

        hyp = [ones(dim, 1)*0.1 ; 1];          % Setup initial hyper-parameters.
        hyp = log(hyp);

        % Initialize model.
        model = init_model(dim, bounds, init_pt, init_f, hyp, 1e-10, 'ard');
        % Do optimization.
        model = sparse_opt(obj_fct, total_iter-1, model);
    end

%% Helper functions.
    function ditance_log_plot(total_iter, fvalues)
        maximazer = trueMaximizer();
        figure;
        dis = zeros(total_iter,1);
        for i =1:total_iter
            dis(i) = -branin(maximazer(:, 1)') - max(fvalues(1:i));
        end
        loglog(1:total_iter, dis);
    end

    function [in_bounds, A] = test_fall_in_bound(used_dim, rm)
        maximizers = trueMaximizer();
        test_bounds = stardardBounds(2);
        scale = max(1.5*log(used_dim));
        test_bounds = test_bounds*scale;

        [prct, A] = success_prctg(high_dim, 1, 2, used_dim, test_bounds,...
            maximizers, rm);
        in_bounds =  prct;

        if in_bounds
            fprintf('At least one maximizer in bounds.\n');
        else
            fprintf('NO maximizer in bounds.\n');
        end
    end

    function [maximizers] = trueMaximizer()
        bounds_branin = [-5,10; 0, 15];
        maximizers = [pi, -pi, 9.42478; 2.275, 12.275, 2.475];
        maximizers = bsxfun(@minus, maximizers, bounds_branin(:, 1));
        maximizers = bsxfun(@rdivide, maximizers, bounds_branin(:, 2) - ...
            bounds_branin(:, 1))*2-1;
    end

    function [prct, A] = success_prctg(high_dim, num_trial, dim,...
        used_dim, bounds, maximizers, rm)

        total = 0;
        cmbnts = combntns(1:used_dim,dim);
        num_maximizers = size(maximizers, 2);

        for i = 1:num_trial
            indices = 1:high_dim;
            A = randn(high_dim, used_dim);

            if nargin > 6
                A = rm*A;
            end
            fail = 1;

            for j = 1:size(cmbnts, 1)
                for k = 1:num_maximizers
                    true_maximizer = inv(A(indices(1:dim), cmbnts(j, :))) * ...
                        maximizers(:, k);
                    if ~(sum(true_maximizer <= bounds(:,2)) < dim || ...
                        sum(true_maximizer >= bounds(:,1)) < dim)
                        fail = 0;
                    end
                end
            end
            total = total + fail;
        end
        prct = 1 - total/num_trial;
    end
end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
function syncPos(RECT1,RECT2,RECT3,ref)
switch ref
    case 1
        pos2 = RECT2.getPosition; pos1 = RECT1.getPosition; pos2([2 4]) = pos1([1 3]); RECT2.newPosition(pos2,1);
        pos3 = RECT3.getPosition; pos1 = RECT1.getPosition; pos3([2 4]) = pos1([2 4]); RECT3.newPosition(pos3,1);
    case 2
        pos2 = RECT2.getPosition; pos1 = RECT1.getPosition; pos1([1 3]) = pos2([2 4]); RECT1.newPosition(pos1,1);
        pos3 = RECT3.getPosition; pos2 = RECT2.getPosition; pos3([1 3]) = pos2([1 3]); RECT3.newPosition(pos3,1);
    case 3
        pos3 = RECT3.getPosition; pos1 = RECT1.getPosition; pos1([2 4]) = pos3([2 4]); RECT1.newPosition(pos1,1);
        pos3 = RECT3.getPosition; pos2 = RECT2.getPosition; pos2([1 3]) = pos3([1 3]); RECT2.newPosition(pos2,1);
end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%% patch2tri
% Below is a demonstration of the features of the |patch2tri| function

%%
clear; close all; clc;

%% Syntax
% |[varargout]=patch2tri(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function draw_mandelbrot(x0, x1, y0, y1)

n_grid = 1000;

[x,y] = meshgrid(linspace(x0,x1,n_grid), linspace(y0,y1,n_grid));

c = x + 1i*y;
z = c;
n_iter = 100;
M = zeros(size(c));

for i = 1:n_iter
    z = z.^2 + c;
    M(M==0 & abs(z)>2) = n_iter-i;
end

imagesc(x(1,:), y(:,1),M)
colormap jet

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function set_tick_timestamps(axes, use_milliseconds)

%% Find nice round numbers to use
number_of_ticks = 6;
x_lim = xlim(axes);
xrange = x_lim(2) - x_lim(1);
xrange = xrange / number_of_ticks;
nice = floor(log10(xrange));
frac = 1;
if nice > 1
    frac = 4;
end
nice = round(xrange * frac / 10^nice) / frac * 10^nice;

x_tick = floor(x_lim(1)):nice:ceil(x_lim(2));

x_tick_minutes = floor(x_tick / 60);
x_tick_seconds = mod(x_tick,60);

if use_milliseconds
        labels = compose('%.0f:%05.2f', x_tick_minutes', x_tick_seconds');
else
        labels = compose('%.0f:%04.1f', x_tick_minutes', x_tick_seconds');
end

x_tick=x_tick(1,2:end-1);
labels=labels(2:end-1,1);

xticks(axes, x_tick)
xticklabels(axes,labels);
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
classdef S2FunTri < S2Fun
% a class represeneting a function on the sphere

  properties
    tri          % S2Triangulation
    values = []  % function values
  end

  properties (Dependent = true)
    vertices
    antipodal
  end

  methods

    function sF = S2FunTri(nodes,values)
      % initialize a spherical function

      if isa(nodes,'function_handle')
        n = equispacedS2Grid('resolution',1.5*degree);
        values = nodes(n);
        nodes = n;
      end

      if isa(nodes,'S2Triangulation')
        sF.tri = nodes;
      else
        sF.tri = S2Triangulation(nodes);
      end

      sF.values = values;
    end

    function v = get.vertices(S2F)
      v = S2F.tri.vertices;
    end

    function v = get.antipodal(S2F)
      v = S2F.tri.antipodal;
    end

    function S2F = set.vertices(S2F,v)
      if ~isempty(S2F.values), S2F.values = S2F.eval(v); end
      S2F.tri.vertices = v;
      S2F.tri.update;
    end

  end

  methods (Static = true)

    function sF = demo

      %mtexdata dubna;

      %sF = S2Fun(pf({1}).r,pf({1}).intensities);

      %plot(sF,'upper');

      odf = SantaFe;

      v = equispacedS2Grid('points',1000);

      values = odf.calcPDF(Miller(1,0,0,odf.CS),v);

      sF = S2FunTri(v,values)

      plot(sF,'upper')

    end
  end
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function varargout = xexpintinv(varargin)
%XEXPINTINV EXPINT(1/Z)/Z

switch class(varargin{1})

    case 'double'
        z = varargin{1};
        varargout{1} = (1./z).*expint(1./z);

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        varargout{1} = [];
        varargout{2} = createOperator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/EXPINT called with CHAR argument?');
end

function operator = createOperator

operator = struct('convexity','none','monotonicity','none','definiteness','positive','model','callback');
operator.derivative = @derivative;
operator.range = [0 1];
operator.domain = [1e-8 inf];

function d = derivative(z);
d = (-1./z.^2).*expint(1./z)+(1./z).*(-z.*exp(-1./z)).*(-1./z.^2);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function rs = norm1(s, low, upp)

%tstoolbox/@signal/norm1
%   Syntax:
%     * rs=norm1(s) => low=0 , upp=1
%     * rs=norm1(s, low) => upp=1
%     * rs=norm1(s, low, upp)
%
%   Scale and move signal values to be within [low,upp].
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,3);

if nargin < 2
	low = 0;
end
if nargin < 3
	upp = 1;
end

c = norm1(s.core, low, upp);
rs = signal(c,s);

rs = addhistory(rs, ['(norm1) Transformed signal to be within [' ...
     num2str(low) ',' num2str(upp) ']']);
rs = addcommandlines(rs, 's = norm1(s');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa