text:
function pixeliZedScratchPad = plotCoordinatesToImagePixels(x_coordinates, y_coordinates, imageSize)
%PLOTCOORDINATESTOIMAGEPIXELS

numberOfPoints = length(x_coordinates);
% make into an image by pixelization the coordinates 1st
x_min = min(x_coordinates);
x_max = max(x_coordinates);
y_min = min(y_coordinates);
y_max = max(y_coordinates);

% [use linspace]
% numberOfGratings = 200;
numberOfGratings = imageSize;
unpixelizedRange_x = linspace(x_min, x_max, numberOfGratings);
unpixelizedRange_y = linspace(y_min, y_max, numberOfGratings);

pixelized_x_range = 1:numberOfGratings;
pixelized_y_range = 1:numberOfGratings;

pixeliZedScratchPad = zeros(numberOfGratings, numberOfGratings);

% pixelized the delineating cooridinates
for i = 1:numberOfPoints
        distances = abs(x_coordinates(i) - unpixelizedRange_x); % find closest distance
        min_distanceIndex = find(min(distances) == distances);
        min_distanceIndexFor_x(i) = min_distanceIndex(1,1); % if there is more than 1, choose 1

        distances = abs(y_coordinates(i) - unpixelizedRange_y); % find closest distance
        min_distanceIndex = find(min(distances) == distances);
        min_distanceIndexFor_y(i) = min_distanceIndex(1,1); % if there is more than 1, choose 1
end

% % pixelized the centroid coordinates
% distances = abs(centreOfShape(1,1) - unpixelizedRange_x);
% min_distanceIndex = find(min(distances) == distances);
% min_distanceIndexFor_Centroid_x0 = min_distanceIndex(1,1); % if there is more than 1, choose 1
%
% distances = abs(centreOfShape(1,2) - unpixelizedRange_y);
% min_distanceIndex = find(min(distances) == distances);
% min_distanceIndexFor_Centroid_y0 = min_distanceIndex(1,1); % if there is more than 1, choose 1

% % pixelized the centre coordinates
% distances = abs(0 - unpixelizedRange_x);
% min_distanceIndex = find(min(distances) == distances);
% min_distanceIndexFor_x0 = min_distanceIndex(1,1); % if there is more than 1, choose 1
%
% distances = abs(0 - unpixelizedRange_y);
% min_distanceIndex = find(min(distances) == distances);
% min_distanceIndexFor_y0 = min_distanceIndex(1,1); % if there is more than 1, choose 1

% connect the dots [in the case that the points are insufficient]
%{
increasedNumberOfPoints = 800;  % [TO DO]: automate this param
methodsOfInterpolation = {'spline', 'pchip', 'linear'};
methodOfInterpolationChosen = char(methodsOfInterpolation(1));
pt = interparc(increasedNumberOfPoints,min_distanceIndexFor_x, ...
    min_distanceIndexFor_y, methodOfInterpolationChosen);
min_distanceIndexFor_x = round(pt(:,1));
min_distanceIndexFor_y = round(pt(:,2));
%}
increasedNumberOfPoints = numberOfPoints; % in the case that the points are sufficient

% mark the delineating pixels
for i = 1:increasedNumberOfPoints
    pixeliZedScratchPad(min_distanceIndexFor_x(i), min_distanceIndexFor_y(i)) = 1;
end

% pixeliZedScratchPad(min_distanceIndexFor_Centroid_x0, min_distanceIndexFor_Centroid_y0) = 1; % mark the centroid pixel
% pixeliZedScratchPad(min_distanceIndexFor_x0, min_distanceIndexFor_y0) = 1; % mark the centre pixel relative to the plot

% invert it
pixeliZedScratchPad = ~pixeliZedScratchPad;

% plot right side up
pixeliZedScratchPad = pixeliZedScratchPad';
pixeliZedScratchPad = flipud(pixeliZedScratchPad);

end

author: Sable

repo: mcbench-benchmarks


text:
function [n,names] = plx_adchan_names(filename)
% plx_adchan_names(filename): gets the names of a/d channels in a .plx or .pl2 file
%
% [n,names] = plx_adchan_names(filename)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%
% OUTPUT:
%   names - array of a/d channel name strings
%   n - number of channels

if nargin ~= 1
    error 'expected 1 input argument';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    pl2 = PL2GetFileIndex(filename);
    n = numel(pl2.AnalogChannels);
    for i=1:numel(pl2.AnalogChannels)
        if i == 1
            names = pl2.AnalogChannels{i}.Name;
        else
            names = char(names, pl2.AnalogChannels{i}.Name);
        end
    end
    names = internalPL2TrimNames(names);
    return;
end

[n,names] = mexPlex(15, filename);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
clear all
close all

disp('Data from multiple subjects');

disp('Estimate mixed effects using Langevin Monte Carlo');
disp('Using LDS model with constrained connectivity');

lds.model='forward';

% Number of dynamical states e.g. brain areas
d=4;

% Observation noise
lds.sd=0.1;

% Prior over initial states
lds.R.pE=linspace(3,1.5,d)';
lds.R.pC=0.5^2*eye(d);

% Number of subjects
lds.Nsub=3;

% Number of observations per subject
lds.Nobs=5;

lds.init_par='random';
lds.flow_par='fixed';

% Generate group data
[lds.pinit,lds.pflow,lds.names,M,U,Y] = mci_lds_group_data (lds);

% Assign init/flow as random/fixed effects
assign.init_par='random';
assign.flow_par='fixed';
assign.out_par='known';

MCI.assign=assign;

MCI.fixed.pE=M{1}.pE;
MCI.fixed.pC=M{1}.pC;

% Initialisation
i0=mci_interp_init(Y,M{1});
a0=[];
if strcmp(assign.init_par,'random')
    MCI.pinit0=i0;
else
    MCI.pinit0=mean(i0,2);
end
if strcmp(assign.flow_par,'random')
    MCI.pflow0=spm_vec(M{1}.pE)*ones(1,lds.Nsub);
    if ~isempty(a0), MCI.pflow0(1:d,:)=a0; end
else
    MCI.pflow0=spm_vec(M{1}.pE);
    if ~isempty(a0), MCI.pflow0(1:d,:)=mean(a0,2); end
end
MCI.pout0=[];

MCI.M=M; MCI.U=U; MCI.Y=Y;
MCI.update_obs_noise=1;
MCI.verbose=1;
MCI.total_its=16;
MCI.rinit=0.25;

tic;
MCI = spm_mci_mfx_dynamic (MCI);
toc

rmse=mci_lds_plot_params (MCI,lds);

for n=1:lds.Nsub,
    mci_lds_plot_fit (MCI,lds,n,1);
end

disp('True dynamics:');
[f,Atrue] = mci_lds_fx (lds.pinit(:,1),U{1},lds.pflow,M{1});
disp(Atrue);

disp('Estimated dynamics:');
[f,Aest] = mci_lds_fx (MCI.pinit,U{1},MCI.pflow,M{1});
disp(Aest);

mci_plot_noiseSD (MCI.Ce,MCI.post_ind);

author: spm


text:
function msm_to_mm_test19 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST19 tests MSM_TO_MM_COORDINATE_REAL_SKEW_SYMMETRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST19\n' );
  fprintf ( 1, '  Convert an MSM to MM coordinate real skew-symmetric format.\n' );

  output_filename = 'msm_to_mm_test19.mm';

  a = r8mat_indicator ( 4, 4 );
  a = a - a';
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'coordinate', 'real', 'skew-symmetric' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [xi] = logSE3(chi)
C = chi(1:3,1:3);
r = chi(1:3,4);
[vecs,eigs] = eig(C);
[~,idx] = min(abs(diag(eigs)-1));
a = vecs(:,idx);

phi = acos(1/2*(trace(C)-1));
if phi == 0
    a = zeros(3,1);
    iJ = eye(3);
else
    C1 = cos(phi)*eye(3) + (1-cos(phi))*a*transpose(a) + sin(phi)*[[0,-a(3),a(2)];[a(3),0,-a(1)];[-a(2),a(1),0]];
    C2 = cos(-phi)*eye(3) + (1-cos(-phi))*a*transpose(a) + sin(-phi)*[[0,-a(3),a(2)];[a(3),0,-a(1)];[-a(2),a(1),0]];
    if norm(C2-C)<norm(C1-C)
        phi = -phi;
    end
    iJ = phi/2*cot(phi/2)*eye(3) + (1-phi/2*cot(phi/2))*a*a' - phi/2*[[0,-a(3),a(2)];[a(3),0,-a(1)];[-a(2),a(1),0]];
end
rho = iJ*r;
xi = real([phi*a;rho]);
end

author: mbrossar


