text:
function gpcf = gpcf_constant(varargin)
%GPCF_CONSTANT  Create a constant covariance function
%
%  Description
%    GPCF = GPCF_CONSTANT('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a constant covariance function structure in which the
%    named parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    GPCF = GPCF_CONSTANT(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for constant covariance function [default]
%      constSigma2       - magnitude (squared) [.1]
%      constSigma2_prior - prior for constSigma2 [prior_sqrtt]
%
%    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.
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*, MEAN_*
%
% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2010 Jaakko Riihimaki, Aki Vehtari
% Copyright (c) 2014 Arno Solin

% 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 = 'GPCF_CONSTANT';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('constSigma2',.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('constSigma2_prior',prior_sqrtt(), @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    init=true;
    gpcf.type = 'gpcf_constant';
  else
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_constant')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  % Initialize parameter
  if init || ~ismember('constSigma2',ip.UsingDefaults)
    gpcf.constSigma2=ip.Results.constSigma2;
  end

  % Initialize prior structure
  if init
    gpcf.p=[];
  end
  if init || ~ismember('constSigma2_prior',ip.UsingDefaults)
    gpcf.p.constSigma2=ip.Results.constSigma2_prior;
  end

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_constant_pak;
    gpcf.fh.unpak = @gpcf_constant_unpak;
    gpcf.fh.lp = @gpcf_constant_lp;
    gpcf.fh.lpg = @gpcf_constant_lpg;
    gpcf.fh.cfg = @gpcf_constant_cfg;
    gpcf.fh.cfdg = @gpcf_constant_cfdg;
    gpcf.fh.cfdg2 = @gpcf_constant_cfdg2;
    gpcf.fh.ginput = @gpcf_constant_ginput;
    gpcf.fh.ginput2 = @gpcf_constant_ginput2;
    gpcf.fh.ginput3 = @gpcf_constant_ginput3;
    gpcf.fh.ginput4 = @gpcf_constant_ginput4;
    gpcf.fh.cov = @gpcf_constant_cov;
    gpcf.fh.trcov  = @gpcf_constant_trcov;
    gpcf.fh.trvar  = @gpcf_constant_trvar;
    gpcf.fh.recappend = @gpcf_constant_recappend;
    gpcf.fh.cf2ss = @gpcf_constant_cf2ss;
  end

end

function [w, s, h] = gpcf_constant_pak(gpcf, w)
%GPCF_CONSTANT_PAK  Combine GP covariance function parameters into
%                   one vector.
%
%  Description
%    W = GPCF_CONSTANT_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W. This is a mandatory subfunction used for example
%    in energy and gradient computations.
%
%       w = [ log(gpcf.constSigma2)
%             (hyperparameters of gpcf.constSigma2)]'
%
%  See also
%    GPCF_CONSTANT_UNPAK

  w = []; s = {}; h=[];

  if ~isempty(gpcf.p.constSigma2)
    w = log(gpcf.constSigma2);
    s = [s 'log(constant.constSigma2)'];
    h = 1;
    % Hyperparameters of constSigma2
    [wh, sh, hh] = gpcf.p.constSigma2.fh.pak(gpcf.p.constSigma2);
    sh=strcat(repmat('prior-', size(sh,1),1),sh);
    w = [w wh];
    s = [s sh];
    h = [h 1+hh];
  end
end

function [gpcf, w] = gpcf_constant_unpak(gpcf, w)
%GPCF_CONSTANT_UNPAK  Sets the covariance function parameters
%                     into the structure
%
%  Description
%    [GPCF, W] = GPCF_CONSTANT_UNPAK(GPCF, W) takes a covariance
%    function structure GPCF and a parameter vector W, and
%    returns a covariance function structure identical to the
%    input, except that the covariance parameters have been set
%    to the values in W. Deletes the values set to GPCF from W
%    and returns the modified W. This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%    Assignment is inverse of
%       w = [ log(gpcf.constSigma2)
%             (hyperparameters of gpcf.constSigma2)]'
%
%  See also
%   GPCF_CONSTANT_PAK

  gpp=gpcf.p;
  if ~isempty(gpp.constSigma2)
    gpcf.constSigma2 = exp(w(1));
    w = w(2:end);
    % Hyperparameters of magnSigma2
    [p, w] = gpcf.p.constSigma2.fh.unpak(gpcf.p.constSigma2, w);
    gpcf.p.constSigma2 = p;
  end
end

function lp = gpcf_constant_lp(gpcf)
%GPCF_CONSTANT_LP  Evaluate the log prior of covariance function parameters
%
%  Description
%    LP = GPCF_CONSTANT_LP(GPCF) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LPG, GP_E

% Evaluate the prior contribution to the error. The parameters that
% are sampled are from space W = log(w) where w is all the
% "real" samples. On the other hand errors are evaluated in the
% W-space so we need take into account also the Jacobian of
% transformation W -> w = exp(W). See Gelman et al. (2013),
% Bayesian Data Analysis, third edition, p. 21.

  lp = 0;
  gpp=gpcf.p;
  if ~isempty(gpp.constSigma2)
    lp = gpp.constSigma2.fh.lp(gpcf.constSigma2, gpp.constSigma2) +log(gpcf.constSigma2);
  end
end

function lpg = gpcf_constant_lpg(gpcf)
%GPCF_CONSTANT_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_CONSTANT_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LP, GP_G

  lpg = [];
  gpp=gpcf.p;

  if ~isempty(gpcf.p.constSigma2)
    lpgs = gpp.constSigma2.fh.lpg(gpcf.constSigma2, gpp.constSigma2);
    lpg = [lpg lpgs(1).*gpcf.constSigma2+1 lpgs(2:end)];
  end
end

function DKff = gpcf_constant_cfg(gpcf, x, x2, mask, i1)
%GPCF_CONSTANT_CFG  Evaluate gradient of covariance function
%                   with respect to the parameters
%
%  Description
%    DKff = GPCF_CONSTANT_CFG(GPCF, X) takes a
%    covariance function structure GPCF, a matrix X of input
%    vectors and returns DKff, the gradients of covariance matrix
%    Kff = k(X,X) with respect to th (cell array with matrix
%    elements). This is a mandatory subfunction used in gradient
%    computations.
%
%    DKff = GPCF_CONSTANT_CFG(GPCF, X, X2) takes a
%    covariance function structure GPCF, a matrix X of input
%    vectors and returns DKff, the gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_CONSTANT_CFG(GPCF, X, [], MASK)
%    takes a covariance function structure GPCF, a matrix X of
%    input vectors and returns DKff, the diagonal of gradients of
%    covariance matrix Kff = k(X,X2) with respect to th (cell
%    array with matrix elements). This subfunction is needed when
%    using sparse approximations (e.g. FIC).
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LP, GP_G

  [n, m] =size(x);

  DKff = {};

  if nargin==5
    % Use memory save option
    if i1==0
      % Return number of hyperparameters
      if ~isempty(gpcf.p.constSigma2)
        DKff=1;
      else
        DKff=0;
      end
      return
    end
  end

  % Evaluate: DKff{1} = d Kff / d constSigma2
  %           DKff{2} = d Kff / d coeffSigma2
  % NOTE! Here we have already taken into account that the parameters are transformed
  % through log() and thus dK/dlog(p) = p * dK/dp

  % evaluate the gradient for training covariance
  if nargin == 2 || (isempty(x2) && isempty(mask))

    if ~isempty(gpcf.p.constSigma2)
      DKff{1}=ones(n)*gpcf.constSigma2;
    end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 || isempty(mask)
    if size(x,2) ~= size(x2,2)
      error('gpcf_constant -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end

    if ~isempty(gpcf.p.constSigma2)
      DKff{1}=ones([n size(x2,1)])*gpcf.constSigma2;
    end

    % Evaluate: DKff{1}    = d mask(Kff,I) / d constSigma2
    %           DKff{2...} = d mask(Kff,I) / d coeffSigma2
  elseif nargin == 4 || nargin == 5

    if ~isempty(gpcf.p.constSigma2)
      DKff{1}=ones(n,1)*gpcf.constSigma2; % d mask(Kff,I) / d constSigma2
    end
  end
  if nargin==5
    DKff=DKff{1};
  end

end

function DKff = gpcf_constant_cfdg(gpcf, x, x2, dims)
%GPCF_CONSTANT_CFDG  Evaluate gradient of covariance function, of
%                which has been taken partial derivative with
%                respect to x, with respect to parameters.
%
%  Description
%    DKff = GPCF_CONSTANT_CFDG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of derivatived covariance matrix
%    dK(df,f)/dhyp = d(d k(X,X)/dx)/dhyp, with respect to the
%    parameters
%
%    Evaluate: DKff{1:m} = d Kff / d coeffSigma2
%    m is the dimension of inputs. If ARD is used, then multiple
%    coefficients. This subfunction is needed when using derivative
%    observations.
%
%  See also
%    GPCF_CONSTANT_GINPUT

[n,m]=size(x);
if nargin<3
    x2=x;
end
if nargin < 4 || isempty(dims)
    dims = 1:m;
end
ii1=0;
DKff={};
if ~isempty(gpcf.p.constSigma2)
    dd=zeros(size(x,1),size(x2,1));
    for i=dims
        ii1=ii1+1;
        DKff{ii1}=dd;
    end
end
end

function DKff = gpcf_constant_cfdg2(gpcf, x, x2, dims1, dims2)
%GPCF_CONSTANT_CFDG2  Evaluate gradient of covariance function, of
%                 which has been taken partial derivatives with
%                 respect to both input variables x, with respect
%                 to parameters.
%
%  Description
%    DKff = GPCF_CONSTANT_CFDG2(GPCF, X) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of derivative covariance matrix
%    dK(df,df)/dhyp = d(d^2 k(X1,X2)/dX1dX2)/dhyp with respect to
%    the parameters
%
%    Evaluate: DKff{1:m} = d Kff / d coeffSigma
%    m is the dimension of inputs. If ARD is used, then multiple
%    lengthScales. This subfunction is needed when using derivative
%    observations.
%
%  See also
%   GPCF_CONSTANT_GINPUT, GPCF_CONSTANT_GINPUT2

if nargin < 4 || isempty(dims1)
    %dims1 = 1:m;
    error('dims1 needs to be given')
end
if nargin < 5 || isempty(dims2)
    %dims2 = 1:m;
    error('dims2 needs to be given')
end

[n,m]=size(x);
ii1=0;
%dd=zeros(size(x,1),size(x2,1));
dd=0;
DKff={};
if ~isempty(gpcf.p.constSigma2)
  ii1=ii1+1;
  DKff{ii1}=dd;
end
end

function DKff = gpcf_constant_ginput(gpcf, x, x2, i1)
%GPCF_CONSTANT_GINPUT  Evaluate gradient of covariance function with
%                      respect to x.
%
%  Description
%    DKff = GPCF_CONSTANT_GINPUT(GPCF, X) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_CONSTANT_GINPUT(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_CONSTANT_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    covariate in X. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LP, GP_G

  [n, m] =size(x);
  if nargin==4
    % Use memory save option
    if i1==0
      % Return number of covariates
      if isfield(gpcf,'selectedVariables')
        DKff=length(gpcf.selectedVariables);
      else
        DKff=m;
      end
      return
    end
  end

  if nargin == 2 || isempty(x2)
    ii1 = 0;
    for j = 1:n
      for i=1:m
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n);
      end
    end

  elseif nargin == 3 || nargin == 4
    %K = feval(gpcf.fh.cov, gpcf, x, x2);

    ii1 = 0;
    for j = 1:n
      for i=1:m
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n, size(x2,1));
        gprior(ii1) = 0;
      end
    end
  end
  if nargin==5
    DKff=DKff{1};
  end
end

function DKff = gpcf_constant_ginput2(gpcf, x, x2, dims, takeOnlyDiag)
%GPCF_CONSTANT_GINPUT2  Evaluate gradient of covariance function with
%                   respect to both input variables x and x2 (in
%                   same dimension).
%
%  Description
%    DKff = GPCF_CONSTANT_GINPUT2(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of twice derivatived covariance
%    matrix K(df,df) = dk(X1,X2)/dX1dX2 (cell array with matrix
%    elements). Input variable's dimensions are expected to be
%    same. The function returns also DKff1 and DKff2 which are
%    parts of DKff and needed with CFDG2. DKff = DKff1 -
%    DKff2. This subfunction is needed when using derivative
%    observations.
%
%  See also
%    GPCF_CONSTANT_GINPUT, GPCF_CONSTANT_GINPUT2, GPCF_CONSTANT_CFDG2

[n,m]=size(x);
if nargin<4 || isempty(dims)
    dims=1:m;
end
ii1=0;
if nargin==5 && isequal(takeOnlyDiag,'takeOnlyDiag')
    for i=dims
        ii1=ii1+1;
        DKff{ii1} = kron(0,zeros(n,1));
    end
    %DKff = kron(zeros(m,1),zeros(n,1));
else
    DK=zeros(size(x,1),size(x2,1));
    for i=dims
        ii1=ii1+1;
        DKff{ii1}=DK;
    end
end
end

function DKff = gpcf_constant_ginput3(gpcf, x, x2, dims1, dims2)
%GPCF_CONSTANT_GINPUT3  Evaluate gradient of covariance function with
%                   respect to both input variables x and x2 (in
%                   different dimensions).
%
%  Description
%    DKff = GPCF_CONSTANT_GINPUT3(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of twice derivatived covariance
%    matrix K(df,df) = dk(X1,X2)/dX1dX2 (cell array with matrix
%    elements). The derivative is calculated in multidimensional
%    problem between input's observation dimensions which are not
%    same. This subfunction is needed when using derivative
%    observations.
%
%  See also
%    GPCF_CONSTANT_GINPUT, GPCF_CONSTANT_GINPUT2, GPCF_CONSTANT_CFDG2

if nargin<4 || isempty(dims1)
    dims1=1:m;
end
if nargin<5 || isempty(dims2)
    dims2=1:m;
end

[n,m]=size(x);
ii1=0;
DK=zeros(size(x,1),size(x2,1));
for i=dims1
  for j=dims2
    ii1=ii1+1;
    DKff{ii1}=DK;
  end
end
end

function DKff = gpcf_constant_ginput4(gpcf, x, x2, dims)
%GPCF_CONSTANT_GINPUT  Evaluate gradient of covariance function with
%                  respect to x. Simplified and faster version of
%                  constant_ginput, returns full matrices.
%
%  Description
%    DKff = GPCF_CONSTANT_GINPUT4(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (whole matrix). This subfunction is needed when
%    using derivative observations.
%
%    DKff = GPCF_CONSTANT_GINPUT4(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (whole matrix). This subfunction
%    is needed when using derivative observations.
%
%  See also
%    GPCF_CONSTANT_PAK, GPCF_CONSTANT_UNPAK, GPCF_CONSTANT_LP, GP_G

[n,m]=size(x);
if nargin<4
    dims=1:m;
end
ii1=0;
if nargin==2
  x2=x;
end
DK=zeros(size(x,1),size(x2,1));
for i=dims
  ii1=ii1+1;
  DKff{ii1}=DK;
end
end

function C = gpcf_constant_cov(gpcf, x1, x2, varargin)
%GP_CONSTANT_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_CONSTANT_COV(GP, TX, X) takes in covariance function
%    of a Gaussian process GP and two matrixes TX and X that
%    contain input vectors to GP. Returns covariance matrix C.
%    Every element ij of C contains covariance between inputs i
%    in TX and j in X. This is a mandatory subfunction used for
%    example in prediction and energy computations.
%
%  See also
%    GPCF_CONSTANT_TRCOV, GPCF_CONSTANT_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  C = ones(n1,n2)*gpcf.constSigma2;
end

function C = gpcf_constant_trcov(gpcf, x)
%GP_CONSTANT_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_CONSTANT_TRCOV(GP, TX) takes in covariance function
%    of a Gaussian process GP and matrix TX that contains
%    training input vectors. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i and j
%    in TX. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%  See also
%    GPCF_CONSTANT_COV, GPCF_CONSTANT_TRVAR, GP_COV, GP_TRCOV

  n =size(x,1);
  C = ones(n,n)*gpcf.constSigma2;

end

function C = gpcf_constant_trvar(gpcf, x)
%GP_CONSTANT_TRVAR  Evaluate training variance vector
%
%  Description
%    C = GP_CONSTANT_TRVAR(GPCF, TX) takes in covariance function
%    of a Gaussian process GPCF and matrix TX that contains
%    training inputs. Returns variance vector C. Every
%    element i of C contains variance of input i in TX. This is
%    a mandatory subfunction used for example in prediction and
%    energy computations.
%
%  See also
%    GPCF_CONSTANT_COV, GP_COV, GP_TRCOV

  n =size(x,1);
  C = ones(n,1)*gpcf.constSigma2;

end

function reccf = gpcf_constant_recappend(reccf, ri, gpcf)
%RECAPPEND Record append
%
%  Description
%    RECCF = GPCF_CONSTANT_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains all
%    the old samples and the current samples from GPCF. This
%    subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_constant';

    % Initialize parameters
    reccf.constSigma2 = [];

    % Set the function handles
    reccf.fh.pak = @gpcf_constant_pak;
    reccf.fh.unpak = @gpcf_constant_unpak;
    reccf.fh.lp = @gpcf_constant_lp;
    reccf.fh.lpg = @gpcf_constant_lpg;
    reccf.fh.cfg = @gpcf_constant_cfg;
    reccf.fh.cfdg = @gpcf_constant_cfdg;
    reccf.fh.cfdg2 = @gpcf_constant_cfdg2;
    reccf.fh.ginput = @gpcf_constant_ginput;
    reccf.fh.ginput2 = @gpcf_constant_ginput2;
    reccf.fh.ginput3 = @gpcf_constant_ginput3;
    reccf.fh.ginput4 = @gpcf_constant_ginput4;
    reccf.fh.cov = @gpcf_constant_cov;
    reccf.fh.trcov  = @gpcf_constant_trcov;
    reccf.fh.trvar  = @gpcf_constant_trvar;
    reccf.fh.recappend = @gpcf_constant_recappend;
    reccf.p=[];
    reccf.p.constSigma2=[];
    if ~isempty(ri.p.constSigma2)
      reccf.p.constSigma2 = ri.p.constSigma2;
    end

  else
    % Append to the record
    gpp = gpcf.p;

    % record constSigma2
    reccf.constSigma2(ri,:)=gpcf.constSigma2;
    if isfield(gpp,'constSigma2') && ~isempty(gpp.constSigma2)
      reccf.p.constSigma2 = gpp.constSigma2.fh.recappend(reccf.p.constSigma2, ri, gpcf.p.constSigma2);
    end
  end
end

function [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = gpcf_constant_cf2ss(gpcf,x)
%GPCF_CONSTANT_CF2SS Convert the covariance function to state space form
%
%  Description
%    Convert the covariance function to state space form such that
%    the process can be described by the stochastic differential equation
%    of the form:
%      df(t)/dt = F f(t) + L w(t),
%    where w(t) is a white noise process. The observation model now
%    corresponds to y_k = H f(t_k) + r_k, where r_k ~ N(0,sigma2).
%
%

  % Check inputs
  if nargin<2, x=[]; end

  % Define the model
  F      = 0;
  L      = 1;
  Qc     = 0;
  H      = 1;
  Pinf   = gpcf.constSigma2;
  dF     = 0;
  dQc    = 0;
  dPinf  = 1;

  % Set params
  params.stationary = true;

  % Check which parameters are optimized
  if isempty(gpcf.p.constSigma2), ind(1) = false; else ind(1) = true; end

  % Return only those derivatives that are needed
  dF    = dF(:,:,ind);
  dQc   = dQc(:,:,ind);
  dPinf = dPinf(:,:,ind);

end

author: gpstuff-dev

repo: gpstuff

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


text:
function status = isequalfigure(a, b)

% ISEQUALFIGURE compares two figures and returns a boolean
%
% Use as
%  isequalfigure(a, b)
% with a and b being the figure handles.

fna = [tempname, '.png'];
fnb = [tempname, '.png'];

print(a, '-dpng', fna);
print(b, '-dpng', fnb);

mata = imread(fna);
matb = imread(fnb);

delete(fna);
delete(fnb);

status = isequal(mata, matb);

% the following is too sensitive
% md5a = ft_hash(fna);
% md5b = ft_hash(fnb);
% status = isequal(md5a, md5b);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [lp,lj,tr] = lineline(pa,pb,pc,pd,varargin)
%LINELINE intersection between lines in d-dimensional space.
%   [LP,LI] = LINELINE(PA,PB,PC,PD) finds intersections bet-
%   ween line segments in d-dimensions. Lines are specified
%   as a set of endpoints [PA,PB] and [PC,PD] where PA, PB,
%   PC and PD are NL-by-ND arrays of coordinates, where ND
%   is the number of dimensions.
%
%   A set of intersecting lines from [PA,PB] is returned for
%   each query line in [PC,PB], such that the II-th query
%   line is associated with the lines LI(LP(II,1):LP(II,2)).
%   Lines without intersections have LP(II,1) == 0.
%
%   [LP,LI,TR] = LINELINE(PA,PB,PC,PD) additionally returns
%   the supporting aabb-tree used internally to compute the
%   query. If the underlying collection [PA,PB] is static,
%   the  tree TR may be recycled for subsequent calls, using
%   [LP,LI,TR] = FINDLINE(PA,PB,PC,PD,TR). This syntax may
%   lead to improved performance, especially when the number
%   of lines is large w.r.t. the number of query lines. Note
%   that in such cases the distribution of underlying lines
%   is NOT permitted to change between calls, or erroneous
%   results may be returned. Additional parameters used to
%   govern the creation of the underlying aabb-tree may be
%   passed via [...] = LINELINE(...,TR,OP). See MAKETREE for
%   additional information.
%
%   See also MAKETREE, FINDTRIA, FINDBALL, FINDLINE

% Please see the following for additional information:
%
%   Darren Engwirda, "Locally-optimal Delaunay-refinement &
%   optimisation-based mesh generation". Ph.D. Thesis, Scho-
%   ol of Mathematics and Statistics, Univ. of Sydney, 2014:
%   http://hdl.handle.net/2123/13148

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 10/10/2017
%-----------------------------------------------------------

    lp = []; lj = []; tr = []; op = [];

%---------------------------------------------- basic checks
    if (nargin < +4 || nargin > +6)
        error('lineline:incorrectNumInputs', ...
            'Incorrect number of inputs.');
    end

    if (nargin >= +5), tr = varargin{1}; end
    if (nargin >= +6), op = varargin{2}; end

%------------------------------ quick return on empty inputs
    if (isempty(pa)), return ; end
    if (isempty(pb)), return ; end
    if (isempty(pc)), return ; end
    if (isempty(pd)), return ; end

%---------------------------------------------- basic checks
    if (~isnumeric(pa) || ~isnumeric(pb) || ...
        ~isnumeric(pc) || ~isnumeric(pd) )
        error('lineline:incorrectInputClass', ...
            'Incorrect input class.') ;
    end

    if (ndims(pa) ~= +2 || size(pa,2) < +2 || ...
        ndims(pb) ~= +2 || size(pb,2) < +2 || ...
        ndims(pc) ~= +2 || size(pc,2) < +2 || ...
        ndims(pd) ~= +2 || size(pd,2) < +2 )
        error('lineline:incorrectDimensions', ...
            'Incorrect input dimensions.');
    end
    if (size(pa,1) ~= size(pb,1) || ...
        size(pc,1) ~= size(pd,1) )
        error('lineline:incorrectDimensions', ...
            'Incorrect input dimensions.');
    end

    if (~isempty(tr) && ~isstruct(tr) )
        error('lineline:incorrectInputClass', ...
            'Incorrect input class.') ;
    end
    if (~isempty(op) && ~isstruct(op) )
        error('lineline:incorrectInputClass', ...
            'Incorrect input class.') ;
    end

    nd = size(pa,2);
    nl = size(pa,1);
    ml = size(pc,1);

    if (isempty(tr))
%------------------------------ compute aabb-tree for d-line
    ab = zeros(nl,nd*+2) ;
    for ax = +1:nd            % compute aabb's
    ab(:,ax+nd*0) = min(pa(:,ax), ...
                        pb(:,ax)) ;
    ab(:,ax+nd*1) = max(pa(:,ax), ...
                        pb(:,ax)) ;
    end
    tr = maketree(ab,op) ;

    end

%------------------------------ compute tree-to-vert mapping
    ab = zeros(ml,nd*+2) ;
    for ax = +1:nd            % compute aabb's
    ab(:,ax+nd*0) = min(pc(:,ax), ...
                        pd(:,ax)) ;
    ab(:,ax+nd*1) = max(pc(:,ax), ...
                        pd(:,ax)) ;
    end
    tm = maprect (tr,ab) ;

%------------------------------ compute line-to-line queries
   [li,ip,lj] = queryset( ...
    tr,tm,@linekern,pc,pd,pa,pb) ;

%------------------------------ re-index onto full obj. list
    lp = zeros(size(pc,1),2) ;
    lp( :,2) = -1 ;

    if (isempty(li)), return ; end

    lp(li,:) = ip ;

end

function [i1,i2] = linekern(l1,l2,pa,pb,pc,pd)
%LINEKERN d-dim. line//line intersection kernel routine.

        m1 = length(l1) ;
        m2 = length(l2) ;

    %-------------------------- push line/vert onto n*m tile
        l1 = l1.' ;

        l1 = l1(ones(m2,1),:) ;
        l1 = l1(:);
        l2 = l2(:,ones(1,m1)) ;
        l2 = l2(:);

    %-------------------------- compute O(n*m) intersections
       [ok,tp,tq] = linenear( ...
           pa(l1,:),pb(l1,:), ...
           pc(l2,:),pd(l2,:)) ;

        rt = +1.+eps;

        ix = abs(tp) <= +rt & ...
             abs(tq) <= +rt ;

        ix(~ok) = false ;

        i1 = l1(ix) ;
        i2 = l2(ix) ;

end

author: dengwirda


text:
function newCell=reformatElements(inputCell,type,delimiter)
%reformatElements  reformat elements of cell array to desired format
% reformatElements
%   convert cell array element format between string and nested cell
%
% Input:
%   inputCell    the input cell array
%   type         two conversion approaches: cell2str and str2cell
%                (default: str2cell)
%   delimiter    specify the delimiter separating values within the
%                element (default: semicolon)
%
% Output:
%   newCell      the output of cell array with refromatted elements
%
% Usage: newCell=reformatElements(inputCell,type,delimiter)
%

newCell={};

if ~iscell(inputCell)
		EM='Wrong input arguments';
		disp(EM);
end

if nargin < 2
    type = 'str2cell';
end

if nargin < 3
    delimiter = ';';
end

% get the index of non-empty elements
index=find(~cellfun(@isempty,inputCell));

if length(index)>0
		if (iscell(inputCell{index(1)}) && strcmpi(type,'str2cell')) || (~iscell(inputCell{index(1)})  &&  strcmpi(type,'cell2str'))
		fprintf('The conversion parameter is conflic with input data!\n');
		return;
elseif length(index)==0
		fprintf('The input cell is empty!\n');
		return;
end

% initilize output
newCell=cell(numel(inputCell),1);
newCell(:)={''};

if strcmp('str2cell',type)
    % from string to nested cell
    inputCell = regexprep(inputCell, '\s', '');   % remove space from input
	newCell(index)=cellfun(@(s) strsplit(s,delimiter),inputCell(index),'UniformOutput', false);
elseif isequal('cell2str',type)
    % combine elements of each cell to string
    delimiter = [delimiter(~isspace(delimiter)) ' '];  % append one space only
    newCell(index)=cellfun(@(s) strjoin(s,delimiter),inputCell(index),'UniformOutput', false);
end

end

author: SysBioChalmers

repo: Human-GEM


text:
function coeffs = decdemo( im, option )
% DECDEMO  demonstrates nonsubsampled Contourlet decomposition and reconstruction.
%
%   DECDEMO shows how to use the nonsubsampled contourlet toolbox to decompose
%   and reconstruct an image.  It provides a sample script that uses
%   basic functions such as nsctdec, nsctrec, and shownsct.
%
%   It can be modified for applications such as image analysis,
%   image retrieval and image processing.
%
% Input:
%	image:  a double or integer matrix for the input image.
%           The default is the zoneplate image.
%   option: option for the demos. The default value is 'auto'
%       'auto' ------  automtatical demo, no input
%       'user' ------  semi-automatic demo, simple interactive inputs
%       'expert' ----  mannual, complete interactive inputs.
%                      (Not implmented in this version)
%
% Output:
%	coeffs: a cell vector for the contourlet decomposition coefficients.
%
% See also:     NSCTDEC, NSCTREC, SHOWNSCT.

% History:
%   08/08/2004  Created by Jianping Zhou.

disp('Welcome to the nonsubsampled Contourlet decomposition demo! :)');
disp('Type help decdemo for help' ) ;
disp('You can also view decdemo.m for details.') ;
disp(' ');

% Input image
if ~exist('im', 'var')
    % Zoneplate image: good for illustrating multiscale and directional
    % decomposition
    im = imread ('zoneplate.png') ;
elseif isstr(im)
    im = imread ( im ) ;
else
    error('You shall input valid image name!');
end

% Show the input image
disp( 'Displaying the input image...');
clf;
imagesc(im, [0, 255]);
title( 'Input image' ) ;
axis image off;
colormap(gray);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Image decomposition by nonsubsampled contourlet transform (NSCT).
% This is the iterated filter bank that computes the nonsubsampled
% contourlet transform.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameteters:
nlevels = [0, 1, 3] ;        % Decomposition level
pfilter = 'maxflat' ;              % Pyramidal filter
dfilter = 'dmaxflat7' ;              % Directional filter

% Nonsubsampled Contourlet decomposition
coeffs = nsctdec( double(im), nlevels, dfilter, pfilter );

% Display the coefficients
disp('Displaying the contourlet coefficients...') ;
shownsct( coeffs ) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonsubsampled Contourlet transform (NSCT) reconstruction.
% This is the inverse of nsctdec, i.e.
% imrec = nsctrec(coeffs, dfilter, pfilter);
% would reconstruct imrec = im
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Reconstruct image
imrec = nsctrec( coeffs, dfilter, pfilter ) ;

disp('Displaying the reconstructed image...') ;
disp('It should be a perfect reconstruction' ) ;
disp(' ') ;

% Show the reconstruction image and the original image
figure;
subplot(1,2,1), imagesc( im, [0, 255] );
title('Original image' ) ;
colormap(gray);
axis image off;
subplot(1,2,2), imagesc( imrec, [0, 255] );
title('Reconstructed image' ) ;
colormap(gray);
axis image off;

mse = sum( sum( (imrec - double(im)).^2 ) );
mse = mse / prod(size(im));

disp( sprintf('The mean square error is: %f', mse ) );
disp(' ');

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [locs pks]=peakseek(x,minpeakdist,minpeakh)
% Alternative to the findpeaks function.  This thing runs much much faster.
% It really leaves findpeaks in the dust.  It also can handle ties between
% peaks.  Findpeaks just erases both in a tie.  Shame on findpeaks.
%
% x is a vector input (generally a timecourse)
% minpeakdist is the minimum desired distance between peaks (optional, defaults to 1)
% minpeakh is the minimum height of a peak (optional)
%
% (c) 2010
% Peter O'Connor
% peter<dot>ed<dot>oconnor .AT. gmail<dot>com

if size(x,2)==1, x=x'; end

% Find all maxima and ties
locs=find(x(2:end-1)>=x(1:end-2) & x(2:end-1)>=x(3:end))+1;

if nargin<2, minpeakdist=1; end % If no minpeakdist specified, default to 1.

if nargin>2 % If there's a minpeakheight
    locs(x(locs)<=minpeakh)=[];
end

if minpeakdist>1
    while 1

        del=diff(locs)<minpeakdist;

        if ~any(del), break; end

        pks=x(locs);

        [garb mins]=min([pks(del) ; pks([false del])]); %#ok<ASGLU>

        deln=find(del);

        deln=[deln(mins==1) deln(mins==2)+1];

        locs(deln)=[];

    end
end

if nargout>1,
    pks=x(locs);
end

end

author: Sable

repo: mcbench-benchmarks


text:
function stroud_test335 ( )

%*****************************************************************************80
%
%% TEST335 tests SPHERE_SHELL_03_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_ND_INDEX;

  n_max = 3;

  num = function_nd_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST335\n' );
  fprintf ( 1, '  For integrals inside a spherical shell in ND:\n' );
  fprintf ( 1, '  SPHERE_SHELL_03_ND approximates the integral.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We compare these results with those computed by\n' );
  fprintf ( 1, '  from the difference of two ball integrals:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  BALL_F1_ND approximates the integral;\n' );
  fprintf ( 1, '  BALL_F3_ND approximates the integral.\n' );
  fprintf ( 1, '\n' );

  for j = 1 : 2

    if ( j == 1 )
      r1 = 0.0;
      r2 = 1.0;
      xc(1:n_max) = 0.0;
    else
      r1 = 2.0;
      r2 = 3.0;
      xc(1:n_max) = [ 1.0, -1.0, 2.0 ];
    end

    for n = 2 : n_max

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Spatial dimension N = %d\n', n );
      fprintf ( 1, '  Sphere center:\n' );
      for i = 1 : n
        fprintf ( 1, '  %12f', xc(i) );
      end
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Inner sphere radius = %f\n', r1 );
      fprintf ( 1, '  Outer sphere radius = %f\n', r2 );
      fprintf ( 1, '  Spherical shell volume = %f\n', ...
        sphere_shell_volume_nd ( n, r1, r2 ) );
      fprintf ( 1, '\n' );
      fprintf ( 1, '\n' );
      fprintf ( 1, '    Rule:      #3       F1(R2)-F1(R1)  F3(R2)-F3(R1)\n' );
      fprintf ( 1, '    F(X)\n' );
      fprintf ( 1, '\n' );

      for i = 1 : num

        FUNC_ND_INDEX = i;

        result1 = sphere_shell_03_nd ( 'function_nd', n, xc, r1, r2 );

        result3 = ball_f1_nd ( 'function_nd', n, xc, r1 );
        result4 = ball_f1_nd ( 'function_nd', n, xc, r2 );

        result5 = ball_f3_nd ( 'function_nd', n, xc, r1 );
        result6 = ball_f3_nd ( 'function_nd', n, xc, r2 );

        fname = function_nd_name ( i );

        fprintf ( 1, '  %s  %12f  %12f  %12f\n', ...
          fname, result1, result4-result3, result6-result5 );

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function str = cstring(x)
%
%CSTRING Print to a string like c
%  CSTRING(x) takes an array x, of type char, and prints it to
%  a string the way C would.
%  the first null, or non-ascii character is the terminator

i=1;
while (x(i)>0 & x(i)<128 & i<length(x))
  i=i+1;
end

if i==1
  str = char(0);
else
  str = char(x(1:i-1));
end

str = sprintf(str,'%s');

return


text:
function MS = createMesh2D(varargin)
% MeshStructure = createMesh2D(Nx, Ny, Lx, Ly)
% MeshStructure = createMesh2D(facelocationX, facelocationY)
% creates a uniform 2D mesh:
% Nx is the number of cells in x (horizontal) direction
% Ny is the number of cells in y (vertical) direction
% Lx is the domain length in x direction
% Ly is the domain length in y direction
%
% SYNOPSIS:
%   MeshStructure = createMesh2D(Nx, Ny, Lx, Ly)
%
% PARAMETERS:
%   Nx: number of cells in the x direction
%   Lx: domain length in x direction
%   Ny: number of cells in the y direction
%   Ly: domain length in y direction
%
% RETURNS:
%   MeshStructure.
%                 dimensions=2 (2D problem)
%                 numbering: shows the indexes of cellsn from left to right
%                 and top to bottom
%                 cellsize: x and y elements of the cell size =[Lx/Nx,
%                 Ly/Ny]
%                 cellcenters.x: location of each cell in the x direction
%                 cellcenters.y: location of each cell in the y direction
%                 facecenters.x: location of interface between cells in the
%                 x direction
%                 facecenters.y: location of interface between cells in the
%                 y direction
%                 numberofcells: [Nx, Ny]
%
%
% EXAMPLE:
%   Nx = 5;
%   Ny = 7;
%   Lx = 10;
%   Ly = 20;
%   m = createMesh2D(Nx, Ny, Lx, Ly);
%   [X, Y] = ndgrid(m.cellcenters.x, m.cellcenters.y);
%   [Xf,Yf]=ndgrid(m.facecenters.x, m.facecenters.y);
%   plot(X, Y, 'or', ...
%        Xf, Yf, '-b', Xf', Yf', '-b');
%
% SEE ALSO:
%     createMesh1D, createMesh3D, createMeshCylindrical1D, ...
%     createMeshCylindrical2D, createCellVariable, createFaceVariable

% Written by Ali A. Eftekhari
% See the license file

if nargin==4
  % uniform 1D mesh
  Nx=varargin{1};
  Ny=varargin{2};
  Width=varargin{3};
  Height=varargin{4};
  % cell size is dx
  dx = Width/Nx;
  dy = Height/Ny;
  G=reshape(1:(Nx+2)*(Ny+2), Nx+2, Ny+2);
  CellSize.x= dx*ones(Nx+2,1);
  CellSize.y= dy*ones(Ny+2,1);
  CellSize.z= [0.0];
  CellLocation.x= [1:Nx]'*dx-dx/2;
  CellLocation.y= [1:Ny]'*dy-dy/2;
  CellLocation.z= [0.0];
  FaceLocation.x= [0:Nx]'*dx;
  FaceLocation.y= [0:Ny]'*dy;
  FaceLocation.z= [0.0];
elseif nargin==2
  % nonuniform 1D mesh
  facelocationX=varargin{1};
  facelocationY=varargin{2};
  facelocationX=facelocationX(:);
  facelocationY=facelocationY(:);
  Nx = length(facelocationX)-1;
  Ny = length(facelocationY)-1;
  G=reshape(1:(Nx+2)*(Ny+2), Nx+2, Ny+2);
  CellSize.x= [facelocationX(2)-facelocationX(1); ...
    facelocationX(2:end)-facelocationX(1:end-1); ...
    facelocationX(end)-facelocationX(end-1)];
  CellSize.y= [facelocationY(2)-facelocationY(1); ...
    facelocationY(2:end)-facelocationY(1:end-1); ...
    facelocationY(end)-facelocationY(end-1)];
  CellSize.z= [0.0];
  CellLocation.x= 0.5*(facelocationX(2:end)+facelocationX(1:end-1));
  CellLocation.y= 0.5*(facelocationY(2:end)+facelocationY(1:end-1));
  CellLocation.z= [0.0];
  FaceLocation.x= facelocationX;
  FaceLocation.y= facelocationY;
  FaceLocation.z= [0.0];
end
c=G([1,end], [1,end]);
MS=MeshStructure(2, ...
  [Nx,Ny], ...
  CellSize, ...
  CellLocation, ...
  FaceLocation, ...
  c(:), ...
  [1]);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function out=spm_api_bmc(F,N,exp_r,xp,family)
% API to select and compare DCMs using Bayesian model comparison
% FORMAT out=spm_api_bmc(F,N,alpha,exp_r,xp)
%
% INPUT:
% F      - Matrix/Vector of log model evidences
% N      - vector of model names
% alpha  - vector of model probabilities
% exp_r  - expectation of the posterior p(r|y)
% xp     - exceedance probabilities
%
% OUTPUT:
% out    - conditional probability of DCMs (when using fixed effect method)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_api_bmc.m 4832 2012-08-08 14:28:00Z will $

if nargin < 4 || isempty(xp)
    inf_method = 'FFX';
else
    inf_method = 'RFX';
end

if nargin>4
    plot_family = 1;
else
    plot_family = 0;
end

nm = length(N);

Fgraph  = spm_figure('GetWin','Graphics');
spm_clf(Fgraph);

switch inf_method

    %======================================================================
    % Fixed Effect
    %======================================================================
    case ('FFX')

        figure(Fgraph);

        %-Compute conditional probability of DCMs under flat priors.
        %------------------------------------------------------------------
        F    = F - min(F);
        i    = F < (max(F) - 32);
        P    = F;
        P(i) = max(F) - 32;
        P    = P - min(P);
        P    = exp(P);
        P    = P/sum(P);

        %-Display results
        %------------------------------------------------------------------
        subplot(2,1,1)
        bar(1:nm,F)
        set(gca,'XTick',1:nm)
        set(gca,'XTickLabel',1:nm)
        ylabel('Log-evidence (relative)','Fontsize',14)
        xlabel('Models','Fontsize',14)
        title(['Bayesian Model Selection: ',inf_method],'Fontsize',14)
        axis square
        grid on

        subplot(2,1,2)
        bar(1:nm,P)
        set(gca,'XTick',1:nm)
        set(gca,'XTickLabel',1:nm)
        ylabel('Model Posterior Probability','Fontsize',14)
        title(['Bayesian Model Selection: ',inf_method],'Fontsize',14)
        xlabel('Models','Fontsize',14)
        axis square
        grid on

        out = P;

        if plot_family

            %-Display results - families
            %--------------------------------------------------------------
            F  = spm_figure('Create','Graphics','BMS: results');
            figure(F);

            Nfam = length(family.post);
            bar(1:Nfam,family.post)
            set(gca,'XTick',1:Nfam)
            set(gca,'XTickLabel',family.names)
            ylabel('Family Posterior Probability','Fontsize',14)
            xlabel('Families','Fontsize',14)
            title(['Bayesian Model Selection: ',inf_method],'Fontsize',14)
            axis square
            grid on

            out = [];

        end

    %======================================================================
    % Random Effect
    %======================================================================
    case ('RFX')

        figure(Fgraph);

        %-Display results
        %------------------------------------------------------------------
        subplot(2,1,1)
        bar(1:length(N),exp_r)
        set(gca,'XTick',1:length(N))
        set(gca,'XTickLabel',1:nm)
        ylabel('Model Expected Probability','Fontsize',14)
        xlabel('Models','Fontsize',14)
        title(['Bayesian Model Selection: ',inf_method],'Fontsize',14)
        axis square
        grid on

        subplot(2,1,2)
        bar(1:length(N),xp')
        set(gca,'XTick',1:length(N))
        set(gca,'XTickLabel',1:nm)
        ylabel('Model Exceedance Probability','Fontsize',14)
        xlabel('Models','Fontsize',14)
        title(['Bayesian Model Selection: ',inf_method],'Fontsize',14)
        axis square
        grid on

        out = [];

        if plot_family

            %-Display results - families
            %--------------------------------------------------------------
            F  = spm_figure('Create','Graphics','BMS: results');
            figure(F);

            %-Display results - families
            %--------------------------------------------------------------
            subplot(2,1,1)
            Nfam = length(family.exp_r);
            bar(1:Nfam,family.exp_r)
            set(gca,'XTick',1:Nfam)
            set(gca,'XTickLabel',family.names)
            ylabel('Family Expected Probability','Fontsize',14)
            xlabel('Families','Fontsize',14)
            title(['Bayesian Model Selection: ',inf_method],'Fontsize',14)
            axis square
            grid on

            subplot(2,1,2)
            bar(1:Nfam,family.xp')
            set(gca,'XTick',1:Nfam)
            set(gca,'XTickLabel',family.names)
            ylabel('Family Exceedance Probability','Fontsize',14)
            xlabel('Families','Fontsize',14)
            title(['Bayesian Model Selection: ',inf_method],'Fontsize',14)
            axis square
            grid on

            out = [];

        end
end

author: spm


text:
function RecreateTspFile(fbasename)

tspdata=load([fbasename '.tsp']);
tspdata = tspdata(:,1);

if size(tspdata,2)>1
  warning('tsp has data in it! Are You sure you want to continue?')
%   keyboard
end

nFramesTsp = size(tspdata,1);
tspnew = [tspdata -1*ones(size(tspdata,1),6)];

%readerobj  = VideoReader([fbasename '.mpg']);

if ~exist([fbasename '.ogg'],'file') && exist([fbasename '.mpg'],'file')
    %eval(['!ffmpeg -i ' fbasename '.mpg -f ogg ' fbasename '.ogg'])
    eval(['!avconv -i ' fbasename '.mpg -f ogg ' fbasename '.ogg'])
elseif ~exist([fbasename '.ogg'],'file') && ~exist([fbasename '.mpg'],'file')
    error('No video file!!')
end

whl = ApproxMedianFilter_RB_LED([fbasename '.ogg']);
figure(1),clf
plot(whl(:,1),whl(:,2))
hold on
plot(whl(:,3),whl(:,4),'r')

nDiff = nFramesTsp - size(whl,1);
if nDiff ~=0 %in case of frame drop, interpolate the output

    twhl = (1:size(whl,1));
    ttsp = (1:size(tspdata,1));
    whl(whl==-1) = NaN;
    interpData = interp1(twhl,whl,ttsp);
    interpData(isnan(interpData)) = -1;
    tspnew(:,[2 3 6 7]) = interpData;
else
    tspnew(:,[2 3 6 7]) = whl;
end

fid = fopen([fbasename '.tsp'],'w');
fprintf(fid,'%i\t %f\t %f\t %f\t %f\t %f\t %f\t \n',tspnew');
fclose(fid);
%eval(['!rm ' fbasename '.ogg'])

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function lynx ( )

%*****************************************************************************80
%
%% LYNX uses MATLAB to draw a curve with datapoints marked.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LYNX:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Make a curve by connecting a series of data points.\n' );
  fprintf ( 1, '  Mark the data points.\n' );

  xy = load ( 'lynx.txt' );
%
%  Plot the data.
%
  plot ( xy(:,1), xy(:,2), 'LineWidth', 2, 'Color', 'g' );
  hold on
  plot ( xy(:,1), xy(:,2), 'ko' );
  grid on
  xlabel ( '<--- Year --->' );
  ylabel ( '<--- Lynx Harvest --->' );
  title ( 'Lynx Population Records' );
  hold off

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LYNX:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
%CONVECTIONDIFFUSION Convection and diffusion physics mode definition.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function reinit_params(net, varargin)
% REINIT_PARAMS re-initializes the parameters of a DAG.

opts.scale = 1;
opts.weightInitMethod = 'xavierimproved';
opts.initBias = 0.1;
opts = vl_argparse(opts, varargin);

for l = 1:numel(net.layers)
    if isempty(net.layers(l).params)
        continue;
    end
    if isa(net.layers(l).block, 'dagnn.Conv')
        sz = net.layers(l).block.size;
        filters_param = net.layers(l).params{1};
        net.params(net.getParamIndex(filters_param)).value = ...
            init_weight(opts, sz(1), sz(2), sz(3), sz(4), 'single');
        % Optional bias.
        if numel(net.layers(l).params) > 1
            bias_param = net.layers(l).params{2};
            net.params(net.getParamIndex(bias_param)).value = ...
                ones(sz(4), 1, 'single') * opts.initBias;
        end
    elseif isa(net.layers(l).block, 'dagnn.BatchNorm')
        gain_param    = net.layers(l).params{1};
        bias_param    = net.layers(l).params{2};
        moments_param = net.layers(l).params{3};
        out = net.layers(l).block.numChannels;
        if isempty(out) || out < 1
            out = numel(net.params(net.getParamIndex(gain_param)).value);
            if out == 0
                error('output size unknown');
            end
        end
        net.params(net.getParamIndex(gain_param)).value = ...
            ones(out, 1, 'single');
        net.params(net.getParamIndex(bias_param)).value = ...
            zeros(out, 1, 'single');
        net.params(net.getParamIndex(moments_param)).value = ...
            zeros(out, 2, 'single');
    else
        % Every layer with params should be reinitialized.
        error(['unknown layer with params: ' class(net.layers(l).block)]);
    end
end

end

function weights = init_weight(opts, h, w, in, out, type)
% From examples/imagenet/cnn_imagenet_init.m

    % See K. He, X. Zhang, S. Ren, and J. Sun. Delving deep into
    % rectifiers: Surpassing human-level performance on imagenet
    % classification. CoRR, (arXiv:1502.01852v1), 2015.
    switch lower(opts.weightInitMethod)
      case 'gaussian'
        sc = 0.01/opts.scale ;
        weights = randn(h, w, in, out, type)*sc;
      case 'xavier'
        sc = sqrt(3/(h*w*in)) ;
        weights = (rand(h, w, in, out, type)*2 - 1)*sc ;
      case 'xavierimproved'
        sc = sqrt(2/(h*w*out)) ;
        weights = randn(h, w, in, out, type)*sc ;
      otherwise
        error('Unknown weight initialization method''%s''', opts.weightInitMethod) ;
    end
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function mask_normalize()
    global config mem;
    curr_layer_idx = config.misc.current_layer - 1;
    mem.activations{curr_layer_idx} = mem.activations{curr_layer_idx} ./ mem.mask_activations{curr_layer_idx};
end

author: jimmy-ren


text:
function [u, erriter, i, timet] = CMF_Cut
%
%   Performing the continuous max-flow algorithm to solve the
%   continuous min-cut problem in 2D
%
%   Usage: [u, erriter, i, timet] = CMF_Cut;
%
%   Inputs: there is no input since all data and parameters can be
%           adjusted within the program
%
%   Outputs:
%       - u: the final results u(x) in [0,1]. As the following paper,
%           the global binary result can be available by threshholding u
%           by any constant alpha in (0,1):
%
%           Nikolova, M.; Esedoglu, S.; Chan, T. F.
%           Algorithms for Finding Global Minimizers of Image Segmentation and Denoising Models
%           SIAM J. App. Math., 2006, 66, 1632-1648
%
%       - erriter: it returns the error evaluation of each iteration,
%           i.e. it shows the convergence rate. One can check the algorithm
%           performance.
%
%       - i: gives the total number of iterations, when the algorithm converges.
%
%       - timet: gives the total computation time.
%
%   Example:
%       >> [u, erriter, i, timet] = CMF_Cut;
%
%       >> us = max(u, beta);  % where beta in (0,1)
%
%       >> imagesc(us), colormap gray, axis image, axis off;figure(gcf)
%
%       >> figure, loglog(erriter,'DisplayName','erriterN');figure(gcf)
%
%
%
%   The original algorithm was proposed in the following papers:
%
%   [1] Yuan, J.; Bae, E.;  Tai, X.-C.
%       A Study on Continuous Max-Flow and Min-Cut Approaches
%       CVPR, 2010
%
%   [2] Yuan, J.; Bae, E.; Tai, X.-C.; Boycov, Y.
%       A study on continuous max-flow and min-cut approaches. Part I: Binary labeling
%       UCLA CAM, Technical Report 10-61, 2010
%
%   The mimetic finite-difference discretization method was proposed for
%   the total-variation function in the paper:
%
%   [1] Yuan, J.; Schn{\"o}rr, C.; Steidl, G.
%       Simultaneous Optical Flow Estimation and Decomposition
%       SIAM J.~Scientific Computing, 2007, vol. 29, page 2283-2304, number 6
%
%   This software can be used only for research purposes, you should cite ALL of
%   the aforementioned papers in any resulting publication.
%
%   Please email cn.yuanjing@gmail.com for any questions, suggestions and bug reports
%
%   The Software is provided "as is", without warranty of any kind.
%
%
%                       Version 1.0
%           https://sites.google.com/site/wwwjingyuan/
%
%           Copyright 2011 Jing Yuan (cn.yuanjing@gmail.com)
%

ur = double(imread('cameraman.jpg'))/255;

[rows, cols] = size(ur);
imgSize = rows*cols;

% define the required parameters:
%
%   - alpha: the penalty parameter to the total-variation term.
%       For the case without incorporating image-edge weights, alpha is given
%       by the constant everywhere. For the case with image-edge weights,
%       alpha is given by the pixelwise weight function:
%
%       For example, alpha(x) = b/(1 + a*| nabla f(x)|) where b and a are positive
%       constants and |nabla f(x)| gives the strength of the local gradient.
%
%   - cc: gives the step-size of the augmented Lagrangian method.
%       The optimal range of cc is [0.3, 3].
%
%   - errbound: the error bound for convergence.
%
%   - numIter: the maximum iteration number.
%
%   - steps: the step-size for the graident-projection step to the
%       total-variation function. The optimal range of steps is [0.1,
%       0.17].
%

alpha = 0.5*ones(rows,cols);
cc = 0.3;
errbound = 1e-4;
numIter = 300;
steps = 0.16;

% build up the data terms
ulab(1) = 0.15;
ulab(2) = 0.6;
Cs = abs(ur - ulab(1));
Ct = abs(ur - ulab(2));

% set the initial values
%   - the initial value of u is set to be an initial cut, see below.
%   - the initial values of two terminal flows ps and pt are set to be the
%     specified legal flows.
%   - the initial value of the spatial flow fiels p = (pp1, pp2) is set to
%   be zero.

u = double((Cs-Ct) >= 0);
ps = min(Cs, Ct);
pt = ps;

pp1 = zeros(rows, cols+1);
pp2 = zeros(rows+1, cols);
divp = pp1(:,2:cols+1)-pp1(:,1:cols)+pp2(2:rows+1,:)-pp2(1:rows,:);

erriter = zeros(numIter,1);

tic
for i = 1:numIter

	% update the spatial flow field p = (pp1, pp2):
    %   the following steps are the gradient descent step with steps as the
    %   step-size.

    pts = divp - (ps - pt  + u/cc);
    pp1(:,2:cols) = pp1(:,2:cols) + steps*(pts(:,2:cols) - pts(:,1:cols-1));
    pp2(2:rows,:) = pp2(2:rows,:) + steps*(pts(2:rows,:) - pts(1:rows-1,:));

    % the following steps give the projection to make |p(x)| <= alpha(x)

    gk = sqrt((pp1(:,1:cols).^2 + pp1(:,2:cols+1).^2 + pp2(1:rows,:).^2 + pp2(2:rows+1,:).^2)*0.5);
    gk = double(gk <= alpha) + double(~(gk <= alpha)).*(gk ./ alpha);
    gk = 1 ./ gk;

    pp1(:,2:cols) = (0.5*(gk(:,2:cols) + gk(:,1:cols-1))).*pp1(:,2:cols);
    pp2(2:rows,:) = (0.5*(gk(2:rows,:) + gk(1:rows-1,:))).*pp2(2:rows,:);

    divp = pp1(:,2:cols+1)-pp1(:,1:cols)+pp2(2:rows+1,:)-pp2(1:rows,:);

    % updata the source flow ps

    pts = divp + pt - u/cc + 1/cc;
    ps = min(pts, Cs);

    % update the sink flow pt

    pts = - divp + ps + u/cc;
    pt = min(pts, Ct);

	% update the multiplier u

	erru = cc*(divp + pt  - ps);
	u = u - erru;

    % evaluate the avarage error

    erriter(i) = sum(sum(abs(erru)))/imgSize;

    if (erriter(i) < errbound)
        break;
    end
end
toc
timet = toc

msg = sprintf('number of iterations = %u. \n', i);
disp(msg);


text:
function [U,S,V] = pca(A,k,its,l)
%PCA  Low-rank approximation in SVD form.
%
%
%   [U,S,V] = PCA(A)  constructs a nearly optimal rank-6 approximation
%             USV' to A, using 2 full iterations of a block Lanczos method
%             of block size 6+2=8, started with an n x 8 random matrix,
%             when A is m x n; the ref. below explains "nearly optimal."
%             The smallest dimension of A must be >= 6 when A is
%             the only input to PCA.
%
%   [U,S,V] = PCA(A,k)  constructs a nearly optimal rank-k approximation
%             USV' to A, using 2 full iterations of a block Lanczos method
%             of block size k+2, started with an n x (k+2) random matrix,
%             when A is m x n; the ref. below explains "nearly optimal."
%             k must be a positive integer <= the smallest dimension of A.
%
%   [U,S,V] = PCA(A,k,its)  constructs a nearly optimal rank-k approx. USV'
%             to A, using its full iterations of a block Lanczos method
%             of block size k+2, started with an n x (k+2) random matrix,
%             when A is m x n; the ref. below explains "nearly optimal."
%             k must be a positive integer <= the smallest dimension of A,
%             and its must be a nonnegative integer.
%
%   [U,S,V] = PCA(A,k,its,l)  constructs a nearly optimal rank-k approx.
%             USV' to A, using its full iterates of a block Lanczos method
%             of block size l, started with an n x l random matrix,
%             when A is m x n; the ref. below explains "nearly optimal."
%             k must be a positive integer <= the smallest dimension of A,
%             its must be a nonnegative integer,
%             and l must be a positive integer >= k.
%
%
%   The low-rank approximation USV' is in the form of an SVD in the sense
%   that the columns of U are orthonormal, as are the columns of V,
%   the entries of S are all nonnegative, and the only nonzero entries
%   of S appear in non-increasing order on its diagonal.
%   U is m x k, V is n x k, and S is k x k, when A is m x n.
%
%   Increasing its or l improves the accuracy of the approximation USV'
%   to A; the ref. below describes how the accuracy depends on its and l.
%
%
%   Note: PCA invokes RAND. To obtain repeatable results,
%         invoke RAND('seed',j) with a fixed integer j before invoking PCA.
%
%   Note: PCA currently requires the user to center and normalize the rows
%         or columns of the input matrix A before invoking PCA (if such
%         is desired).
%
%   Note: The user may ascertain the accuracy of the approximation USV'
%         to A by invoking DIFFSNORM(A,U,S,V).
%
%
%   inputs (the first is required):
%   A -- matrix being approximated
%   k -- rank of the approximation being constructed;
%        k must be a positive integer <= the smallest dimension of A,
%        and defaults to 6
%   its -- number of full iterations of a block Lanczos method to conduct;
%          its must be a nonnegative integer, and defaults to 2
%   l -- block size of the block Lanczos iterations;
%        l must be a positive integer >= k, and defaults to k+2
%
%   outputs (all three are required):
%   U -- m x k matrix in the rank-k approximation USV' to A,
%        where A is m x n; the columns of U are orthonormal
%   S -- k x k matrix in the rank-k approximation USV' to A,
%        where A is m x n; the entries of S are all nonnegative,
%        and its only nonzero entries appear in nonincreasing order
%        on the diagonal
%   V -- n x k matrix in the rank-k approximation USV' to A,
%        where A is m x n; the columns of V are orthonormal
%
%
%   Example:
%     A = rand(1000,2)*rand(2,1000);
%     A = A/normest(A);
%     [U,S,V] = pca(A,2,0);
%     diffsnorm(A,U,S,V)
%
%     This code snippet produces a rank-2 approximation USV' to A such that
%     the columns of U are orthonormal, as are the columns of V, and
%     the entries of S are all nonnegative and are zero off the diagonal.
%     diffsnorm(A,U,S,V) outputs an estimate of the spectral norm
%     of A-USV', which should be close to the machine precision.
%
%
%   Reference:
%   Nathan Halko, Per-Gunnar Martinsson, and Joel Tropp,
%   Finding structure with randomness: Stochastic algorithms
%   for constructing approximate matrix decompositions,
%   arXiv:0909.4061 [math.NA; math.PR], 2009
%   (available at http://arxiv.org).
%
%
%   See also PCACOV, PRINCOMP, SVDS.
%

%   Copyright 2009 Mark Tygert.

%
% Check the number of inputs.
%
if(nargin < 1)
  error('MATLAB:pca:TooFewIn',...
        'There must be at least 1 input.')
end

if(nargin > 4)
  error('MATLAB:pca:TooManyIn',...
        'There must be at most 4 inputs.')
end

%
% Check the number of outputs.
%
if(nargout ~= 3)
  error('MATLAB:pca:WrongNumOut',...
        'There must be exactly 3 outputs.')
end

%
% Set the inputs k, its, and l to default values, if necessary.
%
if(nargin == 1)
  k = 6;
  its = 2;
  l = k+2;
end

if(nargin == 2)
  its = 2;
  l = k+2;
end

if(nargin == 3)
  l = k+2;
end

%
% Check the first input argument.
%
if(~isfloat(A))
  error('MATLAB:pca:In1NotFloat',...
        'Input 1 must be a floating-point matrix.')
end

if(isempty(A))
  error('MATLAB:pca:In1Empty',...
        'Input 1 must not be empty.')
end

%
% Retrieve the dimensions of A.
%
[m n] = size(A);

%
% Check the remaining input arguments.
%
if(size(k,1) ~= 1 || size(k,2) ~= 1)
  error('MATLAB:pca:In2Not1x1',...
        'Input 2 must be a scalar.')
end

if(size(its,1) ~= 1 || size(its,2) ~= 1)
  error('MATLAB:pca:In3Not1x1',...
        'Input 3 must be a scalar.')
end

if(size(l,1) ~= 1 || size(l,2) ~= 1)
  error('MATLAB:pca:In4Not1x1',...
        'Input 4 must be a scalar.')
end

if(k <= 0)
  error('MATLAB:pca:In2NonPos',...
        'Input 2 must be > 0.')
end

if((k > m) || (k > n))
  error('MATLAB:pca:In2TooBig',...
        'Input 2 must be <= the smallest dimension of Input 1.')
end

if(its < 0)
  error('MATLAB:pca:In3Neg',...
        'Input 3 must be >= 0.')
end

if(l < k)
  error('MATLAB:pca:In4ltIn2',...
        'Input 4 must be >= Input 2.')
end

%
% SVD A directly if (its+1)*l >= m/1.25 or (its+1)*l >= n/1.25.
%
if(((its+1)*l >= m/1.25) || ((its+1)*l >= n/1.25))

  if(~issparse(A))
    [U,S,V] = svd(A,'econ');
  end

  if(issparse(A))
    [U,S,V] = svd(full(A),'econ');
  end
%
% Retain only the leftmost k columns of U, the leftmost k columns of V,
% and the uppermost leftmost k x k block of S.
%
  U = U(:,1:k);
  V = V(:,1:k);
  S = S(1:k,1:k);

  return

end

if(m >= n)

%
% Apply A to a random matrix, obtaining H.
%

  if(isreal(A))
    H = A*(2*rand(n,l)-ones(n,l));
  end

  if(~isreal(A))
    H = A*( (2*rand(n,l)-ones(n,l)) + i*(2*rand(n,l)-ones(n,l)) );
  end

%
% Initialize F to its final size and fill its leftmost block with H.
%
  F = zeros(m,(its+1)*l);
  F(1:m, 1:l) = H;

%
% Apply A*A' to H a total of its times,
% augmenting F with the new H each time.
%
  for it = 1:its
    H = (H'*A)';
    H = A*H;
    F(1:m, (1+it*l):((it+1)*l)) = H;
  end

  clear H;

%
% Form a matrix Q whose columns constitute an orthonormal basis
% for the columns of F.
%
  [Q,R,E] = qr(F,0);

  clear F R E;

%
% SVD Q'*A to obtain approximations to the singular values
% and right singular vectors of A; adjust the left singular vectors
% of Q'*A to approximate the left singular vectors of A.
%
  [U2,S,V] = svd(Q'*A,'econ');
  U = Q*U2;

  clear Q U2;

%
% Retain only the leftmost k columns of U, the leftmost k columns of V,
% and the uppermost leftmost k x k block of S.
%
  U = U(:,1:k);
  V = V(:,1:k);
  S = S(1:k,1:k);

end

if(m < n)

%
% Apply A' to a random matrix, obtaining H.
%

  if(isreal(A))
    H = ((2*rand(l,m)-ones(l,m))*A)';
  end

  if(~isreal(A))
    H = (( (2*rand(l,m)-ones(l,m)) + i*(2*rand(l,m)-ones(l,m)) )*A)';
  end

%
% Initialize F to its final size and fill its leftmost block with H.
%
  F = zeros(n,(its+1)*l);
  F(1:n, 1:l) = H;

%
% Apply A'*A to H a total of its times,
% augmenting F with the new H each time.
%
  for it = 1:its
    H = A*H;
    H = (H'*A)';
    F(1:n, (1+it*l):((it+1)*l)) = H;
  end

  clear H;

%
% Form a matrix Q whose columns constitute an orthonormal basis
% for the columns of F.
%
  [Q,R,E] = qr(F,0);

  clear F R E;

%
% SVD A*Q to obtain approximations to the singular values
% and left singular vectors of A; adjust the right singular vectors
% of A*Q to approximate the right singular vectors of A.
%
  [U,S,V2] = svd(A*Q,'econ');
  V = Q*V2;

  clear Q V2;

%
% Retain only the leftmost k columns of U, the leftmost k columns of V,
% and the uppermost leftmost k x k block of S.
%
  U = U(:,1:k);
  V = V(:,1:k);
  S = S(1:k,1:k);

end

repo: REST


text:
function [n_covered,covering_boxes, covering_ids, covered_boxes, uncovered_boxes, object_ids]=computeNCoveredBoxes(ranked_boxes,ground_truth_boxes,n_windows, unique_solution,overlap_criterion)
n_g=size(ground_truth_boxes,1);
n_covered=zeros(length(n_windows),1);
n_windows=[0 n_windows];
object_ids=[];
nDets=size(ranked_boxes,1);

covering_ids=[];
covered_ids=[];
uncovered_ids=[];
if(unique_solution)
  n_windows(n_windows>nDets)=nDets;
  already_found = zeros(size(ground_truth_boxes,1),1);
  for k=2:length(n_windows)
    for i_r=(n_windows(k-1)+1):n_windows(k)
      for i_g=1:n_g
        if(~already_found(i_g))
          pascalScore = computePascalScore(ranked_boxes(i_r,:),ground_truth_boxes(i_g,:));
          if(pascalScore > overlap_criterion)
            n_covered(k-1)=n_covered(k-1)+1;
            covering_ids=[covering_ids;i_r];
            covered_ids=[covered_ids;i_g];
            already_found(i_g)=1;
          end
        end
      end
    end
  end
  object_ids=(1:size(ground_truth_boxes,1))';
else
  for k=2:length(n_windows)
    for i_r=(n_windows(k-1)+1):n_windows(k)
      for i_g=1:n_g
        pascalScore = computePascalScore(ranked_boxes(i_r,:),ground_truth_boxes(i_g,:));
        if(pascalScore > overlap_criterion)
          n_covered(k-1)=n_covered(k-1)+1;
          covering_ids=[covering_ids;i_r];
          covered_ids=[covered_ids;i_g];
          object_ids=[object_ids;i_g];
        end
      end
    end
  end
end

uncovered_ids = setxor(covered_ids,1:size(ground_truth_boxes,1));

covering_boxes=ranked_boxes(covering_ids,:);
covered_boxes=ground_truth_boxes(covered_ids,:);
uncovered_boxes=ground_truth_boxes(uncovered_ids,:);

n_covered=cumsum(n_covered);
end

author: Cloud-CV

repo: object-proposals


text:
function x = p14_start ( n )

%*****************************************************************************80
%
%% P14_START returns a starting point for optimization for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = zeros ( n, 1 );

  for i = 1 : n
    if ( mod ( i, 2 ) == 1 )
      x(i) = - 1.2;
    else
      x(i) = 1.0;
    end
  end

  return
end


text:
classdef MimContextSetMapping < handle
    % MimContextSetMapping. A class for specifying the relationship between
    % context sets.
    %
    % This class is used by the Framework to store the relationship between
    % a context set (a set of related contexts) and its contexts; and
    % also the hierarchy of contexts. For example, the LeftLung and RightLung
    % contexts are both in the context set SingleLung. The SingleLung set is a
    % child of the LungROI context set.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ContextSet  % The context set enum (of type MimContextSet)

        ContextList % The set of MimContextMappings corresponding to contexts of this type

        Parent      % Parent type (of type MimContextSetMapping)
        Children    % Child type (of type MimContextSetMapping)
    end

    methods
        function obj = MimContextSetMapping(context_set_id, parent_context_set_mapping)
            obj.ContextSet = context_set_id;
            obj.Parent = parent_context_set_mapping;
            obj.ContextList = {};
            obj.Children = {};

            if ~isempty(parent_context_set_mapping)
                obj.Parent.AddChild(obj);
            end
        end

        function AddChild(obj, child)
            obj.Children{end + 1} = child;
        end

        function AddContext(obj, context)
            obj.ContextList{end + 1} = context;
        end

        % Determines if the specified context set is higher in the hierarchy
        % than this set
        function is_higher = IsOtherContextSetHigher(obj, other_context_set)
            next_set = obj.Parent;
            while ~isempty(next_set)
                if next_set == other_context_set
                    is_higher = true;
                    return;
                end
                next_set = next_set.Parent;
            end
            is_higher = false;
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function ok = all_set_item(item)

% function ok = all_set_item(item)
% Perform within-item all_set check. For repeats, this is true, if item.val
% has between item.num(1) and item.num(2) elements.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: all_set_item.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok
ok = (numel(item.cfg_item.val) >= item.num(1)) &&...
    (numel(item.cfg_item.val) <= item.num(2));


text:
function [result] = mc(x)

%MC calculates the medcouple, a robust measure of skewness.
%
% The medcouple is described in:
%    Brys, G., Hubert, M. and Struyf, A. (2004),
%    "A Robust Measure of Skewness",
%    Journal of Computational and Graphical Statistics,
%    13, 996-1017.
%
% Required input arguments:
%    x : Data matrix.
%        Rows of x represent observations, and columns represent variables.
%        Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows)
%        with missing or infinite values will automatically be excluded from the computations.
%
% The output of MC is a vector containing the medcouple for each column
%     of the data matrix x.
%
% Example:
%    result = mc([chi2rnd(5,1000,1) trnd(3,1000,1)]);
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Guy Brys
% Last Update: 31/07/2007

if (nargin<1)
    error('No input arguments')
end
[n,p] = size(x);
if n==1
    x = x';
    p = 1;
end
x(sum(~isfinite(x),2)>0,:)=[];
n = size(x,1);
if (n>50000)
    error('When there are more than 50000 observations, the MC may be uncomputable due to memory limitations.')
elseif (n>100)
    result = mcc(x);
elseif n==0
    error('Due to missing values, the data matrix is empty and the MC can not be computed.')
else
    result = 0.5*(-mcc(repmat(max(x),size(x,1),1)-x)+mcc(x));
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit