text:
function IntProj(ana, hs, proj)
%
%  A function to perform Minimum or Maximum intensity projection of a volume of data
%  ana: Name of AFNI brick file
%  hs: Number of slices to use is each direction. Projection at slice i is from i-hs to i+hs
%      hs stands for half slab
%  proj: (-1) minimum intensity projection
%          1  maximum intensity projection
%  For output, the function creates
% see also insane function Int_Proj

FuncName = 'IntProj';

if (proj == -1) proj_s = 'min';
elseif (proj == 1) proj_s = 'max';
else
   fprintf(2,'Error %s: Bad projection parameter\n', FuncName);
   return;
end

[ana_pref, ana_view] = AfniPrefix(ana);
ana = sprintf('%s%s', ana_pref, ana_view);

[err, V, V_Info] = BrikLoad(ana);
Nvox = V_Info.DATASET_DIMENSIONS;

verb = 1;

%padd V by half slab in all directions to fix edge effects
Vp = zeros(Nvox(1)+3.*hs, Nvox(2)+3.*hs, Nvox(3)+3.*hs);
Vp(hs+1:Nvox(1)+hs, hs+1:Nvox(2)+hs, hs+1:Nvox(3)+hs) = V; clear V;

for (pd = 1:1:3), % projection direction
   if (pd == 1) dc = 'I';
   elseif (pd == 2) dc = 'J';
   elseif (pd == 3) dc = 'K';
   else
      fprintf(2,'Error %s: Bad projection direction\n', FuncName);
      return;
   end

   if (verb) fprintf(1,'%s: Processing direction %s\n', FuncName, dc); end
   %create the output set
   Vmi = zeros(Nvox(1), Nvox(2), Nvox(3));

   if (proj == -1),
      if (pd == 1),
         for (i=hs+1:1:Nvox(1)+hs),
            Vmi(i-hs,:,:) = min(Vp(i-hs:i+hs,hs+1:Nvox(2)+hs,hs+1:Nvox(3)+hs), [], pd);
         end
         Vmall = Vmi;
      elseif (pd == 2),
         for (i=hs+1:1:Nvox(2)+hs),
            Vmi(:,i-hs,:) = min(Vp(hs+1:Nvox(1)+hs,i-hs:i+hs,hs+1:Nvox(3)+hs), [], pd);
         end
         Vmall(:) = min(Vmi(:),Vmall(:));
      elseif (pd == 3),
         for (i=hs+1:1:Nvox(3)+hs),
            Vmi(:,:,i-hs) = min(Vp(hs+1:Nvox(1)+hs,hs+1:Nvox(2)+hs,i-hs:i+hs), [], pd);
         end
         Vmall(:) = min(Vmi(:),Vmall(:));
      end
   elseif (proj == 1),
      if (pd == 1),
         for (i=hs+1:1:Nvox(1)+hs),
            Vmi(i-hs,:,:) = max(Vp(i-hs:i+hs,hs+1:Nvox(2)+hs,hs+1:Nvox(3)+hs), [], pd);
         end
         Vmall = Vmi;
      elseif (pd == 2),
         for (i=hs+1:1:Nvox(2)+hs),
            Vmi(:,i-hs,:) = max(Vp(hs+1:Nvox(1)+hs,i-hs:i+hs,hs+1:Nvox(3)+hs), [], pd);
         end
         Vmall(:) = max(Vmi(:),Vmall(:));
      elseif (pd == 3),
         for (i=hs+1:1:Nvox(3)+hs),
            Vmi(:,:,i-hs) = max(Vp(hs+1:Nvox(1)+hs,hs+1:Nvox(2)+hs,i-hs:i+hs), [], pd);
         end
         Vmall(:) = max(Vmi(:),Vmall(:));
      end

   end
   OptW.Prefix = sprintf('%s_%s%g%s', ana_pref, proj_s, hs, dc);
   OptW.Scale = 1;
   WriteBrik(Vmi, V_Info, OptW);
end

   OptW.Prefix = sprintf('%s_%s%gIJK', ana_pref, proj_s, hs);
   OptW.Scale = 1;
   WriteBrik(Vmall, V_Info, OptW);

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test_suite=BatchExample_qmt_bssfp_test
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
    test_functions=localfunctions();
catch % no problem; early Matlab versions can use initTestSuite fine
end
initTestSuite;

function TestSetup
setenv('ISCITEST','1') % go faster! Fit only 2 voxels in FitData.m

function test_batch
curdir = pwd;

tmpDir = tempdir;
mkdir(tmpDir);
cd(tmpDir);

Modellist = {'b1_afi'};

for iModel = 1:length(Modellist)
    disp('===============================================================')
    disp(['Testing: ' Modellist{iModel} ' BATCH...'])
    disp('===============================================================')

    eval(['Model = ' Modellist{iModel}]);

    qMRgenBatch(Model,pwd)

    % Test if any dataset exist
    isdata = true;
    try
        Model.onlineData_url;
    catch
        isdata = false;
    end

    % Run Batch
    if isdata
        starttime = tic;
        eval([Modellist{iModel} '_batch'])
        toc(starttime)
    end
    close all
    cd ..

end
cd(curdir)

function TestTeardown
setenv('ISCITEST','0')

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%%
% test for resolution of analysis prior denoising
% min_x 1/2*|x-y|^2 + lambda*|K*x|_1
% min_x F(x) + G(K*x)   with   G=lambda*|.|_1,  F=1/2*|.-y|^2

addpath('../');
addpath('../toolbox/');

n = 100;
p = 200;
lambda = 1;
K = randn(p,n);
KS = K';

ProxG = @(u,tau)perform_soft_thresholding(u, lambda*tau);
ProxGS = compute_dual_prox(ProxG);
GradFS = @(x)x-y;
F = @(x)1/2*norm(x-y)^2;
G = @(u)lambda*norm(u,1);
options.FS = F;
options.GS = @(u)norm(u,'inf')/lambda;

L = norm(K)^2;
options.method = 'fista';
options.niter = 2000;
[x, fs, Constraint, Hu] = perform_fb_strongly(zeros(n,1), K, KS, GradFS, ProxGS, L, options);

%%
% Compute the primal error.

E = [];
for i=1:size(Hu,2)
    x1 = GradFS(KS*Hu(:,i));
    E(i) = F(x1) + G(K*x1);
end

sel = 1:round(options.niter/10);
clf;
loglog( E(sel)-min(E(:)) );
axis tight;

author: gpeyre

repo: matlab-toolboxes


text:
function[H,P,M1,M2,N1,N2] = countsig(data1,data2,T1,T2,parametric,p,quiet)
% Give the program two spike data sets and one
% or two time intervals and it will decide if
% the counts are significantly different.
% this is either with a non-parametric method
% or with a sqrt transformation followed by a
% t-test
% Usage: [H,P,M1,M2,N1,N2] = countsig(data1,data2,T1,T2,parametric,p,quiet)
%
% Input:
% Note that all times have to be consistent. If data
% is in seconds, so must be sig and t. If data is in
% samples, so must sig and t. The default is seconds.
%
% data1      - structure array of spike times (required)
% data2      - structure array of spike times (required)
% T1         - time interval (default all)
% T2         - time interval (default T1)
% parametric - 0 = non-parametric (Wilcoxon)
%            - 1 = ttest on sqrt of counts
%            - 2 = Poisson assumption
%              (default = 0)
% p          - significance level (0.05)
% quiet      - 1 = no display 0 = display
%
% Output:
%
% H          - 1 if different 0 if not
% P          - prob of result if same
% M1         - mean count for data1
% M2         - mean count for data2
% N1         - counts for data1
% N2         - counts for data2

if nargin < 2;error('I need 2 sets of spike data');end
data1=padNaN(data1); % create a zero padded data matrix from input structural array
data2=padNaN(data2); % create a zero padded data matrix from input structural array
data1=data1'; data2=data2'; % transpose data to get it into a form acceptable to Murray's routine
if nargin < 3,
   T1 = [min(data1(:,1)) max(max(data1))];
end
if nargin < 4,
   T2 = T1;
end
if nargin < 5,
   parametric = 0;
end
if nargin < 6; p = 0.05;end
if nargin < 7; quiet = 0; end

if isempty(T1),
   T1 = [min(data1(:,1)) max(max(data1))];
end
if isempty(T2)
   T2 = T1;
end
if isempty(parametric),
   parametric = 0;
end
if isempty(p)
   p = 0.05;
end
if isempty(quiet),
    quiet = 0;
end

NT1 = length(data1(:,1));
NT2 = length(data2(:,2));

if (NT1 < 4 || NT2 < 4) && parametric ~= 2,
  disp('Low number of trials : switch to Poisson test')
  parametric = 2;
end

if abs((T1(2)-T1(1)) - (T1(2)-T1(1))) > 10.^-6
  error('Time intervals for analysis are different')
end

% get counts...
N1=zeros(1,NT1);
for n=1:NT1
  N1(n) = length(find(data1(n,:) >=  T1(1) & ...
          data1(n,:) <=  T1(2) & ~isnan(data1(n,:))));
end
N2=zeros(1,NT2);
for n=1:NT2
  N2(n) = length(find(data2(n,:) >=  T2(1) & ...
          data2(n,:) <=  T2(2) & ~isnan(data1(n,:))));
end
M1 = mean(N1);
M2 = mean(N2);

% do non parametric test...

if parametric == 0
  [P H] = ranksum(N1,N2,p);
end

% parametric test (with stabilizing transform)...

%  use sqrt transformation from
%  Cox and Lewis to make data more Gaussian
%  the statistical analysis of series of events pg 44

if parametric == 1
  X = sqrt(N1 +0.25);
  Y = sqrt(N2 +0.25);
  [H,P] = ttest2(X,Y,p,0);
end

%  Poisson test.  Use method from Zar
%  pg 580 Ed 3.  Z bahaves as a normal variate under
%  null of same process mean and Poisson processes

if parametric == 2
  X = sum(N1);
  Y = sum(N2);
  Z = abs(X-Y)./sqrt(X+Y);
  P = 2*(1-normcdf(Z));
  if P < p; H = 1;else H = 0;end
end

if quiet == 0
  if H == 1
    disp('Counts are signifcantly different')
  else
    disp('Counts are not signifcantly different')
  end
  disp(['Mean count for data1 = ' num2str(M1)])
  disp(['Mean count for data2 = ' num2str(M2)])
  disp(['P value = ' num2str(P)])
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
clear

addpath(genpath('helpers/'));
find_SEMAINE;

out_loc = './AU_predictions/out_SEMAINE/';
if(~exist(out_loc, 'file'))
    mkdir(out_loc);
end

if(isunix)
    executable = '"../../build/bin/FeatureExtraction"';
else
    executable = '"../../x64/Release/FeatureExtraction.exe"';
end
%%
parfor f1=1:numel(devel_recs)

    if(isdir([SEMAINE_dir, devel_recs{f1}]))

        vid_file = dir([SEMAINE_dir, devel_recs{f1}, '/*.avi']);

        f1_dir = devel_recs{f1};

        curr_vid = [SEMAINE_dir, f1_dir, '/', vid_file.name];

        command = sprintf('%s -aus -f "%s" -out_dir "%s" ', executable, curr_vid, out_loc);

        if(isunix)
            unix(command, '-echo');
        else
            dos(command);
        end

    end
end

%% Actual model evaluation
[ labels, valid_ids, vid_ids, vid_names ] = extract_SEMAINE_labels(SEMAINE_dir, devel_recs, aus_SEMAINE);

labels_gt = cat(1, labels{:});

%% Identifying which column IDs correspond to which AU
tab = readtable([out_loc, vid_names{1}, '.csv']);
column_names = tab.Properties.VariableNames;

% As there are both classes and intensities list and evaluate both of them
aus_pred_class = [];
inds_class_in_file = [];

for c=1:numel(column_names)
    if(strfind(column_names{c}, '_c') > 0)
        aus_pred_class = cat(1, aus_pred_class, int32(str2num(column_names{c}(3:end-2))));
        inds_class_in_file = cat(1, inds_class_in_file, c);
    end
end

%%
inds_au_class = zeros(size(aus_SEMAINE));

for ind=1:numel(aus_SEMAINE)
    if(~isempty(find(aus_pred_class==aus_SEMAINE(ind), 1)))
        inds_au_class(ind) = inds_class_in_file(aus_pred_class==aus_SEMAINE(ind));
    end
end

preds_all = [];
for i=1:numel(vid_names)

    fname = [out_loc, vid_names{i}, '.csv'];
    preds = dlmread(fname, ',', 1, 0);
    preds_all = cat(1, preds_all, preds(vid_ids(i,1):vid_ids(i,2) - 1, :));
end

%%
f = fopen('results/SEMAINE_valid_res.txt', 'w');
f1s = zeros(1, numel(aus_SEMAINE));
for au = 1:numel(aus_SEMAINE)

    if(inds_au_class(au) ~= 0)
        tp = sum(labels_gt(:,au) == 1 & preds_all(:, inds_au_class(au)) == 1);
        fp = sum(labels_gt(:,au) == 0 & preds_all(:, inds_au_class(au)) == 1);
        fn = sum(labels_gt(:,au) == 1 & preds_all(:, inds_au_class(au)) == 0);
        tn = sum(labels_gt(:,au) == 0 & preds_all(:, inds_au_class(au)) == 0);

        precision = tp./(tp+fp);
        recall = tp./(tp+fn);

        f1 = 2 * precision .* recall ./ (precision + recall);
        f1s(au) = f1;
        fprintf(f, 'AU%d class, Precision - %.3f, Recall - %.3f, F1 - %.3f\n', aus_SEMAINE(au), precision, recall, f1);
    end

end
fclose(f);

author: TadasBaltrusaitis

repo: OpenFace


text:
function y = square( x )

%SQUARE   Internal cvx version.

% 0 : all others
% 1 : constant
% 2 : real affine
% 3 : monomial, posynomial
narginchk(1,1);
persistent remap
if isempty( remap ),
    remap1 = cvx_remap( 'constant' );
    remap2 = cvx_remap( 'affine' ) & ~remap1;
    remap3 = cvx_remap( 'log-valid' ) & ~remap1;
    remap  = remap1 + 2 * remap2 + 3 * remap3;
end
v = remap( cvx_classify( x ) );

%
% Perform the computations for each expression type separately
%

vu = sort( v(:) );
vu = vu([true;diff(vu)~=0]);
nv = length( vu );
if nv ~= 1,
    y = cvx( size( x ), [] );
end
for k = 1 : nv,

    %
    % Select the category of expression to compute
    %

    vk = vu( k );
    if nv == 1,
        xt = x;
    else
        t = v == vk;
        xt = cvx_subsref( x, t );
    end

    %
    % Perform the computations
    %

    switch vk,
        case 0,
            % Invalid
            error( 'Disciplined convex programming error:\n    Illegal operation: square( {%s} ).', cvx_class( xt, true, true ) );
        case 1,
            % Constant
            yt = cvx_constant( xt );
            yt = cvx( yt .* yt );
        case 2,
            % Real affine
            yt = quad_over_lin( xt, 1, 0 );
        case 3,
            % Monomial, posynomial
            yt = exp( 2 * log( xt ) );
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    %
    % Store the results
    %

    if nv == 1,
        y = yt;
    else
        y = cvx_subsasgn( y, t, yt );
    end

end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function Y = ttm(X,V,varargin)
%TTM Tensor times matrix.
%
%   Y = TTM(X,A,N) computes the n-mode product of tensor X with a
%   matrix A; i.e., X x_N A.  The integer N specifies the dimension
%   (or mode) of X along which A should be multiplied.  If size(A) =
%   [J,I], then X must have size(X,N) = I.  The result will be the
%   same order and size as X except that size(Y,N) = J.
%
%   Y = TTM(X,{A,B,C,...}) computes the n-mode product of tensor X
%   with a sequence of matrices in the cell array.  The n-mode
%   products are computed sequentially along all dimensions (or modes)
%   of X. The cell array contains ndims(X) matrices.
%
%   Y = TTM(X,{A,B,C,...},DIMS) computes the sequence tensor-matrix
%   products along the dimensions specified by DIMS.
%
%   Y = TTM(...,'t') performs the same computations as above except
%   the matrices are transposed.
%
%   Examples
%   X = tensor(rand(5,3,4,2));
%   A = rand(4,5); B = rand(4,3); C = rand(3,4); D = rand(3,2);
%   Y = ttm(X, A, 1)         %<-- computes X times A in mode-1
%   Y = ttm(X, {A,B,C,D}, 1) %<-- same as above
%   Y = ttm(X, A', 1, 't')   %<-- same as above
%   Y = ttm(X, {A,B,C,D}, [1 2 3 4]) %<-- 4-way multiply
%   Y = ttm(X, {D,C,B,A}, [4 3 2 1]) %<-- same as above
%   Y = ttm(X, {A,B,C,D})            %<-- same as above
%   Y = ttm(X, {A',B',C',D'}, 't')   %<-- same as above
%   Y = ttm(X, {C,D}, [3 4])     %<-- X times C in mode-3 & D in mode-4
%   Y = ttm(X, {A,B,C,D}, [3 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, [1 2 4])   %<-- 3-way multiply
%   Y = ttm(X, {A,B,C,D}, [1 2 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, -3)        %<-- same as above
%   Y = ttm(X, {A,B,C,D}, -3)      %<-- same as above
%
%   See also TENSOR, TENSOR/TTT, TENSOR/TTV.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Check the number of arguments
if (nargin < 2)
    error('TTM requires at least two arguments.');
end

%% Create 'n' and 'tflag' arguments from varargin
n = 1:ndims(X);
tflag = '';
if numel(varargin) == 1
    if ischar(varargin{1})
        tflag = varargin{1};
    else
        n = varargin{1};
    end
elseif numel(varargin) == 2
    n = varargin{1};
    tflag = varargin{2};
end

%% Handle cell array
if iscell(V)
    % Copy n into dims
    dims = n;
    % Check that the dimensions are valid
    [dims,vidx] = tt_dimscheck(dims,ndims(X),numel(V));
    % Calculate individual products
    Y = ttm(X, V{vidx(1)}, dims(1), tflag);
    for i = 2 : numel(dims)
        Y = ttm(Y, V{vidx(i)}, dims(i), tflag);
    end
    % All done
    return;
end

%% Check the second argument
if ndims(V) ~= 2
    error('tensor/ttm: 2nd argument must be a matrix.');
end

%% Check n
if (numel(n) ~= 1 || (n < 0) || n > ndims(X))
    error('Dimension N must be between 1 and NDIMS(X).');
end

%% COMPUTE THE PRODUCT

N = ndims(X);
sz = size(X);
order = [n,1:n-1,n+1:N];
newdata = double(permute(X,order));
newdata = reshape(newdata,sz(n),prod(sz([1:n-1,n+1:N])));
if tflag == 't'
    newdata = V'*newdata;
    p = size(V,2);
else
    newdata = V*newdata;
    p = size(V,1);
end
newsz = [p,sz(1:n-1),sz(n+1:N)];
Y = tensor(newdata,newsz);
Y = ipermute(Y,order);

return;

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function area = annulus_area_2d ( r1, r2 )

%*****************************************************************************80
%
%% ANNULUS_AREA_2D computes the area of a circular annulus in 2D.
%
%  Discussion:
%
%    A circular annulus with center (XC,YC), inner radius R1 and
%    outer radius R2, is the set of points (X,Y) so that
%
%      R1**2 <= (X-XC)**2 + (Y-YC)**2 <= R2**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the inner and outer radii of the annulus.
%
%    Output, real AREA, the area of the annulus.
%
  area = pi * ( r2 + r1 ) * ( r2 - r1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = r8vec_sct ( n, x )

%*****************************************************************************80
%
%% R8VEC_SCT computes a "slow" cosine transform of an R8VEC.
%
%  Discussion:
%
%    This routine is provided for illustration and testing.  It is inefficient
%    relative to optimized routines that use fast Fourier techniques.
%
%      Y(1) = Sum ( 1 <= J <= N ) X(J)
%
%      For 2 <= I <= N-1:
%
%        Y(I) = 2 * Sum ( 1 <= J <= N ) X(J)
%          * cos ( PI * ( I - 1 ) * ( J - 1 ) / ( N - 1 ) )
%
%      Y(N) = Sum ( X(1:N:2) ) - Sum ( X(2:N:2) )
%
%    Applying the routine twice in succession should yield the original data,
%    multiplied by 2 * ( N + 1 ).  This is a good check for correctness
%    and accuracy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of data values.
%
%    Input, real X(N), the data sequence.
%
%    Output, real Y(N), the transformed data.
%
  for i = 1 : n

    y(i) = x(1) / 2.0;

    for j = 2 : n - 1
      angle = pi * mod ( ( i - 1 ) * ( j - 1 ), 2 * ( n - 1 ) ) / ( n - 1 );
      y(i) = y(i) + x(j) * cos ( angle );
    end

    j = n;

    angle = pi * mod ( ( i - 1 ) * ( j - 1 ), 2 * ( n - 1 ) ) / ( n - 1 );

    y(i) = y(i) + x(n) * cos ( angle ) / 2.0;

  end

  y(1:n) = 2.0 * y(1:n) * sqrt ( n / ( n - 1 ) );

  return
end


text:
function lik = lik_coxph(varargin)
%LIK_COXPH    Create a Cox proportional hazard likelihood structure
%
%  Description
%    LIK = LIK_COXPH('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a proportional hazard model where a piecewise log-constant
%    baseline hazard is assumed.
%
%    The likelihood contribution for the ith observation is
%
%      l_i = h_i(y_i)^(1-z_i)*[exp(-int_0^y_i*h_i dt)],
%
%    where hazard is h_i=h_0(y_i)*exp(f_i). A zero mean Gaussian process
%    prior is placed for f = [f_1, f_2,...,f_n] ~ N(0, C). C is the
%    covariance matrix, whose elements are given as C_ij = c(x_i, x_j |
%    th). The function c(x_i, x_j| th) is covariance function and th its
%    parameters, hyperparameters. We place a hyperprior for
%    hyperparameters, p(th).
%
%    The time axis is partioned into K intervals with equal lengths:
%    0 = s_0 < s_1 < ... < s_K, where s_K > y_i for all i. The baseline
%    hazard rate function h_0 is piecewise constant,
%
%      h_0(t) = la_k,
%
%    when t belongs to the interval (s_{k-1},s_k] and where ft_k=log(la_k).
%    The hazard rate function is smoothed by assuming another Gaussian
%    process prior ft = [ft_1, ft_2,...,ft_K] ~ N(0, C).
%
%    z is a vector of censoring indicators with z = 0 for uncensored event
%    and z = 1 for right censored event.
%
%    When using the Coxph likelihood you need to give the vector z
%    as an extra parameter to each function that requires also y.
%    For example, you should call gpla_e as follows: gpla_e(w, gp,
%    x, y, 'z', z)
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%

% Copyright (c) 2007-2010 Jarno Vanhatalo & Jouni Hartikainen
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2011 Jaakko Riihimäki

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

  ip=inputParser;
  ip.FunctionName = 'LIK_COXPH';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('S', linspace(0,1.001,50), @(x) isvector(x));
  ip.addParamValue('stratificationVariables', [], @(x) isvector(x) && all(rem(x,1)==0));
  ip.addParamValue('removeStratificationVariables', [], @(x) ismember(x, {'on' 'off'}));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Coxph';
    lik.nondiagW=true;
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Coxph')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  if init || ~ismember('S', ip.UsingDefaults)
    s=ip.Results.S;
    xtmp=zeros(length(s)-1,1);
    for i2=1:(length(s)-1)
      xtmp(i2,1)=mean([s(i2) s(i2+1)]);
    end
    lik.xtime=xtmp;
    lik.stime=s;
  end

  if ~ismember('stratificationVariables', ip.UsingDefaults)
    lik.stratificationVariables = ip.Results.stratificationVariables;
    if ~ismember('stratificationVariables', ip.UsingDefaults)
      if isequal(ip.Results.removeStratificationVariables, 'on')
        lik.removeStratificationVariables=true;
      end
    end
  end

  if init
    % Set the function handles to the nested functions
    lik.fh.pak = @lik_coxph_pak;
    lik.fh.unpak = @lik_coxph_unpak;
    lik.fh.ll = @lik_coxph_ll;
    lik.fh.llg = @lik_coxph_llg;
    lik.fh.llg2 = @lik_coxph_llg2;
    lik.fh.llg3 = @lik_coxph_llg3;
    lik.fh.tiltedMoments = @lik_coxph_tiltedMoments;
    lik.fh.predy = @lik_coxph_predy;
    lik.fh.invlink = @lik_coxph_invlink;
    lik.fh.recappend = @lik_coxph_recappend;
    lik.fh.predcdf= @lik_coxph_predcdf;
  end

  function [w,s] = lik_coxph_pak(lik)
  %LIK_COXPH_PAK  Combine likelihood parameters into one vector.
  %
  %  Description
  %    W = LIK_COXPH_PAK(LIK) takes a likelihood structure LIK and
  %    combines the parameters into a single row vector W. This is
  %    a mandatory subfunction used for example in energy and
  %    gradient computations.
  %
  %       w = log(lik.disper)
  %
  %   See also
  %   LIK_COXPH_UNPAK, GP_PAK

    w=[];s={};
  end

  function [lik, w] = lik_coxph_unpak(lik, w)
  %LIK_COXPH_UNPAK  Extract likelihood parameters from the vector.
  %
  %  Description
  %    [LIK, W] = LIK_COXPH_UNPAK(W, LIK) takes a likelihood
  %    structure LIK and extracts the parameters from the vector W
  %    to the LIK structure. This is a mandatory subfunction used
  %    for example in energy and gradient computations.
  %
  %   Assignment is inverse of
  %       w = log(lik.disper)
  %
  %   See also
  %   LIK_COXPH_PAK, GP_UNPAK

    lik=lik;
    w=w;

  end

  function ll = lik_coxph_ll(lik, y, f, z)
  %LIK_COXPH_LL  Log likelihood
  %
  %  Description
  %    LL = LIK_COXPH_LL(LIK, Y, F, Z) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z, and
  %    latent values F. Returns the log likelihood, log p(y|f,z).
  %    This subfunction is needed when using Laplace approximation
  %    or MCMC for inference with non-Gaussian likelihoods. This
  %    subfunction is also used in information criteria (DIC, WAIC)
  %    computations.
  %
  %  See also
  %    LIK_COXPH_LLG, LIK_COXPH_LLG3, LIK_COXPH_LLG2, GPLA_E

    if isempty(z)
      error(['lik_coxph -> lik_coxph_ll: missing z!    '...
             'Coxph likelihood needs the expected number of    '...
             'occurrences as an extra input z. See, for         '...
             'example, lik_coxph and gpla_e.               ']);
    end

    [n,ny]=size(y);
    ntime=size(lik.xtime,1);
    i3v=ones(size(y,1),1);
    if isfield(lik, 'stratificationVariables')
      f1=f(1:ntime*lik.n_u);
      f2=f((ntime*lik.n_u+1):(ntime*lik.n_u+n));
      for i=1:lik.n_u
        i3v(lik.stratind{i})=i;
      end
    else
      f1=f(1:ntime);
      f2=f((ntime+1):(ntime+n));
    end

    la1=exp(f1);
    eta2=exp(f2);

    nu=1-z;
    sd=lik.stime(2)-lik.stime(1);

    if ny==1
      ll=0;
      for i1=1:n
        i3=i3v(i1);
        nft=(i3-1)*ntime;
        si=sum(y(i1)>lik.stime);
        ll=ll + nu(i1).*(f1(si+nft)+f2(i1)) - (y(i1)-lik.stime(si)).*la1(si+nft).*eta2(i1) ...
              - sum(sd.*la1((1+nft):(si-1+nft)).*eta2(i1));
      end
    else

      ll=0;
      sb=sum(bsxfun(@gt,y(:,1),lik.stime),2);
      se=sum(bsxfun(@gt,y(:,2),lik.stime),2);
      for i1=1:n
        i3=i3v(i1);
        nft=(i3-1)*ntime;
        if sb(i1)==0
          ll=ll + nu(i1).*(f1(se(i1)+nft)+f2(i1)) - (y(i1,2)-lik.stime(se(i1))).*la1(se(i1)+nft).*eta2(i1) ...
                - sum(sd.*la1((1+nft):(se(i1)-1+nft)).*eta2(i1));
        else

          if se(i1)==sb(i1)
            ll=ll + nu(i1).*(f1(se(i1)+nft)+f2(i1)) - (y(i1,2)-y(i1,1)).*la1(se(i1)+nft).*eta2(i1);
          else
            ll=ll + nu(i1).*(f1(se(i1)+nft)+f2(i1)) - (y(i1,2)-lik.stime(se(i1))).*la1(se(i1)+nft).*eta2(i1) ...
                - sum(sd.*la1((sb(i1)+1+nft):(se(i1)-1+nft)).*eta2(i1)) - (lik.stime(sb(i1)+1)-y(i1,1)).*la1(sb(i1)+nft).*eta2(i1);
          end
        end
      end
    end
  end

  function llg = lik_coxph_llg(lik, y, f, param, z)
  %LIK_COXPH_LLG  Gradient of the log likelihood
  %
  %  Description
  %    LLG = LIK_COXPH_LLG(LIK, Y, F, PARAM) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z and
  %    latent values F. Returns the gradient of the log likelihood
  %    with respect to PARAM. At the moment PARAM can be 'param' or
  %    'latent'. This subfunction is needed when using Laplace
  %    approximation or MCMC for inference with non-Gaussian
  %    likelihoods.
  %
  %  See also
  %    LIK_COXPH_LL, LIK_COXPH_LLG2, LIK_COXPH_LLG3, GPLA_E

    if isempty(z)
      error(['lik_coxph -> lik_coxph_llg: missing z!    '...
             'Coxph likelihood needs the expected number of    '...
             'occurrences as an extra input z. See, for         '...
             'example, lik_coxph and gpla_e.               ']);
    end

    ntime=size(lik.xtime,1);

    [n,ny]=size(y);
    i3v=ones(size(y));
    if isfield(lik, 'stratificationVariables')
      f1=f(1:ntime*lik.n_u);
      f2=f((ntime*lik.n_u+1):(ntime*lik.n_u+n));
      llg=zeros(ntime*lik.n_u+n,1);
      for i=1:lik.n_u
        i3v(lik.stratind{i})=i;
      end
      nf1=lik.n_u*ntime;
    else
      f1=f(1:ntime);
      f2=f((ntime+1):(ntime+n));
      llg=zeros(ntime+n,1);
      nf1=ntime;
    end

    la1=exp(f1(:));
    eta2=exp(f2(:));

    nu=1-z;
    sd=lik.stime(2)-lik.stime(1);

    switch param
      case 'latent'

      if ny==1
        if ~isfield(lik, 'stratificationVariables')
          for i1=1:ntime
            ind=y>=lik.stime(i1) & y<lik.stime(i1+1);
            llg(i1)= sum(nu(ind) - (y(ind)-lik.stime(i1)).*la1(i1).*eta2(ind)) - sum(sum(sd.*la1(i1)).*eta2(~ind & y>=lik.stime(i1+1)));
          end
        else
          for j=1:lik.n_u
            for i1=1:ntime
              ind=y>=lik.stime(i1) & y<lik.stime(i1+1);
              ind2=lik.stratind{j};
              llg(i1+(j-1)*ntime)= sum(nu(ind & ind2) - (y(ind & ind2)-lik.stime(i1)).*la1(i1+(j-1)*ntime).*eta2(ind & ind2)) - sum(sum(sd.*la1(i1+(j-1)*ntime)).*eta2((~ind & ind2) & y>=lik.stime(i1+1)));
            end
          end
        end
        for i1=1:n
          i3=i3v(i1);
          nft=(i3-1)*ntime;
          si=sum(y(i1)>lik.stime);
          llg(i1+nf1)= nu(i1) - (y(i1)-lik.stime(si)).*la1(si+nft).*eta2(i1) - sum(sd.*la1((1+nft):(si-1+nft)).*eta2(i1));
        end
      else
        for i1=1:ntime

          % left truncated + follow-up entry: (1)
          ind_vkst = y(:,1)>=lik.stime(i1) & y(:,1)<lik.stime(i1+1) & y(:,2)>=lik.stime(i1+1);
          % follow-up entry + follow-up exit: (2)
          ind_stsp = y(:,1)>=lik.stime(i1) & y(:,1)<lik.stime(i1+1) & y(:,2)>=lik.stime(i1) & y(:,2)<lik.stime(i1+1);
          % follow-up: (3)
          ind_s = y(:,1)<lik.stime(i1) & y(:,2)>=lik.stime(i1+1);
          % follow-up exit: (4)
          ind_sp = y(:,1)<lik.stime(i1) & y(:,2)>=lik.stime(i1) & y(:,2)<lik.stime(i1+1);

          if ~isfield(lik, 'stratificationVariables')
            % (1)
            s2b=sum(-(lik.stime(i1+1)-y(ind_vkst,1)).*la1(i1).*eta2(ind_vkst));
            % (2)
            s3b=sum(nu(ind_stsp)) - sum((y(ind_stsp,2)-y(ind_stsp,1)).*la1(i1).*eta2(ind_stsp));
            % (3)
            s4= - sum(sd.*la1(i1).*eta2(ind_s));
            % (4)
            s5=sum(nu(ind_sp)) - sum((y(ind_sp,2)-lik.stime(i1)).*la1(i1).*eta2(ind_sp));

            llg(i1) = s2b+s3b+s4+s5;
          else
            for j=1:lik.n_u
              ind2 = lik.stratind{j};
              nft=(j-1)*ntime;
              % (1)
              s2b=sum(-(lik.stime(i1+1)-y(ind_vkst & ind2,1)).*la1(i1+nft).*eta2(ind_vkst & ind2));
              % (2)
              s3b=sum(nu(ind_stsp & ind2)) - sum((y(ind_stsp & ind2,2)-y(ind_stsp & ind2,1)).*la1(i1+nft).*eta2(ind_stsp & ind2));
              % (3)
              s4= - sum(sd.*la1(i1+nft).*eta2(ind_s & ind2));
              % (4)
              s5=sum(nu(ind_sp & ind2)) - sum((y(ind_sp & ind2,2)-lik.stime(i1)).*la1(i1+nft).*eta2(ind_sp & ind2));

              llg(i1+nft) = s2b+s3b+s4+s5;
            end
          end

        end

        sb=sum(bsxfun(@gt,y(:,1),lik.stime),2);
        se=sum(bsxfun(@gt,y(:,2),lik.stime),2);
        for i1=1:n
          i3=i3v(i1);
          nft=(i3-1)*ntime;
          if sb(i1)==0
            llg(i1+nf1)= nu(i1) - (y(i1,2)-lik.stime(se(i1))).*la1(se(i1)+nft).*eta2(i1) - sum(sd.*la1((1+nft):(se(i1)-1+nft)).*eta2(i1));
          else
            if se(i1)==sb(i1)
              llg(i1+nf1) = nu(i1) - (y(i1,2)-y(i1,1)).*la1(se(i1)+nft).*eta2(i1);
            else
              llg(i1+nf1) = nu(i1) - (y(i1,2)-lik.stime(se(i1))).*la1(se(i1)+nft).*eta2(i1) - sum(sd.*la1((sb(i1)+1+nft):(se(i1)-1+nft)).*eta2(i1)) - (lik.stime(sb(i1)+1)-y(i1,1)).*la1(sb(i1)+nft).*eta2(i1);
            end
          end
        end

      end

    end
  end

  function [llg2,llg2mat] = lik_coxph_llg2(lik, y, f, param, z)
  %function [pi_vec, pi_mat] = lik_coxph_llg2(lik, y, ff, param, z)
  %LIK_COXPH_LLG2  Second gradients of the log likelihood
  %
  %  Description
  %    LLG2 = LIK_COXPH_LLG2(LIK, Y, F, PARAM) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z, and
  %    latent values F. Returns the Hessian of the log likelihood
  %    with respect to PARAM. At the moment PARAM can be only
  %    'latent'. LLG2 is a vector with diagonal elements of the
  %    Hessian matrix (off diagonals are zero). This subfunction
  %    is needed when using Laplace approximation or EP for inference
  %    with non-Gaussian likelihoods.
  %
  %  See also
  %    LIK_COXPH_LL, LIK_COXPH_LLG, LIK_COXPH_LLG3, GPLA_E

    if isempty(z)
      error(['lik_coxph -> lik_coxph_llg2: missing z!   '...
             'Coxph likelihood needs the expected number of    '...
             'occurrences as an extra input z. See, for         '...
             'example, lik_coxph and gpla_e.               ']);
    end

    [n,ny]=size(y);
    ntime=size(lik.xtime,1);
    i3v=ones(size(y,1),1);
    if isfield(lik, 'stratificationVariables')
      f1=f(1:ntime*lik.n_u);
      f2=f((ntime*lik.n_u+1):(ntime*lik.n_u+n));
      llg2=zeros(ntime*lik.n_u+n,1);
      llg2mat=zeros(ntime*lik.n_u,n);
      for i=1:lik.n_u
        i3v(lik.stratind{i})=i;
      end
      nf1=ntime*lik.n_u;
    else
      f1=f(1:ntime);
      f2=f((ntime+1):(ntime+n));
      llg2=zeros(ntime+n,1);
      llg2mat=zeros(ntime,n);
      nf1=ntime;
    end

    la1=exp(f1);
    eta2=exp(f2);

    %nu=1-z;
    sd=lik.stime(2)-lik.stime(1);
    switch param
      case 'latent'

        if ny==1
          % 11
          if ~isfield(lik, 'stratificationVariables')
            for i1=1:ntime
              ind=y>=lik.stime(i1) & y<lik.stime(i1+1);
              llg2(i1)= sum(-(y(ind)-lik.stime(i1)).*la1(i1).*eta2(ind)) - sum(sum(sd.*la1(i1)).*eta2(~ind & y>=lik.stime(i1+1)));
              %llg2(i1,i1)= sum(-(y(ind)-lik.stime(i1)).*la1(i1).*eta2(ind)) - sum(sum(sd.*la1(i1)).*eta2(~ind & y>=lik.stime(i1+1)));
            end
          else
            for j=1:lik.n_u
              for i1=1:ntime
                ind=y>=lik.stime(i1) & y<lik.stime(i1+1);
                ind2=lik.stratind{j};
                llg2(i1+(j-1)*ntime)= sum(-(y(ind & ind2)-lik.stime(i1)).*la1(i1+ntime*(j-1)).*eta2(ind & ind2)) - sum(sum(sd.*la1(i1+ntime*(j-1))).*eta2((~ind & ind2) & y>=lik.stime(i1+1)));
                %llg2(i1,i1)= sum(-(y(ind)-lik.stime(i1)).*la1(i1).*eta2(ind)) - sum(sum(sd.*la1(i1)).*eta2(~ind & y>=lik.stime(i1+1)));
              end
            end
          end

          % 22
          for i1=1:n
            si=sum(y(i1)>lik.stime);
            i3=i3v(i1);
            llg2(i1+nf1)= -(y(i1)-lik.stime(si)).*la1(si+(i3-1)*ntime).*eta2(i1) - sum(sd.*la1((1+(i3-1)*ntime):(si-1+(i3-1)*ntime)).*eta2(i1));
            %llg2(i1+ntime,i1+ntime)= -(y(i1)-lik.stime(si)).*la1(si).*eta2(i1) - sum(sd.*la1(1:(si-1)).*eta2(i1));
          end

          % derivative wrt f1 and f2:
          if ~isfield(lik, 'stratificationVariables')
            for i1=1:ntime
              ind=y>=lik.stime(i1) & y<lik.stime(i1+1);
              llg2mat(i1,find(ind))= -(y(ind)-lik.stime(i1)).*la1(i1).*eta2(ind);
              llg2mat(i1,find((~ind & y>=lik.stime(i1+1)))) = - sd.*la1(i1).*eta2((~ind & y>=lik.stime(i1+1)));
              %llg2(i1,ntime+find(ind))= -(y(ind)-lik.stime(i1)).*la1(i1).*eta2(ind);
              %llg2(i1,ntime+find((~ind & y>=lik.stime(i1+1)))) = - sd.*la1(i1).*eta2((~ind & y>=lik.stime(i1+1)));
              %llg2(ntime+find(ind),i1)=llg2(i1,ntime+find(ind));
              %llg2(ntime+find((~ind & y>=lik.stime(i1+1))),i1)=llg2(i1,ntime+find((~ind & y>=lik.stime(i1+1))));
            end
          else
            for j=1:lik.n_u
              ind2=lik.stratind{j};
              for i1=1:ntime
                ind=y>=lik.stime(i1) & y<lik.stime(i1+1);
                llg2mat(i1+(j-1)*ntime,find(ind & ind2))= -(y(ind & ind2)-lik.stime(i1)).*la1(i1+(j-1)*ntime).*eta2(ind & ind2);
                llg2mat(i1+(j-1)*ntime,find(((~ind & ind2) & y>=lik.stime(i1+1)))) = - sd.*la1(i1+(j-1)*ntime).*eta2(((~ind & ind2) & y>=lik.stime(i1+1)));
                %llg2(i1,ntime+find(ind))= -(y(ind)-lik.stime(i1)).*la1(i1).*eta2(ind);
                %llg2(i1,ntime+find((~ind & y>=lik.stime(i1+1)))) = - sd.*la1(i1).*eta2((~ind & y>=lik.stime(i1+1)));
                %llg2(ntime+find(ind),i1)=llg2(i1,ntime+find(ind));
                %llg2(ntime+find((~ind & y>=lik.stime(i1+1))),i1)=llg2(i1,ntime+find((~ind & y>=lik.stime(i1+1))));
              end
            end
          end

        else

          % 11
          for i1=1:ntime

            % left truncated + follow-up entry: (1)
            ind_vkst = y(:,1)>=lik.stime(i1) & y(:,1)<lik.stime(i1+1) & y(:,2)>=lik.stime(i1+1);
            % follow-up entry + follow-up exit: (2)
            ind_stsp = y(:,1)>=lik.stime(i1) & y(:,1)<lik.stime(i1+1) & y(:,2)>=lik.stime(i1) & y(:,2)<lik.stime(i1+1);
            % follow-up: (3)
            ind_s = y(:,1)<lik.stime(i1) & y(:,2)>=lik.stime(i1+1);
            % follow-up exit: (4)
            ind_sp = y(:,1)<lik.stime(i1) & y(:,2)>=lik.stime(i1) & y(:,2)<lik.stime(i1+1);

            if ~isfield(lik, 'stratificationVariables')
              % (1)
              s2b=sum(-(lik.stime(i1+1)-y(ind_vkst,1)).*la1(i1).*eta2(ind_vkst));
              % (2)
              s3b=-sum((y(ind_stsp,2)-y(ind_stsp,1)).*la1(i1).*eta2(ind_stsp));
              % (3)
              s4=-sum(sd.*la1(i1).*eta2(ind_s));
              % (4)
              s5=-sum((y(ind_sp,2)-lik.stime(i1)).*la1(i1).*eta2(ind_sp));
              llg2(i1) = s2b+s3b+s4+s5;
            else
              for j=1:lik.n_u
                ind2=lik.stratind{j};
                nft=(j-1)*ntime;
                % (1)
                s2b=sum(-(lik.stime(i1+1)-y(ind_vkst & ind2,1)).*la1(i1+nft).*eta2(ind_vkst & ind2));
                % (2)
                s3b=-sum((y(ind_stsp & ind2,2)-y(ind_stsp & ind2,1)).*la1(i1+nft).*eta2(ind_stsp & ind2));
                % (3)
                s4=-sum(sd.*la1(i1+nft).*eta2(ind_s & ind2));
                % (4)
                s5=-sum((y(ind_sp & ind2,2)-lik.stime(i1)).*la1(i1+nft).*eta2(ind_sp & ind2));
                llg2(i1+nft) = s2b+s3b+s4+s5;
              end
            end
          end

          % 22
          sb=sum(bsxfun(@gt,y(:,1),lik.stime),2);
          se=sum(bsxfun(@gt,y(:,2),lik.stime),2);
          for i1=1:n
            i3=i3v(i1);
            nft=(i3-1)*ntime;
            if sb(i1)==0
              llg2(i1+nf1)= -(y(i1,2)-lik.stime(se(i1))).*la1(se(i1)+nft).*eta2(i1) - sum(sd.*la1((1+nft):(se(i1)-1+nft)).*eta2(i1));
              %llg2(i1+ntime,i1+ntime)= -(y(i1,2)-lik.stime(se)).*la1(se).*eta2(i1) - sum(sd.*la1(1:(se-1)).*eta2(i1));
            else
              if se(i1)==sb(i1)
                llg2(i1+nf1) = -(y(i1,2)-y(i1,1)).*la1(se(i1)+nft).*eta2(i1);
                %llg2(i1+ntime,i1+ntime) = -(y(i1,2)-y(i1,1)).*la1(se).*eta2(i1);
              else
                llg2(i1+nf1) = -(y(i1,2)-lik.stime(se(i1))).*la1(se(i1)+nft).*eta2(i1) - sum(sd.*la1((sb(i1)+1+nft):(se(i1)-1+nft)).*eta2(i1)) - (lik.stime(sb(i1)+1)-y(i1,1)).*la1(sb(i1)+nft).*eta2(i1);
                %llg2(i1+ntime,i1+ntime) = -(y(i1,2)-lik.stime(se)).*la1(se).*eta2(i1) - sum(sd.*la1((sb+1):(se-1)).*eta2(i1)) - (lik.stime(sb+1)-y(i1,1)).*la1(sb).*eta2(i1);
              end
            end
          end

          % derivative wrt f1 and f2:
          for i1=1:ntime

            % left truncated + follow-up entry: (1)
            ind_vkst = y(:,1)>=lik.stime(i1) & y(:,1)<lik.stime(i1+1) & y(:,2)>=lik.stime(i1+1);
            % follow-up entry + follow-up exit: (2)
            ind_stsp = y(:,1)>=lik.stime(i1) & y(:,1)<lik.stime(i1+1) & y(:,2)>=lik.stime(i1) & y(:,2)<lik.stime(i1+1);
            % follow-up: (3)
            ind_s = y(:,1)<lik.stime(i1) & y(:,2)>=lik.stime(i1+1);
            % follow-up exit: (4)
            ind_sp = y(:,1)<lik.stime(i1) & y(:,2)>=lik.stime(i1) & y(:,2)<lik.stime(i1+1);

            if ~isfield(lik, 'stratificationVariables')
              % (1)
              llg2mat(i1,ind_vkst)=-(lik.stime(i1+1)-y(ind_vkst,1)).*la1(i1).*eta2(ind_vkst);
              % (2)
              llg2mat(i1,ind_stsp)=-(y(ind_stsp,2)-y(ind_stsp,1)).*la1(i1).*eta2(ind_stsp);
              % (3)
              llg2mat(i1,ind_s)= -sd.*la1(i1).*eta2(ind_s);
              % (4)
              llg2mat(i1,ind_sp)=-(y(ind_sp,2)-lik.stime(i1)).*la1(i1).*eta2(ind_sp);
            else
              for j=1:lik.n_u
                ind2=lik.stratind{j};
                nft=(j-1)*ntime;
                % (1)
                llg2mat(i1+nft,ind_vkst & ind2)=-(lik.stime(i1+1)-y(ind_vkst & ind2,1)).*la1(i1+nft).*eta2(ind_vkst & ind2);
                % (2)
                llg2mat(i1+nft,ind_stsp & ind2)=-(y(ind_stsp & ind2,2)-y(ind_stsp & ind2,1)).*la1(i1+nft).*eta2(ind_stsp & ind2);
                % (3)
                llg2mat(i1+nft,ind_s & ind2)= -sd.*la1(i1+nft).*eta2(ind_s & ind2);
                % (4)
                llg2mat(i1+nft,ind_sp & ind2)=-(y(ind_sp & ind2,2)-lik.stime(i1)).*la1(i1+nft).*eta2(ind_sp & ind2);
              end
            end
          end
      end
    end
  end

  function [llg3,llg3mat] = lik_coxph_llg3(lik, y, f, param, z, j1)
  %LIK_COXPH_LLG3  Third gradients of the log likelihood
  %
  %  Description
  %    LLG3 = LIK_COXPH_LLG3(LIK, Y, F, PARAM) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z and
  %    latent values F and returns the third gradients of the log
  %    likelihood with respect to PARAM. At the moment PARAM can be
  %    only 'latent'. LLG3 is a vector with third gradients. This
  %    subfunction is needed when using Laplace approximation for
  %    inference with non-Gaussian likelihoods.
  %
  %  See also
  %    LIK_COXPH_LL, LIK_COXPH_LLG, LIK_COXPH_LLG2, GPLA_E, GPLA_G

    if isempty(z)
      error(['lik_coxph -> lik_coxph_llg3: missing z!   '...
             'Coxph likelihood needs the expected number of    '...
             'occurrences as an extra input z. See, for         '...
             'example, lik_coxph and gpla_e.               ']);
    end

    ntime=size(lik.xtime,1);

    [n,ny]=size(y);
    if isfield(lik, 'stratificationVariables')
      f1=f(1:ntime*lik.n_u);
      f2=f((ntime*lik.n_u+1):(ntime*lik.n_u+n));
      llg3=zeros(ntime*lik.n_u+n,1);
      llg3mat=zeros(ntime*lik.n_u,n);
      nf1=ntime*lik.n_u;
      if j1>nf1
        for j=1:lik.n_u
          if lik.stratind{j}(j1-nf1)==1;
            i3=j;
            break;
          end
        end
      end
    else
      f1=f(1:ntime);
      f2=f((ntime+1):(ntime+n));
      llg3=zeros(ntime+n,1);
      llg3mat=zeros(ntime,n);
      nf1=ntime;
      i3=1;
    end

    la1=exp(f1);
    eta2=exp(f2);

    %nu=1-z;
    sd=lik.stime(2)-lik.stime(1);

    switch param
      case 'latent'

        if ny==1

          if j1<=nf1

            indt=rem(j1,ntime);
            if indt==0
              indt=ntime;
            end

            % 11
            ind=y>=lik.stime(indt) & y<lik.stime(indt+1);

            if isfield(lik, 'stratificationVariables')
              ind2=lik.stratind{(j1-indt)/ntime+1};
            else
              ind2=ones(size(ind));
            end

            fi=find(ind & ind2);
            fni=find((~ind & ind2) & y>=lik.stime(indt));

            llg3(j1) = sum(-(y(ind & ind2)-lik.stime(indt)).*la1(j1).*eta2(ind & ind2)) - sum(sum(sd.*la1(j1)).*eta2((~ind & ind2) & y>=lik.stime(indt+1)));

            if ~isempty(fi)
              valtmp=(-(y(ind & ind2)-lik.stime(indt)).*la1(j1).*eta2(ind & ind2));
              for m2i=1:length(valtmp)
                llg3( nf1+fi(m2i))  = valtmp(m2i);
              end
            end
            if ~isempty(fni)
              valtmp2=(-sd.*la1(j1).*eta2(((~ind & ind2) & y>=lik.stime(indt+1))));
              for m2i=1:length(valtmp2)
                llg3( nf1+fni(m2i))  = valtmp2(m2i);
              end
            end

            % 12/21
            % derivative wrt f1 and f2:
            val1tmp=-(y(ind & ind2)-lik.stime(indt)).*la1(j1).*eta2(ind & ind2);
            llg3mat(j1,fi)= val1tmp;

            val2tmp = - sd.*la1(j1).*eta2(((~ind & ind2) & y>=lik.stime(indt+1)));
            llg3mat(j1,fni) = val2tmp;

          else

            % 11
            s1=sum(y(j1-nf1)>lik.stime);
            llg3((1+(i3-1)*ntime):(s1-1+(i3-1)*ntime)) = - sd.*la1((1+(i3-1)*ntime):(s1-1+(i3-1)*ntime)).*eta2(j1-nf1);
            llg3(s1+(i3-1)*ntime) = -(y(j1-nf1)-lik.stime(s1)).*la1(s1+(i3-1)*ntime).*eta2(j1-nf1);

            % 22
            llg3(j1) = -(y(j1-nf1)-lik.stime(s1)).*la1(s1+(i3-1)*ntime).*eta2(j1-nf1) ...
                        - sum(sd.*la1((1+(i3-1)*ntime):(s1-1+(i3-1)*ntime)).*eta2(j1-nf1));

            % 12/21
            % derivative wrt f1 and f2:
            val3tmp = - sd.*la1((1+(i3-1)*ntime):(s1-1+(i3-1)*ntime)).*eta2(j1-nf1);
            llg3mat((1+(i3-1)*ntime):(s1-1+(i3-1)*ntime),j1-nf1)= val3tmp;
            llg3mat(s1+(i3-1)*ntime,j1-nf1) = -(y(j1-nf1)-lik.stime(s1)).*la1(s1+(i3-1)*ntime).*eta2(j1-nf1);

          end

        else

          if j1<=nf1

            indt=rem(j1,ntime);
            if indt==0
              indt=ntime;
            end
            if isfield(lik, 'stratificationVariables')
              ind2=lik.stratind{(j1-indt)/ntime+1};
            else
              ind2=ones(size(y,1),1);
            end

            % 11
            % left truncated + follow-up entry: (1)
            ind_vkst = y(:,1)>=lik.stime(indt) & y(:,1)<lik.stime(indt+1) & y(:,2)>=lik.stime(indt+1);
            % follow-up entry + follow-up exit: (2)
            ind_stsp = y(:,1)>=lik.stime(indt) & y(:,1)<lik.stime(indt+1) & y(:,2)>=lik.stime(indt) & y(:,2)<lik.stime(indt+1);
            % follow-up: (3)
            ind_s = y(:,1)<lik.stime(indt) & y(:,2)>=lik.stime(indt+1);
            % follow-up exit: (4)
            ind_sp = y(:,1)<lik.stime(indt) & y(:,2)>=lik.stime(indt) & y(:,2)<lik.stime(indt+1);

            % (1)
            s2b=sum(-(lik.stime(indt+1)-y(ind_vkst & ind2,1)).*la1(j1).*eta2(ind_vkst & ind2));
            % (2)
            s3b=-sum((y(ind_stsp & ind2,2)-y(ind_stsp & ind2,1)).*la1(j1).*eta2(ind_stsp & ind2));
            % (3)
            s4=-sum(sd.*la1(j1).*eta2(ind_s & ind2));
            % (4)
            s5=-sum((y(ind_sp & ind2,2)-lik.stime(indt)).*la1(j1).*eta2(ind_sp & ind2));

            llg3(j1) = s2b+s3b+s4+s5;

            % 22
            % (1)
            llg3(nf1+find(ind_vkst & ind2))=-(lik.stime(indt+1)-y(ind_vkst & ind2,1)).*la1(j1).*eta2(ind_vkst & ind2);
            % (2)
            llg3(nf1+find(ind_stsp & ind2))=-(y(ind_stsp & ind2,2)-y(ind_stsp & ind2,1)).*la1(j1).*eta2(ind_stsp & ind2);
            % (3)
            llg3(nf1+find(ind_s & ind2))= -sd.*la1(j1).*eta2(ind_s & ind2);
            % (4)
            llg3(nf1+find(ind_sp & ind2))=-(y(ind_sp & ind2,2)-lik.stime(indt)).*la1(j1).*eta2(ind_sp & ind2);

            % 12/21
            llg3mat(j1,find(ind_vkst & ind2))=-(lik.stime(indt+1)-y(ind_vkst & ind2,1)).*la1(j1).*eta2(ind_vkst & ind2);
            % (2)
            llg3mat(j1,find(ind_stsp & ind2))=-(y(ind_stsp & ind2,2)-y(ind_stsp & ind2,1)).*la1(j1).*eta2(ind_stsp & ind2);
            % (3)
            llg3mat(j1,find(ind_s & ind2))= -sd.*la1(j1).*eta2(ind_s & ind2);
            % (4)
            llg3mat(j1,find(ind_sp & ind2))=-(y(ind_sp & ind2,2)-lik.stime(indt)).*la1(j1).*eta2(ind_sp & ind2);
          else

            sb=sum(y(j1-nf1,1)>lik.stime); % begin
            se=sum(y(j1-nf1,2)>lik.stime); % end

            nft=(i3-1)*ntime;
            % 11
            if sb==0
              llg3((1+nft):(se-1+nft))= -sd.*la1((1+nft):(se-1+nft)).*eta2(j1-nf1);
              llg3(se+nft)= -(y(j1-nf1,2)-lik.stime(se)).*la1(se+nft).*eta2(j1-nf1);
            else
              if se==sb
                llg3(se+nft) = -(y(j1-nf1,2)-y(j1-nf1,1)).*la1(se+nft).*eta2(j1-nf1);
              else
                llg3(sb+nft) = - (lik.stime(sb+1)-y(j1-nf1,1)).*la1(sb+nft).*eta2(j1-nf1);
                llg3((sb+1+nft):(se-1+nft)) = - sd.*la1((sb+1+nft):(se-1+nft)).*eta2(j1-nf1);
                llg3(se+nft) = -(y(j1-nf1,2)-lik.stime(se)).*la1(se+nft).*eta2(j1-nf1);
              end
            end

            % 12/21
            if sb==0
              llg3mat((1+nft):(se-1+nft),j1-nf1) = -sd.*la1((1+nft):(se-1+nft)).*eta2(j1-nf1);
              llg3mat(se+nft,j1-nf1)= -(y(j1-nf1,2)-lik.stime(se)).*la1(se+nft).*eta2(j1-nf1);
            else
              if se==sb
                llg3mat(se+nft,j1-nf1) = -(y(j1-nf1,2)-y(j1-nf1,1)).*la1(se+nft).*eta2(j1-nf1);
              else
                llg3mat(sb+nft,j1-nf1) = - (lik.stime(sb+1)-y(j1-nf1,1)).*la1(sb+nft).*eta2(j1-nf1);
                llg3mat((sb+1+nft):(se-1+nft),j1-nf1) = - sd.*la1((sb+1+nft):(se-1+nft)).*eta2(j1-nf1);
                llg3mat(se+nft,j1-nf1) = -(y(j1-nf1,2)-lik.stime(se)).*la1(se+nft).*eta2(j1-nf1);
              end
            end

            % 22
            if sb==0
              llg3(j1)= -(y(j1-nf1,2)-lik.stime(se)).*la1(se+nft).*eta2(j1-nf1) - sum(sd.*la1((1+nft):(se-1+nft)).*eta2(j1-nf1));
            else
              if se==sb
                llg3(j1) = -(y(j1-nf1,2)-y(j1-nf1,1)).*la1(se+nft).*eta2(j1-nf1);
              else
                llg3(j1) = -(y(j1-nf1,2)-lik.stime(se)).*la1(se+nft).*eta2(j1-nf1) - sum(sd.*la1((sb+1+nft):(se-1+nft)).*eta2(j1-nf1)) - (lik.stime(sb+1)-y(j1-nf1,1)).*la1(sb+nft).*eta2(j1-nf1);
              end
            end

          end
        end
    end
  end

  function [logM_0, m_1, sigm2hati1] = lik_coxph_tiltedMoments(lik, y, i1, S2_i, M_i, z)
    %LIK_COXPH_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
    %
    %  Description
    %    [M_0, M_1, M2] = LIK_COXPH_TILTEDMOMENTS(LIK, Y, I, S2,
    %    MYY, Z) takes a likelihood structure LIK, incedence counts
    %    Y, expected counts Z, index I and cavity variance S2 and
    %    mean MYY. Returns the zeroth moment M_0, mean M_1 and
    %    variance M_2 of the posterior marginal (see Rasmussen and
    %    Williams (2006): Gaussian processes for Machine Learning,
    %    page 55). This subfunction is needed when using EP for
    %    inference with non-Gaussian likelihoods.
    %
    %  See also
    %    GPEP_E

      [n,ny]=size(y);

      % M_i(end);
      % S2_i(end,end);
      fgrid=M_i(end)+sqrt(S2_i(end,end))*[-6 6];
      fg=linspace(fgrid(1),fgrid(2),15);
      ng=length(fg);

      if isfield(gp.lik, 'stratificationVariables')
        ntime=size(lik.xtime,1)*gp.lik.n_u;
      else
        ntime=size(lik.xtime,1);
      end

      %f11=f(1:ntime);
      %f2=f((ntime+1):(ntime+n));
      %la1=exp(f1);
      %eta2=exp(f2);

      nu=1-z;
      sd=lik.stime(2)-lik.stime(1);
      if ny==1
        sb=1;
        se=sum(bsxfun(@gt,y(i1,1),lik.stime),2);
      end
      indf=sb:se;
      sdvec=[ones(se-1,1)*sd; y(i1)-lik.stime(se)];
      nutmp=zeros(se,1);
      nutmp(se)=nu(i1);

      for j1=1:ng

        % conditional distribution
        myy=M_i(indf)+S2_i(indf,end)*(1./S2_i(end,end))*(fg(j1)-M_i(end));
        myy0=myy;
        Sigm=S2_i(indf,indf)-S2_i(indf,end)*(1./S2_i(end,end))*S2_i(end,indf);
        Sigm0=Sigm;

        nu_prior=Sigm\myy;

        nt=size(myy,1);
        c1=exp(fg(j1));

        % site parameters
        tautilde=zeros(nt,1);
        nutilde=zeros(nt,1);
        ztilde=zeros(nt,1);

        max_small_ep_iter=50;
        tol=1e-9;
        small_ep_iter=1;

        tautilde0=Inf; nutilde0=Inf; ztilde0=Inf;

        logZep_tmp=0; logZep=Inf;
        %while small_ep_iter <= max_small_ep_iter && (sum(abs(tautilde0-tautilde)>tol) || sum(abs(nutilde0-nutilde)>tol) || sum(abs(ztilde0-ztilde)>tol))
        while small_ep_iter<=max_small_ep_iter && abs(logZep_tmp-logZep)>tol
          logZep_tmp=logZep;

          %tautilde0=tautilde; nutilde0=nutilde; ztilde0=ztilde;

          for k1=1:nt

            tau_i=Sigm(k1,k1)^-1-tautilde(k1);
            nu_i = Sigm(k1,k1)^-1*myy(k1)-nutilde(k1);
            myy_i=nu_i/tau_i;
            sigm2_i=tau_i^-1;

            % marginal moments
            [logM0(k1), muhati, sigm2hati] = coxph_tiltedMoments(sigm2_i, myy_i, nutmp(k1), sdvec(k1), c1);
            %[M0, muhati, sigm2hati] = coxph_tiltedMoments(lik, y(i1,:), k1, sigm2_i, myy_i, c1, sd_vec(i1), ztmp);

            deltatautilde=sigm2hati^-1-tau_i-tautilde(k1);
            tautilde(k1)=tautilde(k1)+deltatautilde;
            nutilde(k1)=sigm2hati^-1*muhati-nu_i;

            apu = deltatautilde/(1+deltatautilde*Sigm(k1,k1));
            Sigm = Sigm - apu*(Sigm(:,k1)*Sigm(:,k1)');

            % The below is how Rasmussen and Williams
            % (2006) do the update. The above version is
            % more robust.
            %apu = deltatautilde^-1+Sigm(k1,k1);
            %apu = (Sigm(:,k1)/apu)*Sigm(:,k1)';
            %Sigm = Sigm - apu;
            %Sigm=Sigm-(deltatautilde^-1+Sigm(k1,k1))^-1*(Sigm(:,k1)*Sigm(:,k1)');

            %myy=Sigm*nutilde;
            myy=Sigm*(nutilde+nu_prior);

            muvec_i(k1,1)=myy_i;
            sigm2vec_i(k1,1)=sigm2_i;

          end

          if tautilde > 0
            Stilde=tautilde;
            Stildesqroot=diag(sqrt(tautilde));
            B=eye(nt)+Stildesqroot*Sigm0*Stildesqroot;
            L=chol(B,'lower');

            V=(L\Stildesqroot)*Sigm0;
            Sigm=Sigm0-V'*V;
            %myy=Sigm*nutilde;
            myy=Sigm*(nutilde+nu_prior);

            %Ls = chol(Sigm);

            % Compute the marginal likelihood
            % Direct formula (3.65):
            % Sigmtilde=diag(1./tautilde);
            % mutilde=inv(Stilde)*nutilde;
            %
            % logZep=-0.5*log(det(Sigmtilde+K))-0.5*mutilde'*inv(K+Sigmtilde)*mutilde+
            %         sum(log(normcdf(y.*muvec_i./sqrt(1+sigm2vec_i))))+
            %         0.5*sum(log(sigm2vec_i+1./tautilde))+
            %         sum((muvec_i-mutilde).^2./(2*(sigm2vec_i+1./tautilde)))

            % 4. term & 1. term
            term41=0.5*sum(log(1+tautilde.*sigm2vec_i))-sum(log(diag(L)));

            % 5. term (1/2 element) & 2. term
            T=1./sigm2vec_i;
            Cnutilde = Sigm0*nutilde;
            L2 = V*nutilde;
            term52 = nutilde'*Cnutilde - L2'*L2 - (nutilde'./(T+Stilde)')*nutilde;
            term52 = term52.*0.5;

            % 5. term (2/2 element)
            term5=0.5*muvec_i'.*(T./(Stilde+T))'*(Stilde.*muvec_i-2*nutilde);

            % 3. term
            term3 = sum(logM0);

            V_tmp=(L\Stildesqroot);
            Sigm_inv_tmp=V_tmp'*V_tmp;

            term_add1=-0.5*myy0'*Sigm_inv_tmp*myy0;
            term_add2=myy0'*(eye(nt)-Sigm_inv_tmp*Sigm0)*nutilde;
            logZep = -(term41+term52+term5+term3+term_add1+term_add2);

            %logZep = -(term41+term52+term5+term3);

            small_ep_iter=small_ep_iter+1;
            %iter=iter+1;
          else
            error('tautilde <= 0')
          end
        end

        ZZ(j1,1)=exp(-logZep);
        MM(:,j1)=myy;
        SS2(:,:,j1)=Sigm;

      end

      %m_0=zeros(1,1);
      %m_1=zeros(nt+1,1);
      %sigm2hati1=zeros(nt+1,nt+1);
      % indf

      W=normpdf(fg,M_i(end),sqrt(S2_i(end,end)))*(fg(2)-fg(1));

      qw=W.*ZZ';
      m_0=sum(qw);
      m_1=[sum(bsxfun(@times,qw,MM),2); sum(qw.*fg)]./m_0;

      m_211=zeros(nt,nt);
      for k1=1:ng
        m_211=m_211+qw(k1)*(SS2(:,:,k1)+MM(:,k1)*MM(:,k1)');
      end
      m_212=(qw.*fg)*MM';
      m_222=(qw.*fg)*fg';

      m_2=[m_211 m_212'; m_212 m_222]./m_0;

      sigm2hati1=m_2 - m_1*m_1';
      logM_0 = log(m_0);

      %figure(1),hold on, plot(fg(j1),logZep,'.')
      %figure(2),hold on, plot(fg(j1),exp(-logZep),'.')

  end

  function [lpyt,Ey, Vary] = lik_coxph_predy(lik, Ef, Covf, yt, zt)
  %LIK_COXPH_PREDY  Returns the predictive mean, variance and density of y
  %
  %  Description
  %    [EY, VARY] = LIK_COXPH_PREDY(LIK, EF, VARF) takes a
  %    likelihood structure LIK, posterior mean EF and posterior
  %    Variance VARF of the latent variable and returns the
  %    posterior predictive mean EY and variance VARY of the
  %    observations related to the latent variables. This
  %    subfunction is needed when computing posterior predictive
  %    distributions for future observations.
  %
  %    [Ey, Vary, PY] = LIK_COXPH_PREDY(LIK, EF, VARF YT, ZT)
  %    Returns also the predictive density of YT, that is
  %        p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
  %    This requires also the incedence counts YT, expected counts ZT.
  %    This subfunction is needed when computing posterior predictive
  %    distributions for future observations.
  %
  %  See also
  %    GPLA_PRED, GPEP_PRED, GPMC_PRED

    if isempty(zt)
      error(['lik_coxph -> lik_coxph_predy: missing zt!'...
             'Coxph likelihood needs the expected number of    '...
             'occurrences as an extra input zt. See, for         '...
             'example, lik_coxph and gpla_e.               ']);
    end
    ntime=size(lik.xtime,1);

    ntest=size(zt,1);
    ny=size(yt,2);

    if isfield(lik, 'stratificationVariables')
      nf1=lik.n_u*ntime;
      i3v=zeros(ntest,1);
      for ii=1:length(lik.stratindt)
        i3v(lik.stratindt{ii})=ii;
      end
    else
      nf1=ntime;
      i3v=ones(ntest,1);
    end

    Py = zeros(size(zt));
    %Ey = zeros(size(zt));
    %EVary = zeros(size(zt));
    %VarEy = zeros(size(zt));

    S=10000;
    sd=lik.stime(2)-lik.stime(1);
    nu=1-zt;

    [nn1,nn2,c] =size(Covf);
    if (c>1) || (nn1==nn2)
      mcmc=false;
    else
      mcmc=true;
    end

    for i1=1:ntest
      i3=i3v(i1);
      if mcmc
        Sigm_tmp=([Covf(1:ntime,:); Covf(i1+ntime,:)]);
        f_star=bsxfun(@plus,Ef([(1+(i3-1)*ntime):(i3*ntime) i1+nf1])', ...
          bsxfun(@times,sqrt(Sigm_tmp'),randn(S,nf1+1)));
      else
        Sigm_tmp=Covf([(1+(i3-1)*ntime):(i3*ntime) i1+nf1],[(1+(i3-1)*ntime):(i3*ntime) i1+nf1]);
        Sigm_tmp=(Sigm_tmp+Sigm_tmp')./2;
        f_star=mvnrnd(Ef([(1+(i3-1)*ntime):(i3*ntime) i1+nf1]), Sigm_tmp, S);
      end

      f1=f_star(:,1:ntime);
      f2=f_star(:,(ntime+1):end);

      la1=exp(f1);
      eta2=exp(f2);

      if ny==1
        si=sum(yt(i1)>lik.stime);
        Py(i1)=mean(exp(nu(i1).*(f1(:,si)+f2) - (yt(i1)-lik.stime(si)).*la1(:,si).*eta2 ...
                     - sum(sd.*la1(:,1:(si-1)),2).*eta2));
      else

        sb=sum(bsxfun(@gt,yt(i1,1),lik.stime),2);
        se=sum(bsxfun(@gt,yt(i1,2),lik.stime),2);

        if sb==0
          Py(i1) = mean(exp(nu(i1).*(f1(:,se)+f2) - (yt(i1,2)-lik.stime(se)).*la1(:,se).*eta2 - sum(sd.*la1(:,1:(se-1)),2).*eta2));
        else
          if se==sb
            Py(i1) = mean(exp(nu(i1).*(f1(:,se)+f2) - (yt(i1,2)-yt(i1,1)).*la1(:,se).*eta2));
          else
            Py(i1) = mean(exp(nu(i1).*(f1(:,se)+f2) - (yt(i1,2)-lik.stime(se)).*la1(:,se).*eta2 - sum(sd.*la1(:,(sb+1):(se-1)),2).*eta2 - (lik.stime(sb+1)-yt(i1,1)).*la1(:,sb).*eta2));
          end
        end
      end
    end
    Ey = [];
    Vary = [];
    lpyt=log(Py);

    %     % Evaluate Ey and Vary
%     for i1=1:length(Ef)
%       %%% With quadrature
%       myy_i = Ef(i1);
%       sigm_i = sqrt(Varf(i1));
%       minf=myy_i-6*sigm_i;
%       maxf=myy_i+6*sigm_i;
%
%       F = @(f) exp(log(avgE(i1))+f+norm_lpdf(f,myy_i,sigm_i));
%       Ey(i1) = quadgk(F,minf,maxf);
%
%       F2 = @(f) exp(log(avgE(i1).*exp(f)+((avgE(i1).*exp(f)).^2/r))+norm_lpdf(f,myy_i,sigm_i));
%       EVary(i1) = quadgk(F2,minf,maxf);
%
%       F3 = @(f) exp(2*log(avgE(i1))+2*f+norm_lpdf(f,myy_i,sigm_i));
%       VarEy(i1) = quadgk(F3,minf,maxf) - Ey(i1).^2;
%     end
%     Vary = EVary + VarEy;
%
%     % Evaluate the posterior predictive densities of the given observations
%     if nargout > 2
%       for i1=1:length(Ef)
%         % get a function handle of the likelihood times posterior
%         % (likelihood * posterior = Negative-binomial * Gaussian)
%         % and useful integration limits
%         [pdf,minf,maxf]=init_coxph_norm(...
%           yt(i1),Ef(i1),Varf(i1),avgE(i1),r);
%         % integrate over the f to get posterior predictive distribution
%         Py(i1) = quadgk(pdf, minf, maxf);
%       end
%     end
  end

  function [logM_0, m_1, sigm2hati1] = coxph_tiltedMoments(sigm2_i, myy_i, nutmp, sd, c1)

  integrand = @(f) exp(-c1.*exp(f).*sd + nutmp*(f+log(c1)) - log(sigm2_i)/2 - log(2*pi)/2 - 0.5*(f-myy_i).^2./sigm2_i);
  RTOL = 1.e-6;
  ATOL = 1.e-10;
  minf=myy_i+sqrt(sigm2_i)*(-6);
  maxf=myy_i+sqrt(sigm2_i)*(6);

  [m_0, m_1, m_2] = quad_moments(integrand, minf, maxf, RTOL, ATOL);
  sigm2hati1 = m_2 - m_1.^2;

  % If the second central moment is less than cavity variance
  % integrate more precisely. Theoretically for log-concave
  % likelihood should be sigm2hati1 < sigm2_i.

  if sigm2hati1 >= sigm2_i
    ATOL = ATOL.^2;
    RTOL = RTOL.^2;
    [m_0, m_1, m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
    sigm2hati1 = m_2 - m_1.^2;
    if sigm2hati1 >= sigm2_i
      error('lik_poisson_tilted_moments: sigm2hati1 >= sigm2_i');
    end
  end
  logM_0 = log(m_0);
  end

  function [df,minf,maxf] = init_coxph_norm(yy,myy_i,sigm2_i,avgE,r)
  %INIT_COXPH_NORM
  %
  %  Description
  %    Return function handle to a function evaluating
  %    Negative-Binomial * Gaussian which is used for evaluating
  %    (likelihood * cavity) or (likelihood * posterior) Return
  %    also useful limits for integration. This is private function
  %    for lik_coxph.
  %
  %  See also
  %    LIK_COXPH_TILTEDMOMENTS, LIK_COXPH_SITEDERIV,
  %    LIK_COXPH_PREDY

  % avoid repetitive evaluation of constant part
    ldconst = -gammaln(r)-gammaln(yy+1)+gammaln(r+yy)...
              - log(sigm2_i)/2 - log(2*pi)/2;
    % Create function handle for the function to be integrated
    df = @coxph_norm;
    % use log to avoid underflow, and derivates for faster search
    ld = @log_coxph_norm;
    ldg = @log_coxph_norm_g;
    ldg2 = @log_coxph_norm_g2;

    % Set the limits for integration
    % Negative-binomial likelihood is log-concave so the coxph_norm
    % function is unimodal, which makes things easier
    if yy==0
      % with yy==0, the mode of the likelihood is not defined
      % use the mode of the Gaussian (cavity or posterior) as a first guess
      modef = myy_i;
    else
      % use precision weighted mean of the Gaussian approximation
      % of the Negative-Binomial likelihood and Gaussian
      mu=log(yy/avgE);
      s2=(yy+r)./(yy.*r);
      modef = (myy_i/sigm2_i + mu/s2)/(1/sigm2_i + 1/s2);
    end
    % find the mode of the integrand using Newton iterations
    % few iterations is enough, since the first guess in the right direction
    niter=4;       % number of Newton iterations
    mindelta=1e-6; % tolerance in stopping Newton iterations
    for ni=1:niter
      g=ldg(modef);
      h=ldg2(modef);
      delta=-g/h;
      modef=modef+delta;
      if abs(delta)<mindelta
        break
      end
    end
    % integrand limits based on Gaussian approximation at mode
    modes=sqrt(-1/h);
    minf=modef-8*modes;
    maxf=modef+8*modes;
    modeld=ld(modef);
    iter=0;
    % check that density at end points is low enough
    lddiff=20; % min difference in log-density between mode and end-points
    minld=ld(minf);
    step=1;
    while minld>(modeld-lddiff)
      minf=minf-step*modes;
      minld=ld(minf);
      iter=iter+1;
      step=step*2;
      if iter>100
        error(['lik_coxph -> init_coxph_norm: ' ...
               'integration interval minimun not found ' ...
               'even after looking hard!'])
      end
    end
    maxld=ld(maxf);
    step=1;
    while maxld>(modeld-lddiff)
      maxf=maxf+step*modes;
      maxld=ld(maxf);
      iter=iter+1;
      step=step*2;
      if iter>100
        error(['lik_coxph -> init_coxph_norm: ' ...
               'integration interval maximun not found ' ...
               'even after looking hard!'])
      end
    end

%     while minld>(modeld-lddiff)
%       minf=minf-modes;
%       minld=ld(minf);
%       iter=iter+1;
%       if iter>100
%         error(['lik_coxph -> init_coxph_norm: ' ...
%                'integration interval minimun not found ' ...
%                'even after looking hard!'])
%       end
%     end
%     maxld=ld(maxf);
%     while maxld>(modeld-lddiff)
%       maxf=maxf+modes;
%       maxld=ld(maxf);
%       iter=iter+1;
%       if iter>100
%         error(['lik_coxph -> init_coxph_norm: ' ...
%                'integration interval maximum not found ' ...
%                'even after looking hard!'])
%       end
%
%     end

    function integrand = coxph_norm(f)
    % Negative-binomial * Gaussian
      mu = avgE.*exp(f);
      integrand = exp(ldconst ...
                      +yy.*(log(mu)-log(r+mu))+r.*(log(r)-log(r+mu)) ...
                      -0.5*(f-myy_i).^2./sigm2_i);
    end

    function log_int = log_coxph_norm(f)
    % log(Negative-binomial * Gaussian)
    % log_coxph_norm is used to avoid underflow when searching
    % integration interval
      mu = avgE.*exp(f);
      log_int = ldconst...
                +yy.*(log(mu)-log(r+mu))+r.*(log(r)-log(r+mu))...
                -0.5*(f-myy_i).^2./sigm2_i;
    end

    function g = log_coxph_norm_g(f)
    % d/df log(Negative-binomial * Gaussian)
    % derivative of log_coxph_norm
      mu = avgE.*exp(f);
      g = -(r.*(mu - yy))./(mu.*(mu + r)).*mu ...
          + (myy_i - f)./sigm2_i;
    end

    function g2 = log_coxph_norm_g2(f)
    % d^2/df^2 log(Negative-binomial * Gaussian)
    % second derivate of log_coxph_norm
      mu = avgE.*exp(f);
      g2 = -(r*(r + yy))/(mu + r)^2.*mu ...
           -1/sigm2_i;
    end

  end

  function p = lik_coxph_invlink(lik, f, z)
  %LIK_COXPH_INVLINK  Returns values of inverse link function
  %
  %  Description
  %    P = LIK_COXPH_INVLINK(LIK, F) takes a likelihood structure LIK and
  %    latent values F and returns the values of inverse link function P.
  %    This subfunction is needed when using function gp_predprcty.
  %
  %     See also
  %     LIK_COXPH_LL, LIK_COXPH_PREDY

    p = exp(f);
  end

  function [cdf,Ey,Vary] = lik_coxph_predcdf(lik,Ef,Covf,yt)
  %LIK_LOGLOGISTIC_PREDCDF  Returns the predictive cdf evaluated at yt
  %
  %  Description
  %    CDF = LIK_LOGLOGISTIC_PREDCDF(LIK, EF, VARF, YT)
  %    Returns the predictive cdf evaluated at YT given likelihood
  %    structure LIK, posterior mean EF and posterior Variance VARF
  %    of the latent variable. This subfunction is needed when using
  %    functions gp_predcdf or gp_kfcv_cdf.
  %
  %  See also
  %    GP_PREDCDF

    if isfield(lik, 'stratificationVariables')
      ntime=size(lik.xtime,1)*lik.n_u;
    else
      ntime=size(lik.xtime,1);
    end
    Ef1 = Ef(1:ntime); Ef(1:ntime) = []; Ef2 = Ef;
    nsamps = 10000;
    sd=lik.stime(2)-lik.stime(1);
    Sigm_tmp=Covf;
    [nn1,nn2,cc]=size(Sigm_tmp);
    if cc==1 && nn1~=nn2
      f_star=bsxfun(@plus,[Ef1;Ef2]', ...
        bsxfun(@times,sqrt(Sigm_tmp'),randn(nsamps,ntime+size(Ef2,1))));
    else
      Sigm_tmp=(Sigm_tmp+Sigm_tmp')./2;
      % f_star=mvnrnd(Ef1, Sigm_tmp(1:ntime,1:ntime), nsamps);
      f_star=mvnrnd([Ef1;Ef2], Sigm_tmp, nsamps);
    end

    f1=f_star(:,1:ntime);
    f2=f_star(:,(ntime+1):end);

    la1=exp(f1);
    eta2=exp(f2);

    mST=zeros(size(eta2,2),1);
    if size(yt,2) == 1
      % Integrate from zero to yt
      cumsumtmp=cumsum(la1'*sd)';
      %   t=binsgeq(gp.lik.xtime,yt(i));
      for i1=1:size(eta2,2)
        Stime=exp(-bsxfun(@times,cumsumtmp,eta2(:,i1)));
        mStime=mean(Stime);
        % for i=1:size(yt,1)
        mST(i1)=mStime(binsgeq(lik.xtime,yt(i1)));
        %end
      end
      cdf = 1- mST;

    else
      error('Size(yt,2) ~= 1');
    end

  end

  function reclik = lik_coxph_recappend(reclik, ri, lik)
  %RECAPPEND  Append the parameters to the record
  %
  %  Description
  %    RECLIK = GPCF_COXPH_RECAPPEND(RECLIK, RI, LIK) takes a
  %    likelihood record structure RECLIK, record index RI and
  %    likelihood structure LIK with the current MCMC samples of
  %    the parameters. Returns RECLIK which contains all the old
  %    samples and the current samples from LIK. This subfunction
  %    is needed when using MCMC sampling (gp_mc).
  %
  %  See also
  %    GP_MC

  % Initialize record
    if nargin == 2
      reclik=ri;

      % Set the function handles
      reclik.fh.pak = @lik_coxph_pak;
      reclik.fh.unpak = @lik_coxph_unpak;
      reclik.fh.lp = @lik_coxph_lp;
      reclik.fh.lpg = @lik_coxph_lpg;
      reclik.fh.ll = @lik_coxph_ll;
      reclik.fh.llg = @lik_coxph_llg;
      reclik.fh.llg2 = @lik_coxph_llg2;
      reclik.fh.llg3 = @lik_coxph_llg3;
      reclik.fh.tiltedMoments = @lik_coxph_tiltedMoments;
      reclik.fh.predy = @lik_coxph_predy;
      reclik.fh.invlink = @lik_coxph_invlink;
      reclik.fh.recappend = @lik_coxph_recappend;
      return
    else

      reclik.xtime=lik.xtime;
      reclik.stime=lik.stime;

      if isfield(lik, 'stratificationVariables')
        reclik.stratificationVariables=lik.stratificationVariables;
        if isfield(lik,removeStratificationVariables)
          reclik.removeStratificationVariables=lik.removeStratificationVariables;
        end
      end
    end
  end

end


text:
function J=SimFisherMatrix(obj,Prot,x,variables,sigma)
% Alexander, D.C., 2008. A general framework for experiment design in diffusion MRI and its application in measuring direct tissue-microstructure features. Magn. Reson. Med. 60, 439?448.
obj.Prot.(obj.MRIinputs{1}).Mat = Prot;
if ~exist('variables','var'), variables=1:5; end
if ~exist('sigma','var'), sigma=0.1; end

J=zeros(max(variables));
if ismethod(obj,'equation_x')
    S_nominal = obj.equation_x(x);
else
    S_nominal = obj.equation(x);
end
% Gaussian Noise
for i=variables
    for j=variables
        Xi = x;
        Xistep = max(1e-10,Xi(i) / 100);
        Xi(i) = Xi(i) + Xistep;
        Xj = x;
        Xjstep = max(1e-10,Xj(j) / 100);
        Xj(j) = Xj(j) + Xjstep;
        if ismethod(obj,'equation_x')
            J(i,j) = sum( 1./sigma.^2 .* (S_nominal - obj.equation_x(Xi))/Xistep .* (S_nominal - obj.equation_x(Xj))/Xjstep);
        else
            J(i,j) = sum( 1./sigma.^2 .* (S_nominal - obj.equation(Xi))/Xistep .* (S_nominal - obj.equation(Xj))/Xjstep);
        end
    end
end
J=J(variables,variables);

% % Rician Noise
% for i=1:variables
%     for j=variables
%         Xi=x;
%         Xi(i)=Xi(i)+Xi(i)/100;
%         Xj=x;
%         Xj(j)=Xj(j)+Xj(j)/100;
%
%         A = scd_model_GPD_composite(x,Ax);
%
%         Z=zeros(size(scheme,1),1);
%         for k=1:size(scheme,1)
%             Z(k) = integral(@(a) A(k)^2*besseli(1,A(k)*a/sigma^2)^2*besseli(0,A(k)*a/sigma^2)^-2  *a/sigma^2.*besseli(0,A(k)*a/sigma^2).*exp(-(A(k)^2+a^2)/2*sigma^2)  ,0,inf);
%         end
%
%         dAi = (scd_model_GPD_composite(x,Ax)-scd_model_GPD_composite(Xi,Ax))/(x(i)/100);
%         dAj = (scd_model_GPD_composite(x,Ax)-scd_model_GPD_composite(Xj,Ax))/(x(j)/100);
%
%         J(i,j)= sum(1/sigma^4.*(dAi.*dAj).*(Z-A.^2));
%     end
% end
% J=J(variables,variables);
% if J<0, error('bug'); end

end


text:
function value = r8_invgam_sample ( beta, alpha )

%*****************************************************************************80
%
%% R8_INVGAM_SAMPLE samples an inverse gamma distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    Original FORTRAN90 version by Guannan Zhang.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real BETA, the rate parameter.
%    0.0 < BETA.
%
%    Input, real ALPHA, the shape parameter.
%    0.0 < ALPHA.
%
%    Output, real VALUE, a sample value.
%
  value = r8_gamma_sample ( beta, alpha );

  if ( value ~= 0.0 )
    value = 1.0 / value;
  end

  return
end


text:
function [cRow, cCol, ra, rb, phi] = gauss2ellipse( mu, C, rad )
% Creates an ellipse representing the 2D Gaussian distribution.
%
% Creates an ellipse representing the 2D Gaussian distribution with mean mu
% and covariance matrix C.  Returns 5 parameters that specify the ellipse.
%
% USAGE
%  [cRow, cCol, ra, rb, phi] = gauss2ellipse( mu, C, [rad] )
%
% INPUTS
%  mu      - 1x2 vector representing the center of the ellipse
%  C       - 2x2 cov matrix
%  rad     - [2] Number of std to create the ellipse to
%
% OUTPUTS
%  cRow    - the row location of the center of the ellipse
%  cCol    - the column location of the center of the ellipse
%  ra      - semi-major axis length (in pixels) of the ellipse
%  rb      - semi-minor axis length (in pixels) of the ellipse
%  phi     - rotation angle (radians) of semimajor axis from x-axis
%
% EXAMPLE
%  [cRow, cCol, ra, rb, phi] = gauss2ellipse( [5 5], [1 0; .5 2] )
%  plotEllipse( cRow, cCol, ra, rb, phi );
%
% See also PLOTELLIPSE, PLOTGAUSSELLIPSES, MASKELLIPSE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if (nargin<3 || isempty(rad) ); rad=2; end;

% error check
if (~all(size(mu)==[1,2]) || ~all(size(C)==[2,2]))
  error('Works only for 2D Gaussians'); end

% decompose using SVD
[~,D,R] = svd(C);
normstd = sqrt( diag( D ) );

% get angle of rotation (in row/column format)
phi = acos(R(1,1));
if (R(2,1) < 0); phi = 2*pi - phi; end
phi = pi/2 - phi;

% get ellipse radii
ra = rad*normstd(1);
rb = rad*normstd(2);

% center of ellipse
cRow = mu(1);
cCol = mu(2);

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [annotation, img] = LMimpad(annotation, img, PADSIZE, PADVAL)
% [annotation, img] = LMimpad(annotation, img, PADSIZE, PADVAL)
%
% [annotation, img] = LMimpad(annotation, img, [256 256], 0)
% PADSIZE = nrows x ncols (final image size)
% PADVAL = value for the border

[nrows, ncols, cols]=size(img);

if PADSIZE(1)<nrows || PADSIZE(2)<ncols
    error('ERROR: image is larger than target size. Use LMimcrop instead.')
end

Dy = fix((PADSIZE(1)-nrows)/2);
Dx = fix((PADSIZE(2)-ncols)/2);

img = [repmat(PADVAL, [PADSIZE(1) Dx cols]) ...
    [repmat(PADVAL, [Dy ncols cols]); img; repmat(PADVAL, [PADSIZE(1)-nrows-Dy ncols cols])] ...
    repmat(PADVAL, [PADSIZE(1) PADSIZE(2)-ncols-Dx cols])];

% Change the polygon coordinates
if isfield(annotation, 'object')
    Nobjects = length(annotation.object); n=0;
    for i = 1:Nobjects
        [x,y] = getLMpolygon(annotation.object(i).polygon);
        x = round(x + Dx);
        y = round(y + Dy);
        annotation.object(i).polygon = setLMpolygon(x,y);

%         Npoints = length(annotation.object(i).polygon.pt);
%         for j = 1:Npoints
%             x=str2num(annotation.object(i).polygon.pt(j).x);
%             y=str2num(annotation.object(i).polygon.pt(j).y);
%
%             x = round(x + Dx);
%             y = round(y + Dy);
%
%             annotation.object(i).polygon.pt(j).x = num2str(x);
%             annotation.object(i).polygon.pt(j).y = num2str(y);
%         end
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
% This is a demo for segmentation using local gaussian distribution (LGD)
% fitting energy
%
% Reference: <Li Wang, Lei He, Arabinda Mishra, Chunming Li.
% Active Contours Driven by Local Gaussian Distribution Fitting Energy.
% Signal Processing, 89(12), 2009,p. 2435-2447>
%
% Please DO NOT distribute this code to anybody.
% Copyright (c) by Li Wang
%
% Author:       Li Wang
% E-mail:       li_wang@med.unc.edu
% URL:          http://www.unc.edu/~liwa/
%
% 2010-01-02 PM

clc;clear all;close all;

Img=imread('2.bmp');
Img = double(Img(:,:,1));

NumIter = 1000; %iterations
timestep=0.1; %time step
mu=0.1/timestep;% level set regularization term, please refer to "Chunming Li and et al. Level Set Evolution Without Re-initialization: A New Variational Formulation, CVPR 2005"
sigma = 5;%size of kernel
epsilon = 1;
c0 = 2; % the constant value
lambda1=1.05;%outer weight, please refer to "Chunming Li and et al,  Minimization of Region-Scalable Fitting Energy for Image Segmentation, IEEE Trans. Image Processing, vol. 17 (10), pp. 1940-1949, 2008"
lambda2=1.0;%inner weight
%if lambda1>lambda2; tend to inflate
%if lambda1<lambda2; tend to deflate
nu = 0.001*255*255;%length term
alf = 20;%data term weight

figure,imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal
[Height Wide] = size(Img);
[xx yy] = meshgrid(1:Wide,1:Height);
phi = (sqrt(((xx - 64).^2 + (yy - 65).^2 )) - 17);
phi = sign(phi).*c0;

Ksigma=fspecial('gaussian',round(2*sigma)*2 + 1,sigma); %  kernel
ONE=ones(size(Img));
KONE = imfilter(ONE,Ksigma,'replicate');
KI = imfilter(Img,Ksigma,'replicate');
KI2 = imfilter(Img.^2,Ksigma,'replicate');

figure,imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal,
hold on,[c,h] = contour(phi,[0 0],'r','linewidth',1); hold off
pause(0.5)

tic
for iter = 1:NumIter
    phi =evolution_LGD(Img,phi,epsilon,Ksigma,KONE,KI,KI2,mu,nu,lambda1,lambda2,timestep,alf);

    if(mod(iter,100) == 0)
        figure(2),
        imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal,title(num2str(iter))
        hold on,[c,h] = contour(phi,[0 0],'r','linewidth',1); hold off
        pause(0.02);
    end

end
toc

author: Sable

repo: mcbench-benchmarks


text:
function [w, state] = adam(w, state, grad, opts, lr)
%SGD
%   Example SGD solver, with momentum, for use with CNN_TRAIN and
%   CNN_TRAIN_DAG.
%
%   The convergence of SGD depends heavily on the learning rate (set in the
%   options for CNN_TRAIN and CNN_TRAIN_DAG).
%
%   If called without any input argument, returns the default options
%   structure.
%
%   Solver options: (opts.train.solverOpts)
%
%   `momentum`:: 0.9
%      Parameter for Momentum SGD; set to 0 for standard SGD.
%
%   Note: for backwards compatibility, the parameter can also be set in
%   opts.train.momentum.

% Copyright (C) 2016 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin == 0 % Return the default solver options
  w = struct('momentum', 0.9);
  return;
end
if isempty(state)
  state.m = 0;
  state.v = 0;
  state.momentum1t = 1;
  state.momentum2t = 1;
%   momentum = 0 ;
end

opts.momentum1 = 0.9;
opts.momentum2 = 0.999;

state.momentum1t = state.momentum1t*opts.momentum1;
state.momentum2t = state.momentum2t*opts.momentum2;
state.m = opts.momentum1*state.m + (1-opts.momentum1)*grad;
state.v = opts.momentum2*state.v + (1-opts.momentum2)*grad.^2;

momentum = - state.m/(1-state.momentum1t);
momentum = momentum ./ (sqrt(state.v/(1-state.momentum2t)) + 10^(-8));
w = w + lr * momentum ;

repo: CREST-Release


text:
function value = c8_acosh ( z )

%*****************************************************************************80
%
%% C8_ACOSH evaluates the inverse hyperbolic cosine of a C8.
%
%  Discussion:
%
%    A C8 is a complex value.
%
%    Here we use the relationship:
%
%      C8_ACOSH ( Z ) = i * C8_ACOS ( Z ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex Z, the argument.
%
%    Output, complex VALUE, the function value.
%
  value = i * c8_acos ( z );

  return
end


text:
function OUT = Num2NaN(DATA)
% =======================================================================
% Substitute 123456789 values with NaN
% =======================================================================
% OUT = Num2NaN(DATA)
% -----------------------------------------------------------------------
% INPUT
%	- DATA: a (m x n) matrix with some 123456789
% -----------------------------------------------------------------------
% OUTPUT
%	- OUT: a (m x n) matrix with NaN in place of 123456789
% -----------------------------------------------------------------------
% EXAMPLE
%   x = [123456789 2; 123456789 4; 5 6; 7 8; 9 10];
%   OUT = Num2NaN(x)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

OUT=DATA;
for i=1:size(DATA,1)
    for j=1:size(DATA,2)
        if DATA(i,j) == 123456789
            OUT(i,j) = NaN;
        end
    end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function hewma2_plot(Zpop,mu,cp,cp_ind,ooc_indices,maxtime,tthresh,tt,tm,q,Wm,Zm,sterr,varargin)
% :Usage:
% ::
%
%     hewma2_plot(Zpop,mu,cp,cp_ind,ooc_indices,maxtime,tthresh,tt,tm,q,Wm,Zm,sterr,varargin)
%
% see hewma2.m

cp = round(cp);

dogroupplot = 0;
tagname  = 'hewma2display';

if length(varargin) > 0
    dogroupplot = 1;
    tagname  = 'hewma2diffdisplay';
end

f1 = findobj('Tag', tagname);
if isempty(f1)
    f1 = tor_fig(2,2);
    set(f1, 'Tag', tagname);
else
    figure(f1)
    for i = 1:4, subplot(2, 2, i); cla; end
end

ZZ = Zpop-mu;

% ----------------------------------------------------------
% *
% * PLOT 1: GROUP AVERAGE (OR DIFFERENCE)
% *
% ----------------------------------------------------------
subplot(2, 2, 1);

hh(1) = plot(ZZ,'k','LineWidth',2);

% baseline box
% ------------------------------
    yl = get(gca,'YLim');
    cl = tthresh .* mean(sterr);
    yl(1) = min([yl(1)-1 -cl-1]); yl(2) = max([yl(2)+1 cl+1]);

    fill([1 tt tt 1],[yl(1) yl(1) yl(2) yl(2)],[.85 .85 .85]);
    text(5,min(ZZ),'Baseline','FontSize',16);

% group Z stat
% ------------------------------
hh(1) = plot(ZZ,'k','LineWidth',2);
hold on

 % change points
% ------------------------------
    if ~isnan(cp)
        plot(cp,ZZ(cp),'bo','MarkerSize',8,'MarkerFaceColor','g');
        plot([cp cp],get(gca,'YLim'),'g','LineWidth',2);
    end

    % individual change points
    if ~isempty(cp_ind)
        if ~any(isnan(cp_ind))
            plot(cp_ind,ZZ(cp_ind),'go','MarkerSize',6);
        end
    end

% OOC points in red
% ------------------------------
    % all ooc points
    ZZtmp = NaN * zeros(size(ZZ)); ZZtmp(ooc_indices) = ZZ(ooc_indices);
    plot(ZZtmp,'r','LineWidth',2);

% max t value in blue
% ------------------------------

    if ~isnan(maxtime)
        plot(maxtime,ZZ(maxtime),'go','MarkerSize',8,'MarkerFaceColor','b');
    end

 % standard error
% ------------------------------
    fill_around_line(ZZ,sterr,'k');

    title('Group activation');
    if length(varargin) > 0, title('Difference between groups'); end

 % control limits
% ------------------------------
    cl = tthresh .* mean(sterr);
    hh(2) = plot(get(gca,'XLim'),[cl cl],'k--');
    plot(get(gca,'XLim'),[-cl -cl],'k--');
    legend(hh,'Group mean','Control limit (approximate)');

% ----------------------------------------------------------
% *
% * PLOT 2: NULL HYPOTHESIS DISTRIBUTION
% *
% ----------------------------------------------------------

    subplot(2,2,2);
    hist(q,50); hold on; plot([abs(tm) abs(tm)],get(gca,'YLim'),'k-','LineWidth',2);
    title('Null-hypothesis max t-value and observed max t-value');

% ----------------------------------------------------------
% *
% * PLOT 3: SUBJECT WEIGHTS
% *
% ----------------------------------------------------------

    subplot(2,2,3);
    n = 1:length(Wm);
    plot([n;n],[Wm';Wm'],'o','LineWidth',2);
    title('Case weights');

% ----------------------------------------------------------
% *
% * PLOT 4: INDIVIDUAL SUBJECTS (OR GROUPS)
% *
% ----------------------------------------------------------

    subplot(2,2,4);

    % group contrast, or individuals
    if dogroupplot
        gfits = varargin{1};
        plot(gfits');
        title('Groups');
        legend({'Low' 'High'});
        yval = max(gfits(:)) + .01*max(gfits(:));
        %yval = repmat(yval,1,length(ZZtmp));
        ZZtmp(~isnan(ZZtmp)) = yval;
        hold on;

        plot(ZZtmp,'r','LineWidth',3);

        if ~isnan(cp)
            plot([cp cp],get(gca,'YLim'),'g','LineWidth',2);
        end

        % baseline box
        yl = get(gca,'YLim');
        fill([1 tt tt 1],[yl(1) yl(1) yl(2) yl(2)],[.85 .85 .85]);
        text(5,min(ZZ),'Baseline','FontSize',16);

    else
        plot(Zm');
        title('Individual cases');
        drawnow
    end

    return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function h = gaussianNoise3dPlot(noise, plotType, CX, CY, CZ, CZVar, varargin)

% GAUSSIANNOISE3DPLOT Draws a 3D or contour plot for the GAUSSIAN noise model.
% FORMAT
% DESC draws a 3D or contour plot for the Gaussian noise model.
% ARG noise : the noise structure for which the plot is required.
% ARG plotType : string containing the name of the plotting function (for example mesh, contour).
% ARG X : the input X data in the form of a 'mesh' matrix.
% ARG Y : the input Y data in the form of a 'mesh' matrix.
% ARG mu : the input mean in the form of a 'mesh' matrix.
% ARG varSigma : the input variance in the form of a 'mesh' matrix.
% ARG P1, P2, P3 ... : optional additional arguments for the given plot type.
% RETURN h : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : gaussianNoiseParamInit, noise3dPlot,
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

CZ = (CZ+noise.bias);
fhandle = str2func(plotType);
h = fhandle(CX, CY, CZ, varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef MOEADDE_F4 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing MOEA/D-DE

%------------------------------- Reference --------------------------------
% H. Li and Q. Zhang, Multiobjective optimization problems with complicated
% Pareto sets, MOEA/D and NSGA-II, IEEE Transactions on Evolutionary
% Computation, 2009, 13(2): 284-302.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = [0,-ones(1,obj.D-1)];
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            J1 = 3 : 2 : obj.D;
            J2 = 2 : 2 : obj.D;
            PopObj(:,1) = X(:,1)         + 2*mean((X(:,J1)-0.8*repmat(X(:,1),1,length(J1)).*cos(repmat(2*pi*X(:,1),1,length(J1))+repmat(J1*pi/obj.D/3,size(X,1),1))).^2,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + 2*mean((X(:,J2)-0.8*repmat(X(:,1),1,length(J2)).*sin(repmat(6*pi*X(:,1),1,length(J2))+repmat(J2*pi/obj.D,size(X,1),1))).^2,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [y,deriv] = compose_mv(outer,inner,x)
% COMPOSE_MV is an MV2DF (see MV2DF_API_DEFINITION.readme) which represents
% the combination of two functions. If 'outer' is an MV2DF for a function
% g() and 'inner' for a function f(), then this MV2DF represents g(f(x)).

%feature scopedaccelenablement off

if nargin==0
    test_this();
    return;
end

if isempty(x)
    y = @(w)compose_mv(outer,inner,w);
    return;
end

if isa(x,'function_handle')
    fh = compose_mv(outer,inner,[]);  % fh =@(x) outer(inner(x))
    y = compose_mv(fh,x,[]);          %  y =@(w) outer(inner(x(w)))
    return;
end

% if ~isa(outer,'function_handle')
%     outer = const_mv2df([],outer);
% end
% if ~isa(inner,'function_handle')
%     inner = const_mv2df([],inner);
% end

if nargout==1
    y = outer(inner(x));
    return;
end

[y1,deriv1] = inner(x);
[y,deriv2] = outer(y1);
deriv = @(g3) deriv_this(deriv1,deriv2,g3);

function [g,hess,linear] = deriv_this(deriv1,deriv2,g3)

if nargout==1
    g = deriv1(deriv2(g3));
    return;
end

[g2,hess2,lin2] = deriv2(g3);
[g,hess1,lin1] = deriv1(g2);

hess =@(d) hess_this(deriv1,hess1,hess2,lin1,lin2,d);

linear = lin1 && lin2;

function [h,Jv] = hess_this(deriv1,hess1,hess2,lin1,lin2,d)

if nargout==1
    if ~lin2
        [h1,Jv1] = hess1(d);
        h2 = hess2(Jv1);
        h2 = deriv1(h2);
    elseif ~lin1
        h1 = hess1(d);
    end
else
    [h1,Jv1] = hess1(d);
    [h2,Jv] = hess2(Jv1);
    if ~lin2
        h2 = deriv1(h2);
    end
end

if lin1 && lin2
    h=[];
elseif (~lin1) && (~lin2)
    h = h1+h2;
elseif lin1
    h = h2;
else % if lin2
    h = h1;
end

function test_this()

fprintf('-------------- Test 1 ------------------------\n');
fprintf('Composition g(f(w)): f() is non-linear and g() is non-linear:\n');
A = randn(4,5);
B = randn(5,4);

w = [A(:);B(:)];

f = @(w) gemm(w,4,5,4);
g1 = gemm(f,2,4,2);

test_MV2DF(g1,w);
fprintf('--------------------------------------\n\n');

fprintf('-------------- Test 2 ------------------------\n');
fprintf('Composition g(f(w)): f() is linear and g() is non-linear:\n');
A = randn(4,5);
B = randn(5,4);

w = [A(:);B(:)];

T = randn(40);
f = @(w) linTrans(w,@(x)T*x,@(y)T.'*y);
g2 = gemm(f,4,5,4);

test_MV2DF(g2,w);
fprintf('--------------------------------------\n\n');

fprintf('-------------- Test 3 ------------------------\n');
fprintf('Composition g(f(w)): f() is non-linear and g() is linear:\n');
A = randn(4,5);
B = randn(5,4);

w = [A(:);B(:)];
f = @(w) gemm(w,4,5,4);

T = randn(16);
g3 = linTrans(f,@(x)T*x,@(y)T.'*y);

test_MV2DF(g3,w);
fprintf('--------------------------------------\n\n');

fprintf('-------------- Test 4 ------------------------\n');
fprintf('Composition g(f(w)): f() is linear and g() is linear:\n');
w = randn(10,1);
T1 = randn(11,10);
f = @(w) linTrans(w,@(x)T1*x,@(y)T1.'*y);

T2 = randn(5,11);
g4 = linTrans(f,@(x)T2*x,@(y)T2.'*y);

test_MV2DF(g4,w);
fprintf('--------------------------------------\n\n');

author: nesl


text:
function [val, sts] = cfg_eval_valedit(varargin)
% Helper function to evaluate GUI inputs in MATLAB workspace
% FORMAT [val, sts] = cfg_eval_valedit(str)
% Evaluates GUI inputs in MATLAB 'base' workspace. Results are returned
% in val. Expressions in str can be either a pure rhs argument, or a set
% of commands that assign to a workspace variable named 'val'. If
% unsuccessful, sts is false and a message window is displayed.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_eval_valedit.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $';

val = [];
sts = false;
try
    % 1st, try to convert into numeric matrix without evaluation
    % This converts expressions like '1 -1' into [1 -1] instead of
    % evaluating them
    [val, sts] = str2num(varargin{1}); %#ok<ST2NM>
    if ~sts
        % try to evaluate str as rvalue
        val = evalin('base', varargin{1});
        sts = true;
    end
catch
    everr = lasterror;
    if strcmp(everr.identifier, 'MATLAB:m_invalid_lhs_of_assignment')
        try
            evalin('base', varargin{1});
            val = evalin('base','val');
            % test if val variable exists
            if ~exist('val','var')
                cfg_message('cfg_ui:local_eval_valedit:noval','No variable ''val'' assigned.');
            end;
            sts = true;
        catch
            sts = false;
            val = [];
            everr = lasterror;
            msgbox(everr.message,'Evaluation error','modal');
        end;
    end;
end;

author: spm


text:
% Author: Guanghan Ning, Feb 2017
clear all;
clc;
close all;
addpath('../utils_eval_lsp/');
addpath('../eval_LSP/');
mat_folder_path = '../dataset_lsp/results_cropped/quantitative/';
img_folder_path = '../../dataset/LEEDS/lsp_dataset/images/';

scale_in_cpp = load('../scale_in_cpp_LSP.mat');
scale_in_cpp = scale_in_cpp.output_scale;

pred = zeros(2, 14, 1000);
num_test_imgs = 1000;

for img_id = 1:num_test_imgs

    base_name = strcat('im', num2str(img_id + 1000));
    img_path = strcat(img_folder_path, '/', base_name, '.jpg');
    img_info = imfinfo(img_path);
    org_img_size = [img_info.Width, img_info.Height];

    [x_preds_raw, y_preds_raw] = load_pred_result(mat_folder_path, base_name);

    for joint_id = 1:14 %joint_id in our joint order

        % Convert joint to lsp-defined order for evaluation
        id = convert_joint_order(joint_id);  % id in LSP order

        % Convert predictions from 256 to image-size for evaluation
        pred_point_raw = [x_preds_raw(joint_id), y_preds_raw(joint_id)];
        [x_pred, y_pred]= convert_scale_cropped(pred_point_raw, org_img_size, scale_in_cpp(img_id));

        pred(1, id, img_id) = x_pred;
        pred(2, id, img_id) = y_pred;
    end
end
save('../eval_LSP/pred/ning17iccv/pred_keypoints_lsp.mat','pred');

%% Call Official evaluatePCKh() function to process the evaluation
evaluatePCK(33,'PC', true); % 33 is my model index, using person-centric (PC)

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function s=dec2decbin(d,n)
%DEC2BIN Internal function generate binary matrices

[f,e]=log2(max(d)); % How many digits do we need to represent the numbers?
s=rem(floor(d(:)*pow2(1-max(n,e):0)),2);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP