text:
function [g, gdata, gprior] = gp_g(w, gp, x, y, varargin)
%GP_G  Evaluate the gradient of energy (GP_E) for Gaussian Process
%
%  Description
%    G = GP_G(W, GP, X, Y, OPTIONS) takes a full GP parameter
%    vector W, GP structure GP, a matrix X of input vectors and a
%    matrix Y of target vectors, and evaluates the gradient G of
%    the energy function (gp_e). Each row of X corresponds to one
%    input vector and each row of Y corresponds to one target
%    vector.
%
%    [G, GDATA, GPRIOR] = GP_G(W, GP, X, Y, OPTIONS) also returns
%    separately 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_E, GP_PAK, GP_UNPAK, GPCF_*
%

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

% 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.

if isfield(gp,'latent_method') && ~strcmp(gp.latent_method,'MCMC')
  % use an inference specific method
  fh_g = gp.fh.g;
  switch nargout
    case {0 1}
      [g] = fh_g(w, gp, x, y, varargin{:});
    case 2
      [g, gdata] = fh_g(w, gp, x, y, varargin{:});
    case 3
      [g, gdata, gprior] = fh_g(w, gp, x, y, varargin{:});
  end
  return
end

ip=inputParser;
ip.FunctionName = 'GP_G';
ip.addRequired('w', @(x) isvector(x) && isreal(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;
if ~all(isfinite(w(:)));
  % instead of stopping to error, return NaN
  g=NaN;
  gdata = NaN;
  gprior = NaN;
  return;
end

% unpak the parameters
gp=gp_unpak(gp, w);
ncf = length(gp.cf);
if isfield(gp.lik, 'nondiagW')
  % Likelihoods with non-diagonal Hessian
  switch gp.lik.type
    case {'LGP', 'LGPC'}
      % Do nothing
    case {'Softmax', 'Multinom'}
      n=size(x,1);
      nout=size(y(:),1)./n;
%       [n,nout]=size(y);
      nl=cumsum([0 repmat(n,1,nout)]);
    otherwise
      n=size(x,1);
      nout=length(gp.comp_cf);

      % Help indices for latent processes
      if ~isfield(gp.lik, 'xtime')
        nl=[0 repmat(n,1,nout)];
      else
        xtime=gp.lik.xtime;
        ntime=size(xtime,1);
        nl=[0 ntime n];
      end
      nl=cumsum(nl);
  end
  if isfield(gp, 'comp_cf')  % own covariance for each ouput component
    multicf = true;
    if length(gp.comp_cf) ~= nout
      error('GP2_G: the number of component vectors in gp.comp_cf must be the same as number of outputs.')
    end
  else
    multicf = false;
  end
else
  n=size(x,1);
end

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

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

switch gp.type
  case 'FULL'
    % ============================================================
    % FULL
    % ============================================================

    if ~isfield(gp.lik, 'nondiagW') || ismember(gp.lik.type, {'LGP' 'LGPC'})
      % Evaluate covariance
      [K, C] = gp_trcov(gp,x);

      if issparse(C)
        % evaluate the sparse inverse
        [LD, notpositivedefinite] = ldlchol(C);
        invC = spinv(LD,1);
        if notpositivedefinite
          % instead of stopping to chol error, return NaN
          g=NaN;
          gdata = NaN;
          gprior = NaN;
          return;
        end
        if  ~isfield(gp,'meanf')
          b = ldlsolve(LD,y);
        else
          [invNM invAt HinvC]=mean_gf(gp,x,C,LD,[],[],y,'gaussian');
        end
      else
        % evaluate the full inverse
        ws1=warning('off','MATLAB:nearlySingularMatrix');
        ws2=warning('off','MATLAB:SingularMatrix');
        invC = inv(C);
        if  ~isfield(gp,'meanf')
          b = C\y;
        else
          [invNM invAt HinvC]=mean_gf(gp,x,C,invC,[],[],y,'gaussian');
        end
        warning(ws1);
        warning(ws2);
      end
    else
      b = zeros(nl(end),1);
      y=y(:);

      switch gp.lik.type
        case 'Coxph'
          % In Cox-Ph, latent processes have different inputs so stacking in invC
          % is not possible.
          invC = zeros(nl(end),nl(end));
          if multicf
            [tmp, C] = gp_trcov(gp, xtime, gp.comp_cf{1});
            invC(1+nl(1):nl(2),1+nl(1):nl(2)) = inv(C);
            b(nl(1)+1:nl(2)) = C\y(nl(1)+1:nl(2));
            [tmp, C] = gp_trcov(gp, x, gp.comp_cf{2});
            invC(1+nl(2):nl(3),1+nl(2):nl(3)) = inv(C);
            b(1+nl(2):nl(3)) = C\y(1+nl(2):nl(3));
          else
            error('Specify covariance function for time process and input process, when using Cox-Ph likelihood');
          end
        otherwise
          invC = zeros(n,n,nout);
          if multicf
            for i1=1:nout
              [tmp, C] = gp_trcov(gp, x, gp.comp_cf{i1});
              invC(:,:,i1) = inv(C);
              b(1+nl(i1):nl(i1+1)) = C\y(1+nl(i1):nl(i1+1));
              %             b(:,i1) = C\y(:,i1);
            end
          else
            [tmp, C] = gp_trcov(gp, x);
            invCtmp = inv(C);
            for i1=1:nout
              invC(:,:,i1) = invCtmp;
              b(1+nl(i1):nl(i1+1)) = C\y(1+nl(i1):nl(i1+1));
              %             b(:,i1) = C\y(:,i1);
            end
          end
      end
    end

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

        gpcf = gp.cf{i};

        if isfield(gp.lik, 'nondiagW') && ~ismember(gp.lik.type, {'LGP' 'LGPC'})
          % check in which components the covariance function is present
          % for likelihoods with non-diagonal Hessian
          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
        end

        if ~(isfield(gp,'derivobs') && gp.derivobs)
          % No derivative observations
          if ~savememory
            if isfield(gp.lik, 'nondiagW') && isfield(gp,'comp_cf') && ~isempty(intersect(gp.comp_cf{1},i)) && isfield(gp.lik, 'xtime')
              DKffc = gpcf.fh.cfg(gpcf, xtime);
            else
              DKffc = gpcf.fh.cfg(gpcf, x);
            end
            np=length(DKffc);
          else
            % If savememory option is used, just get the number of
            % hyperparameters and calculate gradients later
            np=gpcf.fh.cfg(gpcf,[],[],[],0);
          end
          gprior_cf = -gpcf.fh.lpg(gpcf);
        else
          [n m]=size(x);
          %Case: input dimension is 1
          if m==1

            DKffa = gpcf.fh.cfg(gpcf, x);
            DKdf = gpcf.fh.cfdg(gpcf, x);
            DKdd = gpcf.fh.cfdg2(gpcf, x);
            gprior_cf = -gpcf.fh.lpg(gpcf);

            % DKff{1} -- d K / d magnSigma2
            % DKff{2} -- d K / d lengthScale
            DKffc{1} = [DKffa{1}, DKdf{1}'; DKdf{1}, DKdd{1}];
            DKffc{2} = [DKffa{2}, DKdf{2}'; DKdf{2}, DKdd{2}];
            np=2;

            %Case: input dimension is >1
          else
            DKffa = gpcf.fh.cfg(gpcf, x);
            DKdf = gpcf.fh.cfdg(gpcf, x);
            DKdd = gpcf.fh.cfdg2(gpcf, x);
            gprior_cf = -gpcf.fh.lpg(gpcf);

            %Check whether ARD method is in use (with gpcf_sexp)
            Ard=length(gpcf.lengthScale);

            % DKff{1} - d K / d magnSigma2
            % DKff{2:end} - d K / d lengthScale(1:end)
            for i=1:2
              DKffc{i}=[DKffa{i} DKdf{i}';DKdf{i} DKdd{i}];
            end

            %If ARD is in use
            if Ard>1
              for i=2+1:2+Ard-1
                DKffc{i}=[DKffa{i} DKdf{i}';DKdf{i} DKdd{i}];
              end
            end
            np=length(DKffc);
          end
        end

        % Are there specified mean functions
        if  ~isfield(gp,'meanf')
          % Evaluate the gradient with respect to covariance function
          % parameters
          for i2 = 1:np
            if savememory
              if isfield(gp.lik, 'nondiagW') && isfield(gp,'comp_cf') && ~isempty(intersect(gp.comp_cf{1},i)) && isfield(gp.lik, 'xtime')
                DKff=gpcf.fh.cfg(gpcf,xtime,[],[],i2);
              else
                DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
              end
            else
              DKff=DKffc{i2};
            end
            i1 = i1+1;
            if ~isfield(gp.lik, 'nondiagW')
              Bdl = b'*(DKff*b);
              Cdl = sum(sum(invC.*DKff)); % help arguments
            else
              % Non-diagonalizable likelihoods
              Bdl=0; Cdl=0;
              if isfield(gp.lik,'xtime');
                if do(1)
                  Bdl = Bdl + b(1:ntime)'*(DKff*b(1:ntime));
                  Cdl = Cdl + sum(sum(invC(1:ntime,1:ntime).*DKff)); % help arguments
                end
                if do(2)
                  Bdl = Bdl + b(ntime+1:end)'*(DKff*b(ntime+1:end));
                  Cdl = Cdl + sum(sum(invC(ntime+1:end,ntime+1:end).*DKff)); % help arguments
                end
              else
                for z1=1:nout
                  if do(z1)
                    Bdl = Bdl + b(1+nl(z1):nl(z1+1))'*(DKff*b(1+nl(z1):nl(z1+1)));
                    Cdl = Cdl + sum(sum(invC(:,:,z1).*DKff)); % help arguments
                  end
                end
              end
            end
            gdata(i1)=0.5.*(Cdl - Bdl);
            gprior(i1) = gprior_cf(i2);
          end
        else
          for i2 = 1:np
            if savememory
              DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
            else
              DKff=DKffc{i2};
            end
            i1=i1+1;
            dA = -1*HinvC*DKff*HinvC';                  % d A / d th
            trA = sum(invAt(:).*dA(:));                 % d log(|A|) / dth
            dMNM = invNM'*(DKff*invNM);           % d M'*N*M / d th
            trK = sum(sum(invC.*DKff));       % d log(Ky�?�) / d th
            gdata(i1)=0.5*(-1*dMNM + trK + trA);
            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 Gaussian likelihood function parameters
    if ~isempty(strfind(gp.infer_params, 'likelihood')) && isfield(gp.lik.fh,'trcov')
      % Evaluate the gradient from Gaussian likelihood
      DCff = gp.lik.fh.cfg(gp.lik, x);
      gprior_lik = -gp.lik.fh.lpg(gp.lik);
      for i2 = 1:length(DCff)
        i1 = i1+1;
        if ~isfield(gp,'meanf')
          if size(DCff{i2}) > 1
            yKy = b'*(DCff{i2}*b);
            trK = sum(sum(invC.*DCff{i2})); % help arguments
            gdata_zeromean(i1)=0.5.*(trK - yKy);
          else
            yKy=DCff{i2}.*(b'*b);
            trK = DCff{i2}.*(trace(invC));
            gdata_zeromean(i1)=0.5.*(trK - yKy);
          end
          gdata(i1)=gdata_zeromean(i1);
        else
          if size(DCff{i2}) > 1
            trK = sum(sum(invC.*DCff{i2})); % help arguments
          else
            trK = DCff{i2}.*(trace(invC));
          end
          dA = -1*HinvC*DCff{i2}*HinvC';            % d A / d th
          trA = sum(invAt(:).*dA(:));               % d log(|A|) / dth
          dMNM = invNM'*(DCff{i2}*invNM);           % d M'*N*M / d th
          gdata(i1)=0.5*(-1*dMNM + trA + trK);
        end
        gprior(i1) = gprior_lik(i2);
      end

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

    if ~isempty(strfind(gp.infer_params, 'mean')) && isfield(gp,'meanf')
        notpositivedefinite2 = 0; notpositivedefinite3 = 0;

        nmf=numel(gp.meanf);
        [H,b,B]=mean_prep(gp,x,[]);
        M = H'*b-y;

        if issparse(C)
            [LD, notpositivedefinite] = ldlchol(C);
            if ~notpositivedefinite
                KH = ldlsolve(LD, H');
            end
            [LB, notpositivedefinite2] = chol(B);
            if ~notpositivedefinite2
                A = LB\(LB'\eye(size(B))) + H*KH;
                LA = chol(A);
                a = ldlsolve(LD, M) - KH*(LA\(LA'\(KH'*M)));
                iNH = ldlsolve(LD, H') - KH*(LA\(LA'\(KH'*H')));
            end
        else
            N = C + H'*B*H;
            [LN, notpositivedefinite3] = chol(N);
            if ~notpositivedefinite3
                a = LN\(LN'\M);
                iNH = LN\(LN'\H');
            end
        end
        if (~notpositivedefinite2 && ~notpositivedefinite3)
            Ha=H*a;
            g_bb = (-H*a)';     % b and B parameters are log transformed in packing
            indB = find(B>0);
            for i=1:length(indB)
                Bt = zeros(size(B)); Bt(indB(i))=1;
                BH = Bt*H;
                g_B(i) = 0.5* ( Ha'*Bt*Ha - sum(sum(iNH.*(BH'))) );
            end
            g_BB = g_B.*B(indB)';
            for i=1:nmf
                gpmf = gp.meanf{i};
                [lpg_b, lpg_B] = gpmf.fh.lpg(gpmf);
                ll=length(lpg_b);
                gdata = [gdata -g_bb((i-1)*ll+1:i*ll)];
                gprior = [gprior -lpg_b];
                ll=length(lpg_B);
                gdata = [gdata -g_B((i-1)*ll+1:i*ll)];
                gprior = [gprior -lpg_B];
            end
        else
          g=NaN;
          gdata = NaN;
          gprior = NaN;
          return
        end
    end

    g = gdata + gprior;

  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;
    DKuu_u = 0;
    DKuf_u = 0;

    % First evaluate the needed covariance matrices
    % v defines that parameter is a vector
    [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  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
    [Luu, notpositivedefinite] = chol(K_uu,'lower');
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    % Evaluate the Lambda (La)
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    % Here we need only the diag(Q_ff), which is evaluated below
    B=Luu\(K_fu');
    Qv_ff=sum(B.^2)';
    Lav = Cv_ff-Qv_ff;   % 1 x f, Vector of diagonal elements
                         % iLaKfu = diag(inv(Lav))*K_fu = inv(La)*K_fu
    iLaKfu = zeros(size(K_fu));  % f x u,
    for i=1:n
      iLaKfu(i,:) = K_fu(i,:)./Lav(i);  % f x u
    end
    % ... then evaluate some help matrices.
    % A = K_uu+K_uf*inv(La)*K_fu
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;               % Ensure symmetry
    [A, notpositivedefinite] = chol(A,'upper');
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    L = iLaKfu/A;
    b = y'./Lav' - (y'*L)*L';
    iKuuKuf = Luu'\(Luu\K_fu');
    La = Lav;
    LL = sum(L.*L,2);

    % =================================================================
    % Gradient with respect to covariance function parameters
    if ~isempty(strfind(gp.infer_params, 'covariance'))
      % Loop over the covariance functions
      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(DKffc);
        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
          KfuiKuuKuu = iKuuKuf'*DKuu;
          gdata(i1) = -0.5.*((2*b*DKuf'-(b*KfuiKuuKuu))*(iKuuKuf*b') + 2.*sum(sum(L'.*(L'*DKuf'*iKuuKuf))) - ...
                             sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));

          gdata(i1) = gdata(i1) - 0.5.*(b.*DKff')*b';
          gdata(i1) = gdata(i1) + 0.5.*(2.*b.*sum(DKuf'.*iKuuKuf',2)'*b'- b.*sum(KfuiKuuKuu.*iKuuKuf',2)'*b');
          gdata(i1) = gdata(i1) + 0.5.*(sum(DKff./La) - sum(LL.*DKff));
          gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*sum(DKuf'.*iKuuKuf',2)) - sum(LL.*sum(KfuiKuuKuu.*iKuuKuf',2)));
          gprior(i1) = gprior_cf(i2);
        end

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

    % =================================================================
    % Gradient with respect to Gaussian likelihood function parameters
    if ~isempty(strfind(gp.infer_params, 'likelihood')) && isfield(gp.lik.fh,'trcov')
      % Evaluate the gradient from Gaussian likelihood
      DCff = gp.lik.fh.cfg(gp.lik, x);
      gprior_lik = -gp.lik.fh.lpg(gp.lik);
      for i2 = 1:length(DCff)
        i1 = i1+1;
        gdata(i1)= -0.5*DCff{i2}.*b*b';
        gdata(i1)= gdata(i1) + 0.5*sum(DCff{i2}./La-sum(L.*L,2).*DCff{i2});
        gprior(i1) = gprior_lik(i2);
      end
      % Set the gradients of hyperparameter
      if length(gprior_lik) > length(DCff)
        for i2=length(DCff)+1:length(gprior_lik)
          i1 = i1+1;
          gdata(i1) = 0;
          gprior(i1) = gprior_lik(i2);
        end
      end
    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
            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;
              KfuiKuuKuu = iKuuKuf'*DKuu{i2};

              gdata(i1) = gdata(i1) - 0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuKuu))*(iKuuKuf*b') + ...
                                          2.*sum(sum(L'.*(L'*DKuf{i2}'*iKuuKuf))) - sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));
              gdata(i1) = gdata(i1) + 0.5.*(2.*b.*sum(DKuf{i2}'.*iKuuKuf',2)'*b'- b.*sum(KfuiKuuKuu.*iKuuKuf',2)'*b');
              gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*sum(DKuf{i2}'.*iKuuKuf',2)) - ...
                                          sum(LL.*sum(KfuiKuuKuu.*iKuuKuf',2)));
            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;
    ind = gp.tr_index;
    DKuu_u = 0;
    DKuf_u = 0;

    % First evaluate the needed covariance matrices
    % if they are not in the memory
    % v defines that parameter is a 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
    [Luu, notpositivedefinite] = chol(K_uu,'lower');
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    % Evaluate the Lambda (La)
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    % Here we need only the diag(Q_ff), which is evaluated below
    %B=K_fu/Luu;
    B=Luu\K_fu';
    iLaKfu = zeros(size(K_fu));  % f x u
    for i=1:length(ind)
      Qbl_ff = B(:,ind{i})'*B(:,ind{i});
      [Kbl_ff, Cbl_ff] = gp_trcov(gp, x(ind{i},:));
      la = Cbl_ff - Qbl_ff;
      La{i} = (la + la')./2;
      iLaKfu(ind{i},:) = La{i}\K_fu(ind{i},:);
    end
    % ... then evaluate some help matrices.
    % A = chol(K_uu+K_uf*inv(La)*K_fu))
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;            % Ensure symmetry

    [LA,notpositivedefinite]=chol(A,'upper');
    if notpositivedefinite
      % instead of stopping to chol error, return NaN
      g=NaN; gdata = NaN; gprior = NaN;
      return;
    end
    L = iLaKfu/LA;
    b = zeros(1,n);
    b_apu=(y'*L)*L';
    for i=1:length(ind)
      b(ind{i}) = y(ind{i})'/La{i} - b_apu(ind{i});
    end
    iKuuKuf = Luu'\(Luu\K_fu');

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

    if ~isempty(strfind(gp.infer_params, 'covariance'))
      % Loop over the  covariance functions
      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;
          %            H = (2*K_uf'- KfuiKuuKuu)*iKuuKuf;
          % Here we evaluate  gdata = -0.5.* (b*H*b' + trace(L*L'H)
          gdata(i1) = -0.5.*((2*b*DKuf'-(b*KfuiKuuKuu))*(iKuuKuf*b') + 2.*sum(sum(L'.*(L'*DKuf'*iKuuKuf))) - ...
                             sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));
          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.*(-b(ind{kk})*DKff*b(ind{kk})' ...
                        + 2.*b(ind{kk})*DKuf(:,ind{kk})'*iKuuKuf(:,ind{kk})*b(ind{kk})'- ...
                        b(ind{kk})*KfuiKuuKuu(ind{kk},:)*iKuuKuf(:,ind{kk})*b(ind{kk})' ...
                        +trace(La{kk}\DKff)...
                        - trace(L(ind{kk},:)*(L(ind{kk},:)'*DKff)) ...
                        + 2.*sum(sum(L(ind{kk},:)'.*(L(ind{kk},:)'*DKuf(:,ind{kk})'*iKuuKuf(:,ind{kk})))) - ...
                        sum(sum(L(ind{kk},:)'.*((L(ind{kk},:)'*KfuiKuuKuu(ind{kk},:))*iKuuKuf(:,ind{kk})))));
          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 Gaussian likelihood function parameters
    if ~isempty(strfind(gp.infer_params, 'likelihood')) && isfield(gp.lik.fh,'trcov')
      % Evaluate the gradient from Gaussian likelihood
      DCff = gp.lik.fh.cfg(gp.lik, x);
      gprior_lik = -gp.lik.fh.lpg(gp.lik);
      for i2 = 1:length(DCff)
        i1 = i1+1;
        gdata(i1)= -0.5*DCff{i2}.*b*b';
        for kk=1:length(ind)
          gdata(i1)= gdata(i1) + 0.5*trace((inv(La{kk})-L(ind{kk},:)*L(ind{kk},:)')).*DCff{i2};
        end
        gprior(i1) = gprior_lik(i2);
      end
      % Set the gradients of hyperparameter
      if length(gprior_lik) > length(DCff)
        for i2=length(DCff)+1:length(gprior_lik)
          i1 = i1+1;
          gdata(i1) = 0;
          gprior(i1) = gprior_lik(i2);
        end
      end
    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
            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;
              KfuiKuuDKuu_u = iKuuKuf'*DKuu{i2};
              gdata(i1) = gdata(i1) -0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuDKuu_u))*(iKuuKuf*b') + 2.*sum(sum(L'.*((L'*DKuf{i2}')*iKuuKuf))) - ...
                                         sum(sum(L'.*((L'*KfuiKuuDKuu_u)*iKuuKuf))));

              for kk=1:length(ind)
                gdata(i1) = gdata(i1) + 0.5.*(2.*b(ind{kk})*DKuf{i2}(:,ind{kk})'*iKuuKuf(:,ind{kk})*b(ind{kk})'- ...
                                            b(ind{kk})*KfuiKuuDKuu_u(ind{kk},:)*iKuuKuf(:,ind{kk})*b(ind{kk})' ...
                                            + 2.*sum(sum(L(ind{kk},:)'.*(L(ind{kk},:)'*DKuf{i2}(:,ind{kk})'*iKuuKuf(:,ind{kk})))) - ...
                                            sum(sum(L(ind{kk},:)'.*((L(ind{kk},:)'*KfuiKuuDKuu_u(ind{kk},:))*iKuuKuf(:,ind{kk})))));
              end
            end
          end
        end
      end
    end
    g = gdata + gprior;

  case 'CS+FIC'
    % ============================================================
    % CS+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;
    DKuu_u = 0;
    DKuf_u = 0;

    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 the needed covariance matrices
    % if they are not in the memory
    % v defines that parameter is a vector
    [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  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
    [Luu, notpositivedefinite] = chol(K_uu,'lower');
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    % Evaluate the Lambda (La)
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    % Here we need only the diag(Q_ff), which is evaluated below
    B=Luu\(K_fu');
    Qv_ff=sum(B.^2)';
    Lav = Cv_ff-Qv_ff;   % 1 x f, Vector of diagonal elements

    gp.cf = cf2;
    K_cs = gp_trcov(gp,x);
    La = sparse(1:n,1:n,Lav,n,n) + K_cs;
    gp.cf = cf_orig;

    LD = ldlchol(La);
    %        iLaKfu = La\K_fu;
    iLaKfu = ldlsolve(LD, K_fu);

    % ... then evaluate some help matrices.
    % A = chol(K_uu+K_uf*inv(La)*K_fu))
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;            % Ensure symmetry
    [LA, notpositivedefinite] = chol(A,'upper');
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    L = iLaKfu/LA;
    %b = y'/La - (y'*L)*L';
    b = ldlsolve(LD,y)' - (y'*L)*L';

    siLa = spinv(La);
    idiagLa = diag(siLa);
    iKuuKuf = K_uu\K_fu';
    LL = sum(L.*L,2);

    % =================================================================
    % Gradient with respect to covariance function parameters
    if ~isempty(strfind(gp.infer_params, 'covariance'))
      % Loop over covariance functions
      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
            KfuiKuuKuu = iKuuKuf'*DKuu;
            gdata(i1) = -0.5.*((2*b*DKuf'-(b*KfuiKuuKuu))*(iKuuKuf*b') + 2.*sum(sum(L'.*(L'*DKuf'*iKuuKuf))) - ...
                               sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));

            temp1 = sum(KfuiKuuKuu.*iKuuKuf',2);
            temp2 = sum(DKuf'.*iKuuKuf',2);
            temp3 = 2.*DKuf' - KfuiKuuKuu;
            gdata(i1) = gdata(i1) - 0.5.*(b.*DKff')*b';
            gdata(i1) = gdata(i1) + 0.5.*(2.*b.*temp2'*b'- b.*temp1'*b');
            gdata(i1) = gdata(i1) + 0.5.*(sum(idiagLa.*DKff - LL.*DKff));   % corrected
            gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*temp2) - sum(LL.*temp1));

            %gdata(i1) = gdata(i1) + 0.5.*sum(sum(La\((2.*K_uf') - KfuiKuuKuu).*iKuuKuf',2));
            gdata(i1) = gdata(i1) + 0.5.*sum(sum(ldlsolve(LD,temp3).*iKuuKuf',2));
            gdata(i1) = gdata(i1) - 0.5.*( idiagLa'*(sum(temp3.*iKuuKuf',2)) ); % corrected
            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
            gdata(i1) = 0.5*(sum(sum(siLa.*DKff',2)) - sum(sum(L.*(L'*DKff')')) - b*DKff*b');
            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 Gaussian likelihood function parameters
    if ~isempty(strfind(gp.infer_params, 'likelihood')) && isfield(gp.lik.fh,'trcov')
      % Evaluate the gradient from Gaussian likelihood
      DCff = gp.lik.fh.cfg(gp.lik, x);
      gprior_lik = -gp.lik.fh.lpg(gp.lik);
      for i2 = 1:length(DCff)
        i1 = i1+1;
        gdata(i1)= -0.5*DCff{i2}.*b*b';
        gdata(i1)= gdata(i1) + 0.5*sum(idiagLa-LL).*DCff{i2};
        gprior(i1) = gprior_lik(i2);
      end

      % Set the gradients of hyperparameter
      if length(gprior_lik) > length(DCff)
        for i2=length(DCff)+1:length(gprior_lik)
          i1 = i1+1;
          gdata(i1) = 0;
          gprior(i1) = gprior_lik(i2);
        end
      end
    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
              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;
                KfuiKuuKuu = iKuuKuf'*DKuu{i2};

                gdata(i1) = gdata(i1) - 0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuKuu))*(iKuuKuf*b') + ...
                                            2.*sum(sum(L'.*(L'*DKuf{i2}'*iKuuKuf))) - sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));
                gdata(i1) = gdata(i1) + 0.5.*(2.*b.*sum(DKuf{i2}'.*iKuuKuf',2)'*b'- b.*sum(KfuiKuuKuu.*iKuuKuf',2)'*b');
                gdata(i1) = gdata(i1) + 0.5.*(2.*sum(sum(L.*L,2).*sum(DKuf{i2}'.*iKuuKuf',2)) - ...
                                            sum(sum(L.*L,2).*sum(KfuiKuuKuu.*iKuuKuf',2)));

                gdata(i1) = gdata(i1) + 0.5.*sum(sum(ldlsolve(LD,(2.*DKuf{i2}') - KfuiKuuKuu).*iKuuKuf',2));
                gdata(i1) = gdata(i1) - 0.5.*( idiagLa'*(sum((2.*DKuf{i2}' - KfuiKuuKuu).*iKuuKuf',2)) ); % corrected
              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;
    DKuu_u = 0;
    DKuf_u = 0;

    % First evaluate the needed covariance matrices
    % v defines that parameter is a vector
    [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  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
    [Luu, notpositivedefinite] = chol(K_uu,'lower');
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    % Evaluate the Lambda (La)
    % Q_ff = K_fu*inv(K_uu)*K_fu'
    % Here we need only the diag(Q_ff), which is evaluated below
    B=Luu\(K_fu');
    Qv_ff=sum(B.^2)';
    Lav = Cv_ff-Kv_ff;   % 1 x f, Vector of diagonal elements
                         % iLaKfu = diag(inv(Lav))*K_fu = inv(La)*K_fu
    iLaKfu = zeros(size(K_fu));  % f x u,
    for i=1:n
      iLaKfu(i,:) = K_fu(i,:)./Lav(i);  % f x u
    end
    % ... then evaluate some help matrices.
    % A = K_uu+K_uf*inv(La)*K_fu
    A = K_uu+K_fu'*iLaKfu;
    A = (A+A')./2;               % Ensure symmetry
    [LA, notpositivedefinite] = chol(A);
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    L = iLaKfu/LA;
    b = y'./Lav' - (y'*L)*L';
    iKuuKuf = Luu'\(Luu\K_fu');
    La = Lav;
    LL = sum(L.*L,2);
    iLav=1./Lav;

    LL1=iLav-LL;

    % =================================================================

    if ~isempty(strfind(gp.infer_params, 'covariance'))
      % Loop over the covariance functions
      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
          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
            % If savememory option is used, just get the number of
            % hyperparameters and calculate gradients later
            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

          KfuiKuuKuu = iKuuKuf'*DKuu;
          gdata(i1) = -0.5.*((2*b*DKuf'-(b*KfuiKuuKuu))*(iKuuKuf*b'));
          gdata(i1) = gdata(i1) + 0.5.*(2.*(sum(iLav'*sum(DKuf'.*iKuuKuf',2))-sum(sum(L'.*(L'*DKuf'*iKuuKuf))))...
                                        - sum(iLav'*sum(KfuiKuuKuu.*iKuuKuf',2))+ sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));

          if strcmp(gp.type, 'VAR')
            gdata(i1) = gdata(i1) + 0.5.*(sum(iLav.*DKff)-2.*sum(iLav'*sum(DKuf'.*iKuuKuf',2)) + ...
                                          sum(iLav'*sum(KfuiKuuKuu.*iKuuKuf',2))); % trace-term derivative
          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 Gaussian likelihood function parameters
    if ~isempty(strfind(gp.infer_params, 'likelihood')) && isfield(gp.lik.fh,'trcov')
      % Evaluate the gradient from Gaussian likelihood
      DCff = gp.lik.fh.cfg(gp.lik, x);
      gprior_lik = -gp.lik.fh.lpg(gp.lik);
      for i2 = 1:length(DCff)
        i1 = i1+1;
        gdata(i1)= -0.5*DCff{i2}.*b*b';
        gdata(i1)= gdata(i1) + 0.5*sum(DCff{i2}./La-sum(L.*L,2).*DCff{i2});
        if strcmp(gp.type, 'VAR')
          gdata(i1)= gdata(i1) + 0.5*(sum((Kv_ff-Qv_ff)./La));
        end

        gprior(i1) = gprior_lik(i2);
      end
      % Set the gradients of hyperparameter
      if length(gprior_lik) > length(DCff)
        for i2=length(DCff)+1:length(gprior_lik)
          i1 = i1+1;
          gdata(i1) = 0;
          gprior(i1) = gprior_lik(i2);
        end
      end
    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
            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;
              KfuiKuuKuu = iKuuKuf'*DKuu{i2};
              gdata(i1) = gdata(i1) - 0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuKuu))*(iKuuKuf*b'));
              gdata(i1) = gdata(i1) + 0.5.*(2.*(sum(iLav'*sum(DKuf{i2}'.*iKuuKuf',2))-sum(sum(L'.*(L'*DKuf{i2}'*iKuuKuf))))...
                                          - sum(iLav'*sum(KfuiKuuKuu.*iKuuKuf',2))+ sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));

              if strcmp(gp.type, 'VAR')
                gdata(i1) = gdata(i1) + 0.5.*(0-2.*sum(iLav'*sum(DKuf{i2}'.*iKuuKuf',2)) + ...
                                            sum(iLav'*sum(KfuiKuuKuu.*iKuuKuf',2)));
              end
            end
          end
        end
      end
    end

    g = gdata + gprior;

  case 'SSGP'
    % ============================================================
    % SSGP
    % ============================================================
    % Predictions with sparse spectral sampling approximation for GP
    % The approximation is proposed by M. Lazaro-Gredilla, J.
    % Quinonero-Candela and A. Figueiras-Vidal in Microsoft
    % Research technical report MSR-TR-2007-152 (November 2007)
    % NOTE! This does not work at the moment.

    % First evaluate the needed covariance matrices
    % v defines that parameter is a vector
    [Phi, S] = gp_trcov(gp, x);        % n x m and nxn sparse matrices
    Sv = diag(S);

    m = size(Phi,2);

    A = eye(m,m) + Phi'*(S\Phi);
    [LA, notpositivedefinite] = chol(A,'lower');
    if notpositivedefinite
      % If not positive definite, return NaN
      g=NaN; gdata=NaN; gprior=NaN;
      return;
    end
    L = (S\Phi)/A';

    b = y'./Sv' - (y'*L)*L';
    iSPhi = S\Phi;

    % =================================================================
    if ~isempty(strfind(gp.infer_params,'covariance'))
      % Loop over the covariance functions
      for i=1:ncf
        i1=0;
        if ~isempty(gprior)
          i1 = length(gprior);
        end

        gpcf = gp.cf{i};

        % Get the gradients of the covariance matrices
        % and gprior from gpcf_* structures
        DKff = gpcf.fh.cfg(gpcf, x);
        gprior_cf = -gpcf.fh.lpg(gpcf);

        % Evaluate the gradient with respect to lengthScale
        for i2 = 1:length(DKff)
          i1 = i1+1;
          iSDPhi = S\DKff{i2};

          gdata(i1) = 0.5*( sum(sum(iSDPhi.*Phi,2)) + sum(sum(iSPhi.*DKff{i2},2)) );
          gdata(i1) = gdata(i1) - 0.5*( sum(sum(L'.*(L'*DKff{i2}*Phi' + L'*Phi*DKff{i2}'),1)) );
          gdata(i1) = gdata(i1) - 0.5*(b*DKff{i2}*Phi' + b*Phi*DKff{i2}')*b';
          gprior(i1) = gprior_cf(i2);
        end

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

    % =================================================================
    % Gradient with respect to Gaussian likelihood function parameters
    if ~isempty(strfind(gp.infer_params, 'likelihood')) && isfield(gp.lik.fh,'trcov')
      % Evaluate the gradient from Gaussian likelihood
      DCff = gp.lik.fh.cfg(gp.lik, x);
      gprior_lik = -gp.lik.fh.lpg(gp.lik);
      for i2 = 1:length(DCff)
        i1 = i1+1;
        gdata(i1)= -0.5*DCff{i2}.*b*b';
        gdata(i1)= gdata(i1) + 0.5*sum(1./Sv-sum(L.*L,2)).*DCff{i2};
        gprior(i1) = gprior_lik(i2);
      end

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

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

        gpcf = gp.cf{i};

        gpcf.GPtype = gp.type;
        [gprior_ind, DKuu, DKuf] = gpcf.fh.gind(gpcf, x, y, g_ind, gdata_ind, gprior_ind);

        for i2 = 1:length(DKuu)
          KfuiKuuKuu = iKuuKuf'*DKuu{i2};

          gdata_ind(i2) = gdata_ind(i2) - 0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuKuu))*(iKuuKuf*b') + ...
                                                2.*sum(sum(L'.*(L'*DKuf{i2}'*iKuuKuf))) - sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));
          gdata_ind(i2) = gdata_ind(i2) + 0.5.*(2.*b.*sum(DKuf{i2}'.*iKuuKuf',2)'*b'- b.*sum(KfuiKuuKuu.*iKuuKuf',2)'*b');
          gdata_ind(i2) = gdata_ind(i2) + 0.5.*(2.*sum(sum(L.*L,2).*sum(DKuf{i2}'.*iKuuKuf',2)) - ...
                                                sum(sum(L.*L,2).*sum(KfuiKuuKuu.*iKuuKuf',2)));
        end
      end
    end

    g = gdata + gprior;

end


text:
function [forecast_record,forecast_estimates]=panel2forecast(N,n,p,k,data_endo_a,data_exo_p,It,Bu,beta_gibbs,sigma_gibbs,Fperiods,const,Fband,Fstartlocation,favar)

% initiate the cell recording the forecast draws
forecast_record={};
forecast_estimates={};

% because the forecasts have to be computed for each unit, loop over units
for ii=1:N
% run the Gibbs sampler for unit ii
forecast_record(:,:,ii)=bear.forecast(data_endo_a(:,:,ii),data_exo_p,It,Bu,beta_gibbs,sigma_gibbs,Fperiods,n,p,k,const,Fstartlocation,favar);
% obtain point estimates and credibility intervals for unit ii
forecast_estimates(:,:,ii)=bear.festimates(forecast_record(:,:,ii),n,Fperiods,Fband);
end

author: european-central-bank

repo: BEAR-toolbox


text:
function uI = faceinterpolate3(u,node,elem,quadOrder)
%% FACEINTERPOLATE3 interpolate to face elements RT0.
%
% uI = faceinterpolate3(u,node,face) interpolates a given function u
% into the lowest order RT0. The coefficient is given by the face integral
% int_f u*n ds.
%
% uI = faceinterpolate3(u,node,elem) when |face| is not given, the function
% will generate the face by |[elem2face,face] = dof3face(elem)| and the
% ascend order is used |elem = sortelem3(elem)|.
%
% uI = faceinterpolate3(u,node,face,6) the last input is the quadrature
% order; see quadpts.
%
% Example
%
%   [node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
%   maxIt = 3;
%   pde = mixBCdata3;
%   err = zeros(maxIt,1);
%   h = zeros(maxIt,1);
%   for i = 1:maxIt
%      [node,elem] = uniformrefine3(node,elem);
%      [elem2face,face] = dof3face(elem);
%      uI = faceinterpolate3(pde.Du,node,face);
%      err(i) = getL2error3RT0(node,elem,pde.Du,uI);
%      h(i) = 2^(-i);
%   end
%   figure;
%   showrateh(h,err,2,'-+','|| u - u_I ||');
%
% See also edgeinterpolate, edgeinterpolate1, edgeinterpolate2
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

% if ~exist('elemType','var'), elemType = 'RT0'; end

%% Construct Data Structure
if size(elem,2) == 3 % the input elem is face
    face = elem;
else
    elem = sortelem3(elem);
    [~,face] = dof3face(elem);
end

%% normal vector of each face
NF = size(face,1);
v12 = node(face(:,2),:) - node(face(:,1),:);
v13 = node(face(:,3),:) - node(face(:,1),:);
nVec = mycross(v12,v13); % |nVec| = 2*area;

%% assemble the right hand side
if ~exist('quadOrder','var'), quadOrder = 3; end
[lambda,weight] = quadpts(quadOrder);
nQuad = size(lambda,1);
uI = zeros(NF,1);
for p = 1:nQuad
    pxyz = lambda(p,1)*node(face(:,1),:) ...
		 + lambda(p,2)*node(face(:,2),:) ...
		 + lambda(p,3)*node(face(:,3),:);
    flux = u(pxyz);
    uI = uI + weight(p)*dot(flux,nVec,2)/2;
end
% if strcmp(elemType,'BDM1')
% end

repo: ifem


text:
%SimulationStateSideEffectItem
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 9/12/2010
classdef SimulationStateSideEffectItem
    properties (SetAccess = protected)
        stateName          %name of simulation state to be updated
        propertyName       %name of state property to be updated
        indexName          %name of index over state property
        componentIndex     %component (row) index within state property/property index to be updated
        compartmentIndex   %compartment (column) index within state property to be updated
        deltaValue         %change in value of state property element
    end

    %constructor
    methods
        function this = SimulationStateSideEffectItem(stateName, propertyName, indexName, componentIndex, compartmentIndex, deltaValue)
            this.stateName        = stateName;
            this.propertyName     = propertyName;
            this.indexName        = indexName;
            this.componentIndex   = componentIndex;
            this.compartmentIndex = compartmentIndex;
            this.deltaValue       = deltaValue;
        end
    end

    %method to update simulation state
    methods
        function simulation = updateSimulationState(this, simulation)
            deltaMass = 0;

            for i = 1:numel(this)
                state = simulation.state(this(i).stateName);

                if this(i).indexName
                    iComponent = state.(this(i).indexName)(this(i).componentIndex);
                else
                    iComponent = this(i).componentIndex;
                end

                iCompartment = this(i).compartmentIndex;
                state.(this(i).propertyName)(iComponent, iCompartment) = ...
                    state.(this(i).propertyName)(iComponent, iCompartment) + ...
                    this(i).deltaValue;
                deltaMass = deltaMass + ...
                    state.molecularWeights(iComponent) * this(i).deltaValue;
            end

            if deltaMass ~= 0
                warning('WholeCell:warning',[...
                    'Simulation state side effects should preserve mass (eg. '...
                    'mass should not be passed from processes state properties to '...
                    'SimulationStateSideEffect objects, and each SimulationStateSideEffectItem '...
                    'array should represent a mass-balanced change to the state of '...
                    'the simulation)'])
            end
        end
    end

    %print to standard output
    methods
        function disp(this, indent)
            if nargin < 2
                indent = 0;
            end

            sizCellArr = cellfun(@num2str, num2cell(size(this)), 'UniformOutput', false);
            fprintf([repmat('  ',1, indent) '%s SimulationStateSideEffectItem\n'], ...
                strjoin('x', sizCellArr{:}));
            stateNameMaxLength = max(cellfun(@length, {this.stateName}));
            propertyNameMaxLength = max(cellfun(@length, {this.propertyName}));
            indexNameMaxLength  = max(cellfun(@length, {this.indexName}));
            for i = 1:numel(this)
                fprintf([repmat('  ',1, indent+1) '%' num2str(stateNameMaxLength) 's\t%' num2str(propertyNameMaxLength) 's\t%' num2str(indexNameMaxLength) 's\t%8d\t%8d\t%4d\n'], ...
                    this(i).stateName, this(i).propertyName, this(i).indexName, this(i).componentIndex, this(i).compartmentIndex, this(i).deltaValue);
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function varargout=lpczz2cc(varargin)
%
% For calling details please see v_lpczz2cc.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpczz2cc.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpczz2cc.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpczz2cc(varargin{:});
else
    v_lpczz2cc(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function net = caffe_load_model( net_def_file, net_weights_file, phase )
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

if ~exist('phase', 'var'),  phase = 'test'; end
assert(exist(net_def_file,'file')>0);
net = caffe.Net(net_def_file, phase); % create net but not load weights

assert(iscell(net_weights_file))
num_weights_files = length(net_weights_file);
for i = 1:num_weights_files
    assert(exist(net_weights_file{i},'file')>0);
    net.copy_from(net_weights_file{i});
end

end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [qq, manips]=path_planning_SCO_SAWYER
close all;
global robot
global parameters
global hfigures

%STOMP PARAMETERS
%conversion from cost to Prob factor
%lambda is a constant to smooth the probability function
%with respect to obstacles and manipulability
parameters.lambda_obstacles = 2;
parameters.lambda_manip = 3;
%cost function starts at this distance
%must be below 0.3 for the 4 DOF robot
parameters.epsilon_distance = 0.1;
%height of the obstacle
%parameters.yo = 2.5;

%parameters.noise_sigma_null_space = 0.01;
parameters.alpha=0.05;
parameters.step_time=0.1;
%Error in XYZ to stop inverse kinematics
parameters.epsilonXYZ=0.01;
%Error in Quaternion to stop inverse kinematics.
parameters.epsilonQ=0.01;
parameters.stop_iterations=500;

%number of waypoints. More waypoints means greater precision in the
%trajectory
parameters.N = 15;
%number of particles. Set the number of hypotheses in the algorithm
parameters.K = 5;
%select the rows of J that should be taken into account when computing
%manipulability
% this considers the ellipsoid in linear speed and angular speed
%parameters.sel_J = [1 2 3 4 5 6];
% this one only considers the ellipsoid in linear speed and discards the
% angular speed component
parameters.sel_J = [1 2 3];
parameters.obstacles = [];

parameters.animate = 0;
close all
hfigures.hpaths = figure;
hfigures.hcosts = figure;
hfigures.hee = figure;
hfigures.htheta = figure;
hfigures.hbest_costs = figure;
hfigures.hdtheta = figure;

%LINE 1
x1 = -0.5;
y1 = 0.3; %m
z1 = 0;
x2 = 0;
y2 = 0.4; %m
z2 = 0;
p0 = [x1 y1 z1]';
pf = [x2 y2 z2]';
p3 = [0 0 0]';

%normal
n = compute_plane([p3 pf p0]);
parameters.obstacles{1}.n = n;
parameters.obstacles{1}.p = p0;
%points defining surface. %define obstacles using point and normal
parameters.obstacles{1}.points = [p3 pf p0];
%define trajectory from plane
parameters.trajectory{1}.line = [p0 pf];
parameters.trajectory{1}.T0 = build_T_from_obstacle(parameters.obstacles{1}, p0);

%LINE 2
x1 = 0;
y1 = 0.4; %m
z1 = 0;
x2 = 0.4;
y2 = 0.4; %m
z2 = 0.4;
p0 = [x1 y1 z1]';
pf = [x2 y2 z2]';
p3 = [0 0 0]';
%normal
n = compute_plane([p3 pf p0]);
parameters.obstacles{2}.n = n;
parameters.obstacles{2}.p = p0;
parameters.obstacles{2}.points = [p3 pf p0];
%define trajectory from plane
parameters.trajectory{2}.line = [p0 pf];
parameters.trajectory{2}.T0 = build_T_from_obstacle(parameters.obstacles{2}, p0);

%
% %LINE 3
x1 = 0.4;
y1 = 0.4; %m
z1 = 0.4;
x2 = 0.4;
y2 = 0.0; %m
z2 = 0.0; %m
p0 = [x1 y1 z1]';
pf = [x2 y2 z2]';
p3 = [0 0 0]';
%normal
n = compute_plane([p3 pf p0]);
parameters.obstacles{3}.n = n;
parameters.obstacles{3}.p = p0;
parameters.obstacles{3}.points = [p3 pf p0];
%define trajectory from plane
parameters.trajectory{3}.line = [p0 pf];
parameters.trajectory{3}.T0 = build_T_from_obstacle(parameters.obstacles{3}, p0);

%LAUNCH SCO given the stored parameters
pk = SCO_null_space(robot);

repo: ARTE


text:
function pass = test_paramODE_nonlin_US(pref)
% Test solving a parameter dependent ODE.
% Nick Hale, August 2011

% NOTE: Taken from V4 test chebop_paramODE.

if ( nargin == 0 )
    pref = cheboppref();
end
tol = 1e3*pref.bvpTol;

% Nonlinear parameter dependent problem (ULTRAS)
pref.discretization = @ultraS;

% Natural setup
x = chebfun('x', [-1 1]);
N = chebop(@(x, u, a) (1-x.^2).*u + .1*diff(u,2) + a.*exp(u));
% N = chebop(@(x, u, a) x.*u + .001*diff(u,2) + a.*diff(u));
N.lbc = @(u, a) [u + a + 1 ; diff(u)];
N.rbc = @(u, a) u - 1;
N.init = [x ; -1];
u = mldivide(N, chebfun(0), pref);
res = N(x, u);
Nlbc = N.lbc(u{1},u{2});
Nrbc = N.rbc(u{1},u{2});
err = norm(res) + Nlbc{1}(-1) + Nlbc{2}(-1) + Nrbc(1);

pass = err < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ n_data, x, fx ] = struve_h0_values ( n_data )

%*****************************************************************************80
%
%% STRUVE_H0_VALUES returns some values of the Struve H0 function.
%
%  Discussion:
%
%    The function is defined by:
%
%      HO(x) = 2/pi * Integral ( 0 <= t <= pi/2 ) sin ( x * cos ( t ) ) dt
%
%    In Mathematica, the function can be evaluated by:
%
%      StruveH[0,x]
%
%    The data was reported by McLeod.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%    special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
      0.12433974658847434366E-02, ...
     -0.49735582423748415045E-02, ...
      0.39771469054536941564E-01, ...
     -0.15805246001653314198E+00, ...
      0.56865662704828795099E+00, ...
      0.66598399314899916605E+00, ...
      0.79085884950809589255E+00, ...
     -0.13501457342248639716E+00, ...
      0.20086479668164503137E+00, ...
     -0.11142097800261991552E+00, ...
     -0.17026804865989885869E+00, ...
     -0.13544931808186467594E+00, ...
      0.94393698081323450897E-01, ...
     -0.10182482016001510271E+00, ...
      0.96098421554162110012E-01, ...
     -0.85337674826118998952E-01, ...
     -0.76882290637052720045E-01, ...
      0.47663833591418256339E-01, ...
     -0.70878751689647343204E-01, ...
      0.65752908073352785368E-01 ];

  x_vec = [ ...
        0.0019531250E+00, ...
       -0.0078125000E+00, ...
        0.0625000000E+00, ...
       -0.2500000000E+00, ...
        1.0000000000E+00, ...
        1.2500000000E+00, ...
        2.0000000000E+00, ...
       -4.0000000000E+00, ...
        7.5000000000E+00, ...
       11.0000000000E+00, ...
       11.5000000000E+00, ...
      -16.0000000000E+00, ...
       20.0000000000E+00, ...
       25.0000000000E+00, ...
      -30.0000000000E+00, ...
       50.0000000000E+00, ...
       75.0000000000E+00, ...
      -80.0000000000E+00, ...
      100.0000000000E+00, ...
     -125.0000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function aistats2012_run_comparison

%
% Theta
%

% Cholesky

nips2011_comparison(3, 1, true, false);
nips2011_comparison(3.5, 1, true, false);
nips2011_comparison(4, 1, true, false);
nips2011_comparison(4.5, 1, true, false);
nips2011_comparison(5, 1, true, false);

% Kronecker

% $$$ nips2011_comparison(3, 1, true, true);
% $$$ nips2011_comparison(3.5, 1, true, true);
% $$$ nips2011_comparison(4, 1, true, true);
% $$$ nips2011_comparison(4.5, 1, true, true);
% $$$ nips2011_comparison(5, 1, true, true);
% $$$ nips2011_comparison(6, 1, true, true);
nips2011_comparison(7, 1, true, true);

%
% F
%

% Kronecker

% $$$ nips2011_comparison(3, 1, false, true);
% $$$ nips2011_comparison(3.5, 1, false, true);
% $$$ nips2011_comparison(4, 1, false, true);
% $$$ nips2011_comparison(4.5, 1, false, true);
% $$$ nips2011_comparison(5, 1, false, true);
% $$$ nips2011_comparison(5.5, 1, false, true);
% $$$ nips2011_comparison(6, 1, false, true);
% $$$ nips2011_comparison(7, 1, false, true);
% $$$ nips2011_comparison(8, 1, false, true);

% Cholesky

nips2011_comparison(3, 1, false, false);
nips2011_comparison(3.5, 1, false, false);
nips2011_comparison(4, 1, false, false);
nips2011_comparison(4.5, 1, false, false);
nips2011_comparison(5, 1, false, false);
nips2011_comparison(5.5, 1, false, false);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [varargout] = istable(varargin)

% ISTABLE is a drop-in replacement for the same function that was
% introduced in MATLAB R2013b.
%
% In all MATLAB versions prior to 2013b this function returns false.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% see https://github.com/fieldtrip/fieldtrip/issues/899

alternatives = which(mfilename, '-all');
if ~iscell(alternatives)
  % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
  alternatives = {alternatives};
end

keep = true(size(alternatives));
for i=1:numel(alternatives)
  keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
  keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
end
alternatives = alternatives(keep);

if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
  % remove this directory from the path
  p = fileparts(mfilename('fullpath'));
  sw = warning('off','backtrace');
  warning('Removing "%s" from your path.\nSee http://www.fieldtriptoolbox.org/faq/should_i_add_fieldtrip_with_all_subdirectories_to_my_matlab_path/', p);
  warning(sw);
  rmpath(p);
  % call the original MATLAB function
  if exist(mfilename, 'builtin')
    [varargout{1:nargout}] = builtin(mfilename, varargin{:});
  else
    [varargout{1:nargout}] = feval(mfilename, varargin{:});
  end
  return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = istable(input)

tf = false;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};

author: spm


text:
function [modelUpdated,modelPruned,Ex_Rxns] = pruneModel(model,minGrowth, biomassRxn)
% This function prunes a model to its most compact subnetwork given some model
% constraints and a minimal growth rate by identifying and removing all blocked reactions.
%
% USAGE:
%
%    [modelUpdated, modelPruned, Ex_Rxns] = pruneModel(model, minGrowth, biomassRxn)
%
% INPUTS:
%    model:         model structure
%    minGrowth:     minimal Growth rate to be set on biomass reaction
%    biomassRxn:    biomass reaction name (default: 'biomass_reaction2')
%
% OUTPUTS:
%    modelUpdated:  same as input model but constraints on blocked reactions
%                   are set to be 0
%    modelPruned:   pruned model, where all blocked reactions are removed
%                   (attention this seems to cause issues with GPRs)
%    Ex_Rxns:       List of exchange reactions in pruned model
%
% .. Authors: - Ines Thiele, 02/2014
%             - Modified by Loic Marx, December 2018

if ~exist('biomassRxn','var')
    biomassRxn = checkObjective(model);
    biomassRxn = biomassRxn{1};
end

modelUpdated = model;
modelForPruning = model;
modelForPruning.lb(modelForPruning.lb < 0) = -1000;
modelForPruning.ub(modelForPruning.ub < 0) = 0; % in case  uptake is enforced
modelForPruning.ub(modelForPruning.ub > 0) = 1000;
modelForPruning.lb(modelForPruning.lb > 0) = 0; % in case  secretion is enforced
%set back biomass constraint
modelForPruning.lb(find(ismember(modelForPruning.rxns,biomassRxn)))=minGrowth;
epsilon =getCobraSolverParams('LP', 'feasTol')*100;
[modelForPruningPruned, BlockedRxns] = identifyBlockedRxns(modelForPruning,epsilon);
cnt =1;
for t=1:length(modelForPruningPruned.rxns)
    if  strfind(modelForPruningPruned.rxns{t}, 'EX_')
        Ex_Rxns(cnt,1) =modelForPruningPruned.rxns(t); %make exchange reaction list
        cnt=cnt+1;
    end
end

modelUpdated.lb(ismember(model.rxns,BlockedRxns.allRxns))=0;
modelUpdated.ub(ismember(model.rxns,BlockedRxns.allRxns))=0;

RxnsInModelMin = setdiff(model.rxns,BlockedRxns.allRxns); % all reactions that are not blocked
modelPruned = extractSubNetwork(model,RxnsInModelMin);

author: opencobra

repo: cobratoolbox


text:
function h = match_plot(img1,img2,points1,points2)
%MATCH_PLOT
%   h = match_plot(img1,img2,points1,points2)
%   plot coresponding points between two images
%
%   points1 = [x1 y1; x2 y2; ...] = coresponding points in img1
%   points2 = [x1 y1; x2 y2; ...] = coresponding points in img2
%
%   origin is the top left of the image
%   x axis pointing right, y axis pointing down.
%   points 2 has to be the same size as points 1
%
%   if 2 images have different size, the smaller one is rescaled.
%
%   coresponding lines are plot in different colors, from red to blue in
%   rainbow order.
%
%   returns the figure handle object h
%
%   written by Gooly 4/21/2011
%   http://computervisionblog.wordpress.com/

h = figure;
colormap = {'b','r','m','y','g','c'};
height = max(size(img1,1),size(img2,1));
img1_ratio = height/size(img1,1);
img2_ratio = height/size(img2,1);
img1 = imresize(img1,img1_ratio);
img2 = imresize(img2,img2_ratio);
points1 = points1 * img1_ratio;
points2 = points2 * img2_ratio;
points1 = [points1(:,2) points1(:,1)];
points2 = [points2(:,2) points2(:,1)];
match_img = zeros(height, size(img1,2)+size(img2,2), size(img2,3));
match_img(1:size(img1,1),1:size(img1,2),:) = img1;
match_img(1:size(img2,1),size(img1,2)+1:end,:) = img2;
imshow(uint8(match_img));
hold on;
for i=1:size(points1,1)
    plot([points1(i,2) points2(i,2)+size(img1,2)],[points1(i,1) points2(i,1)],colormap{mod(i,6)+1});
end

hold off;

end

author: Sable

repo: mcbench-benchmarks


text:
% demos for ch03
clear; close all;
d = 1;
n = 200;
[x,t] = linRnd(d,n);
%% Empirical Bayesian linear regression via EM
[model,llh] = linRegEm(x,t);
plot(llh);
[y,sigma] = linRegPred(model,x,t);
figure
plotCurveBar(x,y,sigma);
hold on;
plot(x,t,'o');
hold off;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function val = getHomer2_paramValue(funcName, paramName, file)

val = [];

DEBUG=0;
if DEBUG
    fprintf('Searching Homer2 output %s for %s: %s value\n', file.name, funcName, paramName);
end
load(file.name);
for jj=1:length(group.subjs)
    for kk=1:length(group.subjs(jj).runs)
        for iFunc=1:length(group.subjs(jj).runs)
            if strcmp(group.subjs(jj).runs(kk).procInput.procFunc.funcName{iFunc}, funcName)
                for iParam=1:group.subjs(jj).runs(kk).procInput.procFunc.nFuncParam(iFunc)
                    if strcmp(group.subjs(jj).runs(kk).procInput.procFunc.funcParam{iFunc}{iParam}, paramName)
                        val = group.subjs(jj).runs(kk).procInput.procFunc.funcParamVal{iFunc}{iParam};
                        break
                    end
                end
            end
            if ~isempty(val)
                break;
            end
        end
        if ~isempty(val)
            break;
        end
    end
    if ~isempty(val)
        break;
    end
end

author: BUNPC


text:
classdef InstallDriverLinux < adi.InstallDriverBase

%   Copyright 2017 The MathWorks, Inc.

  methods
    function obj = InstallDriverLinux(varargin)
      obj@adi.InstallDriverBase(varargin{:});

      obj.Title.Text = 'ADI Board Support Package Installer';

      obj.InfoText1.Position = [20 340 400 30];
      obj.InfoText1.Text = 'Test1';%message('plutoradio:hwsetup:InstallDriverLinux_InfoText1').getString;

      obj.InfoText2.Position = [20 290 400 50];
      obj.InfoText2.Text = 'Test2';%message('plutoradio:hwsetup:InstallDriverLinux_InfoText2').getString;

      obj.HelpText.AboutSelection = '';
      obj.HelpText.WhatToConsider = 'Test3';%message('plutoradio:hwsetup:InstallDriverLinux_WhatToConsider').getString;
    end

    function out = getNextScreenID(obj)
      % HSA infrastructure requires a "getNextScreenID" method implemented
      % in the leaf class. Otherwise, it renders "Finish" button instead of
      % "Next".
      out = getNextScreenIDImpl(obj);
    end
  end

  methods (Access = protected)
    function [status,out] = installDriverImpl(obj)
      udevFileName = '90-plutosdr-mw.rules';
      udevrules = fullfile(tempdir, udevFileName);
      fid = fopen(udevrules, 'wt');
      fprintf(fid, '%s\n%s\n', ...
        '# PlutoSDR', ...
        ['SUBSYSTEMS=="usb", ATTRS{idVendor}=="0456", '...
        'ATTRS{idProduct}=="b673", MODE:="0666"']);
      fclose(fid);

      commandwindow
      [status,out] = run(obj.Workflow.HardwareInterface, 'sudo echo', '-echo');

      if ~status
        [status,out] = run(obj.Workflow.HardwareInterface, ['sudo cp ' udevrules ' /lib/udev/rules.d']);
        if ~status
          [status,out] = run(obj.Workflow.HardwareInterface, 'sudo udevadm control --reload');
          if ~status
            [status,out] = run(obj.Workflow.HardwareInterface, 'sudo udevadm trigger');
          end
        end
      end
      run(obj.Workflow.HardwareInterface, sprintf('rm %s', udevrules));

      if ~status
        obj.Workflow.InstalledUdevRule = ...
          fullfile('/lib/udev/rules.d', udevFileName);
      end
    end
  end
end

author: analogdevicesinc


text:
function [ax_mps2] = calcPathAx(s_m, v_mps)
%_________________________________________________________________
%% Documentation
%
% Authors:      Alexander Wischnewski (alexander.wischnewski@tum.de)
%
% Start Date:   08.02.2018
%
% Description:  calculates the acceleration profile corresponding to a
%               velocity profile for a given path length. Constant
%               acceleration is assumed between the discretization points.
%
% Inputs:
%   s_m         Vector with path length values
%   v_mps       Vector with velocity at these path points
%
% Outputs:
%   ax_mps2     Vector with accelerations valid from the current to the
%               next point.

% initialize variables
ax_mps2 = zeros(length(v_mps), 1);
% calculate acceleration based on the assumption that it is constant
% between two discretization points. Last point can't be calculated as no
% further velocity information is available for the point behind.
dS = diff(s_m);
dV = diff(v_mps);
ax_mps2(1:(end-1)) = (dV.^2 + 2.*dV.*v_mps(1:(end-1)))./(2.*dS);
% copy last point which could be calculated, as no better information is
% available
ax_mps2(end) = ax_mps2(end-1);

author: TUMFTM


text:
% DEMO_VBRFA - A simple comparison of VB PCA and two robust extensions.
%
% The algorithms do not use ARD prior for the loadings W in order to keep
% the visualizations of the subspaces more clear. In addition, it helps
% avoiding pruning out relevant components too easily.
%
% NOTE: The robust algorithms might be less sensitive to the initialization
% if they estimate many components. With only one estimated component, they
% might kill all the components.
%
% This demo doesn't always show good results for robust algorithms. You may
% want to run a few times. In addition, sometimes there are several good
% interpretations of the data and the robust algorithms pick one.

% Last modified 2010-06-10
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function Q = demo_vbrfa(seed)

n = 50;
m = 2;
d = 1;
dh = d + 0; % number of components to estimate

if nargin >= 1
  randn('state', seed);
  rand('state', seed);
end

mu = [10;20];
W = [2; 1];
X = orth(randn(d,n)')' * sqrt(n);
s2 = 0.1;

% Generate normal data
Y = W*X + repmat(mu,1,n);

% Add noise
Yn = Y + sqrt(s2)*randn(m,n);

% Generate some outliers
Yno = Yn;
p = (rand(m,n) < 0.05);
Yno(p) = Yno(p) + (5 + 8 * rand(sum(p(:)),1));

% Generate some missing values
Ynom = Yno;
pmv = (rand(m,n) < 0.0);
Ynom(pmv) = NaN;

% Run different algorithms
options.init.tau = 1e1*ones(m,1); % the initialization of this can be crucial..
options.init.nu = 1;
options.prior.a_alpha = 1e10; % "fix" w to non-informative for W by
options.prior.b_alpha = 1e15; % setting strong prior for w
options.update_nu = 1;
options.update_alpha = 1;
options.update_beta = 1;
options.rotate = false;
options.common_nu = false;
options.common_tau = true;
options.maxiter = 100;

% PCA
disp('Run PCA')
options.robustness = 'none';
results_p = vbrfa(Ynom, dh, options);
W_p = results_p.W;
X_p = results_p.X;
Mu_p = results_p.Mu;
nu_p = results_p.nu;

% Robust PCA with multivariate Student-t
disp('Run robust PCA with multivariate Student-t')
options.robustness = 'multivariate-t';
results_rp = vbrfa(Ynom, dh, options);
W_rp = results_rp.W;
X_rp = results_rp.X;
Mu_rp = results_rp.Mu;
nu_rp = results_rp.nu;

options.maxiter = 1000;

% Robust PCA with independent Student-t
disp('Run robust PCA with independent Student-t')
options.robustness = 'independent-t';
results_rvb = vbrfa(Ynom, dh, options);
W_rvb = results_rvb.W;
X_rvb = results_rvb.X;
Mu_rvb = results_rvb.Mu;
nu_rvb = results_rvb.nu;

% Reconstruct
Y_p = bsxfun(@plus, W_p*X_p, Mu_p);
Y_rp = bsxfun(@plus, W_rp*X_rp, Mu_rp);
Y_rvb = bsxfun(@plus, W_rvb*X_rvb, Mu_rvb);

plot_results(Ynom, []);
plot_results(Ynom, Y_p)
plot_results(Ynom, Y_rp)
plot_results(Ynom, Y_rvb)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
function plot_results(Y, Yh)
%subspace2d(Yh, [], Y)
figure
plot(Y(1,:), Y(2,:), 'x', 'markersize', 5, 'color', [.2 .2 .2]);
if ~isempty(Yh)
  hold on
  plot(Yh(1,:), Yh(2,:), 'ko', 'markersize', 3);
  nans = nan*ones(1,size(Y,2));
  plot([Y(1,:);Yh(1,:);nans], [Y(2,:);Yh(2,:);nans], ':', 'color', [.0 .0 .0]);
end
set(gca, 'DataAspectRatioMode', 'manual');
set(gca, 'DataAspectRatio', [1 1 1]);
set(gca, 'XTick', [], 'YTick', []);
set(gcf, 'units', 'centimeters');
% $$$ pos = get(gcf, 'position');
% $$$ set(gcf, 'position', [pos(1), pos(2), 5 5])
% $$$ set(gcf, 'PaperPositionMode', 'auto', 'paperunits', 'centimeters', 'PaperSize', [5 5]);

xmg = 0.05 * (max(Y(1,:)) - min(Y(1,:)));
ymg = 0.05 * (max(Y(2,:)) - min(Y(2,:)));
xl = [min(Y(1,:))-xmg, max(Y(1,:))+xmg];
yl = [min(Y(2,:))-ymg, max(Y(2,:))+ymg];
set(gca, 'xlim', xl, 'ylim', yl);

% $$$ % Mark outliers
% $$$ hold on
% $$$ indeces = sum(p,1)>0;
% $$$ plot(VYno(1,indeces), VYno(2, indeces), 'go', 'MarkerSize', 8);

% $$$ % Mark observations with missing values
% $$$ indeces = sum(pmv,1)>0;
% $$$ scatter(VYno(1,indeces), VYno(2, indeces), 'yo');

return


text:
classdef nonEllipsoid < referenceEllipsoid
    % used when coordinate system is cartesian
    methods
        function obj = nonEllipsoid(varargin)
            obj@referenceEllipsoid();
            p = inputParser;
            p.addOptional('Name','XYZ');
            p.addOptional('LengthUnit','meter');
            p.parse(varargin{:});
            obj.LengthUnit = p.Results.LengthUnit;
            obj.Name = p.Results.Name;
        end
    end

end

author: CelsoReyes


text:
function ans = calc_minDose(doseBinsMidPtsV, volsHistV, volumeType)
% Calculates the minimum dose for a given DVH
% The last argument 'volumeType' is ignored in this instance
%
%  MODIFICATION ALERT:  THIS FUNCTION IS UTILIZED BY THE DREXLER CODEBASE
%
%  Last modified: AJH 11/05
%
%  Usage:  calc_minDose(doseBinsV, volsHistV)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

ind = find(volsHistV~=0,1,'first');
ans = doseBinsMidPtsV(ind);

return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
 function fit = de_ftab_fit(sl, fm, varargin)
%function fit = de_ftab_fit(sl, fm, [options])
%|
%| Fit to DE table fm(), suitable for subsequent interpolation / extrapolation.
%| Uses either classic polynomial basis functions,
%| or an exponential model: -log(sum_k p_k exp(-m_k . s))
%|
%| in
%|	sl	{L}		sample locations for each of L materials
%|	fm	[s1 ... sL M]	DE tables f_m(s_1,...,s_L), size [(Ns) M]
%|
%| option
%|	'type'	'poly' | 'exp'	type of fit (default 'poly') todo: change?
%|	'wt'	{M}		fit weighting for each of M energies.
%|				(default depends on 'type'; see code)
%|	'show'	1|0		plot? (default false)
%|	'ntype' ''		options for negative sl values (default '')
%|				'abs' - experimental symmetrizing option
%| options for 'poly'
%|	'order'			polynomial order (default 3)
%|	'maxdegree'		argument for poly_string() (default [])
%| options for 'exp'
%|	'kev'	[Ne 1]		kev *for fitting* (default [10:5:200]')
%|	'mac'	[Ne L]		*for fitting* (this or mtype are required)
%|	'mtype' {L}		material types
%|	'macbar' [M L]		if nonempty, use it to match derivative at s=0.
%|	'wls_simplex_arg' {}	arguments for wls_simplex (default {})
%|
%| out
%|	fit	strum		strum object for fitted f_m(s_1,...,s_L)
%|				fit.coef is [nbasis M] for 'poly'
%|	methods:
%|				(sll denotes stacked array [(Ns) L] or [*Ns L])
%|	.fmfun(sll [() L])	fm function evaluation [() L] -> [() M]
%|	.fgrad(sll)		fm gradient evaluation [() L] -> [() L M]
%|	.fhess(sll)		fm hessian evaluation [() L] -> [() L L M]
%|	.show_sp([en, sp])	plot true spectrum vs fitted spectrum
%|	.show_fm(sl, [fm])	mesh plot of fm and its fit
%|	.show_err(sl, fm)	mesh plot of fit error
%|	.mac_eff		effective mass atten coef based on fit
%|				(valid for 'exp' only)
%|
%| Copyright 2006-3-3, Jeff Fessler, University of Michigan

if nargin == 1 && streq(sl, 'test'), de_ftab_fit_test, return, end
if nargin < 2, ir_usage, end

arg.show = false;
arg.type = '';
arg.mtype = {};
arg.ntype = ''; % how to handle negatives
arg.kev = [10:5:200]';
arg.mac = [];
arg.macbar = [];
arg.wls_simplex_arg = {}; % pass to wls_simplex which passes it to lsqlin
arg.order = 3;
arg.maxdegree = [];
arg.wt = {};
arg.dc = false; % default: exclude dc (constant) from (polynomial) fit
arg = vararg_pair(arg, varargin);
if isempty(arg.type), arg.type = 'poly'; end

LL = length(sl);

if ndims(fm) == LL
	MM = 1;
elseif ndims(fm) == LL + 1
	MM = size(fm,ndims(fm)); % last dimension is number of spectra
else
	error 'invalid fm size'
end

for ll=1:LL % make sure arguments match table dimensions
	if length(sl{ll}) ~= size(fm, ll)
		fail('dim mismatch %d', ll)
	end
end

switch arg.type
case 'poly'
	[fit fmfun fgrad fhess] = de_ftab_fit_poly(sl, fm, MM, arg.wt, ...
		arg.order, arg.maxdegree, arg.dc);

case 'exp'
	fit = de_ftab_fit_exp(sl, fm, MM, arg.wt, ...
		arg.mtype, arg.mac, arg.kev, arg.macbar, arg.wls_simplex_arg);
	fmfun = @de_ftab_fit_exp_eval;
	fgrad = @de_ftab_fit_exp_grad;
	fhess = @de_ftab_fit_exp_hess;
	fit.mac_eff = zeros(MM, LL); % [M L] mac effective
	for mm=1:MM
		fit.mac_eff(mm,:) = fit.mac{mm}' * fit.coef{mm};
	end

otherwise
	fail('unknown fit type %s', arg.type)
end

fit.LL = LL;
fit.MM = MM;
fit.type = arg.type;
fit.ntype = arg.ntype;

meth = {'fmfun', fmfun, '(sll [() L]) -> [() M]'; ...
	'fgrad', fgrad, '(sll [() L]) -> [() L M]'; ...
	'fhess', fhess, '(sll [() L]) -> [() L L M]'; ...
	'show_err', @de_ftab_fit_show_err, '(sl, fm)'; ...
	'show_fm', @de_ftab_fit_show_fm, '(sl, [fm])'; ...
	'show_sp', @de_ftab_fit_show_sp, '([en, sp])'; ...
	};
fit = strum(fit, meth);

if arg.show && im
	fit.show_fm(sl, fm);
end

end % de_ftab_fit()

%
% de_ftab_fit_poly()
%
function [fit, ffun, fgrad, fhess] = ...
	 de_ftab_fit_poly(sl, fm, MM, wt, order, maxdegree, dc)

sll = ndgrid_jf('mat', sl{:});

% for fitting, up weight the no-bone part
% because soft tissue is most prevalant
if isempty(wt)
	wt{1} = 1 ./ (stackpick(sll,2) + 1);
	wt{2} = wt{1};
end

fit.basis_order = order;
fit.basis_maxdegree = maxdegree;
fit.basis_dc = dc;

[fit.basis_func fit.basis_d1 fit.basis_d2] = ...
	ir_poly2_fun(order, 'maxdegree', maxdegree, 'dc', dc);

fit.nbasis = numel(fit.basis_func(0,0));

ss1 = col(stackpick(sll,1));
ss2 = col(stackpick(sll,2));
tmp = fit.basis_func(ss1, ss2); % [*Ns nbasis]
fit.coef = zeros(fit.nbasis, MM);
for mm=1:MM
	fit.coef(:,mm) = (diag_sp(wt{mm}(:)) * tmp) ...
		\ (wt{mm}(:) .* col(fm(:,:,mm)));
end

ffun = @(fit,sl) ...
	reshape([fit.basis_func(sl{1}(:), sl{2}(:)) * fit.coef(:,1); ...
		fit.basis_func(sl{1}(:), sl{2}(:)) * fit.coef(:,2)], ...
		[size(sl{1}) 2]);
ffun = @de_ftab_fit_poly_eval;
fgrad = @() error('not done');
fhess = @() error('not done');

end % de_ftab_fit_poly()

%
% de_ftab_fit_poly_eval()
%
function fm = de_ftab_fit_poly_eval(fit, sll)

% if fit.MM ~= 2
s1 = stackpick(sll,1);
s2 = stackpick(sll,2);
fm = [	fit.basis_func(s1(:), s2(:)) * fit.coef(:,1);
	fit.basis_func(s1(:), s2(:)) * fit.coef(:,2)];
fm = reshape(fm, [size(s1) 2]);

end % de_ftab_fit_poly_eval()

%
% de_ftab_fit_exp()
% todo: more documentation
% todo: use a set of exponents, not "mac"
%
function fit = de_ftab_fit_exp(sl, fm, MM, wt, mtype, mac, kev, macbar, ...
	wls_simplex_arg)

sll = ndgrid_jf('mat', sl{:}); % [(Ns)]

if isempty(mac)
	if isempty(mtype), error 'mac or mtype required', end
	mac = xray_read_atten(mtype, kev); % [Ne L]
else
	if size(mac,1) ~= length(kev), fail 'size mismatch', end
end

LL = length(sl);

if isempty(wt), wt = num2cell(ones(MM,1)); end

Ab = exp(-reshapee(sll, [], LL) * mac'); % [*Ns Ne] "over-complete" basis

fit.kev = cell(1,MM);
fit.mac = cell(1,MM);
fit.coef = cell(1,MM);
for mm=1:MM
	warg = wls_simplex_arg;
	if LL == 1 && ~isempty(macbar) % constrain deriv. at 0
		warg = {wls_simplex_arg{:}, 'inprodv', mac' / macbar(mm,1)};
	end

	if MM == 1 % todo: kludgy
		dat = fm;
	else
		dat = stackpick(fm,mm); % [(Ns)]
	end
	y = exp(-dat);
	Wh = spdiag(sqrt(wt{mm}(:)), 'nowarn');
	% initial coefficients for each candidate energy
	x = wls_simplex(Ab, y(:), Wh, [], warg{:}); % [Ne 1]

	ie = x > 1e-6; % find key energies
	fit.kev{mm} = kev(ie);
	fit.mac{mm} = mac(ie,:); % [Ne L] (now Ne may be smaller than before)

	warg = wls_simplex_arg;
	if LL == 1 && ~isempty(macbar) % constrain derivative at 0
		warg = {wls_simplex_arg{:}, 'inprodv', fit.mac{mm}' / macbar(mm,1)};
	end

	A = exp(-reshapee(sll, [], LL) * fit.mac{mm}'); % [*Ns Ne] final basis
	% final coefficients at key enerties:
	fit.coef{mm} = wls_simplex(A, y(:), Wh, [], warg{:}); % [Ne 1]
end

end % de_ftab_fit_exp()

%
% de_ftab_fit_exp_eval()
% evaluate
% in
%	sll	[(Ns) L]	stackup of s1,s2,...,s_L
% out
%	f	[(Ns) M]	stackup of f1,f2,...,f_M
%
function f = de_ftab_fit_exp_eval(fit, sll)
LL = fit.LL;
if LL == 1
	Ns = size(sll);
	if Ns(end) == 1, Ns = Ns(1:end-1); end
else
	Ns = size(sll); if LL ~= Ns(end), fail 'bug', end; Ns = Ns(1:end-1);
end
sll = reshapee(sll, [], LL); % [*Ns L]
MM = fit.MM;

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = 0;
end

f = zeros(prod(Ns),MM);
for mm=1:MM

	switch fit.ntype
	case '' % standard
		A = exp(-sll * fit.mac{mm}'); % [*Ns Ne]
		tmp = -log(A * fit.coef{mm}); % [*Ns 1]
		f(:,mm) = tmp;

	case 'abs' % trick for negatives
		if LL == 1
			A = exp(-abs(sll) * fit.mac{mm}'); % [*Ns Ne]
			tmp = -log(A * fit.coef{mm}); % [*Ns 1]
			f(:,mm) = -sign(sll) .* log(A * fit.coef{mm}); % [*Ns 1]
%			bad = sll < 0;
%			f(bad,mm) = sll(bad) * fit.mac_eff(mm,1);
		else
			if ~warned
				warned = 1;
				warn 'negative trick not done for L>1'
			end
		end

	otherwise
		fail('bad ntype: %s', fit.ntype)
	end

end
f = reshape(f, [Ns MM]);

end % de_ftab_fit_exp_eval()

%
% de_ftab_fit_exp_grad()
% evaluate gradient of f for each of the given s vectors.
% in
%	sll	[(Ns) L]	stackup of s1,s2,...,s_L
% out
%	g	[(Ns) L M]	stackup of gradients of f(s)
%
function g = de_ftab_fit_exp_grad(fit, sll)
LL = fit.LL;
Ns = size(sll);
if LL == 1
	if Ns(end) == 1, Ns = Ns(1:end-1); end
else
	if LL ~= Ns(end); error 'bug', end
	Ns = Ns(1:end-1);
end
sll = reshape(sll, [], LL); % [*Ns L]
MM = fit.MM;

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = 0;
end

g = zeros(prod(Ns), LL, MM);
for mm=1:fit.MM
	alf = fit.coef{mm}; % [Ne 1]
	mac = fit.mac{mm}; % [Ne L]
	Ne = length(alf); % # exponential terms in the fit

	switch fit.ntype
%	case '' % standard
	case 'abs'
		if any(sll(:) < 0) % trick for negatives
			if LL == 1
				sll = abs(sll);
%				bad = sll < 0;
%				g(bad,1,mm) = fit.mac_eff(mm,1);
			else
				if ~warned
					warned = 1;
					warn 'negative trick not done for L>1'
				end
			end
		end
	end

	A = exp(-sll * mac'); % [*Ns Ne]
	Q = A .* repmat(alf', [prod(Ns) 1]); % [*Ns Ne]
	denom = sum(Q,2); % [*Ns 1]
	Q = Q ./ repmat(denom, [1 Ne]);
	g(:,:,mm) = Q * mac;

end
g = reshape(g, [Ns LL MM]);

end % de_ftab_fit_exp_grad()

%
% de_ftab_fit_exp_hess()
% evaluate hessian of f for each of the given s vectors.
% in
%	sll	[(Ns) L]	stackup of s1,s2,...,s_L
% out
%	h	[(Ns) L L M]	stackup of hessians of f(s)
%
function h = de_ftab_fit_exp_hess(fit, sll)
Ns = size(sll); LL = Ns(end); Ns = Ns(1:end-1);
sll = reshape(sll, [], LL); % [*Ns L]
MM = fit.MM;

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = 0;
end

h = zeros(prod(Ns), LL, LL, MM);
for mm=1:fit.MM
	mac = fit.mac{mm}; % [Ne L]
	alf = fit.coef{mm}; % [Ne 1]
	Ne = length(alf); % # exponential terms in the fit

	switch fit.ntype
%	case '' % standard
	case 'abs'
		if any(sll(:) < 0) % trick for negatives
			if LL == 1
				sll_sign = sign(sll);
				sll = abs(sll);
%				bad = sll < 0;
%				hm(bad,1,1) = 0;
			else
				if ~warned
					warned = 1;
					warn 'negative trick not done for L>1'
				end
			end
		end
	end

	A = exp(-sll * fit.mac{mm}'); % [*Ns Ne]

	Q = A .* repmat(alf', [prod(Ns) 1]); % [*Ns Ne]
	denom = sum(Q,2); % [*Ns 1]
	Q = Q ./ repmat(denom, [1 Ne]);
	g = Q * mac; % [*Ns L] gradient

	% loop over LL because it is smaller than *Ns or Ne
	hm = zeros(prod(Ns), LL, LL);
	for l1=1:LL
		for l2=1:LL
			hm(:,l1,l2) = g(:,l1) .* g(:,l2) ...
				- Q * (mac(:,l1) .* mac(:,l2));
		end
	end

	switch fit.ntype
%	case '' % standard
	case 'abs'
		if any(sll(:) < 0) % trick for negatives
			if LL == 1
				hm(:,1,1) = hm(:,1,1) .* sll_sign;
%				bad = sll < 0;
%				hm(bad,1,1) = 0;
			end
		end
	end

	h(:,:,:,mm) = hm;

end
h = reshape(h, [Ns LL LL MM]);

end % de_ftab_fit_exp_hess()

%
% de_ftab_fit_show_sp()
% compare true spectra to fitted spectra
%
function out = de_ftab_fit_show_sp(fit, en, sp)
if nargin < 1, fail 'de_ftab_fit_show_sp(fit, en, sp)', end

if ~streq(fit.type, 'exp'), printm 'show_sp only done for exp', return, end
if im
	clf, pl = (fit.MM+1)*100 + 10 + 1;
	if nargin > 1
		subplot(pl)
		plot(en, sp * diag(1 ./ max(sp)))
		tmp = num2cell(char([1:fit.MM] + '0'));
		legend(tmp{:})
		axisx(minmax(en))
	end

	if isfield(fit, 'kev')
		for mm=1:fit.MM
			subplot(pl+mm)
			bar(fit.kev{mm}, fit.coef{mm})
			xtick(fit.kev{mm}(fit.coef{mm} > 0))
			axis tight
			if nargin > 1, axisx(minmax(en)), end
		end
	else
		warn 'kev unknown'
	end
end

if nargout, out = []; end

end % de_ftab_fit_show_sp()

%
% de_ftab_fit_show_fm()
% compare fit to sampled fm
%
function out = de_ftab_fit_show_fm(fit, sl, fm, sll)
if nargin < 2, error 'de_ftab_fit_show_fm(fit, sl, fm)', end

do_fm = nargin > 2;
if ~isvar('sll') || isempty(sll)
	sll = ndgrid_jf('mat', sl{:});
end

fh = fit.fmfun(sll);

if ~iscell(sl), fail 'sl must be cell', end
switch length(sl) % LL
case 1
	s1 = sl{1};
	im clf
	arg = {};
	leg = {};
	order = 'bgrcm';
	if do_fm
		for mm=1:fit.MM
			arg = {arg{:}, s1, fm(:,mm), [order(mm) '.']};
			leg = {leg{:}, sprintf('true m=%d', mm)};
		end
	end
	for mm=1:fit.MM
		arg = {arg{:}, sll, fh(:,mm), [order(mm) '-']};
		leg = {leg{:}, sprintf('fit m=%d', mm)};
	end
	plot(arg{:})
	xlabel 's1'
	ylabel 'fm(s1)'
	legend(leg{:}, 4)

case 2
	s1 = sl{1};
	s2 = sl{2};
	smax(1) = max(s1);
	smax(2) = max(s2);
	if do_fm
		fmax = max(fm(:));
	else
		fmax = max(fh(:));
	end

	ax = [0 smax(1) 0 smax(2) 0 fmax];
	cax = [0 fmax];

	im clf
	im('pl', 2, fit.MM)
	for mm=1:fit.MM
		if do_fm
		show(mm, s1, s2, fm(:,:,mm), ax, cax, sprintf('f_%d(s)', mm))
		end
		m0 = mm + fit.MM;
		show(m0, s1, s2, fh(:,:,mm), ax, cax, sprintf('f_%d fit', mm))
	end
%	% text(-60, 10, '[cm^2/g]')

case 3
	s1 = sl{1};
	s2 = sl{2};
	s3 = sl{3};
	smax(1) = max(s1);
	smax(2) = max(s2);
	smax(3) = max(s3);
	if do_fm
		fmax = max(fm(:));
	else
		fmax = max(fh(:));
	end

	ax = [0 smax(2) 0 smax(3) 0 fmax];
	cax = [0 fmax];

	im clf
	im('pl', 2, fit.MM)
	for mm=1:fit.MM
		if do_fm
		show(mm, s2, s3, fm(1,:,:,mm), ax, cax, sprintf('f_%d(s)', mm))
		end
		m0 = mm + fit.MM;
		show(m0, s2, s3, fh(1,:,:,mm), ax, cax, sprintf('f_%d fit', mm))
	end
%	% text(-60, 10, '[cm^2/g]')

otherwise
	fail 'not done'
end

if nargout, out = []; end

end % de_ftab_fit_show_fm()

%
% de_ftab_fit_show_err()
% show fit errors, and relate to HU
% f = mac s, mac(70kev,H2O) = 0.2 cm^2 / g = 1000 HU
% mh = mac * 1000 HU / (0.2 g/cm^2)
% so Df/Dmh = Df/Dmac * Dmac/Dmh = (50 g/cm^2) (0.2 cm^2/g) / 1000 HU = 1/100 HU
% so Dmh = 100 HU * Df for 50 cm of H2O.
%
function out = de_ftab_fit_show_err(fit, sl, fm)
if nargin < 3, error 'de_ftab_fit_show_err(fit, sl, fm)', end

sll = ndgrid_jf('mat', sl{:});

fh = fit.fmfun(sll);
err = fh - fm;
printm('worst model error = %g of %g', max(abs(err(:))), max(fm(:)))
printm('worst model error = %g HU over 50cm H2O', 100*max(abs(err(:))))
for mm=1:fit.MM
	ee = stackpick(err,mm);
	printm('worst error (mm=%d): %g', mm, max(abs(ee(:))))
end

if max(abs(err(:))) == 0
	if nargout, out = []; end
return
end

if 0
	e1 = err(:,:,1);
	e2 = err(:,:,2);
	disp([minmax(e1); minmax(e2)]')
	printm('worst error1 %g', max(col(abs(err(:,1,:)))))
	printm('worst error2 %g', max(col(abs(err(:,2,:)))))
end
%err = abs(err);

switch length(sl) % LL
case 1
	s1 = sl{1};
	arg = {};
	leg = {};
	order = 'bgrcm';
	for mm=1:fit.MM
		arg = {arg{:}, s1, err(:,mm), [order(mm) '-']};
		leg = {leg{:}, sprintf('m=%d', mm)};
	end
	plot(arg{:})
	xlabel 's1'
	ylabel 'error'
	legend(leg{:}, 4)

case 2
	s1 = sl{1};
	s2 = sl{2};

	elim = minmax(err(:))';
	%elim = [-1 1] * 0.01; % +/- 1 HU
	ax = [0 max(s1) 0 max(s2) elim];
	im clf, im('pl', 1, fit.MM)
	for mm=1:fit.MM
		show(mm, s1, s2, err(:,:,mm), ax, elim, sprintf('f_%d error', mm))
	end

case 3
	s1 = sl{1};
	s2 = sl{2};
	s3 = sl{3};

	elim = minmax(err(:))';
	ax = [0 max(s2) 0 max(s3) elim];
	im clf, im('pl', 1, fit.MM)
	for mm=1:fit.MM
		show(mm, s2, s3, err(1,:,:,mm), ax, elim, sprintf('f_%d error', mm))
	end

otherwise
	fail 'not done'
end

if nargout, out = []; end

end % de_ftab_fit_show_err()

%
% show()
%
function show(pl, x, y, f, ax, cax, ti)
f = squeeze(f); % for LL=3 case
if ~im, return, end
im('subplot', pl)
if 1
	mesh(x,y,f')
	colormap hsv, caxis(cax), cbar
	axis(ax)
	xtick, ytick, ztick, zwhite
else
	im(x,y,f), cbar
	xtick, ytick
end
xlabel 's_1', ylabel 's_2', title(ti)

end % show()

%
% de_ftab_fit_test()
%
function de_ftab_fit_test

%stype = 'mono,70,100,160';
%stype = 'ps1';
%stype = 'poly1,60';
stype = 'poly1,80,100,160'; % stress test
xrs = xray_read_spectra(stype);

list.sl{1} = linspace(0, 50, 26); % coarse for fast fitting
list.sl{2} = linspace(0, 30, 31);
list.sl{3} = linspace(0, 2, 11);
list.mtype = {'water', 'bone', 'iodine'};
list.wls_simplex_arg = {{'reg', 1e-16}, {'reg', 1e-9}, {'reg', 1e-4}}; % more reg for LL>1 case

if 0 % look at derivatives
	LL = 2;
	sl = {list.sl{1:LL}};
	mtype = {list.mtype{1:LL}};
	warg = list.wls_simplex_arg{LL};
	mas = xray_read_mac(mtype);
	mac = mas.mac(xrs.en);
	sll = ndgrid_jf('mat', sl{:});
	fm = de_ftab_fm(sll, mac, xrs.Ide);
	fit = de_ftab_fit(sl, fm, 'show', 1, 'type', 'exp', 'mtype', mtype, ...
		'wls_simplex_arg', warg);
	g = fit.fgrad(sll);
end

%for LL=1:3
for LL=2
	sl = {list.sl{1:LL}};
	mtype = {list.mtype{1:LL}};
	warg = list.wls_simplex_arg{LL};

	mas = xray_read_mac(mtype);
	mac = mas.mac(xrs.en);
	if im
		clf, semilogy(xrs.en, mac), legend(mtype{:})
	end
	sll = ndgrid_jf('mat', sl{:});
	fm = de_ftab_fm(sll, mac, xrs.Ide);
%	fit = de_ftab_fit(sl, fm, 'show', 1, 'type', 'poly')
	fit = de_ftab_fit(sl, fm, 'show', 1, 'type', 'exp', 'mtype', mtype, ...
		'wls_simplex_arg', warg);

	fh = fit.fmfun(sll);
	g = fit.fgrad(sll);
	h = fit.fhess(sll);

	if 1 && LL == 2 % test gradient (for yong)
		if 1 % for gradient graphically
			list.slfine{1} = linspace(0, 90, 261); % fine for evaluating
			list.slfine{2} = linspace(0, 90, 511);
		else
			list.slfine{1} = linspace(0, 5, 261); % fine for evaluating
			list.slfine{2} = linspace(0, 3, 511);
		end
%		list.slfine{3} = linspace(0, 2, 201);
		slfine = {list.slfine{1:LL}};
		slf = ndgrid_jf('mat', slfine{:});
		fh = fit.fmfun(slf);
		gr = fit.fgrad(slf); % [261 511 L M]
		if im % examine gradient values graphically
			pl = @(mm) subplot(100 + 10 * fit.MM + mm);
			for mm=1:fit.MM
				pl(mm)
				tmp = gr(:,:,:,mm);
				plot(tmp(:,:,1), tmp(:,:,2), '.')
				titlef('%d kvp', xrs.kvp(mm))
				xlabel 'g1', ylabel 'g2'
				axis([0 max(col(tmp(:,:,1))) 0 max(col(tmp(:,:,2)))])
			end
		return
		end
		i1 = 5;
		i2 = 9;
		d1 = slfine{1}(i1+1) - slfine{1}(i1);
		d2 = slfine{2}(i2+1) - slfine{2}(i2);
		gh(:,1) = (fh(i1+1,i2,:) - fh(i1,i2,:)) / d1;
		gh(:,2) = (fh(i1,i2+1,:) - fh(i1,i2,:)) / d2;
		pr transpose(gh)
		gg = squeeze(gr(i1,i2,:,:));
		pr gg
		pr transpose(fit.mac_eff) % at origin
%		gg - fit.mac_eff'
	return
	end

	if im
		fit.show_sp(xrs.en, xrs.sp);
		prompt
		fit.show_fm(sl, fm);
		prompt
		fit.show_err(sl, fm);
		prompt
	end

end

end % de_ftab_fit_test()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt