text:
classdef MultinomialResamplerX < ResamplerX & matlabshared.tracking.internal.MultinomialResampler
% MULTINOMIALRESAMPLERX Class
%
% Summary of MultinomialResamplerX:
% This is a class implementation of a multinomial resampler.
%
% MultinomialResamplerX Properties:
%   None
%
% MultinomialResamplerX Methods:
%    MultinomialResamplerX  - Constructor method
%    resample - Perform multinomial resampling
%
% See also SystematicResamplerX
    properties
    end

    methods
        function this = MultinomialResamplerX(varargin)
        % MULTINOMIALRESAMPLERX Constructor method
        %
        % Usage:
        % * MultinomialResamplerX() returns a "MultinomialResamplerX" object
        %   handle
        %
        % See also MultinomialResamplerX/resample

        end

        function [newSamples, newWeights,idx] = resample(this,samples,weights,Nnew)
        % RESAMPLE Perform multinomial resampling
        %
        % Parameters
        % ----------
        % samples: (NumDims x NumSamples) matrix
        %   A matrix, whose columns correspond to individual samples
        % weights: (1 x NumSamples) row vector
        %   A row vector, whose columns correspond to the weights of the
        %   respective columns/samples of the samples matrix
        % Nnew: scalar, optional
        %   The number of samples to be resampled.
        %   (default = NumSamples, which implies that the same ammount of
        %   samples will be generated, as the number of samples provided)
        %
        % Returns
        % -------
        % newSamples: (NumDims x Nnew) matrix
        %   A matrix, whose columns correspond to resampled samples.
        % newWeights: (1 x Nnew) row vector
        %   A row vector, whose columns correspond to the weights of the
        %   respective columns/samples of the newSamples matrix
        % idx: (1 x Nnew) row vector
        %   A row vector of indices, mapping each sample in newSamples to the
        %   corresponding sample in samples.
        %
        % Usage
        % -----
        % * idx = resample(this,weights) returns a row vector idx containing
        %   the indexes of samples generated through by performing multinomial
        %   resampling based on the vector weights.
        %
        % See also SystematicResamplerX/resample

            N = numel(weights);
            if(nargin<4)
                Nnew = N;
            end
            % idx = randsample(1:N, Nnew, true, weights);
            idx = resample@matlabshared.tracking.internal.MultinomialResampler(this,weights,Nnew);
            newSamples = samples(:,idx);
            newWeights = repmat(1/Nnew, 1, Nnew);
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [x, t] = timeseriesdata(data, nin);

% TIMESERIESDATA make a time series data set with the given window length.

% NDLUTIL

ndata = length(data);
alldata = ones(ndata-nin+1, nin+1); % [y_{t-1], ..., y_{t-m]]

% Creates the delay vectors
for j = 1:ndata-nin
  alldata(j, :) = (data(j:j+nin)');
end

% Training set
x = alldata(:,1:nin);
t = alldata(:,nin+1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function lik = lik_gaussian(varargin)
%LIK_GAUSSIAN  Create a Gaussian likelihood structure
%
%  Description
%    LIK = LIK_GAUSSIAN('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a Gaussian likelihood structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    LIK = LIK_GAUSSIAN(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for Gaussian likelihood function [default]
%      sigma2       - variance [0.1]
%      sigma2_prior - prior for sigma2 [prior_logunif]
%      n            - number of observations per input (See using average
%                     observations below)
%
%    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.
%
%    Using average observations
%    The lik_gaussian can be used to model data where each input vector is
%    attached to an average of varying number of observations. That is, we
%    have input vectors x_i, average observations y_i and sample sizes n_i.
%    Each observation is distributed
%
%        y_i ~ N(f(x_i), sigma2/n_i)
%
%    The model is constructed as lik_gaussian('n', n), where n is the same
%    length as y and collects the sample sizes.
%
%  See also
%    GP_SET, PRIOR_*, LIK_*

% Internal note: Because Gaussian noise can be combined
% analytically to the covariance matrix, lik_gaussian is internally
% little between lik_* and gpcf_* functions.
%
% Copyright (c) 2007-2017 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

% 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_GAUSSIAN';
  ip.addOptional('lik', [], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('sigma2',0.1, @(x) isscalar(x) && x>=0);
  ip.addParamValue('sigma2_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('n',[], @(x) isreal(x) && all(x>0));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Gaussian';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Gaussian')
      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('n',ip.UsingDefaults)
    lik.n = ip.Results.n;
  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
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_gaussian_pak;
    lik.fh.unpak = @lik_gaussian_unpak;
    lik.fh.ll = @lik_gaussian_ll;
    lik.fh.llg = @lik_gaussian_llg;
    lik.fh.llg2 = @lik_gaussian_llg2;
    lik.fh.llg3 = @lik_gaussian_llg3;
    lik.fh.lp = @lik_gaussian_lp;
    lik.fh.lpg = @lik_gaussian_lpg;
    lik.fh.cfg = @lik_gaussian_cfg;
    lik.fh.tiltedMoments = @lik_gaussian_tiltedMoments;
    lik.fh.trcov  = @lik_gaussian_trcov;
    lik.fh.trvar  = @lik_gaussian_trvar;
    lik.fh.predy = @lik_gaussian_predy;
    lik.fh.siteDeriv = @lik_gaussian_siteDeriv;
    lik.fh.recappend = @lik_gaussian_recappend;
  end

end

function [w, s, h] = lik_gaussian_pak(lik)
%LIK_GAUSSIAN_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_GAUSSIAN_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.magnSigma2)]'
%
%  See also
%    LIK_GAUSSIAN_UNPAK

  w = []; s = {}; h=[];
  if ~isempty(lik.p.sigma2)
    w = [w log(lik.sigma2)];
    s = [s; 'log(gaussian.sigma2)'];
    h = [h 0];
    % Hyperparameters of sigma2
    [wh, sh, hh] = lik.p.sigma2.fh.pak(lik.p.sigma2);
    w = [w wh];
    s = [s; sh];
    h = [h hh];
  end

end

function [lik, w] = lik_gaussian_unpak(lik, w)
%LIK_GAUSSIAN_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_GAUSSIAN_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.magnSigma2)]'
%
%  See also
%    LIK_GAUSSIAN_PAK

  if ~isempty(lik.p.sigma2)
    lik.sigma2 = exp(w(1));
    w = w(2:end);

    % Hyperparameters of sigma2
    [p, w] = lik.p.sigma2.fh.unpak(lik.p.sigma2, w);
    lik.p.sigma2 = p;
  end
end

function logLik = lik_gaussian_ll(lik, y, f, ~)
%LIK_GAUSSIAN_LL    Log likelihood
%
%  Description
%    E = LIK_GAUSSIAN_LL(LIK, Y, F, Z) takes a likelihood data
%    structure LIK, incedence counts Y, expected counts Z, and
%    latent values F. Returns the log likelihood, log p(y|f,z).
%    This subfunction is needed when using Laplace approximation
%    or MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria (DIC, WAIC)
%    computations.
%
%  See also
%    LIK_GAUSSIAN_LLG, LIK_GAUSSIAN_LLG3, LIK_GAUSSIAN_LLG2, GPLA_E

  s2 = lik.sigma2;
  r2 = (f-y).^2;
  logLik =  sum(-0.5 * r2./s2 - 0.5*log(s2) - 0.5*log(2*pi));

end

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

switch param
    case 'param'
        % there is also correction due to the log transformation
        s2 = lik.sigma2;
        r2 = (f-y).^2;
        llg =  0.5.* sum( r2./s2 -1 );

    case 'latent'
        s2 = lik.sigma2;
        r = f-y;
        llg =  - r./s2 ;
end
end

function llg2 = lik_gaussian_llg2(lik, y, f, param, ~)
%LIK_GAUSSIAN_LLG2  Second gradients of the log likelihood
%
%  Description
%    G2 = LIK_GAUSSIAN_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z,
%    and latent values F. Returns the Hessian of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. G2 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_GAUSSIAN_LL, LIK_GAUSSIAN_LLG, LIK_GAUSSIAN_LLG3, GPLA_E

switch param
    case 'latent'
        s2 = lik.sigma2;
        llg2 =  - ones(size(y))./s2 ;
    case 'latent+param'
        % there is also correction due to the log transformation
        s2 = lik.sigma2;
        r = f-y;
        llg2 =  r./s2 ;

end
end

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

switch param
    case 'latent'
        llg3 = zeros(size(y));
    case 'latent2+param'
        % there is also correction due to the log transformation
        s2 = lik.sigma2;
        llg3 = ones(size(y))./s2 ;
end
end

function lp = lik_gaussian_lp(lik)
%LIK_GAUSSIAN_LP  Evaluate the log prior of likelihood parameters
%
%  Description
%    LP = LIK_T_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters.
%    This subfunctions is needed when there are likelihood
%    parameters.
%
%  See also
%    LIK_GAUSSIAN_PAK, LIK_GAUSSIAN_UNPAK, LIK_GAUSSIAN_G, GP_E

  lp = 0;

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

function lpg = lik_gaussian_lpg(lik)
%LIK_GAUSSIAN_LPG  Evaluate gradient of the log prior with respect
%                  to the parameters.
%
%  Description
%    LPG = LIK_GAUSSIAN_LPG(LIK) takes a Gaussian likelihood
%    function structure LIK and returns LPG = d log (p(th))/dth,
%    where th is the vector of parameters. This subfunction is
%    needed when there are likelihood parameters.
%
%  See also
%    LIK_GAUSSIAN_PAK, LIK_GAUSSIAN_UNPAK, LIK_GAUSSIAN_E, GP_G

  lpg = [];

  if ~isempty(lik.p.sigma2)
    likp=lik.p;

    lpgs = likp.sigma2.fh.lpg(lik.sigma2, likp.sigma2);
    lpg = lpgs(1).*lik.sigma2 + 1;
    if length(lpgs) > 1
      lpg = [lpg lpgs(2:end)];
    end
  end
end

function DKff = lik_gaussian_cfg(lik, x, x2)
%LIK_GAUSSIAN_CFG  Evaluate gradient of covariance with respect to
%                 Gaussian noise
%
%  Description
%    Gaussian likelihood is a special case since it can be
%    analytically combined with covariance functions and thus we
%    compute gradient of covariance instead of gradient of likelihood.
%
%    DKff = LIK_GAUSSIAN_CFG(LIK, X) takes a Gaussian likelihood
%    function structure LIK, a matrix X of input vectors and
%    returns DKff, the gradients of Gaussian noise covariance
%    matrix Kff = k(X,X) with respect to th (cell array with
%    matrix elements). This subfunction is needed only in Gaussian
%    likelihood.
%
%    DKff = LIK_GAUSSIAN_CFG(LIK, X, X2) takes a Gaussian
%    likelihood function structure LIK, a matrix X of input
%    vectors and returns DKff, the gradients of Gaussian noise
%    covariance matrix Kff = k(X,X) with respect to th (cell
%    array with matrix elements). This subfunction is needed
%    only in Gaussian likelihood.
%
%  See also
%    LIK_GAUSSIAN_PAK, LIK_GAUSSIAN_UNPAK, LIK_GAUSSIAN_E, GP_G

  DKff = {};
  if ~isempty(lik.p.sigma2)
      if isempty(lik.n)
          DKff{1}=lik.sigma2;
      else
          n=size(x,1);
          DKff{1} = sparse(1:n, 1:n, lik.sigma2./lik.n, n, n);
      end
  end
end

function DKff  = lik_gaussian_ginput(lik, x, t, g_ind, gdata_ind, gprior_ind, varargin)
%LIK_GAUSSIAN_GINPUT  Evaluate gradient of likelihood function with
%                     respect to x.
%
%  Description
%    DKff = LIK_GAUSSIAN_GINPUT(LIK, X) takes a likelihood
%    function structure LIK, a matrix X of input vectors and
%    returns DKff, the gradients of likelihood matrix Kff =
%    k(X,X) with respect to X (cell array with matrix elements).
%    This subfunction is needed only in Gaussian likelihood.
%
%    DKff = LIK_GAUSSIAN_GINPUT(LIK, X, X2) takes a likelihood
%    function structure LIK, a matrix X of input vectors and
%    returns DKff, the gradients of likelihood matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed only in Gaussian likelihood.
%
%  See also
%    LIK_GAUSSIAN_PAK, LIK_GAUSSIAN_UNPAK, LIK_GAUSSIAN_E, GP_G

end

function C = lik_gaussian_trcov(lik, x)
%LIK_GAUSSIAN_TRCOV  Evaluate training covariance matrix
%                    corresponding to Gaussian noise
%
%  Description
%    C = LIK_GAUSSIAN_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 subfunction is needed only in Gaussian likelihood.
%
%  See also
%    LIK_GAUSSIAN_COV, LIK_GAUSSIAN_TRVAR, GP_COV, GP_TRCOV

  [n, m] =size(x);
  n1=n+1;

  if isempty(lik.n)
      C = sparse(1:n,1:n,ones(n,1).*lik.sigma2,n,n);
  else
      C = sparse(1:n, 1:n, lik.sigma2./lik.n, n, n);
  end

end

function C = lik_gaussian_trvar(lik, x)
%LIK_GAUSSIAN_TRVAR  Evaluate training variance vector
%                    corresponding to Gaussian noise
%
%  Description
%    C = LIK_GAUSSIAN_TRVAR(LIK, TX) takes in covariance function
%    of a Gaussian process LIK and matrix TX that contains
%    training inputs. Returns variance vector C. Every element i
%    of C contains variance of input i in TX. This subfunction is
%    needed only in Gaussian likelihood.
%
%
%  See also
%    LIK_GAUSSIAN_COV, GP_COV, GP_TRCOV

  [n, m] =size(x);
  if isempty(lik.n)
      C=repmat(lik.sigma2,n,1);
  else
      C=lik.sigma2./lik.n(:);
  end

end

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

  lpy = [];
  Ey = Ef;
  EVary = lik.sigma2;
  VarEy = Varf;
  Vary = EVary + VarEy;
  if numel(yt) ~= 0
      lpy = norm_lpdf(yt,Ey,sqrt(Vary));
  end
end

function [logM_0, m_1, sigm2hati1] = lik_gaussian_tiltedMoments(lik, y, i1, sigm2_i, myy_i, z)
%LIK_PROBIT_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_PROBIT_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY) takes a likelihood structure LIK, class labels Y, 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

%   m_1=myy_i;
%   sigm2hati1=sigm2_i;
%   logM_0=zeros(size(y));

s2 = lik.sigma2;
tau = 1./s2 + 1./sigm2_i;
w = (y(i1)./s2 + myy_i./sigm2_i)./tau;
%Zi = 1 ./( sqrt( 2.*pi.*(s2+sigm2_i) ) ) .* exp( -0.5*(y(i1)-myy_i).^2./(s2+sigm2_i) );

m_1 = w;
sigm2hati1 = 1./tau;
%logM_0 = log(Zi)
logM_0 = -0.5*log( 2.*pi) - 0.5*log( (s2+sigm2_i) )  + ( -0.5*(y(i1)-myy_i).^2./(s2+sigm2_i) );

end

function [g_i] = lik_gaussian_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_NEGBIN_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_NEGBIN_SITEDERIV(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 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

s2 = lik.sigma2;
tau = 1/s2 + 1/sigm2_i;
w = (y(i1)/s2 + myy_i/sigm2_i)/tau;
%Zi = 1/( sqrt(2*pi*(s2+sigm2_i)) )*exp( -0.5*(y(i1)-myy_i)^2/(s2+sigm2_i) );

%g_i = 0.5*( (1/tau + w.^2 -2*w*y(i1) + y(i1).^2 ) /s2 - 1)/s2 * s2;
g_i = 0.5*( (1/tau + w.^2 -2*w*y(i1) + y(i1).^2 ) /s2 - 1);

end

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

  if nargin == 2
    % Initialize the record
    reclik.type = 'Gaussian';

    % Initialize the parameters
    reclik.sigma2 = [];
    reclik.n = [];

    % Set the function handles
    reclik.fh.pak = @lik_gaussian_pak;
    reclik.fh.unpak = @lik_gaussian_unpak;
    reclik.fh.lp = @lik_gaussian_lp;
    reclik.fh.lpg = @lik_gaussian_lpg;
    reclik.fh.ll= @lik_gaussian_ll;
    reclik.fh.llg = @lik_gaussian_llg;
    reclik.fh.llg2 = @lik_gaussian_llg2;
    reclik.fh.llg3 = @lik_gaussian_llg3;
    reclik.fh.tiltedMoments = @lik_gaussian_tiltedMoments;
    reclik.fh.siteDeriv = @lik_gaussian_siteDeriv;
    reclik.fh.cfg = @lik_gaussian_cfg;
    reclik.fh.trcov  = @lik_gaussian_trcov;
    reclik.fh.trvar  = @lik_gaussian_trvar;
    reclik.fh.predy = @lik_gaussian_predy;
    reclik.fh.recappend = @lik_gaussian_recappend;

    reclik.p=[];
    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 n if given
    if isfield(lik,'n') && ~isempty(lik.n)
      reclik.n(ri,:)=lik.n(:)';
    end
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
function phometricParams = SRPhotometricParams(numFrames)
%SRPHOTOMETRICPARAMS Photometric parameters used for super-resolution.

    % The multiplicative part of the photometric model.
    % Default: Affine model with multiplicative factor '1'.
    phometricParams.mult = ones(numFrames, 1);

    % The additive part of the photometric model.
    % Default: Affine model with additive factor '0'.
    phometricParams.add = zeros(numFrames, 1);

author: thomas-koehler

repo: SupER


text:
% Fig. 9.15   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
% script to generate Figure 9.15
sys=tf(1,[1 0.2 1 0]);
rlocus(sys)
title('Figure 9.15  Root locus of 1/(s^2+0.2s + 1)')
axis([-3 1 -1.5 1.5])
hold on
r=roots([1 .2 1 .5]);
 plot(r,'*')
 z=0:.1:.9;
wn= .5:.5:3;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function [ r,f ] = GetReward( x )
% MountainCarGetReward returns the reward at the current state
% x: a vector of position and velocity of the car
% r: the returned reward.
% f: true if the car reached the goal, otherwise f is false

position = x(1);
% bound for position; the goal is to reach position = 0.5
bpright=0.5;
f=false;
r=-1;

if( position >= bpright)
	r = 100;
    f = true;
end

author: epfl-lasa


text:
function y = mod(varargin)
% Allows inputs to be char
if ischar(varargin{1})
    varargin{1} = double(varargin{1});
end
if ischar(varargin{2})
    varargin{2} = double(varargin{2});
end
if ~any(cellfun(@(x) isa(x,'sym'), varargin))
    y = builtin('mod', varargin{:});
else
    y = builtin('@sym/mod', varargin{:});
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [a,t]=v_rotqr2ax(q)
%V_ROTQR2AX converts a real quaternion to the corresponding rotation axis and angle
% Inputs:
%
%     Q(4,1)   real-valued quaternion (with magnitude = 1)
%
% Outputs:
%
%     A(3,1)   Unit vector in the direction of the rotation axis.
%     T        Rotation angle in radians (in range 0 to 2pi)
%
% In the quaternion representation of a rotation, and q(1) = cos(t/2)
% where t is the angle of rotation in the range 0 to 2pi
% and q(2:4)/sin(t/2) is a unit vector lying along the axis of rotation
% a positive rotation about [0 0 1] takes the X axis towards the Y axis.
%
%      Copyright (C) Mike Brookes 2007-2018
%      Version: $Id: v_rotqr2ax.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=q(2:4);
m=sqrt(a'*a);
t=2*atan2(m,q(1));      % avoids problems if unnormalized
if m>0
    a=a/m;
else
    a=[0 0 1]';
end
if ~nargout
    v_rotqr2ro(q); % plot a rotated cube
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% undo(n=1) removes the last n figure items added

function undo(n)

if nargin==0 n=1; end
if ischar(n) n = sscanf(n,'%i'); end

if n<1 error(['Not undoing ' int2str(n) ' items']); end
g=get(gca,'children');
L=length(g);
for i=1:n
	delete(g(i));
end


text:
function ehmm = obsupdate_ehmm(Gamma,ehmm,residuals,XX,Tfactor)
%
% Update observation model
%
% INPUT
% X             observations
% T             length of series
% Gamma         p(state given X)
% ehmm           ehmm data structure
% residuals     in case we train on residuals, the value of those.
%
% OUTPUT
% ehmm           estimated ehmm model
%
% Author: Diego Vidaurre, OHBA, University of Oxford / Aarhus Univ (2022)

% Some stuff that will be later used
if nargin<6, Tfactor = 1; end

K = ehmm.K;

Gamma = [Gamma prod(1-Gamma,2)];
Gamma = rdiv(Gamma,sum(Gamma,2));

XXGXX = cell(K,1);
for k = 1:K+1
    XXGXX{k} = bsxfun(@times, XX, Gamma(:,k))' * XX;
end

% autoregression coefficients
[ehmm,XW] = updateW_ehmm(ehmm,Gamma,residuals,XX,XXGXX,Tfactor);

% Omega
ehmm = updateOmega_ehmm(ehmm,Gamma,residuals,XX,XXGXX,XW,Tfactor);
%
% % autoregression coefficient priors
ehmm = updateSigma_ehmm(ehmm); % sigma - channel x channel coefficients
ehmm = updateAlpha_ehmm(ehmm); % alpha - one per order

end

% function ehmm = updateW_here(ehmm,Gamma,residuals,XX)
%
% K = size(Gamma,2); np = size(XX,2);
% ndim = size(ehmm.state(end).W.Mu_W,2);
% noGamma = prod(1-Gamma,2);
% residuals0 = residuals;
% m0 = zeros(size(residuals0));
%
% for n = 1:ndim
%     m0(:,n) = bsxfun(@times,XX * ehmm.state(end).W.Mu_W(:,n),noGamma);
%     residuals(:,n) = residuals(:,n) - m0(:,n);
% end
%
% %Gamma = [Gamma noGamma];
% X = zeros(size(XX,1),np * K);
% for k = 1:K
%    X(:,(1:np) + (k-1)*np) = bsxfun(@times, XX, Gamma(:,k));
% end
% % estimation
% for n = 1:ndim
%     ehmm.state_shared(n).Mu_W = ...
%         pinv(X) * residuals(:,n);
% end
%
% % meand = zeros(size(XX,1),ndim);
% % for n = 1:ndim
% %      W = ehmm.state_shared(n).Mu_W;
% %     meand(:,n) = meand(:,n) + X * W;
% % end
% % d = residuals - meand;
% % mean(d.^2)
%
%
% end

% function [e,meand] = getError_here(Gamma,hmm,residuals,XX)
%
% C = hmm.Omega.Gam_shape ./ hmm.Omega.Gam_rate;
% meand = computeStateResponses_here(XX,hmm,Gamma);
% d = residuals - meand;
%
% [ mean(d.^2) ]
%
% Cd = bsxfun(@times, C, d)';
% dist = zeros(size(residuals,1),1);
% for n = 1:size(residuals,2)
%     dist = dist + 0.5 * (d(:,n).*Cd(n,:)');
% end
% e = sum(dist);
%
% end
%
%
% function meand = computeStateResponses_here(XX,ehmm,Gamma)
%
% K = size(Gamma,2); np = size(XX,2);
% ndim = size(ehmm.state(end).W.Mu_W,2);
% noGamma = prod(1-Gamma,2);
% meand = zeros(size(XX,1),ndim);
% for n = 1:ndim
%     meand(:,n) = meand(:,n) + ...
%         bsxfun(@times,XX * ehmm.state(end).W.Mu_W(:,n),noGamma);
%
% end
% % meand = zeros(size(XX,1),ndim);
% %Gamma = [Gamma noGamma];
% X = zeros(size(XX,1),np * K);
% for k = 1:K
%    X(:,(1:np) + (k-1)*np) = bsxfun(@times, XX, Gamma(:,k));
% end
% for n = 1:ndim
%     %W = [ehmm.state_shared(n).Mu_W; ehmm.state(end).W.Mu_W(:,n)];
%     W = ehmm.state_shared(n).Mu_W;
%     %meand(:,n) = meand(:,n) + X * ehmm.state_shared(n).Mu_W;
%     meand(:,n) = meand(:,n) + X * W;
% end
%
% end

author: OHBA-analysis

repo: HMM-MAR


text:
%% DEMOSLOPE  Short demonstration of slopes
%

%% Some sample applications of the Hessian toolbox
% Hessians implement second order automatic differentiation in forward mode, which is
% conveniently to implement using the Matlab operator concept.
%
% For a demonstration of the slope toolbox, intervals are wide; therefore, output
% is switched to infimum/supremum representation (see demointval):

setround(0)                           % set rounding to nearest
format compact short infsup

%% Slope expansions
% In order to use slope expansions, the expansion range and expansion
% point (interval) need to be identified and values have to be assigned.
% This is performed by the function "slopeinit", for example

u = slopeinit( [1;2] , infsup([1;2],[1.1;2.1]) )

%%
% The length (number of rows) is the number
% of independent variables, in the example above two. Slopes are always
% of type intval. For definition and use of slopes cf. Neumaier's book
% on interval analysis.

%% Computing with slopes
% Automatic slope expansion is performed by operations with slope variables.
% For example,

x = slopeinit(3.5,midrad(3.5,1e-4));
y = exp(3*x-sqrt(x))

%%
% evaluates the expression using slopes with expansion point 3.5 and expansion
% interval 3.5+/-1e-4.
%
% There is access to the center value, the range and
% the slope of y:

y.c
y.r
y.s

%%
% When evaluating the expression for another argument, e.g. a vector argument,
% use the same statement as before with new values.

xs = [2;3];
x = slopeinit(xs,midrad(xs,1e-4));
y = exp(3*x-sqrt(x))

%% Graphical illustration of slopes
% There is a simple graphic to demonstrate the behaviour of slope expansions for
% one-dimensional functions. The routine "slopeplot" assumes a function
% in one variable which can be evaluated for vectors, the expansion
% point and the expansion interval.
% For example,

format short
slopeplot('(x-2).*exp(x)-sqrt(abs(x))',1,infsup(1,2))

%%
% shows a best possible upper slope, but an underestimated lower slope.

%% Non-differentiable functions
% Functions to be handled by slopes need not be differentiable.
% An example is sqrt(abs(x)). Consider

slopeplot('sqrt(abs(x))',-2,infsup(-1,1))

%%
% The slope is best possible, but the grid is not fine enough to catch the
% extreme point at zero. For some 10000 grid points the picture looks
% as follows:

slopeplot('sqrt(abs(x))',-2,infsup(-1,1),[],10000)

%% Computation of the Gamma function with slopes
% Consider the following example in one unknown.
% According to Stirling's formula it is for u -> inf,
%
%                          1      1       139        571
%    Gamma(u) ~ C * ( 1 + --- + ----- - ------- - --------- + ... )
%                         12u       2         3           4
%                               288u    51840u    2488320u
% with
%
%         -u  u-0.5
%    C = e   u      sqrt(2*pi) .
%
% The following function evaluates Stirling's formula. It is also
% suited for vector input.
%
%
%   function y = g(u)
%       C = exp(-u) .* ( u.^(u-0.5) ) * sqrt(2.0*pi) ;
%       v = (((( -571.0/2488320.0 ./ u - 139.0/51840.0 ) ./ u ...
%                 + 1.0/288.0) ./ u ) + 1.0/12.0 ) ./ u + 1.0;
%       y = C .* v;
%
% A corresponding inline function is

format long e
g = @(u) ( ( exp(-u) .* ( u.^(u-0.5) ) * sqrt(2.0*pi) ) .* ...
           ( (((( -571.0/2488320.0 ./ u - 139.0/51840.0 ) ./ u ...
             + 1.0/288.0) ./ u ) + 1.0/12.0 ) ./ u + 1.0 ) )
u = [ 3.5 61 5 ]
g(u)

%% Computation of the range of the function g
% Next calculate the range of the function g within
% a certain interval. This can be performed, for example, by straightforward
% interval evluation (naive interval arithmetic)

format long
X = infsup(4.1,4.2);
Y = g(X)

%%
% or by slopes:

xs = 4.1;
X = infsup(4.1,4.2);
Ys = g(slopeinit(xs,X))

%%
% The range computed by the slope expansion is better by a factor 2 than
% naive interval arithmetic and less than 20 % overestimation of the true range:

diam(Y)
diam(Ys.r)
g(4.2)-g(4.1)

%%
% Note that Y or Ys.r is an inclusion of the range of g; for an inclusion
% of the range of the Gamma function an error
% term has to be added.

%% Slopes in several unknowns
% Automatic slope expansion with several unknowns works the same way.
% Consider the following example by Broyden.
%
%                                .5*sin(x1*x2) - x2/(4*pi) - x1/2  =  0
%  (1-1/(4*pi))*(exp(2*x1)-exp(1)) + exp(1)*x2/pi - 2*exp(1)*x1 )  =  0
%
% with initial approximation [ .6 ; 3 ] and one solution [ .5 ; pi ].
% The following function evaluates Broyden's function.
%
%  function  y = f1(x)
%     y = x;
%     y(1) = .5*sin(x(1)*x(2)) - x(2)/(4*pi) - x(1)/2;
%     y(2) = (1-1/(4*pi))*(exp(2*x(1))-exp(1)) + exp(1)*x(2)/pi - 2*exp(1)*x(1);
%
% The first statement is for efficiency. It is generally better to fix the size
% of an array before assigning values to the components. An inline function
% is as follows (cf. demogradient):

f = @(x) ( [ .5*sin(x(1)*x(2)) - x(2)/(4*pi) - x(1)/2 ; ...
             (1-1/(4*pi))*(exp(2*x(1))-exp(1)) + exp(1)*x(2)/pi - 2*exp(1)*x(1) ] )

%% Solution of a nonlinear system
% The nonlinear system defined by Broyden's function can solved by Newton's procedure
% as follows (cf. demogradient):

x = gradientinit([ .6 ; 3 ]);
for i=1:5
  y = f1(x);
  x = x - y.dx\y.x;
end
x

%%
% For simplicity, we omitted the stopping criterion.
% Here, y.dx is the Jacobian, y.x the function value at x.x, and -y.dx\y.x
% is the correction obtained by solution of a linear system.

%% Verified solution of a nonlinear system
% For verified solution of the nonlinear system, we need a correct definition
% of the function, see demogradient:
%
%  function  y = f(x)
%    y = x;
%    c1 = typeadj( 1 , typeof(x) );
%    cpi = typeadj( midrad(3.14159265358979323,1e-16) , typeof(x) );
%    y(1) = .5*sin(x(1)*x(2)) - x(2)/(4*cpi) - x(1)/2;
%    y(2) = (1-1/(4*cpi))*(exp(2*x(1))-exp(c1)) + exp(c1)*x(2)/cpi - 2*exp(c1)*x(1);
%
% This code is implemented in the function demotest.m .

%% Verified solution using verifynlss
% The nonlinear system defined by Broyden's function can be
% solved with verification and using slopes by:

y = verifynlss('demotest',[ .6 ; 3 ],'s')

%%
% The first parameter gives the name of the function, in this case "demotest",
% such that "demotest(x)" evaluates the function at "x".
% For the last parameter being 's', slopes are used, otherwise gradients (the default).
%
% We used slopes. However, the inclusion of the error with respect to an approximate
% is computed; therefore, the results for gradient and slope inclusion are identical.
% Note that use of gradients guarantees uniqueness of the zero within the computed
% interval, use of slope does not.

%% Nonlinear systems with uncertain parameters
% Next we artificially introduce an interval parameter of large diameter
% in Broyden's function:
%
%  function  y = demotest(x,radius)
%    y = x;
%    if nargin==1
%      radius = 1e-15;
%    end
%    cPi = typeadj( midrad(3.141592653589793,radius) , typeof(x) );
%    y(1) = .5*sin(x(1)*x(2)) - x(2)/(4*cPi) - x(1)/2;
%    y(2) = (1-1/(4*cPi))*(exp(2*x(1))-exp(1)) + exp(1)*x(2)/cPi - 2*exp(1)*x(1);
%
% This is to show that the range of applicability is larger for slopes than
% for gradients. For radius .04, both gradient expansion and slope expansion
% compute an inclusion (the extra parameter is passed to demotest):

radius = 4e-2; y1 = verifynlss('demotest',[ .6 ; 3 ],'g',[],radius)
radius = 4e-2; y2 = verifynlss('demotest',[ .6 ; 3 ],'s',[],radius)

%%
% The inclusion using slopes is better by 10 to 20 % in radius:

[ rad(y1) rad(y2) ]

%% A demo nonlinear function with large uncertainty
% Finally, we enlarge the radius of the parameter cPi to 0.05 and try to
% calculate an inclusion:

radius = 5e-2; y1 = verifynlss('demotest',[ .6 ; 3 ],'g',[],radius)
radius = 5e-2; y2 = verifynlss('demotest',[ .6 ; 3 ],'s',[],radius)

%%
% Now, only the slope expansion is able to compute an inclusion.

%% Enjoy INTLAB
% INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing,
% Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de

repo: OpenMAS


text:
function [H A b] = affineLS(X1, X2, normalization)

% [H A b] = affineLS(X1, X2, normalization)
%
% DESC:
% computes the affine transformation between the point pairs X1, X2
%
% VERSION:
% 1.0.0
%
% INPUT:
% X1, X2        = point matches (cartesian coordinates)
% normalization = true (default) or false to enable/disable point
%                 normalzation
%
% OUTPUT:
% H             = homography representing the affine transformation
% A             = 4x4 affine coefficient matrix
% b             = affine translation vector

% AUTHOR:
% Marco Zuliani, email: marco.zuliani@gmail.com
% Copyright (C) 2010 by Marco Zuliani
%
% LICENSE:
% This toolbox is distributed under the terms of the GNU GPL.
% Please refer to the files COPYING.txt for more information.

% HISTORY
% 1.0.0         11/22/10 - intial version

if (nargin < 3)
    normalization = true;
end;

N = size(X1, 2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% checks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (size(X2, 2) ~= N)
    error('AffineLS:inputError', ...
        'The set of input points should have the same cardinality')
end;
if N < 2
    error('AffineLS:inputError', ...
        'At least 3 point correspondences are needed')
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% normalize the input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if normalization
    % fprintf('\nNormalizing...')
    [X1, T1] = normalize_points(X1);
    [X2, T2] = normalize_points(X2);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Traditional LS
P = zeros(2*N, 6);
q = zeros(2*N, 1);

ind = 1:2;
for n = 1:N

    P(ind, :) = [X1(1,n) 0 X1(2,n) 0 1 0; 0 X1(1,n) 0 X1(2,n) 0 1];

    q(ind) = X2(1:2, n);

    ind = ind + 2;

end;

% solve the linear system in a least square sense
Theta = P\q;

% compute the corresponding homography
H = [Theta(1) Theta(3) Theta(5); Theta(2) Theta(4) Theta(6); 0 0 1];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% de-normalize the parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if normalization
    H = T2\H*T1;
end;
H = H/H(9);

% prepare the output
if nargout > 1

    A = H(1:2, 1:2);
    b = H(1:2, 3);

end;

return

author: RANSAC

repo: RANSAC-Toolbox


text:
function A10 = synA10max(A11, A00, cmin)
%-----------------------------------------------------------------------------
%
% For each point of colour 10 this function assigns the maximum value at the
% neighbouring gridpoints of colours 11 and 00.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 7, 2001.
% (c) 1998-2002 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
[n00, m00]=size(A00);
[n11, m11]=size(A11);
n10=n11;
m10=m00;
%[n10, m10]=size(A10);
if     m10 == m11
  S=max(A11, stripR(extL(A11, cmin)));
elseif m10 == m11+1
  S=max(extL(A11, cmin), extR(A11, cmin));
else
  disp([' size A11 = ' int2str(size(A11)) ' size A00 = ' int2str(size(A00))]);
  error(' synA10max - A11 and A00 do not match ');
end
if     n10 == n00
  T=max(A00, stripU(extD(A00, cmin)));
elseif n10 == n00-1
  T=max(stripD(A00), stripU(A00));
else
  disp([' size A11 = ' int2str(size(A11)) ' size A00 = ' int2str(size(A00))]);
  error(' synA10max - A11 and A00 do not match ');
end
%Note: all(size(S) == size(T)) & all(size(S) == [n10 m10]) always holds.
A10=max(S, T);
%-----------------------------------------------------------------------------


text:
function x = p36_start ( n )

%*****************************************************************************80
%
%% P36_START returns a starting point for optimization for problem 36.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ 0.5, 1.0 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%*****************************************************************
%% validate: validate data
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function  [blk,At,C,b,dim,nnblk,parbarrier] = ...
              validate(blk,At,C,b,par,parbarrier)

   if (nargin >= 5)
      spdensity = par.spdensity;
   else
      spdensity = 0.4;
   end
%%
   if ~iscell(blk);
      error('validate: blk must be a cell array'); end;
   if (size(blk,2) < 2)
      error('validate: blk must be a cell array with at least 2 columns');
   end
   if ~iscell(At) | ~iscell(C);
      error('validate: At, C must be cell arrays');
   end
   if (size(At,1) ~= size(blk,1))
      if (size(At,2) == size(blk,1));
         At = At';
      else
         error('validate: size of At is not compatible with blk');
      end
   end
   if (size(C,1) ~= size(blk,1))
      if (size(C,2) == size(blk,1))
         C = C';
      else
         error('validate: size of C is not compatible with blk');
      end
   end
   if (min(size(b)) > 1); error('validate: b must be a vector'); end
   if (size(b,1) < size(b,2)); b = b'; end
   m = length(b);
%%
%%-----------------------------------------
%% validate blk, At, C
%%-----------------------------------------
   ntotal = zeros(size(blk,1),1);
   for p = 1:size(blk,1)
      if (size(blk{p,2},1) > size(blk{p,2},2))
         blk{p,2} = blk{p,2}';
      end
      pblk = blk(p,:);
      n = sum(pblk{2});
      numblk = length(pblk{2});
      if strcmp(pblk{1},'s');
         m1 = size(At{p,1},2);
         n2 = sum(pblk{2}.*pblk{2});  n22 = sum(pblk{2}.*(pblk{2}+1))/2;
         ntotal(p) = n22;
         if ~all(size(C{p}) == n)
            error('validate: blk and C are not compatible');
         end
         if (norm(C{p}-C{p}',inf) > 1e-13*norm(C{p},inf));
            error('validate: C is not symmetric');
         end
         if (size(At{p,1}) == [m1, n22] & m1~=n22);
            At{p,1} = At{p,1}';
         end
         if (~isempty(At{p,1})) & (size(At{p,1},1) ~= n22)
            error('validate: blk and At not compatible');
         end
         if (nnz(At{p,1}) < spdensity*n22*m1) || (n <= 20) %%2nd condition: 2017-May-08
            if ~issparse(At{p,1}); At{p,1} = sparse(At{p,1}); end
         else
            if issparse(At{p,1}); At{p,1} = full(At{p,1}); end
         end
         if (length(pblk) > 2) %% for low rank constraints
            len = sum(pblk{3});
            if (size(pblk{1,3},2) < size(pblk{1,3},1))
                blk{p,3} = blk{p,3}';
            end
            if any(size(At{p,2})- [n,len])
               error(' low rank structure specified in blk and At not compatible')
            end
            if (length(At(p,:)) > 2)  & ~isempty(At{p,3})
               if all(size(At{p,3},2)-[1,4])
                  error(' low rank structure in At{p,3} not specified correctly')
               end
               if (size(At{p,3},2) == 1)
                  if (size(At{p,3},1) < size(At{p,3},2)); At{p,3} = At{p,3}'; end
                  lenn = length(At{p,3});
                  constrnum = mexexpand(pblk{3},[1:length(pblk{3})]');
                  At{p,3} = [constrnum, [1:lenn]', [1:lenn]', At{p,3}];
               elseif (size(At{p,3},2) == 4)
                  dd = At{p,3};
                  [dummy,idxsort] = sort(dd(:,1));
                  dd = dd(idxsort,:);
                  lenn = size(dd,1);
                  idxD = [0; find(diff(dd(:,1))); lenn];
                  ii = zeros(lenn,1); jj = zeros(lenn,1);
                  ss = [0,cumsum(pblk{3})];
                  for k = 1:length(pblk{3})
                     idx = [idxD(k)+1 : idxD(k+1)];
                     ii(idx) = dd(idx,2)+ss(k); %% convert to cumulative indexing
                     jj(idx) = dd(idx,3)+ss(k);
                  end
                  At{p,3} = [dd(:,1),ii,jj,dd(:,4)];
               end
            else
               constrnum = mexexpand(pblk{3},[1:length(pblk{3})]');
               At{p,3} = [constrnum,  [1:len]', [1:len]', ones(len,1)];
            end
         end
         if (nnz(C{p}) < spdensity*n2) | (numblk > 1);
            if ~issparse(C{p}); C{p} = sparse(C{p}); end;
         else
            if issparse(C{p}); C{p} = full(C{p}); end; %%added: 2017-May-08
         end
      elseif strcmp(pblk{1},'q') | strcmp(pblk{1},'l') ...
             | strcmp(pblk{1},'u');
         ntotal(p) = n;
         if (min(size(C{p})) ~= 1 | max(size(C{p})) ~= n);
            error(['validate: blk and C are not compatible']);
         end;
         if (size(C{p},1) < size(C{p},2)); C{p} = C{p}'; end
         if (size(At{p,1}) == [m, n] & m~=n);
            At{p,1} = At{p,1}';
         end
         if ~all(size(At{p,1}) == [n,m]);
            error('validate: blk and At not compatible');
         end
         if ~issparse(At{p,1});
            At{p,1} = sparse(At{p,1});
         end
         if (nnz(C{p}) < spdensity*n);
            if ~issparse(C{p}); C{p} = sparse(C{p}); end;
         else
            if issparse(C{p}); C{p} = full(C{p}); end;
         end;
      else
         error(' blk: some fields are not specified correctly');
      end
   end
   if (sum(ntotal) < m)
      error(' total dimension of C should be > length(b)');
   end
%%
%%-----------------------------------------
%% problem dimension
%%-----------------------------------------
%%
   dim = zeros(1,4);  nnblk = zeros(1,2);
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if strcmp(pblk{1},'s')
         dim(1) = dim(1) + sum(pblk{2});
         nnblk(1) = nnblk(1) + length(pblk{2});
         nn(p) = sum(pblk{2});
      elseif strcmp(pblk{1},'q')
         dim(2) = dim(2) + sum(pblk{2});
         nnblk(2) = nnblk(2) + length(pblk{2});
         nn(p) = length(pblk{2});
      elseif strcmp(pblk{1},'l')
         dim(3) = dim(3) + sum(pblk{2});
         nn(p) = sum(pblk{2});
      elseif strcmp(pblk{1},'u')
         dim(4) = dim(4) + sum(pblk{2});
         nn(p) = sum(pblk{2});
      end
   end
%%
%%-----------------------------------------
%% validate parbarrier
%%-----------------------------------------
%%
    if (nargin == 6)
       if ~iscell(parbarrier);
	  if (length(parbarrier) == size(blk,1))
	     tmp = parbarrier;
             clear parbarrier;
	     parbarrier = cell(size(blk,1),1);
	     for p = 1:size(blk,1)
		parbarrier{p} = tmp(p);
             end
          end
       end
       if (size(parbarrier,2) > size(parbarrier,1))
          parbarrier = parbarrier';
       end
       for p = 1:size(blk,1)
          pblk = blk(p,:);
          if (size(parbarrier{p},1) > size(parbarrier{p},2))
             parbarrier{p} = parbarrier{p}';
          end
          len = length(parbarrier{p});
          if strcmp(pblk{1},'s') | strcmp(pblk{1},'q')
             if (len == 1)
                parbarrier{p} = parbarrier{p}*ones(1,length(pblk{2}));
	     elseif (len == 0)
                parbarrier{p} = zeros(1,length(pblk{2}));
	     elseif (len ~= length(pblk{2}))
                error('blk and parbarrier not compatible');
             end
          elseif strcmp(pblk{1},'l')
             if (len == 1)
                parbarrier{p} = parbarrier{p}*ones(1,sum(pblk{2}));
	     elseif (len == 0)
		parbarrier{p} = zeros(1,sum(pblk{2}));
             elseif (len ~= sum(pblk{2}))
                error('blk and parbarrier not compatible');
             end
          elseif strcmp(pblk{1},'u')
             parbarrier{p}= zeros(1,sum(pblk{2}));
          end
       end
    end
%%***********************************************************************

author: shamilmamedov


text:
function [ error_per_image ] = compute_diag_error( ground_truth_all, detected_points_all )
%compute_error
%   compute the average point-to-point Euclidean error normalized by the
%   diag bounding box distance (measured as the Euclidean distance between the
%   outer corners of the eyes)
%
%   Inputs:
%          grounth_truth_all, size: num_of_points x 2 x num_of_images
%          detected_points_all, size: num_of_points x 2 x num_of_images
%   Output:
%          error_per_image, size: num_of_images x 1

num_of_images = size(ground_truth_all,3);
num_of_points = size(ground_truth_all,1);

error_per_image = zeros(num_of_images,1);

for i =1:num_of_images
    detected_points      = detected_points_all(:,:,i);
    ground_truth_points  = ground_truth_all(:,:,i);

    diag_distance = norm([min(ground_truth_points(:,1)) min(ground_truth_points(:,2))]-[max(ground_truth_points(:,1)) max(ground_truth_points(:,2))]);

    sum=0;
    for j=1:num_of_points
        sum = sum+norm(detected_points(j,:)-ground_truth_points(j,:));
    end
    error_per_image(i) = sum/(num_of_points*diag_distance);
end

end

author: jiankangdeng

repo: MenpoBenchmark


text:
function [V,E,F,C] = voronoi_tetgen(P,varargin)
  % VORONOI_TETGEN Compute a voronoi diagram using tetgen
  %
  %
  % [V,E,F,C] = voronoi_tetgen(P,varargin)
  %
  % Inputs:
  %   P  #P by 3 list of points
  % Outputs:
  %   V  #V by 3 list of orthocenters (?)
  %   E  #E by 2 list of edge indices into V
  %   F  #F by 3 list of face indices into ?
  %   C  #C list of cells (not supported)
  %
  % I don't trust the tetgen output...
  %

  % default values
  flags = '';
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Flags'}, {'flags'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  % get a temporary file name prefix
  prefix = tempname;
  prefix = tempname;
  node_filename = [prefix '.node'];
  writeNODE(node_filename,P);

  % call tetgen
  voronoi_flags = '-vNEF';
  command = [path_to_tetgen ' ' voronoi_flags ' ' flags ' ' node_filename];
  [status, result] = system(command);
  if status~=0
    error(result)
  end

  edge_filename = [prefix '.1.v.edge'];
  face_filename = [prefix '.1.v.face'];
  node_filename = [prefix '.1.v.node'];
  cell_filename = [prefix '.1.v.cell'];

  F = readFACE(face_filename,'ForceNoBoundary',true);
  E = readEDGE(edge_filename);
  V = readNODE(node_filename);
  C = [];
  if nargout>3
    error('Cell output not handled yet.');
  end

end

author: alecjacobson

repo: gptoolbox


text:
function [means,stats] = weighted_reg(Y,varargin)
% Calculate weighted average using weighted linear least squares
% See examples below for usage
%
% :Model:
%
%   Y_i = 1*Ypop + noise
%
% :Inputs:
%
%   **Y:**
%        data matrix (nsub x T)
%
%   **w:**
%        weights
%
%   **varY:**
%        variance of data at each time point (nsub x T) + var between
%
% :Outputs:
%
%   **Ymean:**
%        weighted mean of each column of Y
%
%   **dfe:**
%        error degrees of freedom, adjusted for inequality of variance
%        (Sattherwaite) and pooled across data columns
%
% Extended output in stats structure:
%
%   **stats.t:**
%        t-values for weighted t-test
%
%   **stats.p:**
%        2-tailed p-values for weighted t-test
%
%
%   **r:**
%        weighted correlation coeff across columns of Y
%
%   **xy:**
%        weighted covariance matrix
%
%   **v:**
%        weighted variance estimates for each column of Y
%      - sqrt(v) is the standard error of the mean (or grp difference)
%
%
%   **stats.fits:**
%        fits  for each group (Ymean by group), low contrast weight group then high
%        Fastest if no stats are asked for.
%
% Computation time:
% For FULL stats report
%  - Triples from 500 -> 1000 columns of Y, continues to increase
%
% For mean/dfe only, fast for full dataset (many columns of Y)
%
% :Examples:
% ::
%
%    % Basic multivariate stats for 1000 columns of dat, no weighting
%    % Multivariate covariances are meaningful if cols of Y are organized, e.g., timeseries
%    [means,stats] = weighted_reg(dat(:,1:1000));
%
%    % The same, but return univariate stats only (good for large Y)
%    [means,stats] = weighted_reg(dat,'uni');
%
% ..
%    NOTE: TOR CHANGED INPUT TO ASSUME THAT WE SHOULD ENTER VARWI + VARBETWEEN
% ..

% ..
%    Set up arguments
% ..

if nargin == 0, error('Must at least enter data as 1st argument.'); end

domultivariate = 1;     % multivariate covariance est for Y
dobtwn = 0;             % between-subjects contrast

bcon = []; Ydiff = []; w = []; varY = [];

for i = 1:length(varargin)
    arg = varargin{i};
    if ischar(arg)
        switch lower(arg)
            case 'w', w = varargin{i+1};
            case 'btwn', bcon = contrast_code(varargin{i+1});
            case 'vary', varY = varargin{i+1};
            case 'uni', domultivariate = 0;
        end
    end
end

[m,n] = size(Y);

% fill in missing inputs with default values

if ~is_entered(w), w = ones(m,1);  end
if ~is_entered(varY), varY = ones(m,1); end
if is_entered(bcon), dobtwn = 1; end

% --------------------------------------
% * Weights and computational steps
% --------------------------------------

W = diag(w);                    % Weight matrix

X = repmat(1,m,1);              % Design matrix - 1 column of all ones to calculate average
% and, separately, use bcon if that's entered

invxwx = inv(X'*W*X);
hat = invxwx * X'* W;         % hat matrix

% Between-observations, if entered
% ----------------------------------------
if dobtwn
    invxwx_diff = inv(bcon'*W*bcon);
    hatdiff = invxwx_diff*bcon'*W;
else
    hatdiff = [];
end

% --------------------------------------
% * Means and contrast
% --------------------------------------

Ymean = hat*Y;

% Output: weighted population mean
means.Ymean = Ymean;

% Between-observation contrast, if entered
% ----------------------------------------
if dobtwn

    Ydiff = hatdiff*Y;

    % for output; fits for each group; low then high
    grpfits = repmat(Ymean,2,1) + repmat(sort(unique(bcon)),1,n) .* repmat(Ydiff,2,1);

    means.Ydiff = Ydiff;
    means.grpmeans = grpfits;
end

if nargout == 1, return, end

% --------------------------------------
% * Residuals
% --------------------------------------

e = Y - repmat(Ymean,m,1);         % residuals

if dobtwn

    % fitted values depending on group
    fits = repmat(Ymean,m,1) + repmat(bcon,1,n) .* repmat(Ydiff,m,1);
    ediff = Y - fits;
end

% --------------------------------------
% * Degrees of freedom
% --------------------------------------

[dfe,dfediff] = get_dfe(m,n,X,hat,varY,dobtwn,hatdiff,bcon);

if ~domultivariate
    % ======================================
    %
    %
    % Univariate stats: MSE, t, and p-values
    %
    %
    % ======================================

    % --------------------------------------
    % * Mean squared error
    % --------------------------------------
    % Loop version of MSE: avoids out of memory errors for large voxel sets
    MSE = zeros(1,n); for i=1:n, MSE(i) = e(:,i)'*W*e(:,i); end, MSE = MSE/dfe;
    v = invxwx * MSE;       % variances for mean

    if dobtwn
        MSEdiff = zeros(1,n);
        for i=1:n, MSEdiff(i) = ediff(:,i)'*W*ediff(:,i); end, MSEdiff = MSEdiff/dfediff;
        vdiff = invxwx_diff * MSEdiff;       % variances for mean
    end

    % output
    stats.descrip1 = 'Univariate stats for test against zero:';
    stats.v = v;
    stats.v_descrip = 'V = ste^2; variance of mean estimate';
    stats.t = Ymean ./ sqrt(v);
    stats.p = 2 * ( 1 - tcdf(abs(stats.t),dfe) );
    stats.dfe = dfe;

    if dobtwn
        stats.descrip2 = 'Univariate stats for between-case contrast:';
        stats.bcon = bcon;
        stats.vdiff = vdiff;
        stats.tdiff = fits ./ sqrt(v);
        stats.pdiff = 2 * ( 1 - tcdf(abs(stats.tdiff),dfediff) );
        stats.dfediff = dfediff;
    end

else
    % ======================================
    %
    %
    % Multivariate stats: MSE, cov(Y), r(Y)
    % Useful for simulating t-values under dependence
    %
    % ======================================

    % --------------------------------------
    % * Mean squared error
    % --------------------------------------

    % additional output: covariance matrix for Ymean and zdiff across time
    % (columns)
    % and correlation matrix for Ymean and zdiff
    % used in Monte Carlo simulations for controlling false positives
    % across columns

    MSE = (e'*W*e)/dfe;             % Mean square error

    if dobtwn
        MSEdiff = (ediff'*W*ediff)/dfediff;
    end

    % --------------------------------------
    % * Estimated covariance and correlation
    %   Estimated between-subjects variance (v)
    % --------------------------------------

    xy = invxwx * MSE;           % Covariance matrix for Ymean;

    xy = 0.5*(xy+xy');              % Remove rounding error

    if dobtwn
        xydiff = inv(bcon'*W*bcon)*MSEdiff;           % Covariance matrix for Ymean;
        xydiff = 0.5*(xydiff+xydiff');
    end

    v = diag(xy);                   % Variance for Ymean

    if dobtwn
        vdiff = diag(xydiff);
    end

    r = xy./sqrt(v*v');             % Correlation matrix for Ymean

    if dobtwn
        rdiff = xydiff./sqrt(vdiff*vdiff');         % Correlation matrix for Ymean
    end

    stats.descrip1 = 'Multivariate stats for test against zero:';
    stats.r = r;
    stats.v = v;
    stats.xy = xy;

    stats.dfe = dfe;
    stats.t = Ymean ./ sqrt(v');
    stats.p = 2 * ( 1 - tcdf(abs(stats.t),dfe) );

    if dobtwn
        stats.descrip2 = 'Univariate stats for between-case contrast:';
        stats.rdiff = rdiff;
        stats.vdiff = vdiff;
        stats.xydiff = xydiff;

        stats.dfediff = dfediff;

        stats.tdiff = Ydiff ./ sqrt(vdiff');
        stats.p = 2 * ( 1 - tcdf(abs(stats.tdiff),dfediff) );
    end
end

return

function [dfe,dfediff] = get_dfe(m,n,X,hat,varY,dobtwn,hatdiff,bcon)

dfediff = [];

% Set up residual-forming matrix
% --------------------------------------
dfe_v = zeros(n,1);
R = eye(m) - X*hat;    % residual inducing matrix

% contrast, if entered
if dobtwn
    dfe_vdiff = zeros(n,1);
    Rdiff = eye(m) - bcon * hatdiff;
end

% Calculate effective degrees of freedom
% --------------------------------------

have_unique_vars = size(varY,2) == n;

if ~have_unique_vars

    % Only one (pooled?) vector of variance estimates
    % --------------------------------------
    V = diag(varY(:,1));
    dfe = (trace(R*V)^2)/trace(R*V*R*V);       % Satherwaite approximation
    if dobtwn, dfediff = (trace(Rdiff*V)^2)/trace(Rdiff*V*Rdiff*V); end
else
    % Variance estimates for each data vector
    % --------------------------------------
    for i=1:n,

        % make diagonal matrix of variances
        V = diag(varY(:,i));

        dfe_v(i) = (trace(R*V)^2)/trace(R*V*R*V);       % Satherwaite approximation

        if dobtwn
            dfe_vdiff(i) = (trace(Rdiff*V)^2)/trace(Rdiff*V*Rdiff*V);
        end

    end

    dfe = mean(dfe_v);               % Calculate average df over all columns (pool over data vectors)
    if dobtwn
        dfediff = mean(dfe_vdiff);
    end

end

return

function bool = is_entered(x)

bool = exist('x','var') && ~isempty(x);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [string,p_st,p_ed]=position(str_long,str_ID)
% Retrieve the value of an attribute (`str_ID`) in the line (`str_long`) of the
% XML file and identify the starting and ending indices of the value in the
% attribute line of the XML file.
%
% USAGE:
%
%    [string, p_st, p_ed] = position(str_long, str_ID)
%
% INPUTS:
%    str_long:       A string of the line of the XML file
%    str_ID:         The name of the attribute
%
% OUTPUTS:
%    string:         The value of the attribute
%    p_st:           The starting index of the value
%    p_ed:           The ending index of the value
%
% .. Author: - Longfei Mao Oct/2014

        ind_pos=strfind(str_long,str_ID); % name='metaid';
        l=length(str_ID)+2;
        p_st=ind_pos{1}+l
       % end_rem=strfind(str_long{1}(ind_pos{1}:end),'"');
        end_rem=strfind(str_long{1}(p_st:end),'"');
        % disp(str_long{1}(ind_pos{1}:end));
        disp(str_long{1}(p_st:end));

%         if size(end_rem)
%         end_rem{1}=[];
%         end

        p_ed=end_rem(1)+p_st-2;
        disp(p_ed);

       string=str_long{1}(p_st:p_ed);

end

author: opencobra

repo: cobratoolbox


text:
% Codes written by Soheil Bahrampour
% Feb 14, 2014

% To perform online supervised task driven dictionary Learning on the train
% XArr with lable trls. The classifier is the multiclass qudratic classifier.

% Inputs:
%       XArr onsisting of train samples where different modalities are concatinated.
%       trls is the train lables
%       n is a vector consiting of feature dimension for each modality
%       d is the number of columns in dictionary
%       opts contains the parameters for multi-task optimization
%       nu is the regularization parameter for the classifier
%       ro is the constant for computing learnning rate
%       InitW consists of parameters for initializing W

% Output:
%       D is the Learned Dictionary which is cell array with same size as X
%           consiting of n*d dictionaries learned from different sensors
%       modelQuad is a Quadratic classifier consisting of the linear coeficients W
%       and bias term b

%%% sparce codes from different modalities are averaged to form the final
%%% feature vector which will be used for classification

function [D, modelQuad] = OnlineSupTaskDrivDicLeaDecFusJointQuadC(XArr, Y, n, d, opts, nu, ro, varargin)

iter = opts.iterSupDic;
rho = opts.rho; % regularization term for ADMM
computeCost = opts.computeCost; % flag for computing and ploting cost
batchSize = opts.batchSize;
intercept = opts.intercept;

sum_n = sum(n);
S = length(n); % number of sensors
N = size(XArr,2); % number of train samples
number_classes = size(Y,1);

if nargin > 7 % initial D is provided
    D  = varargin{1};
else
    % Initialize D using (randomy taken) train samples + at least NumPos samples the positive lables
    D = zeros(sum_n,d);
    permut = randperm(N); %randomly shuffle data
    tempIndex = find(trls==lable);
    NumPos = length(tempIndex); % May need to be tuned for the application
    permut2 = randperm(length(tempIndex));
    for s = 1: S
        if d > N
            error('Number of dictionary columns for a given class should be smaller than or equal to the number of train samples in that class');
        end
        D(:,1:d-NumPos) = XArr(:,permut(1:d-NumPos)); %initialize to first d train samples
        D(:,d-NumPos+1:end) = XArr(:, tempIndex(permut2(1,1:NumPos))); % at least one sample from the  positive class
    end
end

if nargin > 8 % % initial model parameters are provided
   modelQuad =  varargin{2};
   W = zeros(d, number_classes,S); % initialize W to small random numbers rather than setting all to zeros
   b = zeros(S, number_classes);
   for s=1:S
       W(:,:,s) = modelQuad{1,s}.W;
       b(s,:) = modelQuad{1,s}.b;
   end
else
    % intialize W using unsisupervised learning
    W = 0.01*randn(d, number_classes,S); % initialize W to small random numbers rather than setting all to zeros
    b = 0.01*zeros(S, number_classes);
end

% learning rate of the stochastic gradient descent algorithm
t0 = floor(N/batchSize)*iter/10; % after finalizing N % for setting the learning rate according to the task-driven dic learning paper: we set t0=T/10 where T is total number of updates
permut = randperm(N);
XArr = XArr(:,permut);
Y = Y(:,permut);

% optimization
if computeCost
    costStep = floor(N/batchSize); % Compute cost every costStep over the last costStep batch of train samples. For each batch of train samples, the cost will be computed before updating the dic using that trian samples
    costTemp = zeros(costStep,1); % to store cost over lasr costStep samples
    costTempCount = 0; % to count how many train sample are passed
    cost = zeros(iter*N,1); %cost value at each iteration
    costIter = 0;
end

step = 0;
BetaTemp2 = zeros(d,batchSize,S);
Do = zeros(sum_n, S*d);
Atemp2 = zeros(d, S, N);
gradD = zeros(sum_n,d);
L = zeros(d*S, d);  % for ADMM
U = zeros(d*S, d);
gradW = zeros(d, number_classes, S);
gradb = zeros(S, number_classes);
for iteration = 1: iter % number of iterations over whole training samples
    for t = 1: batchSize: N-rem(N,batchSize) % to loop over all train samples
        step = step + 1;
        % find sparse code A
        Atemp = zeros(S*d,batchSize); % concatinate sparse vector of differen modals
        temp = 1;
        for s = 1:S
            [L((s-1)*d+1:s*d,:), U((s-1)*d+1:s*d,:)] = factor(D(temp:temp+n(s,1)-1,:), rho);  % cash  L U factroziation for solving ADMM for whole batch
            Do(temp: temp+n(s,1)-1, s:S:S*d) = D(temp:temp+n(s,1)-1, :);
            temp = temp+n(s,1);
        end
        DoTDoAll = Do'*Do; % computate Do'Do for all columns
        for j= 1: batchSize  % if not running in paralle
            alpha = JointADMMEigenMex(D, XArr(:,j+t-1), n, opts.lambda, rho, L, U, opts.iterADMM);
            Atemp(:,j) = alpha(:);
            Atemp2(:,:,j+t-1) = alpha;
            act = zeros(d, 1); % to maintain active rows
            Gamma = zeros(d*S,d*S);
            temp = 1;
            temp_norm = sqrt(sum(alpha.^2,2));
            for row = 1:d
                if  temp_norm(row,1) > 10^(-8)
                    act(row) = 1;
                    Gamma(temp:temp+S-1, temp:temp+S-1) = (eye(S)-alpha(row,:)'*alpha(row,:)/temp_norm(row,1)^2)/temp_norm(row,1);
                    temp = temp + S;
                end
            end
            Gamma = Gamma(1:temp-1,1:temp-1);
            num_act = sum(act); % number of active rows
            act = logical(act);
            %%%
            %  flag = isAlphaValid(D, XArr(:,j+t-1), n, alpha, opts.lambda, act); % only to check optimality condition for alpha, should be commented after debug.
            %%%
            Beta = zeros(d*S,1);
            W_act = W(act,:,:);
            gradLs_actAlphaVec = zeros(num_act,S);
            for s= 1:S
                gradLs_actAlphaVec(:,s) = W_act(:,:,s)*(W_act(:,:,s)'*alpha(act,s) + b(s,:)' - Y(:,j+t-1));
            end
            gradLs_actAt =  gradLs_actAlphaVec';
            gradLs_actAt_vec = gradLs_actAt(:);
            acts = act(:, ones(S,1))'; %repmat(act', S, 1); repmat is slow
            DoTDo = DoTDoAll(acts(:),acts(:));  % compute Do'Do for whole columns and then select relevant columns
            Beta(acts(:),1) = (DoTDo + opts.lambda*Gamma + opts.lambda2*eye(num_act*S))\gradLs_actAt_vec;   %use this iff ill-conditioned (lambda2 not equal zero)
            Beta(acts(:),1) = (DoTDo + opts.lambda*Gamma)\gradLs_actAt_vec;
            BetaTemp2(:,j,:) = reshape(Beta,S,d)';
%             temp = 1; % replaceed with a vectoried version for efficiency
%             for s= 1:S
%                 gradD(temp:temp+n(s,1)-1,:) = -(D(temp:temp+n(s,1)-1,:)*Beta(s:S:S*d,1))*alpha(:, s)' + (XArr(temp:temp+n(s,1)-1,j+t-1) - D(temp:temp+n(s,1)-1,:)*alpha(:, s))*Beta(s:S:S*d,1)' + gradD(temp:temp+n(s,1)-1,:);
%                 temp = temp+n(s,1);
%             end
%             temp2 = y(1,j+t-1)*gradLogReg(y(1,j+t-1)*(Atemp(:,j)'*W+b));
%             gradW = gradW + Atemp(:,j)*temp2;
%             if intercept
%                 gradb = gradb + temp2;
%             end
        end
        temp = 1;
        for s= 1:S
            temp3 = BetaTemp2(:,:,s)*squeeze(Atemp2(:, s, t:t+batchSize-1))';
            gradD(temp:temp+n(s,1)-1,:) = - D(temp:temp+n(s,1)-1,:)*temp3/batchSize + (XArr(temp:temp+n(s,1)-1,t:t+batchSize-1)*BetaTemp2(:,:,s)' - D(temp:temp+n(s,1)-1,:)*temp3')/batchSize;
%             gradD(temp:temp+n(s,1)-1,:) = -(D(temp:temp+n(s,1)-1,:)*BetaTemp2(:,:,s))*squeeze(Atemp2(:, s, t:t+batchSize-1))'/batchSize + (XArr(temp:temp+n(s,1)-1,t:t+batchSize-1) - D(temp:temp+n(s,1)-1,:)*squeeze(Atemp2(:, s, t:t+batchSize-1)))*BetaTemp2(:,:,s)'/batchSize;
            temp = temp+n(s,1);
        end

        temp = 1;
        for s= 1:S
            temp2 = W(:,:,s)'*Atemp(temp:temp+d-1,:) + repmat(b(s,:)', 1, batchSize) - Y(:,t:t+batchSize-1);
            gradW(:,:,s)  = Atemp(temp:temp+d-1,:)*temp2'/batchSize + nu*W(:,:,s);
            if intercept
                gradb(s,:) = sum(temp2,2)'/batchSize;
            end
            temp = temp + d;
        end

%          error = checkGradientW(W,b,gradW, Atemp, Y(:,t), nu); % set batch size to 1 for checking
 %         error = checkGradientD(D,gradD, W,b, XArr(:,t), Y(:,t), nu, n, S, d, opts, rho, opts.iterADMM); % set batch size to 1 for checking

        % compute cost
        if computeCost
            if costTempCount == costStep
                costIter = costIter + 1;
                cost(costIter,1) = mean(costTemp);
                costTempCount = 0;
                costTemp = zeros(costStep,1);
            end
            costTempCount = costTempCount + 1;
             temp = 1;
             for s = 1 : S
                costTemp(costTempCount,1) = costTemp(costTempCount,1) + 0.5*sum(sum((Y(:,t:t+batchSize-1) - W(:,:,s)'*Atemp(temp:temp+d-1,:) - repmat(b(s,:)', 1, batchSize)).^2));
                temp = temp + d;
             end
             costTemp(costTempCount,1)= costTemp(costTempCount,1)/batchSize + nu/2*sum(sum(sum(W.^2))); % note that b is not regularized
        end

        % update W, b, D
        learnRate = min(ro, ro*t0/step);
        W = W - learnRate*gradW;
        if intercept
            b = b - learnRate*gradb;
        end
        D = D - learnRate*gradD;
        temp = 1;
        for s=1:S
            D(temp:temp+n(s,1)-1,:) = projectionDic(D(temp:temp+n(s,1)-1,:));
            temp = temp + n(s,1);
        end
    end
end
if computeCost
    cost = cost(1:costIter,1); % remove extra elements
    figure;plot(cost);
end

for s= 1:S
    modelQuad{1,s}.W = W(:,:,s);
    modelQuad{1,s}.b = b(s,:);
end
end

% function DoTDo = computeDotDo(D, S, n, sum_n, num_act, act)
% Do = zeros(sum_n, S*num_act);
% ind = 1;
% for s= 1:S
%     Do(ind: ind+n(s,1)-1, s:S:S*num_act) = D(ind:ind+n(s,1)-1, act);
%     ind = ind+n(s,1);
% end
% DoTDo = Do'*Do;
% end
%
% function DoTDo = computeDotDo2(D, S, n, sum_n, d) % compute Do'Do for whole columns, not just active and save that
% Do = zeros(sum_n, S*d);
% ind = 1;
% for s= 1:S
%     Do(ind: ind+n(s,1)-1, s:S:S*d) = D(ind:ind+n(s,1)-1, :);
%     ind = ind+n(s,1);
% end
% DoTDo = Do'*Do;
% end

%%%% only with batchSize equal to 1
function error = checkGradientW(W, b, gradW, Atemp, y, nu)
epsil = 0.0001;
error = zeros(size(W,1), 1);
S = size(b,1);
d = size(W,1);
for i=1:size(W,1)
    W2 = W;
    W2(i,1) = W2(i,1) + epsil;
    temp = 1;
    costp = 0;
    for s = 1 : S
        costp = costp + 0.5*sum(sum((y - W2(:,:,s)'*Atemp(temp:temp+d-1,:) - b(s,:)').^2));
        temp = temp + d;
    end
    costp = costp + nu/2*sum(sum(sum(W2.^2))); % Note that b is not regularized
    W2(i,1) = W2(i,1) - 2*epsil;
    temp = 1;
    costn = 0;
    for s = 1 : S
        costn = costn + 0.5*sum(sum((y - W2(:,:,s)'*Atemp(temp:temp+d-1,:) - b(s,:)').^2));
        temp = temp + d;
    end
    costn = costn + nu/2*sum(sum(sum(W2.^2))); % Note that b is not regularized
    error(i,1)=(costp-costn)/2/epsil-gradW(i,1); % this should be close to zero
end
end

%%%% only with batchSize equal to 1
function error = checkGradientD(D,gradD, W, b, X, y, nu, n, S, d, opts, rho, iterADMM)
epsil = 0.0001;
L = zeros(d*S, d);
U = zeros(d*S, d);
error = zeros(d*S, d);

for s= 1:S%S
    temp = 1;
    for s1 = 1:S
        [L((s1-1)*d+1:s1*d,:), U((s1-1)*d+1:s1*d,:)] = factor(D(temp:temp+n(s1,1)-1,:), rho);  % cash  L U factroziation for solving ADMM for whole batch
        temp = temp + n(s1,1);
    end
    temp = sum(n(1:s-1))+1;
    for j =1 : 20%size(D,2)
        for i=1:10%n(s,1)
            D2 = D;
            D2(temp+i-1,j) = D2(temp+i-1,j) + epsil;
            [L((s-1)*d+1:s*d,:), U((s-1)*d+1:s*d,:)] = factor(D2(temp:temp+n(s,1)-1,:), rho);
            alpha = JointADMMEigenMex(D2, X, n, opts.lambda, rho, L, U, iterADMM);
            costp = 0;
            for s1=1:S
                costp = costp + 0.5*sum(sum((y - W(:,:,s1)'*alpha(:,s1) - b(s1,:)').^2));
            end
            costp = costp + nu/2*sum(sum(sum(W.^2))); % Note that b is not regularized
            D2(temp+i-1,j) = D2(temp+i-1,j) - 2*epsil;
            [L((s-1)*d+1:s*d,:), U((s-1)*d+1:s*d,:)] = factor(D2(temp:temp+n(s,1)-1,:), rho);
            alpha = JointADMMEigenMex(D2, X, n, opts.lambda, rho, L, U, iterADMM);
            costn = 0;
            for s1=1:S
                costn = costn + 0.5*sum(sum((y - W(:,:,s1)'*alpha(:,s1) - b(s1,:)').^2));
            end
            costn = costn + nu/2*sum(sum(sum(W.^2))); % Note that b is not regularized
            error(temp+i-1,j)=(costp-costn)/2/epsil-gradD(temp+i-1,j); % this should be close to zero
        end
    end
end
end

function flag = isAlphaValid(D, X, n, alpha, lambda, act) % check optimality condition for joint sparsity given a candidate solution
flag = 0;
checkMat = zeros(size(alpha));
S = size(alpha,2);
temp = 1;
for s=1:S
    checkMat(:,s) = D(temp:temp+n(s,1)-1,:)'*(X(temp:temp+n(s,1)-1,:)-D(temp:temp+n(s,1)-1,:)*alpha(:,s)) ;
    temp = temp + n(s,1);
end
RHS = lambda*alpha(act,:)./repmat(sqrt(sum(alpha(act,:).^2,2)), 1,S);
if norm(checkMat(act,:)-RHS)< 1e-10 % if condition is valid for active set (approximately)
    if isempty(find(sqrt(sum(alpha(~act,:).^2,2)) > lambda, 1))
        flag = 1;
    end
end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%{
load('dataset/trafficdb/traffic_patches.mat');
[M,m,n,p] = convert_video3d_to_2d(im2double(imgdb{100}));
out = run_algorithm('MC', 'MC_logdet', M, [])
show_results(M.*out.Omega,out.L,out.S,out.O,p,m,n);
%}

mu = 0.1; % 0.1; 0.01;
kappa = 1.2; % 1.2; 1.0;
toler = 1e-10;
maxiter = 1000;
L = MC_LogDet_v3(M,Omega,mu,kappa,toler,maxiter);
S = M - L;

author: andrewssobral

repo: lrslibrary


text:
% Copyright (C) 2000 Paul Kienzle  <pkienzle@users.sf.net>
% Copyright (C) 2007 Peter L. Soendergaard
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

% -*- texinfo -*-
% @deftypefn {Function File} {@var{h} =} hilbert (@var{f}, @var{N}, @var{dim})
% Analytic extension of real valued signal.
%
% @code{@var{h} = hilbert (@var{f})} computes the extension of the real
% valued signal @var{f} to an analytic signal. If @var{f} is a matrix,
% the transformation is applied to each column. For N-D arrays,
% the transformation is applied to the first non-singleton dimension.
%
% @code{real (@var{h})} contains the original signal @var{f}.
% @code{imag (@var{h})} contains the Hilbert transform of @var{f}.
%
% @code{hilbert (@var{f}, @var{N})} does the same using a length @var{N}
% Hilbert transform. The result will also have length @var{N}.
%
% @code{hilbert (@var{f}, [], @var{dim})} or
% @code{hilbert (@var{f}, @var{N}, @var{dim})} does the same along
% dimension @var{dim}.
% @end deftypefn

function f = oc_hilbert(f, N, dim)

% ------ PRE: initialization and dimension shifting ---------

if (nargin<1 || nargin>3)
    error('Invalid call');
end
if (nargin < 3)
    dim = [];
end
if (nargin < 2)
    N = [];
end

if ~isreal(f)
    warning ('HILBERT: ignoring imaginary part of signal');
    f = real (f);
end

D=ndims(f);

% Dummy assignment.
order=1;

if isempty(dim)
    dim=1;

    if sum(size(f)>1)==1
        % We have a vector, find the dimension where it lives.
        dim=find(size(f)>1);
    end

else
    if (numel(dim)~=1 || ~isnumeric(dim))
        error('HILBERT: dim must be a scalar.');
    end
    if rem(dim,1)~=0
        error('HILBERT: dim must be an integer.');
    end
    if (dim<1) || (dim>D)
        error('HILBERT: dim must be in the range from 1 to %d.',D);
    end

end

if (numel(N)>1 || ~isnumeric(N))
    error('N must be a scalar.');
elseif (~isempty(N) && rem(N,1)~=0)
    error('N must be an integer.');
end

if dim>1
    order=[dim, 1:dim-1,dim+1:D];

    % Put the desired dimension first.
    f=permute(f,order);

end

Ls=size(f,1);

% If N is empty it is set to be the length of the transform.
if isempty(N)
    N=Ls;
end

% Remember the exact size for later and modify it for the new length
permutedsize=size(f);
permutedsize(1)=N;

% Reshape f to a matrix.
f=reshape(f,size(f,1),numel(f)/size(f,1));
W=size(f,2);

if ~isempty(N)
    f = oc_postpad(f,N);
end

% ------- actual computation -----------------
if N>2
    f=fft(f);

    if rem(N,2)==0
        f=[f(1,:);
            2*f(2:N/2,:);
            f(N/2+1,:);
            zeros(N/2-1,W)];
    else
        f=[f(1,:);
            2*f(2:(N+1)/2,:);
            zeros((N-1)/2,W)];
    end

    f=ifft(f);
end

% ------- POST: Restoration of dimensions ------------

% Restore the original, permuted shape.
f=reshape(f,permutedsize);

if dim>1
    % Undo the permutation.
    f=ipermute(f,order);
end

end

%!demo
%! % notice that the imaginary signal is phase-shifted 90 degrees
%! t=linspace(0,10,256);
%! z = hilbert(sin(2*pi*0.5*t));
%! grid on; plot(t,real(z),';real;',t,imag(z),';imag;');

%!demo
%! % the magnitude of the hilbert transform eliminates the carrier
%! t=linspace(0,10,1024);
%! x=5*cos(0.2*t).*sin(100*t);
%! grid on; plot(t,x,'g;z;',t,abs(hilbert(x)),'b;|hilbert(z)|;');

author: brainstorm-tools


text:
%% LP1
clc
%Objective & Constraints
f = -[6 5]';
A = ([1,4; 6,4; 2, -5]);
b = [16;28;6];
lb = [0;0]; ub = [10;10];
%Setup Options
opts = optiset('solver','clp');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'bounds',lb,ub,'options',opts)
[x,fval,exitflag,info] = solve(Opt)
%Plot
plot(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% LP1 ROW
rl = -Inf(3,1);
ru = b;
Opt = opti('grad',f,'lin',A,rl,ru,'bounds',lb,ub,'options',opts)
[x,fval,exitflag,info] = solve(Opt)
%Plot
plot(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% Raw Test
opts.debug = 1;
opts.display = 2;
opts.maxiter = 5;
% rl = zeros(3,1);
n = length(f);
H = spalloc(n,n,0);
Aeq = spalloc(0,n,0);
beq = zeros(0,1);

[x,fv,ef,iter,l] = ooqp([],f,sparse(A),rl,ru,[],[],lb,ub,opts)

%% LP2
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
Aeq = [1 1 1]; beq = 40;
lb = [0 0 0]'; ub = [40 inf inf]';
%Setup Options
opts = optiset('solver','ooqp');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'options',opts)
[x,fval,exitflag,info] = solve(Opt)

%% LP2 ROW
Ar = [A;Aeq];
rl = [-Inf(2,1); beq];
ru = [b;beq];
Opt = opti('grad',f,'lin',Ar,rl,ru,'bounds',lb,ub,'options',opts)
[x,fval,exitflag,info] = solve(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% LP2 ROW + EQ (OOQP FORMAT) (NOT SUPPORTED CURRENTLY)
% rl = -Inf(2,1);
% ru = b;
% prob = optiprob('grad',f,'lin',A,rl,ru,'eq',Aeq,beq,'bounds',lb,ub);
% Opt = opti(prob,optiset('solver','clp'))
% [x,fval,exitflag,info] = solve(Opt)
% %Check Solution
% [ok,msg] = checkSol(Opt)

%% LP3
clc
%Objective & Constraints
f = [8,1]';
A = [-1,-2;1,-4;3,-1;1,5;-1,1;-1,0;0,-1];
b = [-4,2,21,39,3,0,0]';
%Setup Options
opts = optiset('solver','clp');
%Build & Solve
x0 = [9,6]';
Opt = opti('grad',f,'ineq',A,b,'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0)
plot(Opt,10)

%% LP3 ROW
rl = -Inf(7,1);
ru = b;
Opt = opti('grad',f,'lin',A,rl,ru,'options',opts)
[x,fval,exitflag,info] = solve(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% MILP1
clc
%Objective & Constraints
f = -[6 5]';
A = [1,4; 6,4; 2, -5];
b = [16;28;6];
lb = [0;0]; ub=[10;10];
%Setup Options
opts = optiset('solver','glpk');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'bounds',lb,ub,'int','II','options',opts)
[x,fval,exitflag,info] = solve(Opt)
plot(Opt)

%% MILP1 ROW
rl = -Inf(3,1);
ru = b;
Opt = opti('grad',f,'lin',A,rl,ru,'bounds',lb,ub,'int','II','options',opts)
[x,fval,exitflag,info] = solve(Opt)
plot(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% MILP2
clc
%Objective & Constraints
f = -[1 2 3 1]';
A = [-1 1 1 10; 1 -3 1 0];
b = [20;30];
Aeq = [0 1 0 -3.5];
beq = 0;
lb = [0 0 0 2]';
ub = [40 inf inf 3]';
int = 'CCCI';
%Setup Options
opts = optiset('solver','cbc');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'int',int,'options',opts)
[x,fval,exitflag,info] = solve(Opt)

%% MILP2 ROW
Ar = [A;Aeq];
rl = [-Inf(2,1); beq];
ru = [b;beq];
Opt = opti('grad',f,'lin',Ar,rl,ru,'bounds',lb,ub,'int',int,'options',opts)
[x,fval,exitflag,info] = solve(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% BILP1
clc
%Objective & Constraints
f = -[6 5]';
A = [-3,5; 6,4; 3, -5; -6, -4];
b = [6;9;1;3];
%Setup Options
opts = optiset('solver','cbc');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'int','BB','options',opts)
[x,fval,exitflag,info] = solve(Opt)
plot(Opt,3)

%% BILP1 ROW
rl = -Inf(4,1);
ru = b;
Opt = opti('grad',f,'lin',A,rl,ru,'int','BB','options',opts)
[x,fval,exitflag,info] = solve(Opt)
plot(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% BILP2
clc
%Objective & Constraints
f = -[9 5 6 4]';
A = [6 3 5 2; 0 0 1 1; -1 0 1 0; 0 -1 0 1];
b = [9; 1; 0; 0];
%Setup Options
opts = optiset('solver','cbc');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'int','BBBB','options',opts)
[x,fval,exitflag,info] = solve(Opt)

%% NLP1 Hock & Schittkowski #71
clc
%Objective & Gradient
obj = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
grad = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3));
              x(1)*x(4);
              x(1)*x(4) + 1;
              x(1)*sum(x(1:3)) ];
%Linear Constraints
lb = ones(4,1);
ub = 5*ones(4,1);
%Nonlinear Constraints
nlcon = @(x) [ prod(x);
               sum(x.^2)];
nljac = @(x) [ prod(x)./x';
                2*x' ];
nlrhs = [25 40]';
nle = [1 0]'; % (>=, ==)
%Setup Options
opts = optiset('solver','matlab','warnings','on','display','iter');
%Build & Solve
x0 = [1 5 5 1]';
Opt = opti('obj',obj,'grad',grad,'nlmix',nlcon,nlrhs,nle,'nljac',nljac,'bounds',lb,ub,'options',opts)
[x,fval,exitflag,info]= solve(Opt,x0)

%% NLP3 Hock & Schittkowski #51
clc
%Objective & Gradient
obj = @(x) (x(1) - x(2))^2 + (x(2) + x(3) - 2)^2 + (x(4) - 1)^2 + (x(5) - 1)^2;
grad = @(x) 2*[ x(1) - x(2);
                x(2) + x(3) - 2 - x(1) + x(2);
                x(2) + x(3) - 2;
                x(4) - 1;
                x(5) - 1 ];

%Nonlinear Constraints & Jacobian Structure
nlcon = @(x) [ x(1) + 3*x(2);
               x(3) + x(4) - 2*x(5)];
%                x(2) - x(5)];
nljac = @(x) sparse([ 1  3  0  0  0;
	                  0  0  1  1 -2]);
% 	                  0  1  0  0 -1]);
nljacstr = @() sparse([1 1 0 0 0;
                       0 0 1 1 1]);
%                        0 1 0 0 1]);

cl = [4;0];
cu = [4;0];
nlrhs = [4 0 1 15]';
nle = [0 0 1 -1]';

A = [0 1 0 0 -1];
rl = 1;
ru = 15;

%Setup Options
opts = optiset('solver','ipopt','display','iter','maxiter',1e4);
%Build & Solve
x0 = [ 2.5 0.5 2 -1 2.5 ];
Opt = opti('obj',obj,'grad',grad,'lin',A,rl,ru,'nl',nlcon,cl,cu,'nljac',nljac,'nljacstr',nljacstr,'options',opts)
[x,fval,exitflag,info]= solve(Opt,x0)

info.Lambda
info.Lambda.ineqlin
info.Lambda.ineqnonlin
info.Lambda.eqnonlin

%% Test Conversion Function

Opt.prob.nlcon(x0)
full(Opt.prob.nljac(x0))
full(Opt.prob.nljacstr())

%% SINGLE BOUNDS NO JAC
clc
nlcon = @(x) [ x(1) + 3*x(2);
               x(3) + x(4) - 2*x(5);
               x(2) - x(5)];

cl = [-Inf 1 2];
cu = [4 1 Inf];

x0 = [ 2.5 0.5 2 -1 1.5 ];

prob = optiprob('nl',nlcon,cl,cu,'x0',x0)
prob = nrow2mix(prob)

prob.nlrhs
prob.nle
prob.nlcon(x0)

%% SINGLE BOUNDS W JAC
clc
nlcon = @(x) [ x(1) + 3*x(2);
               x(3) + x(4) - 2*x(5);
               x(2) - x(5)];

nljac = @(x) sparse([ 1 3 0 0 0;
                      0 0 1 1 -2;
                      0 1 0 0 -1]);

cl = [-Inf 1 2];
cu = [4 1 Inf];

x0 = [ 2.5 0.5 2 -1 1.5 ];

prob = optiprob('nl',nlcon,cl,cu,'nljac',nljac,'x0',x0)
prob = nrow2mix(prob)

prob.nlrhs
prob.nle
prob.nlcon(x0)

%% DUAL BOUNDS W JAC
clc
nlcon = @(x) [ x(1) + 3*x(2);
               x(3) + x(4) - 2*x(5);
               x(2) - x(5);
               x(1) + 3*x(2);
               x(3) + x(4) - 5*x(5)];

nljac = @(x) sparse([ 1 3 0 0 0;
                      0 0 1 1 -2;
                      0 1 0 0 -1;
                      1 3 0 0 0;
                      0 0 1 1 -5;]);

cl = [-Inf 1 -Inf 2 5];
cu = [4 1 12 16 Inf];

x0 = [ 2.5 0.5 2 -1 1.5 ];

prob = optiprob('nl',nlcon,cl,cu,'nljac',nljac,'x0',x0)
prob = nrow2mix(prob)

prob.nlrhs
prob.nle
prob.nlcon(x0)
full(prob.nljac(x0))

%% DUAL BOUNDS W JAC & STR
clc
nlcon = @(x) [ x(1) + 3*x(2);
               x(3) + x(4) - 2*x(5);
               x(2) - x(5);
               x(1) + 3*x(2);
               x(3) + x(4) - 5*x(5)];

nljac = @(x) sparse([ 1 3 0 0 0;
                      0 0 1 1 -2;
                      0 1 0 0 -1;
                      1 3 0 0 0;
                      0 0 1 1 -5;]);

nljacstr = @() sparse([ 1 1 0 0 0;
                        0 0 1 1 1;
                        0 1 0 0 1;
                        1 1 0 0 0;
                        0 0 1 1 1]);

cl = [-Inf 1 -Inf 2 5];
cu = [4 1 12 16 Inf];

x0 = [ 2.5 0.5 2 -1 1.5 ];

prob = optiprob('nl',nlcon,cl,cu,'nljac',nljac,'nljacstr',nljacstr,'x0',x0)
prob = nrow2mix(prob)

prob.nlrhs
prob.nle
prob.nlcon(x0)
full(prob.nljac(x0))
full(prob.nljacstr())

%% Ali Constraint Linearity Index Problem
clc
obj=@(x) (x)^2;
A=1; ru=1; rl=-1; xtype='I'; x0=1;
opt=opti('obj',obj,'xtype',xtype,'x0',x0,'lin',A,rl,ru)
opt.prob.sizes
solve(opt,x0)

%% Identical form of above
clc
obj=@(x) (x)^2;
A=[1;-1];  b=[1;1];  xtype='I'; x0=1;
opt=opti('obj',obj,'xtype',xtype,'x0',x0,'ineq',A,b);
opt.prob.sizes
solve(opt,x0)

%% Nonlinear Constraints Double Sided Row, Linear
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
A = [-1 1; 1 1];
rl = [-3;5];
ru = [-1;5];
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ivars',1,'lin',A,rl,ru)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
plot(Opt,[],1)

%% Nonlinear Constraints Double Sided Ineq, Linear
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
A = [-1 1; 1 -1];
b = [-1;3];
Aeq = [1 1];
beq = 5;
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ivars',1,'ineq',A,b,'eq',Aeq,beq)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
plot(Opt,[],1)

%% Nonlinear Constraints Double Sided Row, NonLinear
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
nlcon = @(x) [-x(1) + x(2); x(1) + x(2)];
cl = [-3;5];
cu = [-1;5];
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ivars',1,'nl',nlcon,cl,cu)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
plot(Opt,[],1)

%% Nonlinear Constraints Double Sided Ineq, Nonlinear
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
nlcon = @(x) [-x(1) + x(2); x(1) - x(2); x(1) + x(2)];
nlrhs = [-1;3;5];
nle = [-1;-1;0];
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ivars',1,'nlmix',nlcon,nlrhs,nle)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
plot(Opt,[],1)

%% Nonlinear Constraints Double Sided Ineq, Nonlinear + Linear
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
A = [-1 1; 1 -1];
b = [0;2];
nlcon = @(x) [-x(1) + x(2); x(1) - x(2); x(1) + x(2)];
nlrhs = [-1;3;5];
nle = [-1;-1;0];
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ivars',1,'nlmix',nlcon,nlrhs,nle,'ineq',A,b)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
plot(Opt,[],1)

%% Nonlinear Constraints Double Sided Row, NonLinear + Linear
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
A = [-1 1];
rl = -2; ru = 0;
nlcon = @(x) [-x(1) + x(2); x(1) + x(2)];
cl = [-3;5];
cu = [-1;5];
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ivars',1,'nl',nlcon,cl,cu,'lin',A,rl,ru)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
plot(Opt,[],1)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function trx = SetLandmarkParameters_px(trx,arenatype,arenacenterx,arenacentery,...
  arenaradius)

switch lower(arenatype),
  case 'none',
    return;
  case 'circle',
    for i = 1:numel(trx),
      trx(i).arena.r = arenaradius;
      trx(i).arena.x = arenacenterx;
      trx(i).arena.y = arenacentery;
    end
  case 'rectangle',
end

author: kristinbranson

repo: JAABA