text:
function f = floor(f)
%FLOOR   Pointwise floor function of a CHEBTECH.
%   G = FLOOR(F) returns the CHEBTECH G such that G(X) = FLOOR(F(x)) for each x
%   in [-1, 1].
%
%   If F is complex, then the G = FLOOR(REAL(F)) + 1i*FLOOR(IMAG(F)).
%
%   Note that FLOOR() assumes the output G(X) is a constant. If it is not, then
%   garbage is returned with no warning.
%
% See also CEIL, ROUND, FIX.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Evaluate at the two end points, and an arbitrary interior point:
arbitraryPoint = 0.1273881594;
fx = feval(f, [-1 ; arbitraryPoint ; 1]);
% Take the mean:
meanfx = mean(fx, 1);
% Compute the floor:
f.coeffs = floor(meanfx);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%  Figure 10.05      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_05.m is a script to generate Fig. 10.5
% root locus, for the PD design for the satellite
clf;
% parameters of two-mass spring model
m=[1, 0.1]; k0=[0, 0.091] ; d0=[0, 0.0036]; k1=[0, 0.4];

% call two mass-spring model function
[f,g,h,j] = twomass(m,k0,d0);
nc1=0.25*[2, 1];
dc1=[1/40, 1];

% convert controller to state-space
[ac,bc,cc,dc]=tf2ss(nc1, dc1);

% series of controller and plant
[aol,bol,col,dol]= series(ac, bc,cc,dc,f,g,h,j);
[acl]=aol-bol*col;

hold off ; clf

rlocus(aol,bol,col,dol)
grid;
v =[-2.0000,    2.0000,   -1.5000,    1.5000];
axis(v);
title('Fig. 10.5 Root locus for the PD design of the satellite')
%


text:
function Z = uplus(Y)
%UPLUS (overloaded)

if isa(Y,'blkvar')
    Y = sdpvar(Y);
end

switch Y.typeflag
case {0,5}
	Z = Y;
case {1,3,4}
	error('Relational objects canot be manipulated')
otherwise
	error('Please report internal bug in sdpvar/uplus')
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [] = main(training_file, test_file, layers, units_per_layer, rounds)
    obj = neural_network(training_file, test_file, layers, units_per_layer, rounds);
    obj = obj.initialise(obj);
    for i = 1:obj.rounds
        obj = obj.feed_forward(obj, i-1);
    end
    obj = obj.testing(obj);
end

repo: Machine-Learning-Models


text:
function covmats_unflipped = getAllCovMats(data,T,options)

if ~isfield(options,'maxlag'), options.maxlag = 10; end
if ~isfield(options,'partial'), options.partial = 0; end
if ~isfield(options,'standardise'), options.standardise = 1; end

N = length(T);
if iscell(data)
    for j = 1:N
        if ischar(data{j})
            fsub = data{j};
            loadfile_sub;
        else
            X = data{j};
        end
        Tj = T{j}; Nj = length(Tj);
        if options.standardise
            for jj = 1:Nj
                ind = (1:Tj(jj)) + sum(Tj(1:jj-1));
                X(ind,:) = bsxfun(@minus,X(ind,:),mean(X(ind,:)));
                sd = std(X(ind,:));
                if any(sd==0)
                    error('At least one channel in at least one trial has variance equal to 0')
                end
                X(ind,:) = X(ind,:) ./ repmat(sd,Tj(jj),1);
            end
        end
        covmats_unflipped_j = getCovMats(X,sum(Tj),options.maxlag,options.partial);
        if j==1
            covmats_unflipped = covmats_unflipped_j;
        else % do some kind of weighting here according to the number of samples?
            covmats_unflipped = cat(4,covmats_unflipped,covmats_unflipped_j);
        end
    end

else
    if isstruct(data), data = data.X; end
    if options.standardise
        for j = 1:N
            ind = (1:T(j)) + sum(T(1:j-1));
            data(ind,:) = bsxfun(@minus,data(ind,:),mean(data(ind,:)));
            sd = std(data(ind,:));
            if any(sd==0)
                error('At least one channel in at least one trial has variance equal to 0')
            end
            data(ind,:) = data(ind,:) ./ repmat(sd,T(j),1);
        end
    end
    covmats_unflipped = getCovMats(data,T,options.maxlag,options.partial);
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function i4vec_aminz_test ( )

%*****************************************************************************80
%
%% I4VEC_AMINZ_TEST tests I4VEC_AMINZ;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_AMINZ_TEST\n' );
  fprintf ( 1, '  For an integer vector:\n' );
  fprintf ( 1, '  I4VEC_AMINZ:  minimum nonzero absolute entry;\n' );

  seed = 123456789;

  b = -n;
  c = n;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Input vector:' );

  fprintf ( 1, '\n' );

  aval = i4vec_aminz ( n, a );

  fprintf ( 1, '  Minimum abs nonzero:       %f\n', aval );

  return
end


text:
function pix_hdr = GE_readHeaderPixel(fid, byte_align)
%
% pix_hdr = read_pixel_header(fid, byte_align)
%
% Loads the pixel header from the file with filed id fid
% and returns it as a structure.
% if byte_align = 1 then 32-bit alignment (SGI, LX2 format)
% if byte_align = 0 then 16-bit alignment (Sun, 5.X format)
%
% Souheil J. Inati
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%

% define the structure and read
pix_hdr = struct('img_magic', fread(fid,1,'int32'));
pix_hdr = setfield(pix_hdr, 'img_hdr_length', fread(fid,1,'int32'));% length of total header in bytes and
                                                                    % a byte displacement to the 'pixel data area'
pix_hdr = setfield(pix_hdr, 'img_width', fread(fid,1,'int32'));  % width (pixels) of image
pix_hdr = setfield(pix_hdr, 'img_height', fread(fid,1,'int32')); % height (pixels) of image
pix_hdr = setfield(pix_hdr, 'img_depth', fread(fid,1,'int32'));  % depth (1, 8, 16, or 24 bits) of pixel
pix_hdr = setfield(pix_hdr, 'img_compress', fread(fid,1,'int32'));   % type of compression; see IC_* below
pix_hdr = setfield(pix_hdr, 'img_dwindow', fread(fid,1,'int32'));    % default window setting
pix_hdr = setfield(pix_hdr, 'img_dlevel', fread(fid,1,'int32')); % default level setting
pix_hdr = setfield(pix_hdr, 'img_bgshade', fread(fid,1,'int32'));    % background shade to use for non-image
pix_hdr = setfield(pix_hdr, 'img_ovrflow', fread(fid,1,'int32'));    % overflow value
pix_hdr = setfield(pix_hdr, 'img_undflow', fread(fid,1,'int32'));    % underflow value
pix_hdr = setfield(pix_hdr, 'img_top_offset', fread(fid,1,'int32')); % number of blank lines at image top
pix_hdr = setfield(pix_hdr, 'img_bot_offset', fread(fid,1,'int32')); % number of blank lines at image bottom
pix_hdr = setfield(pix_hdr, 'img_version', fread(fid,1,'int16'));    % version of the header structure
if byte_align, fseek(fid,2,0); end % 32-bit alignment
pix_hdr = setfield(pix_hdr, 'img_checksum', fread(fid,1,'uint16')); % 16 bit end_around_carry sum of pixels
pix_hdr = setfield(pix_hdr, 'img_p_id', fread(fid,1,'int32'));   % a byte disp to unique image identifier
pix_hdr = setfield(pix_hdr, 'img_l_id', fread(fid,1,'int32'));   % byte length of unique image identifier
pix_hdr = setfield(pix_hdr, 'img_p_unpack', fread(fid,1,'int32'));   % a byte disp to 'unpack control'
pix_hdr = setfield(pix_hdr, 'img_l_unpack', fread(fid,1,'int32'));   % byte length of 'unpack control'
pix_hdr = setfield(pix_hdr, 'img_p_compress', fread(fid,1,'int32')); % a byte disp to 'compression control'
pix_hdr = setfield(pix_hdr, 'img_l_compress', fread(fid,1,'int32')); % byte length of 'compression control'
pix_hdr = setfield(pix_hdr, 'img_p_histo', fread(fid,1,'int32'));    % a byte disp to 'histogram control'
pix_hdr = setfield(pix_hdr, 'img_l_histo', fread(fid,1,'int32'));    % byte length of 'histogram control'
pix_hdr = setfield(pix_hdr, 'img_p_text', fread(fid,1,'int32')); % a byte disp to 'text plane data'
pix_hdr = setfield(pix_hdr, 'img_l_text', fread(fid,1,'int32')); % byte length of 'text plane data'
pix_hdr = setfield(pix_hdr, 'img_p_graphics', fread(fid,1,'int32')); % a byte disp to 'graphics plane data'
pix_hdr = setfield(pix_hdr, 'img_l_graphics', fread(fid,1,'int32')); % byte length of 'graphics plane data'
pix_hdr = setfield(pix_hdr, 'img_p_dbHdr', fread(fid,1,'int32'));    % a byte disp to 'data base header data'
pix_hdr = setfield(pix_hdr, 'img_l_dbHdr', fread(fid,1,'int32'));    % byte length of 'data base header data'
pix_hdr = setfield(pix_hdr, 'img_levelOffset', fread(fid,1,'int32'));% value to add to stored Pixel Data values
                                                                     % to get the correct presentation value
pix_hdr = setfield(pix_hdr, 'img_p_user', fread(fid,1,'int32')); % byte displacement to user defined data
pix_hdr = setfield(pix_hdr, 'img_l_user', fread(fid,1,'int32')); % byte length of user defined data

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [Channel] = elec_3dspace2brainstorm(filename,N,bsfile)

% elec_3dspace2brainstorm - Convert NeuroScan 3Dspace ascii to brainstorm file
%
% The ascii file format is that of NeuroScan 3Dspace export files.  Each row
% of the file comprises an electrode label, an electrode type code, and
% the x,y,z coordinates (cm).  Each field is separated by spaces.  See
% ELEC_LOAD for more information.
%
% Useage: Channel = elec_3dspace2brainstorm(filename,[N],[brainstormfile])
%
% where:    file = 'path\filename' with format described in ELEC_LOAD
%
%           N is how many electrodes to load (rows of 3Dspace file, 129 default).
%
% Result:   Channel is an array of structures.  The fields are:
%
%           Loc     - a 3x2 matrix of electrode coordinates (x,y,z in rows).
%                     BrainStorm (x,y,z meters) = 3Dspace (x,y,z cm) / 100.
%           Orient  - a corresponding matrix of sensor orientations (MEG);
%                     all zero for EEG.
%           Weight  - a vector of relative or absolute weights (eg, amplification);
%                     all ones for this routine.
%           Type    - a character string, 'EEG' in this case.
%           Name    - a charater string indicating the electrode name.
%           Comment - a charater string indicating the reference electrode.  Empty
%                     for active electrodes and 'EEG REF' for reference electrode.
%
% See brainstorm website at http://neuroimage.usc.edu/, including a
% download pdf file describing the brainstorm database formats.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Licence:  GNU GPL, no express or implied warranties
% History:  20/05/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('N', 'var'), N = 129;
elseif isempty(N),     N = 129;
end

file = char(filename);

[elec,type,X,Y,Z] = elec_load(file,[],[],[],[],N);

tic;

fprintf('\nELEC_3DSPACE2BRAINSTORM...\n');
fprintf('...Converting to brainstorm structure.\n');

elecindex = find(type == 69);
Ee = elec(elecindex);
Xe = X(elecindex) ./ 100;   % 3Dspace is cm, BrainStorm is m
Ye = Y(elecindex) ./ 100;
Ze = Z(elecindex) ./ 100;

reftype = ones(size(type)) * 120;
refindex = find(type == reftype);
Eref = elec(refindex);
Xref = X(refindex) ./ 100;   % 3Dspace is cm, BrainStorm is m
Yref = Y(refindex) ./ 100;
Zref = Z(refindex) ./ 100;

for i=1:length(elecindex),
    Channel(i).Loc = [[Xe(i) Ye(i) Ze(i)]',[Xref Yref Zref]'];
    Channel(i).Orient = [];     % used for MEG rather than EEG
    Channel(i).Weight = 1;      % Like Amplification
    Channel(i).Type = 'EEG';
    Channel(i).Name = char(Ee(i));
    Channel(i).Comment = '';
end
Channel(i+1).Loc = [[Xref Yref Zref]',[Xref Yref Zref]'];
Channel(i+1).Orient = [];
Channel(i+1).Weight = 1;
Channel(i+1).Type = 'EEG';
Channel(i+1).Name = char(Eref);
Channel(i+1).Comment = 'EEG REF';

if ~exist('bsfile', 'var'),
    bsfile = 'channel';
elseif isempty(bsfile),
    bsfile = 'channel';
end

if findstr('.mat',bsfile),
    bsfile = strrep(bsfile,'.mat','');
end

fprintf('...saving BrainStorm channel data to:\n\t%s.mat\n',bsfile);
save(bsfile, 'Channel');

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
function [g, gdata, gprior] = gp_g(w, gp, x, y, varargin)
%GP_G  Evaluate the gradient of energy (GP_E) for Gaussian Process
%
%  Description
%    G = GP_G(W, GP, X, Y, OPTIONS) takes a full GP parameter
%    vector W, GP structure GP, a matrix X of input vectors and a
%    matrix Y of target vectors, and evaluates the gradient G of
%    the energy function (gp_e). Each row of X corresponds to one
%    input vector and each row of Y corresponds to one target
%    vector.
%
%    [G, GDATA, GPRIOR] = GP_G(W, GP, X, Y, OPTIONS) also returns
%    separately the data and prior contributions to the gradient.
%
%    OPTIONS is optional parameter-value pair
%      z - optional observed quantity in triplet (x_i,y_i,z_i)
%          Some likelihoods may use this. For example, in case of
%          Poisson likelihood we have z_i=E_i, that is, expected
%          value for ith case.
%
%  See also
%    GP_E, GP_PAK, GP_UNPAK, GPCF_*
%
% Copyright (c) 2007-2011 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2010 Heikki Peura
% Copyright (c) 2014 Arno Solin and Jukka Koskenranta

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

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

ip=inputParser;
ip.FunctionName = 'GP_G';
ip.addRequired('w', @(x) isvector(x) && isreal(x));
ip.addRequired('gp',@isstruct);
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.parse(w, gp, x, y, varargin{:});
z=ip.Results.z;
if ~all(isfinite(w(:)));
    % instead of stopping to error, return NaN
    g=NaN;
    gdata = NaN;
    gprior = NaN;
    return;
end

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

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

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

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

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

        if ~isfield(gp.lik, 'nondiagW') || ismember(gp.lik.type, {'LGP' 'LGPC'})
            % Evaluate covariance
            if isfield(gp, 'lik_mono')
                [K,C] = gp_dtrcov(gp, x, gp.xv);
                %         if isequal(gp.lik.type, 'Gaussian')
                %           % Condition the derivatives on the observations
                %           cc=C(size(x,1)+1:end,size(x,1)+1:end);
                %           cy=C(size(x,1)+1:end,1:size(x,1));
                %           cyy=C(1:size(x,1),1:size(x,1));
                %           C=cc - cy*(cyy\cy');
                %           %C=0.5.*(C+C')+1e-10.*eye(size(C));
                %           meany=cy*(cyy\y(1:n));
                %           y=(y(n+1:end)-meany);
                %           n=length(y);
                %         end
            else
                [~, C] = gp_trcov(gp, x);
            end

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

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

        % =================================================================
        % Gradient with respect to covariance function parameters
        i1=0;
        if ~isempty(strfind(gp.infer_params, 'covariance'))

            if isfield(gp,'deriv') && gp.deriv  % derivative observations in use
                % initialize few variables related to derivative observations
                [n,m] = size(x);
                ind_Ddim = x(:,gp.deriv);
                ind_Ddim_derivs = ind_Ddim(ind_Ddim>0);
                uDdim = unique(ind_Ddim_derivs);
                x1 = x(:,setdiff(1:m,gp.deriv));   % Take only the non-index columns
            end

            for i=1:ncf

                gpcf = gp.cf{i};

                if isfield(gp.lik, 'nondiagW') && ~ismember(gp.lik.type, {'LGP' 'LGPC'})
                    % check in which components the covariance function is present
                    % for likelihoods with non-diagonal Hessian
                    do = false(nout,1);
                    if multicf
                        for z1=1:nout
                            if any(gp.comp_cf{z1}==i)
                                do(z1) = true;
                            end
                        end
                    else
                        do = true(nout,1);
                    end
                end

                if ~(isfield(gp,'deriv') && gp.deriv)
                    % No derivative observations
                    if ~savememory
                        if isfield(gp.lik, 'nondiagW') && isfield(gp,'comp_cf') && ~isempty(intersect(gp.comp_cf{1},i)) && isfield(gp.lik, 'xtime')
                            DKffc = gpcf.fh.cfg(gpcf, xtime);
                        else
                            DKffc = gpcf.fh.cfg(gpcf, x);
                        end
                        np=length(DKffc);
                    else
                        % If savememory option is used, just get the number of
                        % hyperparameters and calculate gradients later
                        np=gpcf.fh.cfg(gpcf,[],[],[],0);
                    end
                    gprior_cf = -gpcf.fh.lpg(gpcf);
                else
                    [n,m]=size(x);
                    gprior_cf = -gpcf.fh.lpg(gpcf);
                    if (~isfield(gpcf, 'selectedVariables') || any(ismember(gpcf.selectedVariables,uDdim)))
                        % !!! Note. The check whether to calculate derivative
                        % matrices for a covariance function could/should be made
                        % nicer
                            DKffa = gpcf.fh.cfg(gpcf, x1(ind_Ddim==0,:));
                            np=length(DKffa);
                            for inp=1:np
                                % the block of covariance matrix
                                DKffc{inp}(ind_Ddim==0,ind_Ddim==0) = DKffa{inp};
                            end
                            for u1 = 1:length(uDdim)
                                % the blocks on the left side, below Kff
                                Kdf = gpcf.fh.cfdg(gpcf, x1(ind_Ddim==uDdim(u1),:), x1(ind_Ddim==0,:), uDdim(u1));
                                D = gpcf.fh.cfdg2(gpcf, x1(ind_Ddim==uDdim(u1),:), x1(ind_Ddim==uDdim(u1),:), uDdim(u1), uDdim(u1));
                                for inp=1:np
                                    DKffc{inp}(ind_Ddim==uDdim(u1),ind_Ddim==0) = Kdf{inp};
                                    DKffc{inp}(ind_Ddim==0,ind_Ddim==uDdim(u1)) = Kdf{inp}';
                                    DKffc{inp}(ind_Ddim==uDdim(u1),ind_Ddim==uDdim(u1)) = D{inp};
                                end

                                uDdim2 = uDdim(u1+1:end);
                                for u2=1:length(uDdim2)
                                    Kdf2 = gpcf.fh.cfdg2(gpcf, x1(ind_Ddim==uDdim(u1),:) ,x1(ind_Ddim==uDdim2(u2),:), uDdim(u1), uDdim2(u2));
                                    for inp=1:np
                                        DKffc{inp}(ind_Ddim==uDdim(u1),ind_Ddim==uDdim2(u2)) = Kdf2{inp};
                                        DKffc{inp}(ind_Ddim==uDdim2(u2),ind_Ddim==uDdim(u1)) = Kdf2{inp}';
                                    end
                                end

                            end

                            %DKffc{inp} = Ktemp;

                    else
                        % A covariance function that does not use any
                        % of the input dimensions with respect to which
                        % derivatives are observed
                        %warning('not tested yet')
                        DKffa = gpcf.fh.cfg(gpcf, x(ind_Ddim==0,:));
                        np=length(DKffa);
                        for inp=1:np
                            Ktemp = sparse(n,n);
                            Ktemp(ind_Ddim==0,ind_Ddim==0) = DKffa{inp};
                            DKffc{inp} = Ktemp;
                        end
                    end

                end

                % Are there specified mean functions
                if  ~isfield(gp,'meanf')
                    % Evaluate the gradient with respect to covariance function
                    % parameters
                    for i2 = 1:np
                        if savememory
                            if isfield(gp.lik, 'nondiagW') && isfield(gp,'comp_cf') && ~isempty(intersect(gp.comp_cf{1},i)) && isfield(gp.lik, 'xtime')
                                DKff=gpcf.fh.cfg(gpcf,xtime,[],[],i2);
                            else
                                DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
                            end
                        else
                            DKff=DKffc{i2};
                        end
                        i1 = i1+1;
                        if ~isfield(gp.lik, 'nondiagW') || ismember(gp.lik.type, {'LGP' 'LGPC'})
                            Bdl = b'*(DKff*b);
                            Cdl = sum(sum(invC.*DKff)); % help arguments
                        else
                            % Non-diagonalizable likelihoods
                            Bdl=0; Cdl=0;
                            if isfield(gp.lik,'xtime');
                                if do(1)
                                    Bdl = Bdl + b(1:ntime)'*(DKff*b(1:ntime));
                                    Cdl = Cdl + sum(sum(invC(1:ntime,1:ntime).*DKff)); % help arguments
                                end
                                if do(2)
                                    Bdl = Bdl + b(ntime+1:end)'*(DKff*b(ntime+1:end));
                                    Cdl = Cdl + sum(sum(invC(ntime+1:end,ntime+1:end).*DKff)); % help arguments
                                end
                            else
                                for z1=1:nout
                                    if do(z1)
                                        Bdl = Bdl + b(1+nl(z1):nl(z1+1))'*(DKff*b(1+nl(z1):nl(z1+1)));
                                        Cdl = Cdl + sum(sum(invC(:,:,z1).*DKff)); % help arguments
                                    end
                                end
                            end
                        end
                        gdata(i1)=0.5.*(Cdl - Bdl);
                    end
                else
                    for i2 = 1:np
                        if savememory
                            DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
                        else
                            DKff=DKffc{i2};
                        end
                        i1=i1+1;
                        dA = -1*HinvC*DKff*HinvC';                  % d A / d th
                        trA = sum(invAt(:).*dA(:));                 % d log(|A|) / dth
                        dMNM = invNM'*(DKff*invNM);           % d M'*N*M / d th
                        trK = sum(sum(invC.*DKff));       % d log(Ky�?�) / d th
                        gdata(i1)=0.5*(-1*dMNM + trK + trA);
                    end
                end

                gprior = [gprior gprior_cf];
            end
        end

        % =================================================================
        % Gradient with respect to Gaussian likelihood function parameters
        if ~isempty(strfind(gp.infer_params, 'likelihood')) && isfield(gp.lik.fh,'trcov')
            % Evaluate the gradient from Gaussian likelihood
            gprior_lik = -gp.lik.fh.lpg(gp.lik);
            if ~isempty(gprior_lik)
                DCff = gp.lik.fh.cfg(gp.lik, x);
                if isfield(gp, 'lik_mono')
                    DCff{1}=diag(DCff{1}*[ones(size(x,1),1); zeros(size(gp.xv,1)*length(nvd),1)]);
                end
                for i2 = 1:length(DCff)
                    i1 = i1+1;
                    if ~isfield(gp,'meanf')
                        if size(DCff{i2}) > 1
                            yKy = b'*(DCff{i2}*b);
                            trK = sum(sum(invC.*DCff{i2})); % help arguments
                            gdata_zeromean(i1)=0.5.*(trK - yKy);
                        else
                            yKy=DCff{i2}.*(b'*b);
                            trK = DCff{i2}.*(trace(invC));
                            gdata_zeromean(i1)=0.5.*(trK - yKy);
                        end
                        gdata(i1)=gdata_zeromean(i1);
                    else
                        if size(DCff{i2}) > 1
                            trK = sum(sum(invC.*DCff{i2})); % help arguments
                        else
                            trK = DCff{i2}.*(trace(invC));
                        end
                        dA = -1*HinvC*DCff{i2}*HinvC';            % d A / d th
                        trA = sum(invAt(:).*dA(:));               % d log(|A|) / dth
                        dMNM = invNM'*(DCff{i2}*invNM);           % d M'*N*M / d th
                        gdata(i1)=0.5*(-1*dMNM + trA + trK);
                    end
                end
            end
            gprior = [gprior gprior_lik];
        end

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

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

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

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

        u = gp.X_u;
        DKuu_u = 0;
        DKuf_u = 0;

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

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

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

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

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

                gprior = [gprior gprior_cf];
            end
        end

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

        % =================================================================
        % Gradient with respect to inducing inputs
        if ~isempty(strfind(gp.infer_params, 'inducing'))
            if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
                m = size(gp.X_u,2);
                st=0;
                if ~isempty(gdata)
                    st = length(gdata);
                end

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

                % Loop over the covariance functions
                for i=1:ncf
                    i1 = st;
                    gpcf = gp.cf{i};

                    if savememory
                        % If savememory option is used, just get the number of
                        % covariates in X and calculate gradients later
                        np=gpcf.fh.ginput(gpcf,u,[],0);
                    else
                        np=1;
                        DKuu = gpcf.fh.ginput(gpcf, u);
                        DKuf = gpcf.fh.ginput(gpcf, u, x);
                    end
                    for i3=1:np
                        if savememory
                            DKuu=gpcf.fh.ginput(gpcf,u,[],i3);
                            DKuf=gpcf.fh.ginput(gpcf,u,x,i3);
                        end
                        for i2 = 1:length(DKuu)
                            if savememory
                                i1 = st + (i2-1)*np+i3;
                                %i1 = st + 2*i2 - (i3==1);
                            else
                                i1 = i1+1;
                            end
                            KfuiKuuKuu = iKuuKuf'*DKuu{i2};
                            gdata(i1) = gdata(i1) - 0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuKuu))*(iKuuKuf*b') + ...
                                2.*sum(sum(L'.*(L'*DKuf{i2}'*iKuuKuf))) - sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));
                            gdata(i1) = gdata(i1) + 0.5.*(2.*b.*sum(DKuf{i2}'.*iKuuKuf',2)'*b'- b.*sum(KfuiKuuKuu.*iKuuKuf',2)'*b');
                            gdata(i1) = gdata(i1) + 0.5.*(2.*sum(LL.*sum(DKuf{i2}'.*iKuuKuf',2)) - ...
                                sum(LL.*sum(KfuiKuuKuu.*iKuuKuf',2)));
                        end
                    end
                end
            end
        end

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

        u = gp.X_u;
        ind = gp.tr_index;
        DKuu_u = 0;
        DKuf_u = 0;

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

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

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

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

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

                for i2 = 1:np
                    i1 = i1+1;
                    if savememory
                        DKuu=gpcf.fh.cfg(gpcf,u,[],[],i2);
                        DKuf=gpcf.fh.cfg(gpcf,u,x,[],i2);
                    else
                        DKuu=DKuuc{i2};
                        DKuf=DKufc{i2};
                    end
                    KfuiKuuKuu = iKuuKuf'*DKuu;
                    %            H = (2*K_uf'- KfuiKuuKuu)*iKuuKuf;
                    % Here we evaluate  gdata = -0.5.* (b*H*b' + trace(L*L'H)
                    gdata(i1) = -0.5.*((2*b*DKuf'-(b*KfuiKuuKuu))*(iKuuKuf*b') + 2.*sum(sum(L'.*(L'*DKuf'*iKuuKuf))) - ...
                        sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));
                    for kk=1:length(ind)
                        if savememory
                            DKff=gpcf.fh.cfg(gpcf,x(ind{kk},:),[],[],i2);
                        else
                            DKff=DKffc{kk}{i2};
                        end
                        gdata(i1) = gdata(i1) ...
                            + 0.5.*(-b(ind{kk})*DKff*b(ind{kk})' ...
                            + 2.*b(ind{kk})*DKuf(:,ind{kk})'*iKuuKuf(:,ind{kk})*b(ind{kk})'- ...
                            b(ind{kk})*KfuiKuuKuu(ind{kk},:)*iKuuKuf(:,ind{kk})*b(ind{kk})' ...
                            +trace(La{kk}\DKff)...
                            - trace(L(ind{kk},:)*(L(ind{kk},:)'*DKff)) ...
                            + 2.*sum(sum(L(ind{kk},:)'.*(L(ind{kk},:)'*DKuf(:,ind{kk})'*iKuuKuf(:,ind{kk})))) - ...
                            sum(sum(L(ind{kk},:)'.*((L(ind{kk},:)'*KfuiKuuKuu(ind{kk},:))*iKuuKuf(:,ind{kk})))));
                    end
                end
                gprior = [gprior gprior_cf];
            end
        end

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

        % =================================================================
        % Gradient with respect to inducing inputs
        if ~isempty(strfind(gp.infer_params, 'inducing'))
            if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
                m = size(gp.X_u,2);

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

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

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

                    for i3 = 1:np
                        if savememory
                            DKuu=gpcf.fh.ginput(gpcf, u, [], i3);
                            DKuf=gpcf.fh.ginput(gpcf, u, x, i3);
                        end
                        for i2 = 1:length(DKuu)
                            if savememory
                                i1 = st + (i2-1)*np+i3;
                                %i1 = st + 2*i2 - (i3==1);
                            else
                                i1 = i1+1;
                            end
                            KfuiKuuDKuu_u = iKuuKuf'*DKuu{i2};
                            gdata(i1) = gdata(i1) -0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuDKuu_u))*(iKuuKuf*b') + 2.*sum(sum(L'.*((L'*DKuf{i2}')*iKuuKuf))) - ...
                                sum(sum(L'.*((L'*KfuiKuuDKuu_u)*iKuuKuf))));

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

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

        u = gp.X_u;
        DKuu_u = 0;
        DKuf_u = 0;

        cf_orig = gp.cf;

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

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

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

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

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

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

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

                gpcf = gp.cf{i};

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

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

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

                        %gdata(i1) = gdata(i1) + 0.5.*sum(sum(La\((2.*K_uf') - KfuiKuuKuu).*iKuuKuf',2));
                        gdata(i1) = gdata(i1) + 0.5.*sum(sum(ldlsolve(LD,temp3).*iKuuKuf',2));
                        gdata(i1) = gdata(i1) - 0.5.*( idiagLa'*(sum(temp3.*iKuuKuf',2)) ); % corrected
                    end

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

                    for i2 = 1:np
                        i1 = i1+1;
                        if savememory
                            DKff=gpcf.fh.cfg(gpcf,x,[],[],i2);
                        else
                            DKff=DKffc{i2};
                        end
                        gdata(i1) = 0.5*(sum(sum(siLa.*DKff',2)) - sum(sum(L.*(L'*DKff')')) - b*DKff*b');
                    end
                end
                gprior = [gprior gprior_cf];
            end
        end

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

        % =================================================================
        % Gradient with respect to inducing inputs
        if ~isempty(strfind(gp.infer_params, 'inducing'))
            if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
                m = size(gp.X_u,2);
                st=0;
                if ~isempty(gdata)
                    st = length(gdata);
                end

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

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

                        for i3 = 1:np
                            if savememory
                                DKuu = gpcf.fh.ginput(gpcf,u,[],i3);
                                DKuf = gpcf.fh.ginput(gpcf,u,x,i3);
                            end
                            for i2 = 1:length(DKuu)
                                if savememory
                                    i1 = st + (i2-1)*np+i3;
                                    %i1 = st + 2*i2 - (i3==1);
                                else
                                    i1 = i1+1;
                                end
                                KfuiKuuKuu = iKuuKuf'*DKuu{i2};

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

                                gdata(i1) = gdata(i1) + 0.5.*sum(sum(ldlsolve(LD,(2.*DKuf{i2}') - KfuiKuuKuu).*iKuuKuf',2));
                                gdata(i1) = gdata(i1) - 0.5.*( idiagLa'*(sum((2.*DKuf{i2}' - KfuiKuuKuu).*iKuuKuf',2)) ); % corrected
                            end
                        end
                    end
                end
            end
        end

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

        u = gp.X_u;
        DKuu_u = 0;
        DKuf_u = 0;

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

        LL1=iLav-LL;

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

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

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

                for i2 = 1:np
                    i1 = i1+1;
                    if savememory
                        % If savememory option is used, just get the number of
                        % hyperparameters and calculate gradients later
                        DKff=gpcf.fh.cfg(gpcf,x,[],1,i2);
                        DKuu=gpcf.fh.cfg(gpcf,u,[],[],i2);
                        DKuf=gpcf.fh.cfg(gpcf,u,x,[],i2);
                    else
                        DKff=DKffc{i2};
                        DKuu=DKuuc{i2};
                        DKuf=DKufc{i2};
                    end

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

                    if strcmp(gp.type, 'VAR')
                        gdata(i1) = gdata(i1) + 0.5.*(sum(iLav.*DKff)-2.*sum(iLav'*sum(DKuf'.*iKuuKuf',2)) + ...
                            sum(iLav'*sum(KfuiKuuKuu.*iKuuKuf',2))); % trace-term derivative
                    end
                end
                gprior = [gprior gprior_cf];
            end
        end

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

        % =================================================================
        % Gradient with respect to inducing inputs
        if ~isempty(strfind(gp.infer_params, 'inducing'))
            if isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
                m = size(gp.X_u,1);
                st=0;
                if ~isempty(gdata)
                    st = length(gdata);
                end

                gdata(st+1:st+length(gp.X_u(:))) = 0;
                i1 = st+1;
                gprior_ind=[];
                if iscell(gp.p.X_u) % Own prior for each inducing input
                    for i = 1:size(gp.X_u,1)
                        pr = gp.p.X_u{i};
                        gprior_ind =[gprior_ind -pr.fh.lpg(gp.X_u(i,:), pr)];
                    end
                else % One prior for all inducing inputs
                    gprior_ind = -gp.p.X_u.fh.lpg(gp.X_u(:)', gp.p.X_u);
                end
                gprior = [gprior gprior_ind];
                i1 = i1 + m;
                % Loop over the covariance functions
                for i=1:ncf
                    i1 = st;
                    gpcf = gp.cf{i};
                    if savememory
                        % If savememory option is used, just get the number of
                        % covariates in X and calculate gradients later
                        np=gpcf.fh.ginput(gpcf,u,[],0);
                    else
                        np=1;
                        DKuu = gpcf.fh.ginput(gpcf, u);
                        DKuf = gpcf.fh.ginput(gpcf, u, x);
                    end

                    for i3 = 1:np
                        if savememory
                            DKuu=gpcf.fh.ginput(gpcf,u,[],i3);
                            DKuf=gpcf.fh.ginput(gpcf,u,x,i3);
                        end
                        for i2 = 1:length(DKuu)
                            if savememory
                                i1 = st + (i2-1)*np+i3;
                                %i1 = st + 2*i2 - (i3==1);
                            else
                                i1 = i1+1;
                            end
                            KfuiKuuKuu = iKuuKuf'*DKuu{i2};
                            gdata(i1) = gdata(i1) - 0.5.*((2*b*DKuf{i2}'-(b*KfuiKuuKuu))*(iKuuKuf*b'));
                            gdata(i1) = gdata(i1) + 0.5.*(2.*(sum(iLav'*sum(DKuf{i2}'.*iKuuKuf',2))-sum(sum(L'.*(L'*DKuf{i2}'*iKuuKuf))))...
                                - sum(iLav'*sum(KfuiKuuKuu.*iKuuKuf',2))+ sum(sum(L'.*((L'*KfuiKuuKuu)*iKuuKuf))));

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

    case {'KALMAN'}
        % ============================================================
        % Kalman filtering and smoothing
        % ============================================================
        %
        % The implementation below is primarily based on the methods
        % presented in the following publication. If you find this
        % useful as a part of your own research, please cite the papers.
        %
        %  [1] Simo Sarkka, Arno Solin, Jouni Hartikainen (2013).
        %      Spatiotemporal learning via infinite-dimensional Bayesian
        %      filtering and smoothing. IEEE Signal Processing Magazine,
        %      30(4):51-61.
        %
        %  [2] Simo Sarkka (2013). Bayesian filtering and smoothing.
        %      Cambridge University Press.
        %
        %  [3] Simo Sarkka (2006). Recursive Bayesian inference on stochastic
        %      differential equations. Doctoral dissertation, Helsinki
        %      University of Technology, Finland.
        %

        % Ensure that this is a purely temporal problem
        if size(x,2) > 1,
            error('The ''KALMAN'' option only supports one-dimensional data.')
        end

        % Extract the noise magnitude from the GP likelihood model
        R = gp.lik.sigma2;

        % Initialize model matrices
        F   = []; L     = []; Qc = [];
        H   = []; Pinf  = []; dF = [];
        dQc = []; dPinf = []; isstable = true;

        % For each covariance function
        for j=1:length(gp.cf)

            % Form correpsonding state space model for this covariance function
            try
                [jF,jL,jQc,jH,jPinf,jdF,jdQc,jdPinf,p] = gp.cf{j}.fh.cf2ss(gp.cf{j},x);
            catch
                gdata = nan*w; gprior = nan*w; g = nan*w;
                return
            end

            % Stack model
            F  = blkdiag(F,jF);
            L  = blkdiag(L,jL);
            Qc = blkdiag(Qc,jQc);
            H  = [H jH];
            Pinf = blkdiag(Pinf,jPinf);

            % Should the covariance parameters be inferred?
            if ~isempty(strfind(gp.infer_params, 'covariance'))

                % Add derivatives
                dF      = mblk(dF,jdF);
                dQc     = mblk(dQc,jdQc);
                dPinf   = mblk(dPinf, jdPinf);

            end

            % Set options
            isstable = isfield(p,'stationary') && (isstable && p.stationary);

        end

        % Number of partial derivatives (not including R)
        nparam = min(size(dF,3),numel(dF));

        % Gradient with respect to Gaussian likelihood function parameters
        if ~isempty(gp.lik.p.('sigma2')) && ...
                ~isempty(strfind(gp.infer_params, 'likelihood')) && ...
                isfield(gp.lik.fh,'trcov')

            % Include noise magnitude R into parameters
            nparam = nparam + 1;

            % Derivative of noise magnitude w.r.t. itself (1)
            dR = zeros(1,1,nparam);
            dR(end) = 1;

            % Derivatives of model matrices w.r.t. noise magnitude (0)
            dF(:,:,nparam)    = zeros(size(F));
            dQc(:,:,nparam)   = zeros(size(Qc));
            dPinf(:,:,nparam) = zeros(size(Pinf));

        else

            % Noise magnitude is not optimized
            dR = zeros(1,1,nparam);

        end

        % Sort values
        [x,ind] = sort(x(:));
        y = y(ind);

        % State dimension, number of data points and number of parameters
        n      = size(F,1);
        steps  = numel(y);

        % Allocate for results
        % edata  = 0;
        gdata = zeros(1,nparam);

        % Set up
        Z  = zeros(n);
        m  = zeros(n,1);
        P  = Pinf;
        dm = zeros(n,nparam);
        dP = dPinf;
        dt = -inf;
        QC = L*Qc*L';

        % Allocate space for expm results
        AA  = zeros(2*n,2*n,nparam);
        AAA = zeros(4*n,4*n,nparam);

        % Loop over all observations
        for k=1:steps

            % The previous time step
            dt_old = dt;

            % The time discretization step length
            if (k>1)
                dt = x(k)-x(k-1);
            else
                dt = 0;
            end

            % Loop through all parameters (Kalman filter prediction step)
            for j=1:nparam

                % Should we recalculate the matrix exponential?
                if abs(dt-dt_old) > 1e-9

                    % The first matrix for the matrix factor decomposition
                    FF = [ F        Z;
                        dF(:,:,j) F];

                    % Solve the matrix exponential
                    AA(:,:,j) = expm2(FF*dt);

                end

                % Solve the differential equation
                foo     = AA(:,:,j)*[m; dm(:,j)];
                mm      = foo(1:n,:);
                dm(:,j) = foo(n+(1:n),:);

                if isstable

                    % For stable systems we can use the method by Davison,
                    % which simplifies everything (and speeds things up).

                    % The discrete-time dynamical model
                    if (j==1)
                        A  = AA(1:n,1:n,j);
                        Q  = Pinf - A*Pinf*A';
                        PP = A*P*A' + Q;
                    end

                    % The derivatives of A and Q
                    dA = AA(n+1:end,1:n,j);
                    %dQ = dPinf(:,:,j) - dA*Pinf*A' - A*dPinf(:,:,j)*A' - A*Pinf*dA';
                    dAPinfAt = dA*Pinf*A';
                    dQ = dPinf(:,:,j) - dAPinfAt - A*dPinf(:,:,j)*A' - dAPinfAt';

                    % The derivatives of P
                    %dP(:,:,j) = dA*P*A' + A*dP(:,:,j)*A' + A*P*dA' + dQ;
                    dAPAt = dA*P*A';
                    dP(:,:,j) = dAPAt + A*dP(:,:,j)*A' + dAPAt' + dQ;

                else

                    % The more general way for closed-form integration of
                    % the covariance by matrix fraction decomposition.

                    % Should we recalculate the matrix exponential?
                    if abs(dt-dt_old) > 1e-9

                        % Define W and G
                        W = L*dQc(:,:,j)*L';
                        G = dF(:,:,j);

                        % The second matrix for the matrix factor decomposition
                        FFF = [F  QC   Z   Z;
                            Z  -F'  Z   Z;
                            G  W    F   QC;
                            Z  -G'  Z   -F'];

                        % Solve the matrix exponential
                        AAA(:,:,j) = expm2(FFF*dt);

                    end

                    % Solve using matrix fraction decomposition
                    foo = AAA(:,:,j)*[P; eye(size(P)); dP(:,:,j); Z];

                    % Pick the parts
                    C  = foo(    (1:n),:);
                    D  = foo(  n+(1:n),:);
                    dC = foo(2*n+(1:n),:);
                    dD = foo(3*n+(1:n),:);

                    % The prediction step covariance
                    if j==1, PP = C/D; end

                    % Sove dP for j (C/D == P_{k|k-1})
                    dP(:,:,j) = (dC - PP*dD)/D;

                end
            end

            % Set predicted m and P
            m = mm;
            P = PP;

            % Start the Kalman filter update step and precalculate variables
            S = H*P*H' + R;
            [LS,notposdef] = chol(S,'lower');

            % If matrix is not positive definite, add jitter
            if notposdef>0,
                jitter = gp.jitterSigma2*diag(rand(size(S,1),1));
                [LS,notposdef] = chol(S+jitter,'lower');

                if notposdef>0,
                    gdata = nan; gprior = nan; g = nan;
                    return;
                end
            end

            % Continue update
            HtiS = H'/LS/LS';
            iS   = eye(size(S))/LS/LS';
            K    = P*HtiS;
            v    = y(k) - H*m;
            vtiS = v'/LS/LS';

            % Loop through all parameters (Kalman filter update step derivative)
            for j=1:nparam

                % Innovation covariance derivative
                dS = H*dP(:,:,j)*H' + dR(:,:,j);

                % Evaluate the energy derivative for j (optimized from above)
                gdata(j) = gdata(j) ...
                    + .5*sum(iS(:).*dS(:)) ...
                    - .5*(H*dm(:,j))*vtiS' ...
                    - .5*vtiS*dS*vtiS'     ...
                    - .5*vtiS*(H*dm(:,j));

                % Kalman filter update step derivatives
                dK        = dP(:,:,j)*HtiS - P*HtiS*dS/LS/LS';
                dm(:,j)   = dm(:,j) + dK*v - K*H*dm(:,j);
                dKSKt     = dK*S*K';
                dP(:,:,j) = dP(:,:,j) - dKSKt - K*dS*K' - dKSKt';

            end

            % Evaluate the energy (but only gradient required)
            %edata = edata + .5*size(S,1)*log(2*pi) + sum(log(diag(LS))) + .5*vtiS*v;

            % Finish Kalman filter update step
            m = m + K*v;
            P = P - K*S*K';

        end

        % Take all priors into account in the gradient
        gprior = [];
        if ~isempty(strfind(gp.infer_params, 'covariance'))
            for i=1:length(gp.cf)
                gprior = [gprior, -gp.cf{i}.fh.lpg(gp.cf{i})];
            end
        end

        gprior_lik = -gp.lik.fh.lpg(gp.lik);
        if ~isempty(gp.lik.p.('sigma2')) && ...
                ~isempty(strfind(gp.infer_params, 'likelihood')) && ...
                isfield(gp.lik.fh,'trcov')
            gprior = [gprior gprior_lik];
        end

        % Account for log transform
        w = gp_pak(gp);

        % Take correct values from w
        if length(gdata) < length(w)
            if any(hier==0)
                w = w([hier(1:find(hier==0,1)-1)==1 ... % Covariance function
                    hier(find(hier==0,1):find(hier==0,1)+...
                    length(gprior_lik)-1)==0]);      % Likelihood function
            else
                w = w(hier==1);
            end
        end

        % Log-transformation
        [w,ws] = gp_pak(gp);
        ind = strncmpi(ws,'log(',4) & ~strncmpi(ws,'log(log(',8);
        gdata(ind) = gdata(ind).*exp(w(ind));

        % Log-log-transformation
        ind = strncmpi(ws,'log(log(',8);
        gdata(ind) = gdata(ind).*exp(w(ind)+exp(w(ind)));

    otherwise
        error('Unknown type of Gaussian process!')
end

% If ther parameters of the model (covariance function parameters,
% likelihood function parameters, inducing inputs, mean function
% parameters) have additional hyperparameters that are not fixed,
% set the gradients in correct order
if length(gprior) > length(gdata)
    %gdata(gdata==0)=[];
    tmp=gdata;
    gdata = zeros(size(gprior));
    % Set the gradients to right place
    if any(hier==0)
        gdata([hier(1:find(hier==0,1)-1)==1 ...  % Covariance function
            hier(find(hier==0,1):find(hier==0,1)+length(gprior_lik)-1)==0 ... % Likelihood function
            hier(find(hier==0,1)+length(gprior_lik):end)==1 | ... % Inducing inputs or ...
            hier(find(hier==0,1)+length(gprior_lik):end)==-1]) = tmp; % Mean function parameters
    else
        if any(hier<0)
            hier(hier==-1)=1;
        end
        gdata(hier==1) = tmp;
    end
end
g = gdata + gprior;

function C = mblk(A,B)
%% mblk - Stack matrices by third dimension

% Get sizes
sA=size(A); sB=size(B);

% Check if A or B is empty
Ae = ~any(sA==0); Be = ~any(sB==0);

% Numel of sizes to 3
sA = [sA Ae]; sA = sA(1:3);
sB = [sB Be]; sB = sB(1:3);

% Assign space for C
C = zeros(sA+sB);

% Set values of A if there is any
if Ae
    C(1:sA(1), 1:sA(2), 1:sA(3)) = A;
end

% Set values of B if there is any
if Be
    C(sA(1)+(1:sB(1)), ...
        sA(2)+(1:sB(2)), ...
        sA(3)+(1:sB(3))) = B;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = sprice_4_2_fast(a, b, r, n, f, k, t,m, mu, nu, l, u, cp)
% sabr prices using the risk neutral density psabr_4
% and integrating this density with respect to the payoff
% fast version since cumsum is applied

eps = 0.001;
nl = length(k);
y = ones(1,nl);

rr = eps/10:eps:.5;

F1 = @(x) x .* psabr_4_2(a, b, r, n, f, x, t, m, mu, nu, l, u);
F2 = @(x) psabr_4_2(a, b, r, n, f, x, t, m, mu, nu, l, u);

y1 = F1(rr);
y2 = F2(rr);

y1 = cumsum(y1(end:-1:1)); y1 = y1(end:-1:1);
y2 = cumsum(y2(end:-1:1)); y2 = y2(end:-1:1);

    % call
for j = 1:nl;
    index = find(rr>=k(j),1,'first');
    y(j) = eps*(y1(index) - k(j) * y2(index));
end

if (cp ~= 1)
     y = k - f + y;
end

end


text:
function [ bestThresh ] = ThresholdSelection(  trainDir, image, s, widthOfBins, thresh, p)
% ThresholdSelection - This function is used for selecting the optimal
% threshold by comparing the
% fraction of the non-thresholded content which lies in the principal
% component.  The variable "p" should represent the fraction of the image
% inputed that contains the desired object... example: p=0.04.
%--------------------------------------------------------------------------
%   Params: trainDir - directory of training images.  Note there should be
%               a subdirectory in trainDir which contains .jpg images
%           s - the window size that each frame will be split up in to form
%               histograms
%           widthOfBins - the width of the bins for the RGB color
%               histograms
%           thresh - the cutoff distance threshold used to measure whether
%               or not window histograms are close enough to the training
%               histograms.
%           p - proportion of pixels containing object of interest
%
%   Returns: bestThresh - a guess for threshold value to choose
%--------------------------------------------------------------------------

    display(strcat(datestr(now,'HH:MM:SS'),' [INFO] Processing training images...'));
    trainingHistograms = BuildTrainingHistograms(trainDir, widthOfBins);

    display(strcat(datestr(now,'HH:MM:SS'),' [INFO] Reading image...'));
    image = double(image);
    bestThresh = 0;
    bestRatio = 0;
    for t = thresh
        disp('-------------------------------------------------------------');
        disp(strcat('Testing Threshold Value: ',num2str(t)))
        scoreImage = ImageToScoreArray( image, trainingHistograms, s, widthOfBins, t );
        pixels = prod(size(scoreImage)); %#ok<PSIZE>
        [L,num] = bwlabeln(scoreImage);
        max = 0;
        total = 0;
        for i = 1:num
            temp = sum(sum(L==i));
            total = total+temp;
            if (temp>max)
                max = temp;
            end
        end
        display(strcat('Fraction of pixels above threshold:',num2str(total/pixels)));
        display(strcat('Fraction of these pixels in principal component:',num2str(max/total)));
        if (0.9>total/pixels)
            if (total/pixels>p)
                if (max/total>bestRatio)
                    bestRatio = max/total;
                    bestThresh = t;
                end
            end
        end
    end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function genSimTests()

% Read template line by line into cell array
fid = fopen('simTemp.qmr');
allScript = textscan(fid,'%s','Delimiter','\n');
fclose(fid);
allScript = allScript{1}; % This is a cell aray that contains template

MethodList = list_models;
disp('Copy below section and paste them into Travis config.')
disp('------------------------------------------------------')
for ii = 1:length(MethodList)

    Model = str2func(MethodList{ii}); Model = Model();
    if Model.voxelwise

        newScript = replaceJoker('*-modelName-*',MethodList{ii},allScript,1);

        writeName = ['SimTest_' MethodList{ii} '_test' '.m'];

        fileID = fopen(writeName,'w');
        formatSpec = '%s\n';
        [nrows,~] = size(newScript);
        for row = 1:nrows
            fprintf(fileID,formatSpec,newScript{row,:});
        end
        fclose(fileID);

        disp('- stage: Test');
        disp(['  script: travis_wait 30 octave --no-gui --eval "bootstrapTest;cd(' char(39) '/home/travis/build/neuropoly/qMRLab/Test/MoxUnitCompatible' char(39) ');res=moxunit_runtests(' char(39) writeName char(39) ');exit(~res);"']);
    end

end

end

function newScript = replaceJoker(thisJoker,replaceWith,inScript,type)

% thisJoker: Joker definition *- joker_def_here -*
% replaceWith: Cell Aray or cell to be repaced with thisJoker
% inScript : Input script (thisJoker will be replaced in inSCript)
% type: 1 or 2: In-line replacements or block replacements

IndexC = strfind(inScript, char({thisJoker}));
idx = find(not(cellfun('isempty', IndexC)));

if type == 1 % In-line replacements

    for i=1:length(idx)

        inScript(idx(i)) = strrep(inScript(idx(i)),{thisJoker},{replaceWith});

    end

    newScript = inScript;

elseif type ==2 % Code blocks

    lenNew = length(replaceWith);
    newScript = cell(length(inScript)+lenNew-1,1);
    try
        newScript(idx:(idx+lenNew-1)) = replaceWith;
    catch
        newScript(idx:(idx+lenNew-1)) = replaceWith';
    end
    newScript(1:(idx-1)) = inScript(1:(idx-1));
    newScript((idx+lenNew):(length(inScript)+lenNew-1)) = inScript((idx+1):end);

else

end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% function allAcc = LRSearchLambdaSGD(Xtrain, Ytrain, Xvalidation, Yvalidation, lambdas)
% For each value of lambda provided, fit parameters to the training data and return
% the accuracy in the validation data in the corresponding entry of allAcc.
% For instance, allAcc(i) = accuracy in the validation set using lambdas(i).
%
% Inputs:
% Xtrain        training data features   (numTrainInstances x numFeatures)
% Ytrain        training set labels      (numTrainInstances x 1)
% Xvalidation   validation data features (numValidInstances x num features)
% Yvalidation   validation set labels    (numValidInstances x 1)
% lambdas       values of lambda to try  (numLambdas x 1)
%
% Output:
% allAcc        vector of accuracies in validation set  (numLambdas x 1)
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function allAcc = LRSearchLambdaSGD(Xtrain, Ytrain, Xvalidation, Yvalidation, lambdas)

  % You may use the functions we have provided such as LRTrainSGD, LRPredict, and LRAccuracy.

  allAcc = zeros(size(lambdas));

  %%%%%%%%%%%%%%
  %%% Student code
 for i = 1: length(lambdas)
	 theta = LRTrainSGD(Xtrain,Ytrain,lambdas(i));
	 pred = LRPredict(Xvalidation,theta);
	 allAcc(i) = LRAccuracy(pred,Yvalidation);
 end

  %%%%%%%%%%%

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [m] = region_sub(map,A,afterest,arr,col)
% this function is estimate initial guess for pixel independet method
% at field map(arr,col) by region growing method
% in
%   map         [n,n]   initial guess field map
%   A           [3*n,3] matrix of [1 arr col]
%   afterest    [n,n]   afterest(arr,col)=1 means we estimate
%                       field map(arr,col) by region growing method
%   arr,col             the position of the pixel now we try to estimate

%% initialize
window = zeros(512,512);
% make 41*41 window whose center is (arr,col)th pixel for region growing
% method
sa = arr-20;    % start array
ea = arr+20;    % end array
sc = col-20;    % start col
ec = col+20;    % end col

if sa<1
    sa = 1;
    ea = 41;
elseif ea>512
    ea = 512;
    sa = 472;
end
if sc<1
    sc = 1;
    ec = 41;
elseif ec>512
    ec = 512;
    sc = 472;
end

window(sa:ea,sc:ec) = 1;
maskgre = afterest.*window;
% we use only pixels which is already estimated and near(41*41) from
% (arr,col)th pixel

choose = find(maskgre==1);
map2 = map(choose);
w= map2.^2; % wieght

B = zeros(size(choose,1),3);

for i = 1:3
    B(:,i) = A(choose,i);
    C(:,i) = B(:,i).*w;
end

theta = inv(C'*B)*C'*map2; % the result of MSE

m = theta(1) + theta(2)*arr+theta(3)*col;
% m is new initla guess value for fieldmap(arr,col)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function x = bt3(p, adata)
  n=5;

  t1=p(1)-p(2);
  t2=p(2)+p(3)-2.0;
  t3=p(4)-1.0;
  t4=p(5)-1.0;

  for i=1:n
    x(i)=t1*t1 + t2*t2 + t3*t3 + t4*t4;
  end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function color_plot(x,colors,varargin)
% COLOR_PLOT    Scatterplot with colored points.
% color_plot(x) makes a scatterplot of x(:,1) versus x(:,2) with points colored
% according to quantiles of x(:,3).
% color_plot(x,n) specifies the number of color quantiles (default 4).
% color_plot(x,colors) specifies an RGB matrix of colors (the number of rows
% determines the number of quantiles).  The default is YlGnBu_colors.
% color_plot(...,'ColorBar',1) adds a color bar with tick marks from the
% quantile values.
%
% Example:
%   xy = ndgridmat(linspace(-12,12,20),linspace(-12,12,20));
%   z = sin(sqrt(xy(:,1).^2 + xy(:,2).^2));
%   color_plot([xy z]);
%
% See also YlGnBu_colors.

% Written by Tom Minka and Charles Sutton

args = makestruct(varargin);
default_args = struct('ColorBar',0,'MarkerSize',6);
args = setfields(default_args,args);

if nargin < 2
  colors = 4;
end
if length(colors) == 1
  nlevels = colors;
  colors = YlGnBu_colors(nlevels);
else
  nlevels = rows(colors);
end
% color groups
[c,q] = cut_quantile(x(:,3),nlevels);
for lev = 1:nlevels
  i = find(c == lev);
  plot(x(i,1),x(i,2),'o','Color',colors(lev,:),'MarkerFaceColor',colors(lev,:),'MarkerSize',args.MarkerSize);
  hold on
end
hold off

colormap(colors);

if args.ColorBar
  caxis ([0,1]);

  cTickLbls = cell(numel(q), 1);
  for i = 1:length(q)
      cTickLbls{i} = num2str(q(i), '%11.2g');
  end

  colorbar('YTick', linspace(0,1,nlevels+1), 'YTickLabel', cTickLbls);
end

set(gca,'Color','none')

author: tminka

repo: lightspeed

save_path: github-repos/MATLAB/tminka-lightspeed


text:
function [node,elem] = delmesh(node,elem,expr)
%% DELMESH delete part of the mesh

dim = size(node,2); elemdim = size(elem,2);
%% delete element
switch elemdim
    case 3
        center = (node(elem(:,1),:)+node(elem(:,2),:)+node(elem(:,3),:))/3;
    case 4
        center = (node(elem(:,1),:) + node(elem(:,2),:) ...
                + node(elem(:,3),:) + node(elem(:,4),:))/4;
end
x = center(:,1);  y = center(:,2); %#ok<*NASGU>
if dim == 3
	z = center(:,3); %#ok<*NASGU>
end
idx = eval(expr);
elem(idx,:) = [];

%% delete vertices
isValidNode = false(size(node,1),1);
isValidNode(elem(:)) = true;
node = node(isValidNode,:);

%% shift index of element
Nnew = sum(isValidNode);
indexMap(isValidNode) = (1:Nnew)';
elem = indexMap(elem);

repo: Aero-matlab


text:
function imdb=getImdb(Name_batch,theConf,meta,IsTrain)
trainRate=0.9;
posRate=0.75;

objset=readAnnotation(Name_batch,theConf);
objset_neg=getNegObjSet(theConf,Name_batch);
num_pos=length(objset);
num_neg=length(objset_neg);
tarN=round(posRate/(1-posRate)*num_neg);
repN=ceil(tarN/(num_pos*2));
list_train_pos=round(linspace(1,num_pos,round(num_pos*trainRate)));
list_train_pos=[list_train_pos.*2-1,list_train_pos.*2];
if(repN>1)
    list_train_pos=repmat(list_train_pos,[repN,1])+repmat((0:repN-1)'.*(num_pos*2),[1,numel(list_train_pos)]);
    list_train_pos=reshape(sort(list_train_pos(:)),[1,numel(list_train_pos)]);
end
list_train_pos=list_train_pos(list_train_pos<=tarN);
list_train_neg=round(linspace(1,num_neg,round(num_neg*trainRate)));
list_train=sort([list_train_pos,list_train_neg+tarN]);

image_size=meta.normalization.imageSize(1:2);
data=zeros(image_size(1),image_size(2),3,num_pos*2,'single');
data_neg=zeros(image_size(1),image_size(2),3,num_neg,'single');
for i=1:num_pos
    tar=i*2-1;
    IsFlip=false;
    [I_patch,~]=getI(objset(i),theConf,image_size,IsFlip);
    data(:,:,:,tar)=I_patch;
    tar=i*2;
    IsFlip=true;
    [I_patch,~]=getI(objset(i),theConf,image_size,IsFlip);
    data(:,:,:,tar)=I_patch;
end
data=repmat(data,[1,1,1,repN]);
data=data(:,:,:,1:tarN);
theConf_neg=theConf;
if(isfield(theConf_neg.data,'imgdir_neg'))
    theConf_neg.data.imgdir=theConf_neg.data.imgdir_neg;
end
for i=1:num_neg
    tar=i;
    IsFlip=false;
    [I_patch,~]=getI(objset_neg(i),theConf_neg,image_size,IsFlip);
    data_neg(:,:,:,tar)=I_patch;
end
total_images=tarN+num_neg;
imdb.images.data=zeros(image_size(1),image_size(2),3,total_images,'single');
list=1:tarN;
imdb.images.data(:,:,:,list)=data;
imdb.images.data(:,:,:,tarN+(1:num_neg))=data_neg;
clear data data_neg
imdb.images.labels=ones(1,total_images,'single').*(-1);
imdb.images.labels(1,list)=1;
imdb.images.set=zeros(1,total_images,'uint8');
imdb.images.set(1:end)=2;
imdb.images.set(list_train)=1;
dataMean=mean(imdb.images.data(:,:,:,imdb.images.set==1),4);
imdb.images.data=bsxfun(@minus,imdb.images.data,dataMean);
if(IsTrain)
    rng(0);
    list=randperm(total_images);
    imdb.images.data=imdb.images.data(:,:,:,list);
    imdb.images.labels=imdb.images.labels(:,list);
    imdb.images.set=imdb.images.set(:,list);
    [~,imdb.images.order]=sort(list);
end
imdb.meta.classes={Name_batch};
imdb.meta.sets={'train','val','test'} ;
imdb.meta.dataMean=dataMean;
end

repo: interpretableCNN