text:
function prime_plot ( n_max )

%*****************************************************************************80
%
%% PRIME_PLOT makes a box plot of primes.
%
%  Usage:
%
%    prime_plot ( n_max )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_MAX, the maximum number to be considered.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_PLOT\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display primes and composite numbers.\n' );
  fprintf ( 1, '  If N is composite, Box (N,D) is blue if N is divisible by D.\n' );
  fprintf ( 1, '  If N is prime, Column (N,:) is red.\n' );
%
%  First argument is N_MAX.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    n_max = input ( '  Enter N, the maximum number to check: ' );
  end

  clf
%
%  Every box starts out WHITE.
%
  for n = 1 : n_max

    for d = 1 : n_max

      x1 = n - 0.44;
      x2 = n + 0.44;
      y1 = d - 0.44;
      y2 = d + 0.44;

      patch ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], 'w', 'EdgeColor', 'none' );

    end

  end
%
%  By convention, 1 is NOT a prime!
%
  n = 1;
  d = 1;

  x1 = n - 0.44;
  x2 = n + 0.44;
  y1 = d - 0.44;
  y2 = d + 0.44;
  patch ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], 'b', 'EdgeColor', 'none' );

  for n = 2 : n_max

    prime = 1;

    d = 1;

    for d = 2 : min ( n_max, n - 1 )

      if ( mod ( n, d ) == 0 )
        x1 = n - 0.44;
        x2 = n + 0.44;
        y1 = d - 0.44;
        y2 = d + 0.44;
        patch ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], 'b', 'EdgeColor', 'none' );
        prime = 0;
      end

    end

    d = n;

    if ( prime )
      x1 = n - 0.44;
      x2 = n + 0.44;
      y1 = 1 - 0.44;
      y2 = d + 0.44;
      patch ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], 'r', 'EdgeColor', 'none' );
    else
      x1 = n - 0.44;
      x2 = n + 0.44;
      y1 = d - 0.44;
      y2 = d + 0.44;
      patch ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], 'b', 'EdgeColor', 'none' );
      x1 = n - 0.44;
      x2 = n + 0.44;
      y1 = 1 - 0.44;
      y2 = 1 + 0.44;
      patch ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], 'b', 'EdgeColor', 'none' );
    end

  end

  xlabel ( '<---N--->' )
  ylabel ( '<--Divisors-->' )
  title ( 'Primes are red columns; blue squares are factors.' )

  axis ( [ 0.5, n + 0.5, 0.5, n + 0.5 ] );
  axis equal
  axis square
  axis tight

  filename = 'prime_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved in file "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_PLOT\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer


text:
%RevoluteMDH Robot manipulator Revolute link class for MDH convention
%
% A subclass of the Link class: holds all information related to a robot
% link such as kinematics parameters, rigid-body inertial parameters, motor
% and transmission parameters.
%
% Notes::
% - This is reference class object
% - Link class objects can be used in vectors and arrays
% - Modified Denavit-Hartenberg parameters are used
%
% References::
% - Robotics, Vision & Control, Chap 7
%   P. Corke, Springer 2011.
%
% See also Link, Prismatic, SerialLink.
% See also Link, PrismaticMDH, Revolute, SerialLink.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

classdef RevoluteMDH < Link
    methods
        function L = RevoluteMDH(varargin)
            L = L@Link(varargin{:});

            if nargin == 0
                L.theta = [];
            end
            L.sigma = 0;
            if isempty(L.d)
                L.d = 0;
            end
            L.mdh = 1;
            if ~isempty(L.theta)
                error('theta cannot be specified for a revolute link');
            end
        end

    end
end

author: Allopart

repo: rbpf-gmapping


text:
function result = parseBC(N, BC, type)
%PARSEBC  Parse boundary conditions for CHEBOP object.
%   This method is not intended for end users. For information about boundary
%   conditions in CHEBOPs, see CHEBOP.
%
%   This method is invoked by the set methods for LBC, RBC, and BC. The types
%   and meanings of the allowed input are described in the documentation for
%   CHEBOP. The result is either empty or a function handle that represents the
%   given condition as needed internally.

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

numIn = nargin(N);

if ( isempty(BC) )
    result = [];

elseif ( isnumeric(BC) )
    % We reach here if we get assigned BCs like N.lbc = [1; 3].
    %   * In the scalar case, this means setting up conditions so that the
    %     solution equals the first entry of the vector at the endpoint, its
    %     derivative the second entry etc. In the example above, for a
    %     second order problem on [0,1], this amounts to
    %         u(0) = 1, u'(0) = 3.
    %   * In the system case, this means that the value of the first unknown
    %     function equals the first entry of the vector, the second unknown
    %     function the second entry of the vector, etc. In the example above,
    %     for a coupled system on [0,1], this amounts to
    %         u(0) = 1, v(0) = 3.
    %     Higher order coupled systems are not supported, as that would require
    %     inspecting the diffOrder of the linearized chebop, a (relatively)
    %     computionally expensive operation. That information only becomes
    %     available for free once we've called \ for starting to solve the
    %     problem.
    result = @(varargin) setupNumericalConditions(BC, varargin);

elseif ( isa(BC, 'function_handle') )
    % If we are dealing with a scalar problem where the independent variable is
    % not specified in the function handle arguments, allow also passing an
    % input function handle that takes one argument. Otherwise, we request that
    % the number of input to the BC function handle is one less than the number
    % of arguments to the OP part. A special case is when BCs are specified as a
    % vector, this leads to nargin(BC) = -1 (since we assigned the boundary
    % conditions as an anonymous function with a varargin argument), if this is
    % the case and the number of inputs for the BC don't match the operator, a
    % meaningful error will be thrown later on.
    if ( (numIn == 0) || ( (numIn == 1) && (nargin(BC) == 1) ) || ...
            ( strcmp(type,'lrbc') && (nargin(BC) == (numIn - 1)) ) || ...
            ( strcmp(type,'lrbc') && (nargin(BC) == -1) ) || ...
            ( strcmp(type,'bc') && (nargin(BC) == numIn) ) || ...
            ( strcmp(type,'bc') && (nargin(BC) == numIn + 1) ) )
        % We've got an anonymous function we're happy with. Do we need to
        % vectorize it?
        if ( N.vectorize && ~strcmp(func2str(BC), vectorize(BC)) )
            result = N.vectorizeOp(BC);
        else
            result = BC;
        end
    else
        error('CHEBFUN:CHEBOP:parseBC:inputs', ...
            'Number of inputs to BCs do not match operator.');
    end

elseif ( strcmpi(BC, 'neumann') )
    % Homogeneous Neumann.
    if ( numIn <= 2 )
        result = @(u) diff(u);
    else
        error('CHEBFUN:CHEBOP:parseBC:neuman', ...
            'Can only assign scalar BCs to scalar problems.');
    end

elseif ( strcmpi(BC, 'dirichlet') )
    % Homogeneous Dirichlet.
    if ( numIn <= 2 )
        result = @(u) u;
    else
        error('CHEBFUN:CHEBOP:parseBC:dirichlet', ...
            'Can only assign scalar BCs to scalar problems.');
    end

elseif ( iscell(BC) ) && ( length(BC) == 2 ) && ( ischar(BC{2}) ) && ...
        ( isnumeric(BC{1}) )
    % A cell may have a numerical value followed by a keyword of 'dirichlet' or
    % 'neumann'.

    % This behavior is retained only for backward compatibility and only for
    % problems with one variable.
    warning('CHEBFUN:CHEBOP:parseBC:keywordbc',...
        ['Keyword/value specifications of boundary conditions are ', ...
         'deprecated and may be removed in future versions of Chebfun.'])

    if ( strcmpi(BC{2}, 'neumann') && ( numIn <= 2 ) )
        result = @(u) diff(u) - BC{1};
    elseif ( strcmpi(BC{2}, 'dirichlet') && ( numIn <= 2 ) )
        result = @(u) u - BC{1};
    else
         error('CHEBFUN:CHEBOP:parseBC:cell', ...
            'Unable to parse cell input to set BC.');
    end

else
    error('CHEBFUN:CHEBOP:parseBC:unknown', 'Unsupported format of BCs.')

end

end

function out = setupNumericalConditions(BC, varargin)
% SETUPNUMERICALCONDITIONS  Return anonymous function for evaluating BCs.
%
% We need this function so that we can construct a function handle that allows
% specifying boundary conditions with a vector. In the scalar case, we will be
% subtracting the nth entry of the BC vector from the (n-1)st derivative of the
% solution. In the system case, we will be subtracting the nth entry of the
% vector from the nth unknown function at the endpoint:

u = varargin{1};
% Check if we're dealing with scalar or system case:
if (length(u) == 1)
    % Scalar case, subtract entries of BC from U and its derivatives:
    u = u{1};
    % We must always have at least one condition
    out = u - BC(1);
    for bcCounter = 2:length(BC)
        % Add to the boundary conditions vector. It's tricky to initialize the
        % OUT argument to the correct dimensions, as we need to be able to
        % evaluate this both with CHEBFUNs, ADCHEBFUNs and TREEVARs, hence a
        % cheeky growth of the array.
        out = [out; diff(u, bcCounter-1) - BC(bcCounter)]; %#ok<AGROW>
    end
else
    assert(length(u) == length(BC), ...
        'CHEBFUN:CHEBOP:parseBC:numberOfConditions', ...
        ['The number of initial/final conditions appears to be greater than' ...
        ' the number of variables in problem. Specifying boundary ' ...
        'conditions as a vector for systems is only supported for first ' ...
        'order systems.']);
    % We must always have at least one condition
    out = u{1} - BC(1);
    for bcCounter = 2:length(BC)
        % Add to the boundary conditions vector. Like above, it's tricky to
        % initialize the OUT argument to the correct dimensions, as we need to
        % be able to evaluate this both with CHEBFUNs, ADCHEBFUNs and TREEVARs.
        % So, allow ourselves to grow the array, by subtracting elements from
        % the BC vector from the solution components.
        out = [out; u{bcCounter} - BC(bcCounter)]; %#ok<AGROW>
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [  net,res,opts ] = sgd2(  net,res,opts )
% Stochastic gradient descent using second-order information.
%   Ye, C., Yang, Y., Fermuller, C., & Aloimonos, Y. (2017).
%   On the Importance of Consistency in Training Deep Neural Networks. arXiv preprint arXiv:1708.00631.

    if ~isfield(opts.parameters,'weightDecay')
        opts.parameters.weightDecay=1e-4;
    end
    if ~isfield(opts.parameters,'lambda_sgd2')
        opts.parameters.lambda_sgd2=1e0;
    end
    if ~isfield(opts.parameters,'large_matrix_inversion')
        opts.parameters.large_matrix_inversion=0;
    end
    if ~isfield(opts.parameters,'max_inv_size')
        opts.parameters.max_inv_size=500;
    end
    if ~isfield(opts.parameters,'clip')
        opts.parameters.clip=0;
    end
    if ~isfield(opts.parameters,'decorr_bias')
        opts.parameters.decorr_bias=1;
    end
    if ~isfield(net,'iterations')||(isfield(opts,'reset_mom')&&opts.reset_mom==1)
        net.iterations=0;
    end

    net.iterations=net.iterations+1;

    mom_factor=(1-opts.parameters.mom.^net.iterations);
    max_inv_size=opts.parameters.max_inv_size;
    lambda=opts.parameters.lambda_sgd2;

    for layer=1:numel(net.layers)
        if isfield(net.layers{layer},'weights')&&~isempty(net.layers{layer}.weights)
            if ~isfield(net.layers{layer},'momentum')||(isfield(opts,'reset_mom')&&opts.reset_mom==1)
                net.layers{layer}.momentum{1}=zeros(size(net.layers{layer}.weights{1}),'like',net.layers{layer}.weights{1});
                net.layers{layer}.momentum{2}=zeros(size(net.layers{layer}.weights{2}),'like',net.layers{layer}.weights{2});
            end

            dzdw=res(layer).dzdw;
            dzdb=res(layer).dzdb;

            if length(net.layers{layer}.weights)==2
                x=res(layer).x;
                batch_dim=length(size(x));%This assumes the batch size must be >1
                if batch_dim==4%2d cnn
                    x=permute(x,[3,1,2,4]);x=reshape(x,size(x,1),[]);
                    dzdw=permute(dzdw,[1,2,4,3]);new_size=size(dzdw);dzdw=reshape(dzdw,prod(new_size(1:3)),new_size(4));
                    K=size(dzdw,1)/numel(dzdb);dzdb=repelem(dzdb(:),K,1);
                end
                if batch_dim==3%1d cnn
                    x=permute(x,[2,1,3]);x=reshape(x,size(x,1),[]);
                    dzdw=permute(dzdw,[1,3,2]);new_size=size(dzdw);dzdw=reshape(dzdw,prod(new_size(1:2)),new_size(3));
                    K=size(dzdw,1)/numel(dzdb);dzdb=repelem(dzdb(:),K,1);
                end
                subsample=1;batch_size=size(x,2);
                if batch_size>1e4,subsample=ceil(min(50,batch_size/1e4));end
                if subsample>1,x=x(:,1:subsample:end);end
                if opts.parameters.decorr_bias==1
                    %insert bias
                    x=[ones(1,size(x,2),'like',x);x];
                    dzdw=[dzdb,dzdw];
                end
                if size(dzdw,2)<=max_inv_size %small scale inversion
                    dzdw=dzdw/(x*x'./size(x,2)+lambda*eye(size(x,1),'like',x));
                elseif opts.parameters.large_matrix_inversion %divide large scale into smaller scale
                    order=randperm(size(dzdw,2));
                    for i=1:max_inv_size:length(order) %could have been parallelized
                       block_size=min(max_inv_size,length(order)-i+1);
                       idx=order(i:i+block_size-1);x_tmp=x(idx,:);
                       dzdw(:,idx)=dzdw(:,idx)/(x_tmp*x_tmp'./size(x_tmp,2)+lambda*eye(size(x_tmp,1),'like',x));
                    end
                end
                if opts.parameters.decorr_bias==1
                    dzdb=dzdw(:,1);dzdw(:,1)=[];
                end
                if batch_dim==4,dzdw=reshape(dzdw,new_size);dzdw=permute(dzdw,[1,2,4,3]);end
                if batch_dim==3,dzdw=reshape(dzdw,new_size);dzdw=permute(dzdw,[1,3,2]);end
                if batch_dim>2%for cnn:
                    %dzdb is decorrelated with dzdw, take average to smooth the results.
                    dzdb=reshape(mean(reshape(dzdb(:),K,[]),1),size(res(layer).dzdb));
                end
            end

            if opts.parameters.clip>0
                mask=abs(res(layer).dzdw)>opts.parameters.clip;
                res(layer).dzdw(mask)=sign(res(layer).dzdw(mask)).*opts.parameters.clip;%%this type of processing seems to be very helpful
                mask=abs(res(layer).dzdb)>opts.parameters.clip;
                res(layer).dzdb(mask)=sign(res(layer).dzdb(mask)).*opts.parameters.clip;
            end

            %sgd updates
            net.layers{layer}.momentum{1}=opts.parameters.mom.*net.layers{layer}.momentum{1}-(1-opts.parameters.mom).*dzdw - opts.parameters.weightDecay * net.layers{layer}.weights{1};
            net.layers{layer}.weights{1}=net.layers{layer}.weights{1}+opts.parameters.lr*net.layers{layer}.momentum{1}./mom_factor;

            net.layers{layer}.momentum{2}=opts.parameters.mom.*net.layers{layer}.momentum{2}-(1-opts.parameters.mom).*dzdb;
            net.layers{layer}.weights{2}=net.layers{layer}.weights{2}+opts.parameters.lr*net.layers{layer}.momentum{2}./mom_factor;

        end
    end

    if ~isfield(opts,'reset_mom')||opts.reset_mom==1
        opts.reset_mom=0;
    end
end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function monogrid_poisson_1d_test01_multi ( )

%*****************************************************************************80
%
%% MONOGRID_POISSON_1D_TEST01_MULTI tests MULTIGRID_POISSON_1D on test case 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 December 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONOGRID_POISSON_1D_TEST01_MULTI\n' );
  fprintf ( 1, '  MULTIGRID_POISSON_1D solves a 1D Poisson BVP\n' );
  fprintf ( 1, '  using the multigrid method.\n' );

  a = 0.0;
  b = 1.0;
  ua = 0.0;
  ub = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  -u"(x) = 1, for %g < x < %g\n', a, b );
  fprintf ( 1, '  u(%g) = %g, u(%g) = %g.\n', a, ua, b, ub );
  fprintf ( 1, '  Solution is u(x) = ( -x^2 + x ) / 2\n' );

  for k = 5 : 5

    n = 2^k;
    x = ( linspace ( a, b, n + 1 ) )';

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Mesh index K = %d\n', k );
    fprintf ( 1, '  Number of intervals N=2^K = %d\n', n );
    fprintf ( 1, '  Number of nodes = 2^K+1 =   %d\n', n + 1 );

    [ u, it_num ] = multigrid_poisson_1d ( n, a, b, ua, ub, @force1, @exact1 );

    fprintf ( 1, '\n' );
    fprintf ( 1, '     I        X(I)      U(I)         U Exact(X(I))\n' );
    fprintf ( 1, '\n' );
    for i = 1 : n + 1
      fprintf ( 1, '  %4d  %10f  %14g  %14g\n', i, x(i), u(i), exact1 ( x(i) ) );
    end

    fprintf ( 1, '\n' );

    difmax = 0.0;
    for i = 1 : n + 1
      difmax = max ( difmax, abs ( u(i) - exact1 ( x(i) ) ) );
    end
    fprintf ( 1, '  Maximum error = %g\n', difmax );
    fprintf ( 1, '  Number of iterations = %d\n', it_num );

  end

  return
end

repo: jburkardt-m


text:
function value = p05_autonomous ( )

%*****************************************************************************80
%
%% P05_AUTONOMOUS reports whether problem P05 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 0;

  return
end


text:
function [V,F] = engraving(im,w,t,s,l)
  % ENGRAVING Create an engraving of an image
  %
  % [V,F] = engraving(im,w,t,s)
  %
  % Inputs:
  %   im  width by height grayscale image
  %   w   desired width of engraving model (in mm)
  %   t  desired thickness of engraving model (in mm)
  %   s  desired span of thickness devoted to levels (in mm)
  %   l  number of levels
  % Outputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of triangle indices
  %
  % Example:
  %  % load image
  %  im = imresize(rgb2gray(im2double(imread('hans-hass.jpg'))),0.5);
  %  % pad image by 10% of width
  %  im = padarray(im,ceil(0.1*repmat(size(im,2),1,2)));
  %  % engrave: 50mm wide, 5mm thick, 1mm devoted to 4 layers
  %  [V,F] = engraving(im,50,5,1,4);
  %

  assert(isfloat(im));
  [V,F] = create_regular_grid(size(im,2),size(im,1),0,0);
  V(:,1) = V(:,1)*w;
  V(:,2) = (1-V(:,2))*size(im,1)/size(im,2)*w;
  [V,F] = extrude(V,F);
  V(:,3) = V(:,3)*t;
  V(1:numel(im),3) = V(1:numel(im),3)-s*round(im(:)*(l-1))/(l-1);

end

author: alecjacobson

repo: gptoolbox


text:
function mm = minmax(data)
mm = zero_crossing(diff(data));

author: alexandrebarachant


text:
 function st = nufft_init(om, Nd, Jd, Kd, varargin)
%function st = nufft_init(om, Nd, Jd, Kd, [n_shift,] ...)
%|
%| Initialize structure for d-dimension NUFFT using KB interpolator,
%| particularly the interpolation matrix in sparse format.
%| caution: this routine can require a lot of memory!
%| in
%|	om [M d]	"digital" frequencies in radians
%|	Nd [d]		image dimensions (N1,N2,...,Nd)
%|	Jd [d]		# of neighbors used (in each direction)
%|	Kd [d]		FFT sizes (should be >= N1,N2,...)
%|
%| options
%|	n_shift [d]	n = 0-n_shift to N-1-n_shift (must be first)
%|	'minmax:kb'	minmax interpolator with excellent KB scaling!
%|				(minmax:kb is recommended, and used by default)
%|	'minmax:tuned'	minmax interpolator, somewhat numerically tuned
%|	'minmax:user'	minmax interpolator with user ({alpha}, {beta})
%|	'uniform'	uniform scaling factors (not recommended)
%|	'kaiser'	kaiser-bessel (KB) interpolator (minmax best alpha, m)
%|			or 'kaiser', [alpha], [m] to specify parameter (vectors)
%|	'linear'	linear interpolator (a terrible straw man)
%|	kernel		user-provided function handle interpolation kernel(k,J)
%|			(or a cell array of kernels, one for each dimension)
%|	'table'		use table-based interpolation rather than sparse matrix.
%|			this can save a lot of memory for large problems.
%|			example ..., 'table', 2^11, 'minmax:kb'
%|			where 2^11 is the table over-sampling factor.
%|
%| out
%|	st.p		[M *Kd]		sparse interpolation matrix
%|					(or empty if table-based)
%|	st.sn		[(Nd)]		scaling factors
%|	st.Nd,Jd,Kd,om	copies of inputs
%|
%| *Nd is shorthand for prod(Nd).
%| (Nd) is shorthand for (N1,N2,...,Nd)
%|
%| Like fft(), the NUFFT expects the signals to be x(0,0), ...
%| Use n_shift = [N1/2, N2/2, ...] for x(-N1/2,-N2/2,...), ...
%|
%| Copyright 2002-5-30, Jeff Fessler, University of Michigan

if nargin == 1 && streq(om, 'test'), nufft_init_test, return, end
if nargin < 4, ir_usage, end

% dimensionality of input space (usually 2 or 3)
dd = length(Nd);
if dd ~= length(Jd) || dd ~= length(Kd)
	fail 'inconsistent dim'
end
if any(Kd < Nd), warn 'Kd < Nd unlikely to work.  Try Kd=2*Nd', end
if any(Jd > Kd)
	Jd_new = min(Jd, Kd);
	fun = @(x) regexprep(num2str(x), '\s+', ' ');
	warn('Jd = [%s] > Kd = [%s]; changing to Jd = [%s]', ...
		fun(Jd), fun(Kd), fun(Jd_new))
	Jd = Jd_new; clear Jd_new fun
end

% special cases of input sampling pattern
if ischar(om)
	om = nufft_samples(om, Nd);
end
if dd ~= size(om,2), fail('omega needs %d columns', dd), end

%
% process optional arguments
%

% n_shift argument? (must be first)
if length(varargin) > 0 && isnumeric(varargin{1})
	n_shift = varargin{1};
	if dd ~= length(n_shift)
		fail('n_shift needs %d columns', dd)
	end
	varargin = {varargin{2:end}};
else
	n_shift = zeros(size(Nd));
end
st.n_shift = n_shift;

% default/recommended interpolator is minmax with KB scaling factors
if length(varargin) == 0
	varargin = {'minmax:kb'};
end

st.alpha = {};
st.beta = {};
is_kaiser_scale = false;

% table based?
if ischar(varargin{1}) && streq(varargin{1}, 'table')
	st = nufft_table_init(om, Nd, Jd, Kd, n_shift, varargin{2:end});
	return
end

ktype = varargin{1};

% cell array of kernel functions: {kernel1, kernel2, ..., kernelD}
if isa(ktype, 'cell')
	if isa(ktype{1}, 'inline') || isa(ktype{1}, 'function_handle')
		ktype = 'function_handle';
		if length(varargin) > 1, fail 'excess arguments?', end
		if length(varargin{1}) ~= dd, fail 'wrong # of kernels', end
		st.kernel = varargin{1};
	else
		fail 'cell array should be function_handles!?'
	end

% or a single function handle for all dimension
elseif isa(ktype, 'inline') || isa(ktype, 'function_handle')
	ktype = 'function_handle';
	if length(varargin) > 1, fail 'excess arguments?', end
	for id = 1:dd
		st.kernel{id} = varargin{1}; % all same
	end

% or a string that describes the type of interpolator
elseif ~ischar(ktype)
	fail 'non-string kernel type?'

end

st.ktype = ktype;

% set up whatever is needed for each interpolator
switch ktype
case 'function_handle'
	% already did it above

% linear interpolator straw man
case 'linear'
	ktype = 'function_handle';
	kernel = @(k,J) (1 - abs(k/(J/2))) .* (abs(k) < J/2);
	for id = 1:dd
		st.kernel{id} = kernel;
	end

% KB interpolator
case 'kaiser'
	is_kaiser_scale = true;

	% with minmax-optimized parameters
	if length(varargin) == 1
		for id = 1:dd
			[st.kernel{id} st.kb_alf(id) st.kb_m(id)] = ...
				kaiser_bessel('inline', Jd(id));
		end

	% with user-defined parameters
	elseif length(varargin) == 3
		alpha_list = varargin{2};
		m_list = varargin{3};
		if (length(alpha_list) ~= dd) || (length(m_list) ~= dd)
			fail('#alpha=%d #m=%d vs dd=%d', ...
				length(alpha_list), length(m_list), dd)
		end
		for id = 1:dd
			[st.kernel{id} st.kb_alf(id) st.kb_m(id)] = ...
				kaiser_bessel('inline', Jd(id), ...
					alpha_list(id), m_list(id));
		end
	else
		fail 'kaiser should have no arguments, or both alpha and m'
	end

% minmax interpolator with KB scaling factors (recommended default)
case 'minmax:kb'
	for id = 1:dd
		[st.alpha{id}, st.beta{id}] = ...
			nufft_alpha_kb_fit(Nd(id), Jd(id), Kd(id));
	end

% minmax interpolator with numerically "tuned" scaling factors
case 'minmax:tuned'
	for id = 1:dd
		[st.alpha{id}, st.beta{id}, ok] = ...
			nufft_best_alpha(Jd(id), 0, Kd(id)/Nd(id));
		if ~ok, fail 'unknown J,K/N', end
	end

% minmax interpolator with user-provided scaling factors
case 'minmax:user'
	if length(varargin) ~= 3, fail 'user must provide alpha/beta', end
	st.alpha = varargin{2};
	st.beta = varargin{3};
	if length(st.alpha) ~= dd || length(st.beta) ~= dd
		fail 'alpha/beta size mismatch'
	end

case 'uniform'
	for id = 1:dd
		st.alpha{id} = 1;
		st.beta{id} = 0;
	end

otherwise
	fail('unknown kernel type %s', ktype)
end

nufft_check_dft(om, Nd, ktype)

st.tol = 0;

st.Jd = Jd;
st.Nd = Nd;
st.Kd = Kd;

M = size(om,1);
st.M = M;
st.om = om;

% scaling factors: "outer product" of 1D vectors
st.sn = 1;
for id=1:dd
	if is_kaiser_scale
		nc = [0:Nd(id)-1]'-(Nd(id)-1)/2;
		tmp = 1 ./ kaiser_bessel_ft(...
			nc/Kd(id), Jd(id), st.kb_alf(id), st.kb_m(id), 1);
	elseif streq(ktype, 'function_handle')
		tmp = 1 ./ nufft_interp_zn(0, Nd(id), Jd(id), Kd(id), st.kernel{id});
	else
		tmp = nufft_scale(Nd(id), Kd(id), st.alpha{id}, st.beta{id});
	end
	st.sn = st.sn(:) * tmp';
end
if length(Nd) > 1
	st.sn = reshape(st.sn, Nd); % [(Nd)]
else
	st.sn = st.sn(:); % [(Nd)]
end

% [J? M] interpolation coefficient vectors.  will need kron of these later
for id=1:dd
	N = Nd(id);
	J = Jd(id);
	K = Kd(id);
	if isfield(st, 'kernel')
		[c, arg] = ...
		nufft_coef(om(:,id), J, K, st.kernel{id}); % [J? M]
	else
		alpha = st.alpha{id};
		beta = st.beta{id};
		T = nufft_T(N, J, K, st.tol, alpha, beta); % [J? J?]
		[r, arg] = ...
		nufft_r(om(:,id), N, J, K, alpha, beta); % [J? M]
		c = T * r;	clear T r
	end

	gam = 2*pi/K;
	phase_scale = 1i * gam * (N-1)/2;

	phase = exp(phase_scale * arg); % [J? M] linear phase
	ud{id} = phase .* c; % [J? M]

	% indices into oversampled FFT components
	koff = nufft_offset(om(:,id), J, K); % [M 1] to leftmost near nbr
	kd{id} = mod(outer_sum([1:J]', koff'), K) + 1; % [J? M] {1,...,K?}
	if id > 1 % trick: pre-convert these indices into offsets!
		kd{id} = (kd{id}-1) * prod(Kd(1:(id-1)));
	end

end, clear c arg gam phase phase_scale koff N J K

% build sparse matrix that is [M *Kd]
% with *Jd nonzero entries per frequency point
if dd >= 3
	tmp = prod(Jd)*M*8/10^9*2;
	if tmp > 1. % only display if more than 1GB
		printm('Needs at least %g Gbyte RAM', tmp)
	end
end

kk = kd{1}; % [J1 M]
uu = ud{1}; % [J1 M]
for id = 2:dd
	Jprod = prod(Jd(1:id));
	kk = block_outer_sum(kk, kd{id}); % outer sum of indices
	kk = reshape(kk, Jprod, M);
	uu = block_outer_prod(uu, ud{id}); % outer product of coefficients
	uu = reshape(uu, Jprod, M);
end % now kk and uu are [*Jd M]

% apply phase shift
% pre-do Hermitian transpose of interpolation coefficients
phase = exp(1i * (om * n_shift(:))).'; % [1 M]
uu = conj(uu) .* phase(ones(1,prod(Jd)),:); % [*Jd M]

mm = [1:M]; mm = mm(ones(prod(Jd),1),:); % [*Jd M]
% make sparse matrix, ensuring arguments are double for stupid matlab
st.p = sparse(mm(:), double(kk(:)), double(uu(:)), M, prod(Kd));
% sparse object, to better handle single precision operations!
st.p = Gsparse(st.p, 'odim', [M 1], 'idim', [prod(Kd) 1]);

% block_outer_sum()
%
% in
%	x1	[J1 M]
%	x2	[J2 M]
% out
%	y	[J1 J2 M]	y(i1,i2,m) = x1(i1,m) + x2(i2,m)
%
function y = block_outer_sum(x1, x2)
[J1 M] = size(x1);
[J2 M] = size(x2);
xx1 = reshape(x1, [J1 1 M]); % [J1 1 M] from [J1 M]
xx1 = xx1(:,ones(J2,1),:); % [J1 J2 M], emulating ndgrid
xx2 = reshape(x2, [1 J2 M]); % [1 J2 M] from [J2 M]
xx2 = xx2(ones(J1,1),:,:); % [J1 J2 M], emulating ndgrid
y = xx1 + xx2; % [J1 J2 M]

% block_outer_prod()
function y = block_outer_prod(x1, x2)
[J1 M] = size(x1);
[J2 M] = size(x2);
xx1 = reshape(x1, [J1 1 M]); % [J1 1 M] from [J1 M]
xx1 = xx1(:,ones(J2,1),:); % [J1 J2 M], emulating ndgrid
xx2 = reshape(x2, [1 J2 M]); % [1 J2 M] from [J2 M]
xx2 = xx2(ones(J1,1),:,:); % [J1 J2 M], emulating ndgrid
y = xx1 .* xx2; % [J1 J2 M]

% nufft_check_dft()
% see if they are DFT samples; if so, give warning
function nufft_check_dft(om, Nd, ktype)
kk = om / (2*pi) .* repmat(Nd(:)', [nrow(om) 1]);
tol = 1e-6;
if all(col(abs(round(kk) - kk)) < tol) ...
	&& any(om(:)) ... % trick: ignore all-zero om used in nufft_table_init()
	&& (streq(ktype, 'minmax:kb') || streq(ktype, 'kaiser'))
	warn('DFT samples with ktype="%s" has suboptimal accuracy', ktype)
end

% nufft_init_test()
% for a more complete test, use "nufft test"
function nufft_init_test
Nd = [20 10];
st = nufft_init('epi', Nd, [5 5], 2*Nd);
if 0
	clf
	om = st.om;
	plot(om(:,1), om(:,2), 'o-')
	axis_pipi set
end
st;
st.alpha{1};

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [j, dist] = LMhistintersectionquery(Query, h)
%
% Assumes histograms are normalized and sum = 1

N = size(h,1);

dist = 1-sum(min(h, repmat(Query, [N 1])),2); % histogram intersection
[dist, j] = sort(dist);

author: CSAILVision

repo: LabelMeToolbox


text:
  function d = max_percent_diff(s1, s2, varargin)
%|function d = max_percent_diff(s1, s2, [options])
%|
%| compute the "maximum percent difference" between two signals: s1, s2
%| options
%|	1		use both arguments as the normalizer
%|	string		print this
%|
%| Copyright 2000-9-16, Jeff Fessler, University of Michigan

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

base = '';

if ischar(s1)
	t1 = s1;
	s1 = evalin('caller', t1);
	base = [caller_name ': '];
else
	t1 = inputname(1);
end
if ischar(s2)
	t2 = s2;
	s2 = evalin('caller', t2);
else
	t2 = inputname(2);
end

use_both = 0;
doprint = (nargout == 0);
while (length(varargin))
	arg1 = varargin{1};
	if isnumeric(varargin{1})
		use_both = 1;
		varargin = {varargin{2:end}};
		continue
	end

	if ischar(arg1)
		doprint = 1;
		base = [arg1 ': '];
		varargin = {varargin{2:end}};
		continue
	end
end

% first check that we have comparable signals!
if ~isequal(size(s1), size(s2))
	printm(['size(%s) = %s'], inputname(1), mat2str(size(s1)))
	printm(['size(%s) = %s'], inputname(2), mat2str(size(s2)))
	error 'dimension mismatch'
end

if any(isnan(s1(:)))
	warn([mfilename ': NaN values in %s in %s!?'], t1, caller_name)
end
if any(isnan(s2(:)))
	warn([mfilename ': NaN values in %s in %s!?'], t2, caller_name)
end
s1 = doubles(s1);
s2 = doubles(s2);

if use_both
	denom = max(abs([s1(:); s2(:)]));
	if ~denom
		d = 0;
	else
		d = max(abs(s1(:)-s2(:))) / denom;
	end
else
	denom = max(abs(s1(:)));
	if ~denom
		denom = max(abs(s2(:)));
	end
	if ~denom
		d = 0;
	else
		d = max(abs(s1(:)-s2(:))) / denom;
	end
end

d = full(d) * 100;

if doprint
	printf([base 'max_percent_diff(%s, %s) = %g%%'], t1, t2, d)
	if ~nargout
		clear d
	end
end

%
% max_percent_diff_test
%
function max_percent_diff_test
v1 = [0 1000];
v2 = [0 1001];
max_percent_diff([0 1000], [1 1000])
max_percent_diff(v1, v2, 'numeric')
max_percent_diff '[0 1000]' '[1 1000]'


text:
function J=calcRuvConvJacob(zRUV,useHalfRange,lTx,lRx,M)
%%CALCRUVCONVJACOB Calculate the Jacobian for a monostatic or bistatic
%            range and direction cosines measurement in 3D, ignoring
%            atmospheric effects, with respect to Cartesian position.
%            This type of Jacobian is useful when performing tracking using
%            Cartesian-converted measurements where the clutter density is
%            specified in the measurement coordinate system, not the
%            converted measurement coordinate system.
%
%INPUTS: zRUV A 3X1 point in bistatic range and u and v coordinates in the
%           format [bistatic range;u;v].
% useHalfRange A boolean va0lue specifying whether the bistatic range value
%           should be divided by two. This normally comes up when operating
%           in monostatic mode, so that the range reported is a one-way
%           range. The default if this parameter is not provided, or an
%           empty matrix is passed, is false.
%       lTx The 3X1 [x;y;z] location vector of the transmitter in global
%           Cartesian coordinates. If this parameter is omitted or an
%           empty matrix is passed, then the transmitter is assumed to be
%           at the origin.
%       lRx The 3X1 [x;y;z] location vector of the receiver in Cartesian
%           coordinates. If this parameter is omitted or an empty matrix
%           is passed, then the receiver is assumed to be at the origin.
%         M A 3X3 rotation matrix to go from the alignment of the global
%           coordinate system to that at the receiver. The z-axis of the
%           local coordinate system of the receiver is the pointing
%           direction of the receiver. If omitted or an empty matrix is
%           passed, then it is assumed that the local coordinate system is
%           aligned with the global and M=eye(3) --the identity matrix is
%           used.
%
%OUTPUTS: J The 3X3 Jacobian matrix. Each row is a components of
%           [range;u;v] in that order with derivatives taken with respect
%           to [x,y,z] by column.
%
%This function converts the measurement into Cartesian coordinates and then
%calls rangeGradient and uvGradient.
%
%February 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(M))
    M=eye(3,3);
end

if(nargin<4||isempty(lRx))
    lRx=zeros(3,1);
end

if(nargin<3||isempty(lTx))
    lTx=zeros(3,1);
end

if(nargin<2||isempty(useHalfRange))
    useHalfRange=false;
end

x=ruv2Cart(zRUV,useHalfRange,lTx,lRx,M);

J=zeros(3,3);
J(1,:)=rangeGradient(x,useHalfRange,lTx,lRx);
J(2:3,:)=uvGradient(x,lRx,M);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Y = VEC(x)  Given an m x n matrix x, this produces the vector Y of length
%   m*n that contains the columns of the matrix x, stacked below each other.
%
% See also mat.

function x = vec(X)

%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

[m n] = size(X);
x = reshape(X,m*n,1);

author: andrewssobral

repo: lrslibrary


text:
function drawArrowhead(point,dir,rad,color)
%drawArrowhead Adds arrowhead to curve
%
% drawArrowhead(point,dir,rad,color)
%
%INPUTS
% point     Coordinates
% dir       Direction of arrowhead
% rad       Controls width of arrowhead
% color     Color of arrowhead
%
    global mapHandle
    global CB_MAP_OUTPUT
    angle = atan2(dir(2),dir(1));
    l = rad*.9;
    spread = pi/6*.9;
    x = [(point(1)+l*cos(angle+spread)) point(1) (point(1)+l*cos(angle-spread))];
    y = [(point(2)+l*sin(angle+spread)) point(2) (point(2)+l*sin(angle-spread))];
    if strcmp(CB_MAP_OUTPUT, 'matlab')
        if find(color>1)
            color = color/255;
        end
        fill(x,-y,color);
    elseif strcmp(CB_MAP_OUTPUT, 'java')
        %missing code
    elseif strcmp(CB_MAP_OUTPUT, 'svg')
        %determine type of color input
        if ischar(color)
            colorFill = color;
        else if isvector(color)
                colorFill = strcat('rgb(',num2str(color(1)),',',num2str(color(2)), ',',num2str(color(3)),')');
            end
        end
        colorStroke = colorFill;
        fprintf(mapHandle,'<path d="M %8.2f %8.2f L %8.2f %8.2f L %8.2f %8.2f z" fill="%s" stroke="%s"/>\n',x(1),y(1),x(2),y(2),x(3),y(3),colorFill,colorStroke);
    else
        display('no render found');
    end
end

author: opencobra

repo: cobratoolbox


text:
classdef CAMOEA < ALGORITHM
% <multi> <real/integer/label/binary/permutation>
% Clustering based adaptive multi-objective evolutionary algorithm

%------------------------------- Reference --------------------------------
% Y. Hua, Y. Jin, K. Hao, A clustering-based adaptive evolutionary
% algorithm for multiobjective optimization with irregular Pareto fronts,
% IEEE Transactions on Cybernetics, 2019, 49(7): 2758-2770.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Yicun Hua

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % Generate offspring randomly
                MatingPool = randperm(Problem.N);
                Offspring  = OperatorGA(Problem,Population(MatingPool));

                % Elitism strategy
                UniPop = [Population,Offspring];
                PopObj = UniPop.objs;
                [FrontNo,MaxFNo] = NDSort(PopObj,Problem.N);

                % The number of individuals to be selected in the last
                % non-dominated front
                K = Problem.N - sum(FrontNo<MaxFNo);

                if K ~= 0
                    % Normalization
                    pareto_population = find(FrontNo<MaxFNo);
                    last_population   = find(FrontNo == MaxFNo);
                    Zmin = min(PopObj(FrontNo == MaxFNo,:));
                    Zmax = max(PopObj(FrontNo == MaxFNo,:));
                    S = sum(FrontNo == MaxFNo);
                    MaxFnorm = (PopObj(last_population,:)-repmat(Zmin,S,1))./repmat(Zmax-Zmin,S,1);

                    % Clustering-based reference points generation
                    [Ref] = Reference_Generation( MaxFnorm,Problem.M,K);

                    % Clustering-based environmental selection
                    [reference_population] = Reference_Point_Selection(MaxFnorm,last_population,Ref,K,Problem.M);
                else
                    pareto_population    = find(FrontNo<=MaxFNo);
                    reference_population = [];
                end
                Population = UniPop([pareto_population,reference_population]);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function displayTrackingResult(sceneInfo, res, seqName)
% Display Tracking Result
%
% Take scene information sceneInfo and
% the tracking result from stateInfo
%
%
% (C) Anton Andriyenko, 2012
%
% The code may be used free of charge for non-commercial and
% educational purposes, the only requirement is that this text is
% preserved within the derivative work. For any other purpose you
% must contact the authors for permission. This code may not be
% redistributed without written permission from the authors.

%% convert res to stateInfo
stateInfo = [];
stateInfo.F = max(sceneInfo.frameNums);
stateInfo.frameNums = sceneInfo.frameNums;
index = 0;
cur_id = -1;

res = sortrows(res,2);

for i = 1:size(res,1)
    if (cur_id ~= res(i,2))
        cur_id = res(i,2);
        index = index + 1;
        stateInfo.X(:,index) = zeros(stateInfo.F,1);
        stateInfo.Y(:,index) = zeros(stateInfo.F,1);
        stateInfo.Xi(:,index) = zeros(stateInfo.F,1);
        stateInfo.Yi(:,index) = zeros(stateInfo.F,1);
        stateInfo.W(:,index) = zeros(stateInfo.F,1);
        stateInfo.H(:,index) = zeros(stateInfo.F,1);
    end
    bbox = res(i,:);
    n = bbox(1);
    stateInfo.X(n,index) = bbox(3)+0.5*bbox(5);
    stateInfo.Y(n,index) = bbox(4)+bbox(6);
    stateInfo.Xi(n,index) = stateInfo.X(n,index);
    stateInfo.Yi(n,index) = stateInfo.Y(n,index);
    stateInfo.W(n,index) = bbox(5);
    stateInfo.H(n,index) = bbox(6);
end

reopenFig(['Tracking Results of Sequence ' seqName]);
displayBBoxes(sceneInfo, stateInfo);

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function [vardistx, model] = vargplvmPartExpand(model, x, update)
% VARGPLVMPARTEXPAND
% This function is used when the model is using jointly the test and
% training data (possibly with missing values from the first) to optimise
% and in that case only part of the parameters are being optimised. This
% function, thus, expands the model by taking into account only the
% optimised parameters (the old ones are being kept fixed).

% In this version of the function the params that are considered to be
% optimised again are the var. distribution, dyn. kernel hyperparams and
% inducing points. The rest (base kernel and beta) are being kept fixed.

% x = [mu_bar lambda theta_t X_u] where mu_bar and lambda are augmented to
% include the vardist. of the test points.

% Note on the order of the parameters:
% [(N*Q*2)              |theta_t| Q*k |theta_f| 1] -> vargplvmExtractParam(model)
% [(N*Q*2)+(Nstar*Q*2)  |theta_t| Q*k]             -> x
%
% COPYRIGHT: Andreas C. Damianou, Michalis K. Titsias 2011
% VARGPLVM

% The model params before optimisation
paramsOrig = vargplvmExtractParam(model);

% The optimised parameters x here are not only the vardist. but also
% gDynKern and inducing points.
paramsNew = x;

% x will now hold only the vardist. params.
vardistxNparams=length(model.dynamics.t_star) * model.q * 2;
x = paramsNew(1:(model.dynamics.vardist.nParams + vardistxNparams));

% now separate the variational disribution into the training part and the
% testing part and update the original training model (only with the new training
% variational distribution) and the test variational distribution
% this is doing the expand
vardistNumData = size(x,2)/(model.q*2);
x = reshape(x, vardistNumData, model.dynamics.q*2);
xtrain = x(1:model.N,:);
xtest = x(model.N+1:end,:);
%model.dynamics.vardist = vardistExpandParam(model.dynamics.vardist, xtrain);
vardistx = vardistExpandParam(model.vardistx, xtest(:)');

% The model will be expanded with the origParams vector where all the
% new values (optimised) will overwrite the old ones.
params = paramsOrig;

startVal = 1;
endVal = model.dynamics.vardist.nParams;

% This will overwrite the vardist original params with the newly optimised
% ones
params(startVal:endVal) = xtrain(:)';

% This will overwrite the dyn.kern params and the inducing points.
startVal = endVal+1;
endVal = endVal + model.dynamics.kern.nParams + model.q * model.k;
params(startVal:endVal) = paramsNew((startVal+vardistxNparams):(endVal+vardistxNparams));

% Update the model
if exist('update') && update
    model = vargplvmExpandParam(model, params);
else
    % Update the model withoug calling update stats (this is used in the
    % optimisation likelihood and wrapper functions which do all the
    % precomputations themselves)
    model = vargplvmExpandParamNoUpdate(model, params);
end

%-----------

% This is copy-paste of vargplvmExpandParam without the last call to
% vargplvmUpdateStats (because this is done inside the objective and gradient functions)
function model = vargplvmExpandParamNoUpdate(model,params)
startVal = 1;
if isfield(model, 'dynamics') & ~isempty(model.dynamics)
    % variational parameters (reparametrized) AND dyn.kernel's parameters
    endVal = model.dynamics.nParams;
    model.dynamics = modelExpandParam(model.dynamics, params(startVal:endVal));
else
    % variational parameters (means and covariances), original ones
    endVal = model.vardist.nParams;
    model.vardist = modelExpandParam(model.vardist, params(startVal:endVal));
end

% inducing inputs
startVal = endVal+1;
if model.fixInducing
    if isfield(model, 'dynamics') & ~isempty(model.dynamics)
        Kt = kernCompute(model.dynamics.kern, model.dynamics.t);%%%%%%%%%%%%5
        model.X_u = Kt*model.dynamics.vardist.means; %dynamics
        model.X_u = model.X_u(model.inducingIndices,:);
    else
        model.X_u = model.vardist.means(model.inducingIndices, :); % static
    end
    % X_u values are taken from X values.
    % model.X_u = model.X(model.inducingIndices, :);
else
    % Parameters include inducing variables.
    endVal = endVal + model.q*model.k;
    model.X_u = reshape(params(startVal:endVal),model.k,model.q);
end

% kernel hyperparameters
startVal = endVal+1;
endVal = endVal + model.kern.nParams;
model.kern = kernExpandParam(model.kern, params(startVal:endVal));

% likelihood beta parameters
if model.optimiseBeta
    startVal = endVal + 1;
    endVal = endVal + prod(size(model.beta));
    fhandle = str2func([model.betaTransform 'Transform']);
    model.beta = fhandle(params(startVal:endVal), 'atox');
end

model.nParams = endVal;

author: SheffieldML

repo: vargplvm


text:
function d = setlabel(d, n)

%tstoolbox/@description/setlabel
%   Syntax:
%     * d = setlabel(d, label)
%
%   the label field of a description is used to give a signal some 'tag'
%   which remains constant through various processing steps
%   e.g. which topic this signal belongs to
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(2,2);

d.label = n;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ d2q ] = gen_DirectDynamics(q,Pcii,Icii,mcii,dq,taw)
%% This function is used to calculate the direct dynamics of the KUKA iiwa 7  R 800
% This function proposes that the robot is mounted with the base in the
% horizontal poistion.

% Arreguments:
%--------------------
% q: is 1x7 vector, joint nagles vector of the manipulator
% dq: is 1x7 vector, joint angular velocity vector
% taw: 1x7 vector, the torques on the joints.
% Pcii: is 3X7 matrix while each column represents the local coordinates
% of the center of mass of each link.
% Icii: is (3x3x7) matrix, each 3x3 matrix of which represnets the
% associated link inertial tensor represented in its local inertial frame
% mcii: is (1x7) vector, each element of which specifies a mass of one of
% the links

% Return value:
%--------------------
% d2q: is 7x1 vector, joint angular acceleration vector

% Copyright: Mohammad SAFEEA, 9th-April-2018

[M]=gen_MassMatrix(q,Pcii,Icii,mcii);
[B]=gen_CoriolisMatrix(q,Pcii,Icii,mcii,dq);
[ G ] = gen_GravityVector(q,Pcii,mcii);
% convert angular velocity/acceleration into column vectors
dq=columnVec(dq);
taw=columnVec(taw);
d2q=M\(taw-B*dq-G);
end

function y=columnVec(x)
if(size(x,2)==1)
    y=x;
else
    y=x';
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
%WEIGHT_CALC    MATLAB file that performs weight calculations
%   Calls a file of ship data.  Prints an output file called weights.out.
%   For details see Subsection 7.2.2 in the book.
%   Companion file to Biran, A. (2003), Ship Hydrostatics and Stability,
%   Oxford: Butterworth-Heinemann.

!rename weights.out weights.old      % prepare space for new data file
disp('Enter name of data file, then write RETURN and press ENTER ')
keyboard
fname = 'weights.out';
fid = fopen(fname, 'w');
title = [ sname ', weight calculations' ];
fprintf(fid, '%40s\n', title)
fprintf(fid, '---------------------------------------------------------------------------\n')
fprintf(fid, '     Weight item          Mass      vcg      z-Moment    vcg      x-Moment\n')
fprintf(fid, '----------------------------------------------------------------------------\n')

Displ = sum(wdata(:, 1));
vmom  = wdata(:, 1).*wdata(:, 2);
KG    = sum(vmom)/Displ;
lmom  = wdata(:, 1).*wdata(:, 3);
LCG   = sum(lmom)/Displ;

head = '%16s %11.2f %7.2f %12.2f %8.2f %13.2f\n';
[ m, n ] = size(wdata);
for k = 1:m
    name = names(k, :);
    fprintf(fid, head, name, wdata(k, 1), wdata(k, 2), vmom(k), wdata(k, 3), lmom(k))
end
fprintf(fid, '---------------------------------------------------------------------------\n')
subtitle = 'Total           ';
fprintf(fid, head, subtitle, Displ, KG, sum(vmom), LCG, sum(lmom))

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function lossLayer = addLossLayer(name, loss, scaleFactor)

switch lower(loss)
    case 'mse'
        lossLayer = mseLossLayer(name, scaleFactor);
    case 'mae'
        lossLayer = maeLossLayer(name, scaleFactor);
    case 'mse + cosine'
        lossLayer = mse_cos_LossLayer(name, scaleFactor);
    case 'mae + cosine'
        lossLayer = mae_cos_LossLayer(name, scaleFactor);
    otherwise
        error('Wrong value in inSpace argument');
end

author: mahmoudnafifi


text:
function lik = lik_t(varargin)
%LIK_T  Create a Student-t likelihood structure
%
%  Description
%    LIK = LIK_T('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates Student-t likelihood structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    LIK = LIK_T(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood structure with the named parameters
%    altered with the specified values.
%
%    Parameters for Student-t likelihood [default]
%      sigma2       - scale squared [1]
%      nu           - degrees of freedom [4]
%      sigma2_prior - prior for sigma2 [prior_logunif]
%      nu_prior     - prior for nu [prior_fixed]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                  __ n
%      p(y|f, z) = || i=1 C(nu,s2) * (1 + 1/nu * (y_i - f_i)^2/s2 )^(-(nu+1)/2)
%
%      where nu is the degrees of freedom, s2 the scale and f_i the
%      latent variable defining the mean. C(nu,s2) is constant
%      depending on nu and s2.
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%

% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2011 Pasi Jyl�nki

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

  ip=inputParser;
  ip.FunctionName = 'LIK_T';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('sigma2',0.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('sigma2_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('nu',4, @(x) isscalar(x) && x>0);
  ip.addParamValue('nu_prior',prior_fixed, @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

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

  % Initialize parameters
  if init || ~ismember('sigma2',ip.UsingDefaults)
    lik.sigma2 = ip.Results.sigma2;
  end
  if init || ~ismember('nu',ip.UsingDefaults)
    lik.nu = ip.Results.nu;
  end
  % Initialize prior structure
  if init
    lik.p=[];
  end
  if init || ~ismember('sigma2_prior',ip.UsingDefaults)
    lik.p.sigma2=ip.Results.sigma2_prior;
  end
  if init || ~ismember('nu_prior',ip.UsingDefaults)
    lik.p.nu=ip.Results.nu_prior;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_t_pak;
    lik.fh.unpak = @lik_t_unpak;
    lik.fh.lp = @lik_t_lp;
    lik.fh.lpg = @lik_t_lpg;
    lik.fh.ll = @lik_t_ll;
    lik.fh.llg = @lik_t_llg;
    lik.fh.llg2 = @lik_t_llg2;
    lik.fh.llg3 = @lik_t_llg3;
    lik.fh.tiltedMoments = @lik_t_tiltedMoments;
    lik.fh.tiltedMoments2 = @lik_t_tiltedMoments2;
    lik.fh.siteDeriv = @lik_t_siteDeriv;
    lik.fh.siteDeriv2 = @lik_t_siteDeriv2;
    lik.fh.optimizef = @lik_t_optimizef;
    lik.fh.upfact = @lik_t_upfact;
    lik.fh.invlink = @lik_t_invlink;
    lik.fh.predy = @lik_t_predy;
    lik.fh.predprcty = @lik_t_predprcty;
    lik.fh.recappend = @lik_t_recappend;
  end

end

function [w, s] = lik_t_pak(lik)
%LIK_T_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_T_PAK(LIK) takes a likelihood structure LIK and
%    combines the parameters into a single row vector W. This
%    is a mandatory subfunction used for example in energy and
%    gradient computations.
%
%       w = [ log(lik.sigma2)
%             (hyperparameters of lik.sigma2)
%             log(log(lik.nu))
%             (hyperparameters of lik.nu)]'
%
%  See also
%    LIK_T_UNPAK, GP_PAK

  w = []; s = {};
  if ~isempty(lik.p.sigma2)
    w = [w log(lik.sigma2)];
    s = [s; 'log(lik.sigma2)'];
    [wh sh] = lik.p.sigma2.fh.pak(lik.p.sigma2);
    w = [w wh];
    s = [s; sh];
  end
  if ~isempty(lik.p.nu)
    w = [w log(log(lik.nu))];
    s = [s; 'loglog(lik.nu)'];
    [wh sh] = lik.p.nu.fh.pak(lik.p.nu);
    w = [w wh];
    s = [s; sh];
  end
end

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

  if ~isempty(lik.p.sigma2)
    lik.sigma2 = exp(w(1));
    w = w(2:end);
    [p, w] = lik.p.sigma2.fh.unpak(lik.p.sigma2, w);
    lik.p.sigma2 = p;
  end
  if ~isempty(lik.p.nu)
    lik.nu = exp(exp(w(1)));
    w = w(2:end);
    [p, w] = lik.p.nu.fh.unpak(lik.p.nu, w);
    lik.p.nu = p;
  end
end

function lp = lik_t_lp(lik)
%LIK_T_LP  log(prior) of the likelihood parameters
%
%  Description
%    LP = LIK_T_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters.
%    This subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_T_LLG, LIK_T_LLG3, LIK_T_LLG2, GPLA_E

  v = lik.nu;
  sigma2 = lik.sigma2;
  lp = 0;

  if ~isempty(lik.p.sigma2)
    lp = lp + lik.p.sigma2.fh.lp(sigma2, lik.p.sigma2) +log(sigma2);
  end
  if ~isempty(lik.p.nu)
    lp = lp + lik.p.nu.fh.lp(lik.nu, lik.p.nu)  +log(v) +log(log(v));
  end
end

function lpg = lik_t_lpg(lik)
%LIK_T_LPG  d log(prior)/dth of the likelihood parameters th
%
%  Description
%    LPG = LIK_T_LPG(LIK) takes a likelihood structure LIK
%    and returns d log(p(th))/dth, where th collects the
%    parameters. This subfunction is needed when there are
%    likelihood parameters.
%
%  See also
%    LIK_T_LLG, LIK_T_LLG3, LIK_T_LLG2, GPLA_G

% Evaluate the gradients of log(prior)

  v = lik.nu;
  sigma2 = lik.sigma2;
  lpg = [];
  i1 = 0;

  if ~isempty(lik.p.sigma2)
    i1 = i1+1;
    lpg(i1) = lik.p.sigma2.fh.lpg(lik.sigma2, lik.p.sigma2).*sigma2 + 1;
  end
  if ~isempty(lik.p.nu)
    i1 = i1+1;
    lpg(i1) = lik.p.nu.fh.lpg(lik.nu, lik.p.nu).*v.*log(v) +log(v) + 1;
  end
end

function ll = lik_t_ll(lik, y, f, z)
%LIK_T_LL  Log likelihood
%
%  Description
%    LL = LIK_T_LL(LIK, Y, F) takes a likelihood structure LIK,
%    observations Y, and latent values F. Returns the log
%    likelihood, log p(y|f,z). This subfunction is needed when
%    using Laplace approximation or MCMC for inference with
%    non-Gaussian likelihoods. This subfunction is also used in
%    information criteria (DIC, WAIC) computations.
%
%  See also
%    LIK_T_LLG, LIK_T_LLG3, LIK_T_LLG2, GPLA_E

  r = y-f;
  v = lik.nu;
  sigma2 = lik.sigma2;

  term = gammaln((v + 1) / 2) - gammaln(v/2) -log(v.*pi.*sigma2)/2;
  ll = term + log(1 + (r.^2)./v./sigma2) .* (-(v+1)/2);
  ll = sum(ll);
end

function llg = lik_t_llg(lik, y, f, param, z)
%LIK_T_LLG  Gradient of the log likelihood
%
%  Description
%    LOKLIKG = LIK_T_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y, and latent values F. Returns
%    the gradient of log likelihood with respect to PARAM. At the
%    moment PARAM can be 'param' or 'latent'. This subfunction is
%    needed when using Laplace approximation or MCMC for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_T_LL, LIK_T_LLG2, LIK_T_LLG3, GPLA_E

  r = y-f;
  v = lik.nu;
  sigma2 = lik.sigma2;

  switch param
    case 'param'
      n = length(y);

      i1=0;
      if ~isempty(lik.p.sigma2)
        i1=i1+1;
        % Derivative with respect to sigma2
        llg(i1) = -n./sigma2/2 + (v+1)./2.*sum(r.^2./(v.*sigma2.^2+r.^2*sigma2));
        % correction for the log transformation
        llg(i1) = llg(i1).*sigma2;
      end
      if ~isempty(lik.p.nu)
        i1=i1+1;
        % Derivative with respect to nu
        llg(i1) = 0.5.* sum(psi((v+1)./2) - psi(v./2) - 1./v - log(1+r.^2./(v.*sigma2)) + (v+1).*r.^2./(v.^2.*sigma2 + v.*r.^2));

        % correction for the log transformation
        llg(i1) = llg(i1).*v.*log(v);
      end
    case 'latent'
      llg  = (v+1).*r ./ (v.*sigma2 + r.^2);
  end

end

function llg2 = lik_t_llg2(lik, y, f, param, z)
%LIK_T_LLG2  Second gradients of log likelihood
%
%  Description
%    LLG2 = LIK_T_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y, and latent values F. Returns
%    the Hessian of log likelihood with respect to PARAM. At the
%    moment PARAM can be only 'latent'. LLG2 is a vector with
%    diagonal elements of the Hessian matrix (off diagonals are
%    zero). This subfunction is needed when using Laplace
%    approximation or EP for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_T_LL, LIK_T_LLG, LIK_T_LLG3, GPLA_E

  r = y-f;
  v = lik.nu;
  sigma2 = lik.sigma2;

  switch param
    case 'param'

    case 'latent'
      % The Hessian d^2 /(dfdf)
      llg2 =  (v+1).*(r.^2 - v.*sigma2) ./ (v.*sigma2 + r.^2).^2;
    case 'latent+param'
      % gradient d^2 / (dfds2)
      llg2 = -v.*(v+1).*r ./ (v.*sigma2 + r.^2).^2;

      % Correction for the log transformation
      llg2 = llg2.*sigma2;
      if ~isempty(lik.p.nu)
        % gradient d^2 / (dfdnu)
        llg2(:,2) = r./(v.*sigma2 + r.^2) - sigma2.*(v+1).*r./(v.*sigma2 + r.^2).^2;

        % Correction for the log transformation
        llg2(:,2) = llg2(:,2).*v.*log(v);
      end
  end
end

function llg3 = lik_t_llg3(lik, y, f, param, z)
%LIK_T_LLG3  Third gradients of log likelihood (energy)
%
%  Description
%    LLG3 = LIK_T_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y and latent values F and
%    returns the third gradients of log likelihood with respect
%    to PARAM. At the moment PARAM can be only 'latent'. G3 is a
%    vector with third gradients. This subfunction is needed when
%    using Laplace approximation for inference with non-Gaussian
%    likelihoods.
%
%  See also
%    LIK_T_LL, LIK_T_LLG, LIK_T_LLG2, GPLA_E, GPLA_G

  r = y-f;
  v = lik.nu;
  sigma2 = lik.sigma2;

  switch param
    case 'param'

    case 'latent'
      % Return the diagonal of W differentiated with respect to latent values / dfdfdf
      llg3 = (v+1).*(2.*r.^3 - 6.*v.*sigma2.*r) ./ (v.*sigma2 + r.^2).^3;
    case 'latent2+param'
      % Return the diagonal of W differentiated with respect to
      % likelihood parameters / dfdfds2
      llg3 = (v+1).*v.*( v.*sigma2 - 3.*r.^2) ./ (v.*sigma2 + r.^2).^3;
      llg3 = llg3.*sigma2;
      if ~isempty(lik.p.nu)
        % dfdfdnu
        llg3(:,2) = (r.^2-2.*v.*sigma2-sigma2)./(v.*sigma2 + r.^2).^2 - 2.*sigma2.*(r.^2-v.*sigma2).*(v+1)./(v.*sigma2 + r.^2).^3;
        llg3(:,2) = llg3(:,2).*v.*log(v);
      end
  end
end

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

  zm = @zeroth_moment;

  tol = 1e-8;
  yy = y(i1);
  nu = lik.nu;
  sigma2 = lik.sigma2;

  % Set the limits for integration and integrate with quad
  % -----------------------------------------------------
  mean_app = myy_i;
  sigm_app = sqrt(sigm2_i);

  lambdaconf(1) = mean_app - 8.*sigm_app; lambdaconf(2) = mean_app + 8.*sigm_app;
  test1 = zm((lambdaconf(2)+lambdaconf(1))/2) > zm(lambdaconf(1));
  test2 = zm((lambdaconf(2)+lambdaconf(1))/2) > zm(lambdaconf(2));
  testiter = 1;
  if test1 == 0
    lambdaconf(1) = lambdaconf(1) - 3*sigm_app;
    test1 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(1));
    if test1 == 0
      go=true;
      while testiter<10 & go
        lambdaconf(1) = lambdaconf(1) - 2*sigm_app;
        lambdaconf(2) = lambdaconf(2) - 2*sigm_app;
        test1 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(1));
        test2 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(2));
        if test1==1&test2==1
          go=false;
        end
        testiter=testiter+1;
      end
    end
    mean_app = (lambdaconf(2)+lambdaconf(1))/2;
  elseif test2 == 0
    lambdaconf(2) = lambdaconf(2) + 3*sigm_app;
    test2 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(2));
    if test2 == 0
      go=true;
      while testiter<10 & go
        lambdaconf(1) = lambdaconf(1) + 2*sigm_app;
        lambdaconf(2) = lambdaconf(2) + 2*sigm_app;
        test1 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(1));
        test2 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(2));
        if test1==1&test2==1
          go=false;
        end
        testiter=testiter+1;
      end
    end
    mean_app = (lambdaconf(2)+lambdaconf(1))/2;
  end
  RTOL = 1.e-6;
  ATOL = 1.e-10;

  % Integrate with quadrature
  [m_0, m_1, m_2] = quad_moments(zm,lambdaconf(1), lambdaconf(2), RTOL, ATOL);

  sigm2hati1 = m_2 - m_1.^2;
  logM_0 = log(m_0);
  function integrand = zeroth_moment(f)
    r = yy-f;
    term = gammaln((nu + 1) / 2) - gammaln(nu/2) -log(nu.*pi.*sigma2)/2;
    integrand = exp(term + log(1 + r.^2./nu./sigma2) .* (-(nu+1)/2));
    integrand = integrand.*exp(- 0.5 * (f-myy_i).^2./sigm2_i - log(sigm2_i)/2 - log(2*pi)/2); %
  end
end

function [g_i] = lik_t_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_T_SITEDERIV  Evaluate the expectation of the gradient
%                 of the log likelihood term with respect
%                 to the likelihood parameters for EP
%
%  Description
%    [M_0, M_1, M2] = LIK_T_TILTEDMOMENTS(LIK, Y, I, S2, MYY)
%    takes a likelihood structure LIK, observations Y, index I
%    and cavity variance S2 and mean MYY. Returns E_f [d log
%    p(y_i|f_i) /d a], where a is the likelihood parameter and
%    the expectation is over the marginal posterior. This term is
%    needed when evaluating the gradients of the marginal
%    likelihood estimate Z_EP with respect to the likelihood
%    parameters (see Seeger (2008): Expectation propagation for
%    exponential families).  This subfunction is needed when using
%    EP for inference with non-Gaussian likelihoods and there are
%    likelihood parameters.
%
%  See also
%    GPEP_G

  zm = @zeroth_moment;
  znu = @deriv_nu;
  zsigma2 = @deriv_sigma2;

  tol = 1e-8;
  yy = y(i1);
  nu = lik.nu;
  sigma2 = lik.sigma2;

  % Set the limits for integration and integrate with quad
  mean_app = myy_i;
  sigm_app = sqrt(sigm2_i);

  lambdaconf(1) = mean_app - 6.*sigm_app; lambdaconf(2) = mean_app + 6.*sigm_app;
  test1 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(1));
  test2 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(2));
  testiter = 1;
  if test1 == 0
    lambdaconf(1) = lambdaconf(1) - 3*sigm_app;
    test1 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(1));
    if test1 == 0
      go=true;
      while testiter<10 & go
        lambdaconf(1) = lambdaconf(1) - 2*sigm_app;
        lambdaconf(2) = lambdaconf(2) - 2*sigm_app;
        test1 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(1));
        test2 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(2));
        if test1==1&test2==1
          go=false;
        end
        testiter=testiter+1;
      end
    end
    mean_app = (lambdaconf(2)+lambdaconf(1))/2;
  elseif test2 == 0
    lambdaconf(2) = lambdaconf(2) + 3*sigm_app;
    test2 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(2));
    if test2 == 0
      go=true;
      while testiter<10 & go
        lambdaconf(1) = lambdaconf(1) + 2*sigm_app;
        lambdaconf(2) = lambdaconf(2) + 2*sigm_app;
        test1 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(1));
        test2 = zm((lambdaconf(2)+lambdaconf(1))/2)>zm(lambdaconf(2));
        if test1==1&test2==1
          go=false;
        end
        testiter=testiter+1;
      end
    end
    mean_app = (lambdaconf(2)+lambdaconf(1))/2;
  end

  % Integrate with quad
  [m_0, fhncnt] = quadgk(zm, lambdaconf(1), lambdaconf(2));

  %         t=linspace(lambdaconf(1),lambdaconf(2),100);
  %         plot(t,zm(t))
  %         keyboard

  [g_i(1), fhncnt] = quadgk( @(f) zsigma2(f).*zm(f) , lambdaconf(1), lambdaconf(2));
  g_i(1) = g_i(1)/m_0*sigma2;

  if ~isempty(lik.p.nu)
    [g_i(2), fhncnt] = quadgk(@(f) znu(f).*zm(f) , lambdaconf(1), lambdaconf(2));
    g_i(2) = g_i(2)/m_0.*nu.*log(nu);
  end

  function integrand = zeroth_moment(f)
    r = yy-f;
    term = gammaln((nu + 1) / 2) - gammaln(nu/2) -log(nu.*pi.*sigma2)/2;
    integrand = exp(term + log(1 + r.^2./nu./sigma2) .* (-(nu+1)/2));
    integrand = integrand.*exp(- 0.5 * (f-myy_i).^2./sigm2_i - log(sigm2_i)/2 - log(2*pi)/2);
  end

  function g = deriv_nu(f)
    r = yy-f;
    temp = 1 + r.^2./nu./sigma2;
    g = psi((nu+1)/2)./2 - psi(nu/2)./2 - 1./(2.*nu) - log(temp)./2 + (nu+1)./(2.*temp).*(r./nu).^2./sigma2;
  end

  function g = deriv_sigma2(f)
    r = yy-f;
    g  = -1/sigma2/2 + (nu+1)./2.*r.^2./(nu.*sigma2.^2 + r.^2.*sigma2);
  end

end

function [lnZhat, muhat, sigm2hat] = lik_t_tiltedMoments2(likelih, y, yi, sigm2_i, myy_i, z, eta)
%LIKELIH_T_TILTEDMOMENTS    Returns the marginal moments for EP algorithm
%
%   Description
%   [M_0, M_1, M2] = LIKELIH_T_TILTEDMOMENTS(LIKELIH, Y, I, S2, MYY, Z)
%   takes a likelihood data structure LIKELIH, incedence counts Y,
%   expected counts Z, index I and cavity variance S2 and mean
%   MYY. Returns the zeroth moment M_0, mean M_1 and variance M_2
%   of the posterior marginal (see Rasmussen and Williams (2006):
%   Gaussian processes for Machine Learning, page 55). This subfunction
%   is needed when using robust-EP for inference with non-Gaussian
%   likelihoods.
%
%   See also
%   GPEP_E

  if nargin<7
    eta=1;
  end

  yy = y(yi);
  nu = likelih.nu;
  sigma2 = likelih.sigma2;
  sigma = sqrt(sigma2);

  nuprime = eta*nu+eta-1;
  a=nuprime/2; %a=nu/2;

  u=linspace(log(1e-8),5,200);
  du=u(2)-u(1);
  lnpu=(a-1)*u -a*exp(u)+u;

  %   sigma2 t-likelihood parameter, scale squared
  %   sigm2_i cavity variance
  %   myy_i cavity mean

  sigma2prime = sigma2*nu/nuprime;
  Vu = sigm2_i + (sigma2prime)./exp(u);
  lnZu = 0.5*(-log(2*pi*Vu)) -0.5 * (yy-myy_i)^2 ./Vu;
  lnZt = eta*gammaln((nu+1)/2) - eta/2*log(nu*pi*sigma2) - eta*gammaln(nu/2) - gammaln((nuprime+1)/2) + 0.5*log(nuprime*pi*sigma2prime) + gammaln(nuprime/2);

  ptu=exp(lnpu+lnZu+lnZt);

  Z_0=sum(ptu)*du;
  lnZhat=log(Z_0) + a*log(a)-gammaln(a);

  Vtu=1./(1/sigm2_i +(1/sigma2prime)*exp(u));
  mtu=Vtu.*(myy_i/sigm2_i + (yy/sigma2prime)*exp(u));

  muhat=sum(mtu.*ptu)*du/Z_0;
  sigm2hat=sum((Vtu+mtu.^2).*ptu)*du/Z_0-muhat^2;

  % limiting distribution (nu -> infinity)
%   Vg=1/(1/sigm2_i +eta/sigma2);
%   mg=Vg*(myy_i/sigm2_i +yy*eta/sigma2);
%   sigm_i=sqrt(sigm2_i);
%   sg=sqrt(Vg);
%
%   % set integration limits and scaling
%   nu_lim=1e10;
%   if nu<nu_lim
%
%     if sqrt(sigma2/sigm2_i)<0.05
%       % set the integration limits when the likelihood is very narrow
%
%       % grid resolution
%       dd=10;
%       df = [12*sigm_i/100 2*dd*sigma/100];
%
%       if yy>=myy_i
%         % grid break points
%         bp=[min(myy_i-6*sigm_i,yy-dd*sigma) myy_i-6*sigm_i, ...
%           min(myy_i+6*sigm_i,yy-dd*sigma), yy-dd*sigma, yy+dd*sigma,...
%           max(myy_i+6*sigm_i,yy+dd*sigma)];
%
%         % grid values
%         a=1e-6;
%         fvec =[ bp(1):df(2):bp(2)-a, bp(2):df(1):bp(3)-a, bp(3):max(df):bp(4)-a, ...
%           bp(4):df(2):bp(5)-a, bp(5):df(1):bp(6)];
%       else
%         % grid break points
%         bp=[min(myy_i-6*sigm_i,yy-dd*sigma), yy-dd*sigma, yy+dd*sigma,...
%           max(myy_i-6*sigm_i,yy+dd*sigma), myy_i+6*sigm_i, ...
%           max(myy_i+6*sigm_i,yy+dd*sigma)];
%
%         % grid values
%         a=1e-6;
%         fvec =[ bp(1):df(1):bp(2)-a, bp(2):df(2):bp(3)-a, bp(3):max(df):bp(4)-a, ...
%           bp(4):df(1):bp(5)-a, bp(5):df(2):bp(6)];
%       end
%
%       np=numel(fvec);
%       logpt = lpt(fvec,0);
%       lpt_max = max([logpt lpt([myy_i mg],0)]);
%       lambdaconf=[fvec(1), fvec(end)];
%       for i1=2:np-1
%         if logpt(i1) < lpt_max+log(1e-7) %(exp(logpt(i1))/exp(lpt_max) < 1e-7)
%           lambdaconf(1) = fvec(i1);
%         else
%           break;
%         end
%       end
%       for i1=1:np-2
%         if logpt(end-i1) < lpt_max+log(1e-7) %(exp(logpt(end-i1))/exp(lpt_max) < 1e-7)
%           lambdaconf(2) = fvec(end-i1);
%         else
%           break;
%         end
%       end
%     else
%       % set the integration limits in easier cases
%       np=20;
%       if mg>myy_i
%         lambdaconf=[myy_i-6*sigm_i,max(mg+6*sg,myy_i+6*sigm_i)];
%         fvec=linspace(myy_i,mg,np);
%       else
%         lambdaconf=[min(mg-6*sg,myy_i-6*sigm_i),myy_i+6*sigm_i];
%         fvec=linspace(mg,myy_i,np);
%       end
%       lpt_max=max(lpt(fvec,0));
%     end
%     C=log(1)-lpt_max; % scale the log-density for the quadrature tolerance
%   else
%     lambdaconf=[mg-6*sg,mg+6*sg];
%     C=log(1)-lpt(mg,0);
%   end
%
%   if nu>nu_lim
%     % the limiting Gaussian case
%     Vz=sigm2_i+sigma2/eta;
%     lnZhat = 0.5*(-log(eta) +(1-eta)*log(2*pi*sigma2) -log(2*pi*Vz)) -(0.5/Vz)*(yy-myy_i)^2;
%     muhat = mg;
%     sigm2hat = Vg;
%   else
%     % Integrate with quadrature
%     RTOL = 1.e-6;
%     ATOL = 1.e-7;
%     tic
%     [m_0, m_1, m_2] = quad_moments(@(f) exp(lpt(f,C)),lambdaconf(1), lambdaconf(2), RTOL, ATOL);toc
%     muhat = m_1;
%     sigm2hat = m_2 - m_1.^2;
%     lnZhat = log(m_0) -C;
%   end

  function lpdf = lpt(f,C)
    % logarithm of the tilted distribution
    r = yy-f;
    lpdf = gammaln((nu + 1) / 2) - gammaln(nu/2) -log(nu.*pi.*sigma2)/2;
    lpdf = lpdf + log(1 + r.^2./nu./sigma2) .* (-(nu+1)/2);
    lpdf = lpdf*eta - (0.5/sigm2_i) * (f-myy_i).^2 + (C-log(2*pi*sigm2_i)/2);
  end
end

function [g_i] = lik_t_siteDeriv2(likelih, y, yi, sigm2_i, myy_i, z, eta, lnZhat)
%LIKELIH_T_SITEDERIV   Evaluate the expectation of the gradient
%                           of the log likelihood term with respect
%                           to the likelihood parameters for EP
%
%   Description
%   [M_0, M_1, M2] = LIKELIH_T_TILTEDMOMENTS(LIKELIH, Y, I, S2, MYY)
%   takes a likelihood data structure LIKELIH, observations Y, index I
%   and cavity variance S2 and mean MYY. Returns E_f [d log
%   p(y_i|f_i) /d a], where a is the likelihood parameter and the
%   expectation is over the marginal posterior. This term is
%   needed when evaluating the gradients of the marginal
%   likelihood estimate Z_EP with respect to the likelihood
%   parameters (see Seeger (2008): Expectation propagation for
%   exponential families).  This subfunction is needed when using
%   robust-EP for inference with non-Gaussian likelihoods and there
%   are likelihood parameters.
%
%   See also
%   GPEP_G

  if nargin<7
    eta=1;
  end

  yy = y(yi);
  nu = likelih.nu;
  sigma2 = likelih.sigma2;
  sigma = sqrt(sigma2);

  % limiting distribution (nu -> infinity)
  Vg=1/(1/sigm2_i +eta/sigma2);
  mg=Vg*(myy_i/sigm2_i +yy*eta/sigma2);
  sigm_i=sqrt(sigm2_i);
  sg=sqrt(Vg);

  % set integration limits and scaling
  nu_lim=1e10;
  if nu<nu_lim

    if sqrt(sigma2/sigm2_i)<0.05
      % set the integration limits when the likelihood is very narrow

      % grid resolution
      dd=10;
      df = [12*sigm_i/100 2*dd*sigma/100];

      if yy>=myy_i
        % grid break points
        bp=[min(myy_i-6*sigm_i,yy-dd*sigma) myy_i-6*sigm_i, ...
          min(myy_i+6*sigm_i,yy-dd*sigma), yy-dd*sigma, yy+dd*sigma,...
          max(myy_i+6*sigm_i,yy+dd*sigma)];

        % grid values
        a=1e-6;
        fvec =[ bp(1):df(2):bp(2)-a, bp(2):df(1):bp(3)-a, bp(3):max(df):bp(4)-a, ...
          bp(4):df(2):bp(5)-a, bp(5):df(1):bp(6)];
      else
        % grid break points
        bp=[min(myy_i-6*sigm_i,yy-dd*sigma), yy-dd*sigma, yy+dd*sigma,...
          max(myy_i-6*sigm_i,yy+dd*sigma), myy_i+6*sigm_i, ...
          max(myy_i+6*sigm_i,yy+dd*sigma)];

        % grid values
        a=1e-6;
        fvec =[ bp(1):df(1):bp(2)-a, bp(2):df(2):bp(3)-a, bp(3):max(df):bp(4)-a, ...
          bp(4):df(1):bp(5)-a, bp(5):df(2):bp(6)];
      end

      np=numel(fvec);
      logpt = lpt(fvec,0);
      lpt_max = max([logpt lpt([myy_i mg],0)]);
      lambdaconf=[fvec(1), fvec(end)];
      for i1=2:np-1
        if logpt(i1) < lpt_max+log(1e-7) %(exp(logpt(i1))/exp(lpt_max) < 1e-7)
          lambdaconf(1) = fvec(i1);
        else
          break;
        end
      end
      for i1=1:np-2
        if logpt(end-i1) < lpt_max+log(1e-7) %(exp(logpt(end-i1))/exp(lpt_max) < 1e-7)
          lambdaconf(2) = fvec(end-i1);
        else
          break;
        end
      end
    else
      % set the integration limits in easier cases
      np=20;
      if mg>myy_i
        lambdaconf=[myy_i-6*sigm_i,max(mg+6*sg,myy_i+6*sigm_i)];
        fvec=linspace(myy_i,mg,np);
      else
        lambdaconf=[min(mg-6*sg,myy_i-6*sigm_i),myy_i+6*sigm_i];
        fvec=linspace(mg,myy_i,np);
      end
      lpt_max=max(lpt(fvec,0));
    end
    C=log(1)-lpt_max; % scale the log-density for the quadrature tolerance
  else
    lambdaconf=[mg-6*sg,mg+6*sg];
    C=log(1)-lpt(mg,0);
  end

  if nu>nu_lim
    % the limiting normal observation model
    Vz=sigm2_i+sigma2/eta;
    g_i(1) = 0.5*( (1-eta)/sigma2 -1/Vz/eta  + (yy-myy_i)^2 /Vz^2 /eta ) *sigma2/eta;

    if (isfield(likelih,'p') && ~isempty(likelih.p.nu))
      g_i(2) = 0;
    end
  else

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1e-7;

    % Integrate with quad
    %zm=@(f) exp(lpt(f,C));
    %[m_0, fhncnt] = quadgk(zm, lambdaconf(1), lambdaconf(2),'AbsTol',ATOL,'RelTol',RTOL)

    % Use the normalization determined in the lik_t_tiltedMoments2
    m_0=exp(lnZhat+C);

    zm=@(f) deriv_sigma2(f).*exp(lpt(f,C))*sigma2;
    [g_i(1), fhncnt] = quadgk( zm, lambdaconf(1), lambdaconf(2),'AbsTol',ATOL,'RelTol',RTOL);
    g_i(1) = g_i(1)/m_0;

    if (isfield(likelih,'p') && ~isempty(likelih.p.nu))
      zm=@(f) deriv_nu(f).*exp(lpt(f,C));
      [g_i(2), fhncnt] = quadgk( zm, lambdaconf(1), lambdaconf(2),'AbsTol',ATOL,'RelTol',RTOL);
      g_i(2) = g_i(2)/m_0.*nu.*log(nu);
    end

  end

  function lpdf = lpt(f,C)
    % logarithm of the tilted distribution
    r = yy-f;
    lpdf = gammaln((nu + 1) / 2) - gammaln(nu/2) -log(nu.*pi.*sigma2)/2;
    lpdf = lpdf + log(1 + r.^2./nu./sigma2) .* (-(nu+1)/2);
    lpdf = lpdf*eta - (0.5/sigm2_i) * (f-myy_i).^2 + (C-log(2*pi*sigm2_i)/2);
  end

  function g = deriv_nu(f)
    % derivative of the log-likelihood wrt nu
    r = yy-f;
    temp = r.^2 ./(nu*sigma2);
    g = psi((nu+1)/2) - psi(nu/2) - 1/nu;
    g = g + (1+1/nu).*temp./(1+temp);

    % for small values use a more accurate method for log(1+x)
    ii = temp<1e3;
    g(ii) = g(ii)  - log1p(temp(ii));
    g(~ii) = g(~ii) - log(1+temp(~ii));
    g = g*0.5;

  end

  function g = deriv_sigma2(f)
    % derivative of the log-likelihood wrt sigma2
    r = yy-f;
    temp = r.^2 /sigma2;
    g  = -1/sigma2/2 + ((1+1/nu)/2) * temp ./ (1 + temp/nu) /sigma2;
  end

end

function [f, a] = lik_t_optimizef(gp, y, K, Lav, K_fu)
%LIK_T_OPTIMIZEF  function to optimize the latent variables
%                 with EM algorithm
%
%  Description:
%    [F, A] = LIK_T_OPTIMIZEF(GP, Y, K, Lav, K_fu) Takes Gaussian
%    process structure GP, observations Y and the covariance
%    matrix K. Solves the posterior mode of F using EM algorithm
%    and evaluates A = (K + W)\Y as a sideproduct. Lav and K_fu
%    are needed for sparse approximations. For details, see
%    Vanhatalo, Jyl�nki and Vehtari (2009): Gaussian process
%    regression with Student-t likelihood. This subfunction is
%    needed when using lik_specific optimization method for mode
%    finding in Laplace algorithm.
%

  iter = 1;
  sigma2 = gp.lik.sigma2;
%  if sigma2==0
%    f=NaN;a=NaN;
%    return
%  end
  nu = gp.lik.nu;
  n = length(y);

  switch gp.type
    case 'FULL'
      iV = ones(n,1)./sigma2;
      siV = sqrt(iV);
      B = eye(n) + siV*siV'.*K;
      [L,notpositivedefinite] = chol(B);
      if notpositivedefinite
        f=NaN;a=NaN;
        return
      end
      B=B';
      b = iV.*y;
      a = b - siV.*(L'\(L\(siV.*(K*b))));
      f = K*a;
      while iter < 200
        fold = f;
        iV = (nu+1) ./ (nu.*sigma2 + (y-f).^2);
        siV = sqrt(iV);
        B = eye(n) + siV*siV'.*K;
        L = chol(B)';
        b = iV.*y;
        ws=warning('off','MATLAB:nearlySingularMatrix');
        a = b - siV.*(L'\(L\(siV.*(K*b))));
        warning(ws);
        f = K*a;

        if max(abs(f-fold)) < 1e-8
          break
        end
        iter = iter + 1;
      end
    case 'FIC'
      K_uu = K;

      Luu = chol(K_uu)';
      B=Luu\(K_fu');       % u x f

      K = diag(Lav) + B'*B;

      iV = ones(n,1)./sigma2;
      siV = sqrt(iV);
      B = eye(n) + siV*siV'.*K;
      L = chol(B)';
      b = iV.*y;
      a = b - siV.*(L'\(L\(siV.*(K*b))));
      f = K*a;
      while iter < 200
        fold = f;
        iV = (nu+1) ./ (nu.*sigma2 + (y-f).^2);
        siV = sqrt(iV);
        B = eye(n) + siV*siV'.*K;
        L = chol(B)';
        b = iV.*y;
        a = b - siV.*(L'\(L\(siV.*(K*b))));
        f = K*a;

        if max(abs(f-fold)) < 1e-8
          break
        end
        iter = iter + 1;
      end
  end

end

function upfact = lik_t_upfact(gp, y, mu, ll, z)
  nu = gp.lik.nu;
  sigma = sqrt(gp.lik.sigma2);
  sll = sqrt(ll);

  fh_e = @(f) t_pdf(f, nu, y, sigma).*norm_pdf(f, mu, sll);
  EE = quadgk(fh_e, -40, 40);

  fm = @(f) f.*t_pdf(f, nu, y, sigma).*norm_pdf(f, mu, sll)./EE;
  mm  = quadgk(fm, -40, 40);

  fV = @(f) (f - mm).^2.*t_pdf(f, nu, y, sigma).*norm_pdf(f, mu, sll)./EE;
  Varp = quadgk(fV, -40, 40);

  upfact = -(Varp - ll)./ll^2;
end

function [lpy, Ey, Vary] = lik_t_predy(lik, Ef, Varf, y, z)
%LIK_T_PREDY    Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_T_PREDY(LIK, EF, VARF YT)
%    Returns logarithm of the predictive density PY of YT, that is
%       p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires also the observations YT. This subfunction is
%    needed when computing posterior preditive distributions for
%    future observations.
%
%    [LPY, EY, VARY] = LIK_T_PREDY(LIK, EF, VARF) takes a likelihood
%    structure LIK, posterior mean EF and posterior Variance
%    VARF of the latent variable and returns the posterior
%    predictive mean EY and variance VARY of the observations
%    related to the latent variables. This subfunction is needed when
%    computing posterior preditive distributions for future observations.
%

%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  nu = lik.nu;
  sigma2 = lik.sigma2;
  sigma = sqrt(sigma2);

  Ey = zeros(size(Ef));
  EVary = zeros(size(Ef));
  VarEy = zeros(size(Ef));
  lpy = zeros(size(Ef));
  if nargout > 1
%       for i1=1:length(Ef)
%         %%% With quadrature
%         ci = sqrt(Varf(i1));
%
%         F = @(x) x.*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
%         Ey(i1) = quadgk(F,Ef(i1)-6*ci,Ef(i1)+6*ci);
%
%         F2 = @(x) (nu./(nu-2).*sigma2).*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
%         EVary(i1) = quadgk(F2,Ef(i1)-6*ci,Ef(i1)+6*ci);
%
%         F3 = @(x) x.^2.*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
%         VarEy(i1) = quadgk(F3,Ef(i1)-6*ci,Ef(i1)+6*ci) - Ey(i1).^2;
%       end
%       Vary = EVary + VarEy;

      Ey = Ef;
      if nu>2
        Vary=nu./(nu-2).*sigma2 +Varf;
      else
        warning('Variance of Student''s t-distribution is not defined for nu<=2')
        Vary=NaN+Varf;
      end
  end

  lpy = zeros(length(y),1);
  for i2 = 1:length(y)
    mean_app = Ef(i2);
    sigm_app = sqrt(Varf(i2));

    pd = @(f) t_pdf(y(i2), nu, f, sigma).*norm_pdf(f,Ef(i2),sqrt(Varf(i2)));
    lpy(i2) = log(quadgk(pd, mean_app - 12*sigm_app, mean_app + 12*sigm_app));
  end

end

function prctys = lik_t_predprcty(lik, Ef, Varf, zt, prcty)
%LIK_T_PREDPRCTY  Returns the percentiles of predictive density of y
%
%  Description
%    PRCTY = LIK_T_PREDPRCTY(LIK, EF, VARF YT, ZT)
%    Returns percentiles of the predictive density PY of YT. This
%    subfunction is needed when using function gp_predprcty.
%
%  See also
%    GP_PREDPCTY

  opt=optimset('TolX',1e-5,'Display','off');
  nt=size(Ef,1);
  prctys = zeros(nt,numel(prcty));
  prcty=prcty/100;
  nu = lik.nu;
  nu_p=max(2.5,nu);
  sigma2 = lik.sigma2;
  Vary=nu_p./(nu_p-2).*sigma2 +Varf;
  for i1=1:nt
    ci = sqrt(Varf(i1));
    for i2=1:numel(prcty)
      minf=sqrt(Vary(i1))*tinv(prcty(i2),nu)+(Ef(i1)-2.5*sqrt(Vary(i1)));
      maxf=sqrt(Vary(i1))*tinv(prcty(i2),nu)+(Ef(i1)+2.5*sqrt(Vary(i1)));
      a=(fminbnd(@(a) (quadgk(@(f) tcdf((a-f)/sqrt(Vary(i1)),nu).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)-prcty(i2)).^2,minf,maxf,opt));
%       a=(fminbnd(@(a) (quadgk(@(f) quadgk(@(y) t_pdf(y,nu,Ef(i1),sqrt(Vary(i1))),Ef(i1)-12*sqrt(Vary(i1)),a).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)-prcty(i2)).^2,minf,maxf,opt));
      prctys(i1,i2)=a;
      close all;
    end
  end
end

function mu = lik_t_invlink(lik, f, z)
%LIK_T_INVLINK  Returns values of inverse link function
%
%  Description
%    P = LIK_T_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values MU of inverse link function.
%    This subfunction is needed when using gp_predprctmu.
%
%     See also
%     LIK_T_LL, LIK_T_PREDY

  mu = f;
end

function reclik = lik_t_recappend(reclik, ri, lik)
%RECAPPEND  Record append
%  Description
%    RECCF = GPCF_SEXP_RECAPPEND(RECCF, RI, GPCF) takes old
%    covariance function record RECCF, record index RI, RECAPPEND
%    returns a structure RECCF. This subfunction is needed when
%    using MCMC sampling (gp_mc).

  if nargin == 2
    % Initialize the record
    reclik.type = 'Student-t';

    % Initialize parameters
    reclik.nu = [];
    reclik.sigma2 = [];

    % Set the function handles
    reclik.fh.pak = @lik_t_pak;
    reclik.fh.unpak = @lik_t_unpak;
    reclik.fh.lp = @lik_t_lp;
    reclik.fh.lpg = @lik_t_lpg;
    reclik.fh.ll = @lik_t_ll;
    reclik.fh.llg = @lik_t_llg;
    reclik.fh.llg2 = @lik_t_llg2;
    reclik.fh.llg3 = @lik_t_llg3;
    reclik.fh.tiltedMoments = @lik_t_tiltedMoments;
    reclik.fh.tiltedMoments2 = @lik_t_tiltedMoments2;
    reclik.fh.siteDeriv = @lik_t_siteDeriv;
    reclik.fh.siteDeriv2 = @lik_t_siteDeriv2;
    reclik.fh.optimizef = @lik_t_optimizef;
    reclik.fh.upfact = @lik_t_upfact;
    reclik.fh.invlink = @lik_t_invlink;
    reclik.fh.predy = @lik_t_predy;
    reclik.fh.predprcty = @lik_t_predprcty;
    reclik.fh.recappend = @lik_t_recappend;
    reclik.p.nu=[];
    if ~isempty(ri.p.nu)
      reclik.p.nu = ri.p.nu;
    end
    reclik.p.sigma2=[];
    if ~isempty(ri.p.sigma2)
      reclik.p.sigma2 = ri.p.sigma2;
    end
  else
    % Append to the record
    likp = lik.p;

    % record sigma2
    reclik.sigma2(ri,:) = lik.sigma2;
    if isfield(likp,'sigma2') && ~isempty(likp.sigma2)
      reclik.p.sigma2 = likp.sigma2.fh.recappend(reclik.p.sigma2, ri, likp.sigma2);
    end
    % record nu
    reclik.nu(ri,:) = lik.nu;
    if isfield(likp,'nu') && ~isempty(likp.nu)
      reclik.p.nu = likp.nu.fh.recappend(reclik.p.nu, ri, likp.nu);
    end
  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x, R] = orth_at(x, pos, dir, apply)
    %ORTH_AT Orthogonalize single core.
    %   X = ORTH_AT( X, POS, 'LEFT') left-orthogonalizes the core at position POS
    %   and multiplies the corresponding R-factor with core POS+1. All other cores
    %   are untouched. The modified tensor is returned.
    %
    %   X = ORTH_AT( X, POS, 'RIGHT') right-orthogonalizes the core at position POS
    %   and multiplies the corresponding R-factor with core POS-1. All other cores
    %   are untouched. The modified tensor is returned.
    %
    %   See also ORTHOGONALIZE.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    if ~exist('apply', 'var')
        apply = true;
    end

    sz = size(x.U{pos});

    if length(sz) == 2
        sz = [sz, 1];
    end

    if strcmpi(dir, 'left')
        [Q, R] = qr_unique(unfold(x.U{pos}, 'left'));
        x.U{pos} = reshape(Q, [sz(1), sz(2), size(Q, 2)]);

        if apply
            x.U{pos + 1} = tensorprod_ttemps(x.U{pos + 1}, R, 1);
        end

    elseif strcmpi(dir, 'right')
        % mind the transpose as we want to orthonormalize rows
        [Q, R] = qr_unique(unfold(x.U{pos}, 'right')');
        x.U{pos} = reshape(Q', [size(Q, 2), sz(2), sz(3)]);

        if apply
            x.U{pos - 1} = tensorprod_ttemps(x.U{pos - 1}, R, 3);
        end

    else
        error('Unknown direction specified. Choose either LEFT or RIGHT')
    end

end

author: NicolasBoumal

repo: manopt


text:
function a = companion_inverse ( n, x )

%*****************************************************************************80
%
%% COMPANION_INVERSE returns the inverse of the COMPANION matrix.
%
%  Example:
%
%    N = 5, X = ( 1, 2, 3, 4, 5 )
%
%    0    1    0    0    0
%    0    0    1    0    0
%    0    0    0    1    0
%    0    0    0    0    1
%   1/1 -5/1 -4/1 -3/1 -2/1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub, Charles Van Loan,
%    Matrix Computations, second edition,
%    Johns Hopkins University Press, Baltimore, Maryland, 1989,
%    section 7.4.6.
%
%    Charles Kenney, Alan Laub,
%    Controllability and stability radii for companion form systems,
%    Math. Control Signals Systems,
%    Volume 1, 1988, pages 239-256.
%
%    James Wilkinson,
%    The Algebraic Eigenvalue Problem,
%    Oxford University Press,
%    1965, page 12.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real X(N), the coefficients of the polynomial
%    which define the matrix.  X(1) must not be zero.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for j = 1 : n
    for i = 1 : n

      if ( i == n )

        if ( j == 1 )
          a(i,j) = 1.0 / x(1);
        else
          a(i,j) = - x(n+2-j) / x(1);
        end

      elseif ( i == j - 1 )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end

author: johannesgerer


text:
function [data_out]=ofdmsymbol_fft_cp(data_in,G,TxRx)

% We can make Nfft as a parameter as well
 Nfft = 256;

 if TxRx==10
      data=data_in'; %this makes the input column vector to a array

% %% making odfm symbol and taking IFFT
% function  symbol_ofdm = createsymbol (pilots,data)

%% now first generate the pilot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=complex(-1,0);
B=complex(1,0);
pilots = [A B A B B B A A];  %%% here we direct making the pilot,detail procedure given below

% n_symbol = 1;  % At the time of generating the pilots, I need to know what symbol i am simulating
%                % because the seed to do it depends on it.
% % The values of the pilots are to be modulated are defined in the standard as such(pp 443) :
% % Before beginning, it is necessary to consider that the value to calculate depends on 2 factors :
% % the number of symbols and whether we are in the uplink or dwnlink. We will consider that we are in
% % the downlink and we are transmitting the symbol "1".
% % If we want to consider the other uplink connection, the seed would be "10101010101".
% %
% seed = [1 1 1 1 1 1 1 1 1 1 1];
% for i=1:n_symbol+2
%     wk(i) = seed (11);
%     next = xor(seed(9),seed(11));
%     seed = [next seed(1,1:10)];
% end
%
% % Once the value of wk is found(that depends on the number of symbol with wihich it is working),
% % the values of the subcarriers must be found and of the mapping of them with BPSK constellation.
%
% wk = wk(n_symbol+2)
% A = 1 - 2*wk                           % Values defined in the standard.
% B = 1 - 2*(~wk)
% value_carrier = [A B A B B B A A]
%
% % For uplink, the values should be [A B A B A A A A]
%
% pilot_mapping = 2*mapping(value_carrier,1,Tx);
%
% % The factor of "2" is due to the fact that the pilots are transmitted to a
% % double power of the information bits.

% NOW The guard bands are prepared.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

guard1 = complex (0,0) * ones (1,28);
DC = complex (0,0);
guard2 = complex (0,0) * ones (1,27);

% The pilot and guard subcarriers are placed according to the standard.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbol_ofdm = [guard1 data(1:12) pilots(1) data(13:36)...
    pilots(2) data(37:60) pilots(3) data(61:84) pilots(4)...
    data(85:96) DC data(97:108) pilots(5) data(109:132) pilots(6)...
    data(133:156) pilots(7) data(157:180) pilots(8) data(181:192) guard2];

% here the ofdm symbol is completed
%Now taking IFFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 symbol_ofdm = sqrt(Nfft).*ifft(symbol_ofdm,Nfft);

 %Now adding cyclic prefix%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% we generate the cyclic prefix so that the multipaths do not affect our data so much.
margin = length(symbol_ofdm)*G;
data_tx = [symbol_ofdm((end-margin+1):end) symbol_ofdm];
data_out=data_tx;

   % At receiving end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif TxRx==01

data_rx=data_in;
% First, we must remove the CP.
margin = length(data_rx)*G;
margin = margin/(1+G);
symbol_ofdm_rx= data_rx(margin+1:end);

% After removing the CP, we have to inverse the IFFT, logcally by FFT.

 symbol_rx = fft(symbol_ofdm_rx,Nfft) ./ sqrt(Nfft);

 % Here the pilots are indicated, since i need to know where they are.
 % now we will be able to estimate the channel.

 pilots = [symbol_rx(41) symbol_rx(66) symbol_rx(91) symbol_rx(116) symbol_rx(142) symbol_rx(167) symbol_rx(192) symbol_rx(217)];

% After getting the received symbol, the channel is to be estimated using
% the pilot carriers. This is bypassed for an AWGN channel.

% Next, the values of the data and pilot carriers are extracted
 data_total = [symbol_rx(29:40) symbol_rx(42:65) symbol_rx(67:90) symbol_rx(92:115) symbol_rx(117:128)...
           symbol_rx(130:141) symbol_rx(143:166) symbol_rx(168:191) symbol_rx(193:216) symbol_rx(218:229)];
 data_out=data_total'; %this makes it a column vector
else
      disp('error in ofdmsymbol_fft_cp.m function');
end


text:
function f = abs(f, varargin)
%ABS   Absolute value of a DELTAFUN object.
%   ABS(F) returns the absolute value of F, where F is a DELTAFUN object such
%   that the fun part of F has no roots in the domain of F. If
%   ~isempty(roots(F)), then ABS(F) will return garbage with no warning.

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

% Take the absolute values:
f.funPart = abs(f.funPart, varargin{:});
f.deltaMag = abs(f.deltaMag);

end

author: chebfun

repo: chebfun