text:
function lik = lik_loggaussian(varargin)
%LIK_LOGGAUSSIAN  Create a right censored log-Gaussian likelihood structure
%
%  Description
%    LIK = LIK_LOGGAUSSIAN('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a likelihood structure for right censored log-Gaussian
%    survival model in which the named parameters have the
%    specified values. Any unspecified parameters are set to
%    default values.
%
%    LIK = LIK_LOGGAUSSIAN(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood structure with the named parameters
%    altered with the specified values.
%
%    Parameters for log-Gaussian likelihood [default]
%      sigma2       - variance [1]
%      sigma2_prior - prior for sigma2 [prior_logunif]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                  __ n
%      p(y|f, z) = || i=1 [ (2*pi*s^2)^(-(1-z_i)/2)*y_i^-(1-z_i)
%                            *exp(-1/(2*s^2)*(1-z_i)*(log(y_i) - f_i)^2)
%                            *(1-norm_cdf((log(y_i)-f_i)/s))^z_i ]
%
%
%    where s is the standard deviation of loggaussian distribution.
%    z is a vector of censoring indicators with z = 0 for uncensored event
%    and z = 1 for right censored event.
%
%    When using the log-Gaussian likelihood you need to give the vector z
%    as an extra parameter to each function that requires also y.
%    For example, you should call gpla_e as follows: gpla_e(w, gp,
%    x, y, 'z', z)
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%

% Copyright (c) 2012 Ville Tolvanen

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

  if isempty(lik)
    init=true;
    lik.type = 'Log-Gaussian';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Log-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
  % 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_loggaussian_pak;
    lik.fh.unpak = @lik_loggaussian_unpak;
    lik.fh.lp = @lik_loggaussian_lp;
    lik.fh.lpg = @lik_loggaussian_lpg;
    lik.fh.ll = @lik_loggaussian_ll;
    lik.fh.llg = @lik_loggaussian_llg;
    lik.fh.llg2 = @lik_loggaussian_llg2;
    lik.fh.llg3 = @lik_loggaussian_llg3;
    lik.fh.tiltedMoments = @lik_loggaussian_tiltedMoments;
    lik.fh.siteDeriv = @lik_loggaussian_siteDeriv;
    lik.fh.invlink = @lik_loggaussian_invlink;
    lik.fh.predy = @lik_loggaussian_predy;
    lik.fh.recappend = @lik_loggaussian_recappend;
    lik.fh.predcdf = @lik_loggaussian_predcdf;
  end

end

function [w,s] = lik_loggaussian_pak(lik)
%LIK_LOGGAUSSIAN_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_LOGGAUSSIAN_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)
%
%   See also
%   LIK_LOGGAUSSIAN_UNPAK, GP_PAK

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

function [lik, w] = lik_loggaussian_unpak(lik, w)
%LIK_LOGGAUSSIAN_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    [LIK, W] = LIK_LOGGAUSSIAN_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)
%
%   See also
%   LIK_LOGGAUSSIAN_PAK, GP_UNPAK

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

function lp = lik_loggaussian_lp(lik, varargin)
%LIK_LOGGAUSSIAN_LP  log(prior) of the likelihood parameters
%
%  Description
%    LP = LIK_LOGGAUSSIAN_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This subfunction
%    is needed when there are likelihood parameters.
%
%  See also
%    LIK_LOGGAUSSIAN_LLG, LIK_LOGGAUSSIAN_LLG3, LIK_LOGGAUSSIAN_LLG2, GPLA_E

% If prior for sigma2 parameter, add its contribution
  lp=0;
  if ~isempty(lik.p.sigma2)
    lp = lik.p.sigma2.fh.lp(lik.sigma2, lik.p.sigma2) +log(lik.sigma2);
  end

end

function lpg = lik_loggaussian_lpg(lik)
%LIK_LOGGAUSSIAN_LPG  d log(prior)/dth of the likelihood
%                parameters th
%
%  Description
%    E = LIK_LOGGAUSSIAN_LPG(LIK) takes a likelihood structure LIK and
%    returns d log(p(th))/dth, where th collects the parameters. This
%    subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_LOGGAUSSIAN_LLG, LIK_LOGGAUSSIAN_LLG3, LIK_LOGGAUSSIAN_LLG2, GPLA_G

  lpg=[];
  if ~isempty(lik.p.sigma2)
    % Evaluate the gprior with respect to sigma2
    ggs = lik.p.sigma2.fh.lpg(lik.sigma2, lik.p.sigma2);
    lpg = ggs(1).*lik.sigma2 + 1;
    if length(ggs) > 1
      lpg = [lpg ggs(2:end)];
    end
  end
end

function ll = lik_loggaussian_ll(lik, y, f, z)
%LIK_LOGGAUSSIAN_LL  Log likelihood
%
%  Description
%    LL = LIK_LOGGAUSSIAN_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, survival times Y, censoring indicators 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_LOGGAUSSIAN_LLG, LIK_LOGGAUSSIAN_LLG3, LIK_LOGGAUSSIAN_LLG2, GPLA_E

  if isempty(z)
    error(['lik_loggaussian -> lik_loggaussian_ll: missing z!    '...
           'loggaussian likelihood needs the censoring    '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loggaussian and gpla_e.               ']);
  end

  s2 = lik.sigma2;
  ll = sum(-(1-z)./2*log(2*pi*s2) - (1-z).*log(y) - (1-z)./(2*s2).*(log(y)-f).^2 ...
           + z.*log(1-norm_cdf((log(y)-f)./sqrt(s2))));

end

function llg = lik_loggaussian_llg(lik, y, f, param, z)
%LIK_LOGGAUSSIAN_LLG  Gradient of the log likelihood
%
%  Description
%    LLG = LIK_LOGGAUSSIAN_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, survival times Y, censoring indicators 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_LOGGAUSSIAN_LL, LIK_LOGGAUSSIAN_LLG2, LIK_LOGGAUSSIAN_LLG3, GPLA_E

  if isempty(z)
    error(['lik_loggaussian -> lik_loggaussian_llg: missing z!    '...
           'loggaussian likelihood needs the censoring    '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loggaussian and gpla_e.               ']);
  end

  s2 = lik.sigma2;
  r = log(y)-f;
  switch param
    case 'param'
      llg = sum(-(1-z)./(2.*s2) + (1-z).*r.^2./(2.*s2^2) + z./(1-norm_cdf(r/sqrt(s2))) ...
             .* (r./(sqrt(2.*pi).*2.*s2.^(3/2)).*exp(-1/(2.*s2).*r.^2)));
      % correction for the log transformation
      llg = llg.*lik.sigma2;
    case 'latent'
      llg = (1-z)./s2.*r + z./(1-norm_cdf(r/sqrt(s2))).*(1/sqrt(2*pi*s2) .* exp(-1/(2.*s2).*r.^2));
  end
end

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

  if isempty(z)
    error(['lik_loggaussian -> lik_loggaussian_llg2: missing z!   '...
           'loggaussian likelihood needs the censoring   '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loggaussian and gpla_e.               ']);
  end

  s2 = lik.sigma2;
  r = log(y)-f;
  switch param
    case 'param'

    case 'latent'
      llg2 = (z-1)./s2 + z.*(-exp(-r.^2/s2)./(2*pi*s2.*(1-norm_cdf(r/sqrt(s2))).^2) ...
              + r./(sqrt(2*pi).*s2^(3/2).*(1-norm_cdf(r/sqrt(s2)))).*exp(-r.^2./(2*s2)));
    case 'latent+param'
      llg2 = -(1-z)./s2^2.*(log(y)-f) + z.*(-r./(4*pi*s2^2.*(1-norm_cdf(r/sqrt(s2))).^2) ...
              .* exp(-r.^2./s2) + (-1 + r.^2/s2)./(1-norm_cdf(r/sqrt(s2))).*1./(sqrt(2*pi)*2*s2^(3/2)).*exp(-r.^2./(2*s2)));
      % correction due to the log transformation
      llg2 = llg2.*s2;
  end
end

function llg3 = lik_loggaussian_llg3(lik, y, f, param, z)
%LIK_LOGGAUSSIAN_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_LOGGAUSSIAN_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z and
%    latent values F and returns the third gradients of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. LLG3 is a vector with third gradients. This
%    subfunction is needed when using Laplace approximation for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LOGGAUSSIAN_LL, LIK_LOGGAUSSIAN_LLG, LIK_LOGGAUSSIAN_LLG2, GPLA_E, GPLA_G

  if isempty(z)
    error(['lik_loggaussian -> lik_loggaussian_llg3: missing z!   '...
           'loggaussian likelihood needs the censoring    '...
           'indicators as an extra input z. See, for         '...
           'example, lik_loggaussian and gpla_e.               ']);
  end

  s2 = lik.sigma2;
  r = log(y) - f;
  switch param
    case 'param'

    case 'latent'
      llg3 = 2.*z./(1-norm_cdf(r/sqrt(s2))).^3.*1./(2*pi*s2)^(3/2).*exp(-3/(2*s2)*r.^2) ...
              - z./(1-norm_cdf(r/sqrt(s2))).^2.*r./(pi*s2^2).*exp(-r.^2./s2) ...
              - z./(1-norm_cdf(r/sqrt(s2))).^2.*r./(2*pi*s2^2).*exp(-r.^2/s2) ...
              - z./(1-norm_cdf(r/sqrt(s2))).^1.*1./(s2^(3/2)*sqrt(2*pi)).*exp(-r.^2/(2*s2)) ...
              + z./(1-norm_cdf(r/sqrt(s2))).^1.*r.^2./(sqrt(2*pi*s2)*s2^2).*exp(-r.^2/(2*s2));
    case 'latent2+param'
      llg3 = (1-z)./s2^2 + z.*(1./(1-norm_cdf(r/sqrt(s2))).^3.*r./(sqrt(8*pi^3).*s2.^(5/2)).*exp(-3/(2.*s2).*r.^2) ...
              + 1./(1-norm_cdf(r./sqrt(s2))).^2.*1./(4.*pi.*s2^2).*exp(-r.^2./s2) ...
              - 1./(1-norm_cdf(r./sqrt(s2))).^2.*r.^2./(2*pi*s2^3).*exp(-r.^2./s2) ...
              + 1./(1-norm_cdf(r./sqrt(s2))).^2.*1./(4*pi*s2^2).*exp(-r.^2/s2) ...
              - 1./(1-norm_cdf(r./sqrt(s2))).^1.*r./(sqrt(2*pi)*2*s2^(5/2)).*exp(-r.^2/(2*s2)) ...
              - 1./(1-norm_cdf(r./sqrt(s2))).^2.*r.^2./(4*pi*s2^3).*exp(-r.^2/s2) ...
              - 1./(1-norm_cdf(r./sqrt(s2))).^1.*r./(sqrt(2*pi)*s2^(5/2)).*exp(-r.^2/(2*s2)) ...
              + 1./(1-norm_cdf(r./sqrt(s2))).^1.*r.^3./(sqrt(2*pi)*2*s2^(7/2)).*exp(-r.^2/(2*s2)));
      % correction due to the log transformation
      llg3 = llg3.*lik.sigma2;
  end
end

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

 if isempty(z)
   error(['lik_loggaussian -> lik_loggaussian_tiltedMoments: missing z!'...
          'loggaussian likelihood needs the censoring            '...
          'indicators as an extra input z. See, for                 '...
          'example, lik_loggaussian and gpep_e.                       ']);
 end

  yy = y(i1);
  yc = 1-z(i1);
  s2 = lik.sigma2;
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    % get a function handle of an unnormalized tilted distribution
    % (likelihood * cavity = Negative-binomial * Gaussian)
    % and useful integration limits
    [tf,minf,maxf]=init_loggaussian_norm(yy(i),myy_i(i),sigm2_i(i),yc(i),s2);

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1.e-10;
    [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
    sigm2hati1(i) = m_2 - m_1(i).^2;

    % If the second central moment is less than cavity variance
    % integrate more precisely. Theoretically for log-concave
    % likelihood should be sigm2hati1 < sigm2_i.
    if sigm2hati1(i) >= sigm2_i(i)
      ATOL = ATOL.^2;
      RTOL = RTOL.^2;
      [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
      sigm2hati1(i) = m_2 - m_1(i).^2;
      if sigm2hati1(i) >= sigm2_i(i)
        error('lik_loggaussian_tilted_moments: sigm2hati1 >= sigm2_i');
      end
    end
    logM_0(i) = log(m_0);
  end
end

function [g_i] = lik_loggaussian_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_LOGGAUSSIAN_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_LOGGAUSSIAN_SITEDERIV(LIK, Y, I, S2, MYY, Z) takes a
%    likelihood structure LIK, survival times 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

  if isempty(z)
    error(['lik_loggaussian -> lik_loggaussian_siteDeriv: missing z!'...
           'loggaussian likelihood needs the censoring        '...
           'indicators as an extra input z. See, for             '...
           'example, lik_loggaussian and gpla_e.                   ']);
  end

  yy = y(i1);
  yc = 1-z(i1);
  s2 = lik.sigma2;

  % get a function handle of an unnormalized tilted distribution
  % (likelihood * cavity = Log-Gaussian * Gaussian)
  % and useful integration limits
  [tf,minf,maxf]=init_loggaussian_norm(yy,myy_i,sigm2_i,yc,s2);
  % additionally get function handle for the derivative
  td = @deriv;

  % Integrate with quadgk
  [m_0, fhncnt] = quadgk(tf, minf, maxf);
  [g_i, fhncnt] = quadgk(@(f) td(f).*tf(f)./m_0, minf, maxf);
  g_i = g_i.*s2;

  function g = deriv(f)
    r=log(yy)-f;
    g = -yc./(2.*s2) + yc.*r.^2./(2.*s2^2) + (1-yc)./(1-norm_cdf(r/sqrt(s2))) ...
             .* (r./(sqrt(2.*pi).*2.*s2.^(3/2)).*exp(-1/(2.*s2).*r.^2));
  end
end

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

  if isempty(zt)
    error(['lik_loggaussian -> lik_loggaussian_predy: missing zt!'...
           'loggaussian likelihood needs the censoring    '...
           'indicators as an extra input zt. See, for         '...
           'example, lik_loggaussian and gpla_e.               ']);
  end

  yc = 1-zt;
  s2 = lik.sigma2;

  Ey=[];
  Vary=[];

  % Evaluate the posterior predictive densities of the given observations
  lpy = zeros(length(yt),1);
  for i1=1:length(yt)
    if abs(Ef(i1))>700
      lpy(i1) = NaN;
    else
      % get a function handle of the likelihood times posterior
      % (likelihood * posterior = Negative-binomial * Gaussian)
      % and useful integration limits
      [pdf,minf,maxf]=init_loggaussian_norm(...
        yt(i1),Ef(i1),Varf(i1),yc(i1),s2);
      % integrate over the f to get posterior predictive distribution
      lpy(i1) = log(quadgk(pdf, minf, maxf));
    end
  end
end

function [df,minf,maxf] = init_loggaussian_norm(yy,myy_i,sigm2_i,yc,s2)
%INIT_LOGGAUSSIAN_NORM
%
%  Description
%    Return function handle to a function evaluating
%    loggaussian * Gaussian which is used for evaluating
%    (likelihood * cavity) or (likelihood * posterior) Return
%    also useful limits for integration. This is private function
%    for lik_loggaussian. This subfunction is needed by subfunctions
%    tiltedMoments, siteDeriv and predy.
%
%  See also
%    LIK_LOGGAUSSIAN_TILTEDMOMENTS, LIK_LOGGAUSSIAN_SITEDERIV,
%    LIK_LOGGAUSSIAN_PREDY

% avoid repetitive evaluation of constant part
  ldconst = -yc./2.*log(2*pi*s2) -yc.*log(yy) ...
            - log(sigm2_i)/2 - log(2*pi)/2;

  % Create function handle for the function to be integrated
  df = @loggaussian_norm;
  % use log to avoid underflow, and derivates for faster search
  ld = @log_loggaussian_norm;
  ldg = @log_loggaussian_norm_g;
  ldg2 = @log_loggaussian_norm_g2;

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

  function integrand = loggaussian_norm(f)
  % loggaussian * Gaussian
    integrand = exp(ldconst ...
                    - yc./(2*s2).*(log(yy)-f).^2 + (1-yc).*log(1-norm_cdf((log(yy)-f)/sqrt(s2))) ...
                    -0.5*(f-myy_i).^2./sigm2_i);
  end

  function log_int = log_loggaussian_norm(f)
  % log(loggaussian * Gaussian)
  % log_loggaussian_norm is used to avoid underflow when searching
  % integration interval
    log_int = ldconst ...
              -yc./(2*s2).*(log(yy)-f).^2 + (1-yc).*log(1-norm_cdf((log(yy)-f)/sqrt(s2))) ...
              -0.5*(f-myy_i).^2./sigm2_i;
  end

  function g = log_loggaussian_norm_g(f)
  % d/df log(loggaussian * Gaussian)
  % derivative of log_loggaussian_norm
    g = yc./s2.*(log(yy)-f) + (1-yc)./(1-norm_cdf((log(yy)-f)/sqrt(s2))).*1/sqrt(2*pi*s2)*exp(-(log(yy)-f).^2./(2*s2))  ...
        + (myy_i - f)./sigm2_i;
  end

  function g2 = log_loggaussian_norm_g2(f)
  % d^2/df^2 log(loggaussian * Gaussian)
  % second derivate of log_loggaussian_norm
    g2 = -yc./s2 + (1-yc).*(-exp(-(log(yy)-f).^2/s2)./(2*pi*s2.*(1-norm_cdf((log(yy)-f)/sqrt(s2))).^2) ...
              + (log(yy)-f)./(sqrt(2*pi).*s2^(3/2).*(1-norm_cdf((log(yy)-f)/sqrt(s2)))).*exp(-(log(yy)-f).^2./(2*s2))) ...
              -1/sigm2_i;
  end

end

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

  s2 = lik.sigma2;

  % Evaluate the posterior predictive densities of the given observations
  cdf = zeros(length(yt),1);
  for i1=1:length(yt)
    % Get a function handle of the likelihood times posterior
    % (likelihood * posterior = log-Gaussian * Gaussian)
    % and useful integration limits.
    % yc=0 when evaluating predictive cdf
    [pdf,minf,maxf]=init_loggaussian_norm(...
      yt(i1),Ef(i1),Varf(i1),0,s2);
    % integrate over the f to get posterior predictive distribution
    cdf(i1) = 1-quadgk(pdf, minf, maxf);
  end

end

function p = lik_loggaussian_invlink(lik, f)
%LIK_LOGGAUSSIAN Returns values of inverse link function
%
%  Description
%    P = LIK_LOGGAUSSIAN_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values of inverse link function P.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_LOGGAUSSIAN_LL, LIK_LOGGAUSSIAN_PREDY

p = exp(f);
end

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

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

    % Initialize parameter
    reclik.sigma2 = [];

    % Set the function handles
    reclik.fh.pak = @lik_loggaussian_pak;
    reclik.fh.unpak = @lik_loggaussian_unpak;
    reclik.fh.lp = @lik_loggaussian_lp;
    reclik.fh.lpg = @lik_loggaussian_lpg;
    reclik.fh.ll = @lik_loggaussian_ll;
    reclik.fh.llg = @lik_loggaussian_llg;
    reclik.fh.llg2 = @lik_loggaussian_llg2;
    reclik.fh.llg3 = @lik_loggaussian_llg3;
    reclik.fh.tiltedMoments = @lik_loggaussian_tiltedMoments;
    reclik.fh.invlink = @lik_loggaussian_invlink;
    reclik.fh.predy = @lik_loggaussian_predy;
    reclik.fh.predcdf = @lik_loggaussian_predcdf;
    reclik.fh.recappend = @lik_loggaussian_recappend;
    reclik.p=[];
    reclik.p.sigma2=[];
    if ~isempty(ri.p.sigma2)
      reclik.p.sigma2 = ri.p.sigma2;
    end
  else
    % Append to the record
    reclik.sigma2(ri,:)=lik.sigma2;
    if ~isempty(lik.p)
      reclik.p.sigma2 = lik.p.sigma2.fh.recappend(reclik.p.sigma2, ri, lik.p.sigma2);
    end
  end
end


text:
%READ_BF_FILE Reads in a file of beamforming feedback logs.
%   This version uses the *C* version of read_bfee, compiled with
%   MATLAB's MEX utility.
%
% (c) 2008-2011 Daniel Halperin <dhalperi@cs.washington.edu>
%
function ret = read_bf_file(filename)
%% Input check
error(nargchk(1,1,nargin));

%% Open file
f = fopen(filename, 'rb');
if (f < 0)
    error('Couldn''t open file %s', filename);
    return;
end

status = fseek(f, 0, 'eof');
if status ~= 0
    [msg, errno] = ferror(f);
    error('Error %d seeking: %s', errno, msg);
    fclose(f);
    return;
end
len = ftell(f);

status = fseek(f, 0, 'bof');
if status ~= 0
    [msg, errno] = ferror(f);
    error('Error %d seeking: %s', errno, msg);
    fclose(f);
    return;
end

%% Initialize variables
ret = cell(ceil(len/95),1);     % Holds the return values - 1x1 CSI is 95 bytes big, so this should be upper bound
cur = 0;                        % Current offset into file
count = 0;                      % Number of records output
broken_perm = 0;                % Flag marking whether we've encountered a broken CSI yet
triangle = [1 3 6];             % What perm should sum to for 1,2,3 antennas

%% Process all entries in file
% Need 3 bytes -- 2 byte size field and 1 byte code
while cur < (len - 3)
    % Read size and code
    field_len = fread(f, 1, 'uint16', 0, 'ieee-be');
    code = fread(f,1);
    cur = cur+3;

    % If unhandled code, skip (seek over) the record and continue
    if (code == 187) % get beamforming or phy data
        bytes = fread(f, field_len-1, 'uint8=>uint8');
        cur = cur + field_len - 1;
        if (length(bytes) ~= field_len-1)
            fclose(f);
            return;
        end
    else % skip all other info
        fseek(f, field_len - 1, 'cof');
        cur = cur + field_len - 1;
        continue;
    end

    if (code == 187) %hex2dec('bb')) Beamforming matrix -- output a record
        count = count + 1;
        ret{count} = read_bfee(bytes);

        perm = ret{count}.perm;
        Nrx = ret{count}.Nrx;
        if Nrx == 1 % No permuting needed for only 1 antenna
            continue;
        end
        if sum(perm) ~= triangle(Nrx) % matrix does not contain default values
            if broken_perm == 0
                broken_perm = 1;
                fprintf('WARN ONCE: Found CSI (%s) with Nrx=%d and invalid perm=[%s]\n', filename, Nrx, int2str(perm));
            end
        else
            ret{count}.csi(:,perm(1:Nrx),:) = ret{count}.csi(:,1:Nrx,:);
        end
    end
end
ret = ret(1:count);

%% Close file
fclose(f);

end

author: linteresa

repo: WiAR

save_path: github-repos/MATLAB/linteresa-WiAR


text:
function test_issue1618

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_volumewrite ft_write_mri

% this function was related to github issue 1618

%%
mri      = [];
mri.dim  = [2 2 2];
mri.transform = eye(4);
mri.unit = 'mm';
mri.coordsys = 'ras';

mri.dat1 = randn(2,2,2);
mri.dat2 = reshape(1:8,[2 2 2]);
mri.dat3 = rand(2,2,2)>0.5;

datatype = {'logical', 'uint8', 'int8', 'int16', 'int32', 'single', 'double'};
filetype = {'analyze_old', 'nifti', 'nifti_img', 'nifti_spm', 'nifti_gz', 'analyze_spm', 'mgz', 'mgh', 'vmp', 'vmr'};
fname    = fullfile(tempdir, 'issue1618');

ok_table = true(numel(datatype), numel(filetype), 3);
for m = 1:numel(filetype)
  for k = 1:numel(datatype)
    try
      cfg = [];
      cfg.datatype  = datatype{k};
      cfg.filetype  = filetype{m};
      cfg.parameter = 'dat1';
      cfg.filename  = [fname, '_', datatype{k}, '_', filetype{m},'_', cfg.parameter];
      %if k==2, keyboard; end
      ft_volumewrite(cfg, mri);
    catch
      ok_table(k,m,1) = false;
    end
    fnames{k,m,1} = cfg.filename;

    try
      cfg = [];
      cfg.datatype  = datatype{k};
      cfg.filetype  = filetype{m};
      cfg.parameter = 'dat2';
      cfg.scaling   = 'no';
      cfg.filename = [fname, '_', datatype{k}, '_', filetype{m},'_', cfg.parameter];
      ft_volumewrite(cfg, mri);
    catch
      ok_table(k,m,2) = false;
    end
    fnames{k,m,2} = cfg.filename;

    try
      cfg = [];
      cfg.datatype  = datatype{k};
      cfg.filetype  = filetype{m};
      cfg.parameter = 'dat3';
      cfg.scaling   = 'no';
      cfg.filename = [fname, '_', datatype{k}, '_', filetype{m},'_', cfg.parameter];
      ft_volumewrite(cfg, mri);
    catch
      ok_table(k,m,3) = false;
    end
    fnames{k,m,3} = cfg.filename;
  end
end

data = cell(size(fnames));
for k = 1:size(fnames,1)
  for kk = 1:size(fnames,2)
    for kkk = 1:size(fnames,3)
      if ok_table(k,kk,kkk) && ~contains(fnames{k,kk,kkk},'vmr') && ~contains(fnames{k,kk,kkk},'vmp')
        % get the filename + extension but only if file has been written, and
        % if it's not a brainvoyager type, because no reading functionality
        % availablt
        d = dir([fnames{k,kk,kkk} '*']);
        if numel(d)>1
          % could be a hdr/img combo, take the first one
          d = d(1);
        end
        fname = fullfile(d.folder,d.name);
        data{k,kk,kkk} = ft_read_mri(fname);
      end
    end
  end
end

% FIXME this ideally should contain some check as to whether the values
% that have been written to disk make sense and/or are a scaled version of
% the original input data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% How do I avoid the slow Monte Carlo significance test in wtc?
% You can do that by simply specifying the MonteCarloCount to be zero.
% Example:

figure('color',[1 1 1])
t=(0:1:500)';
X=sin(t*2*pi/11)+randn(size(t))*.1;
Y=sin(t*2*pi/11+.4)+randn(size(t))*.1;
wtc([t X],[t Y],'mcc',0); %MCC:MonteCarloCount

%%
% Note that the significance contour can not be trusted with out running the
% Monte Carlo test.

author: grinsted

repo: wavelet-coherence


text:
% Author: Guosheng Lin (guosheng.lin@gmail.com)

function exp_info = my_net_tool(opts, imdb, net_config, input_exp_info)

net_run_opts=opts.net_run_opts;

exp_info=[];
assert(isa(net_config, 'ref_obj'));

work_infos=cell(0);

work_info_trn=gen_work_info_basic();
work_info_trn.ref.run_type='train';
work_info_trn.ref.split_name='train';
work_infos{end+1, 1}=work_info_trn;

work_info_eva_val=gen_work_info_basic();
work_info_eva_val.ref.run_type='evaluate';
work_info_eva_val.ref.split_name='val';
work_infos{end+1, 1}=work_info_eva_val;

for w_idx=1:length(work_infos)
    one_work_info=work_infos{w_idx};
    one_work_info.ref.run_trn=strcmp(one_work_info.ref.run_type, 'train');
    one_work_info.ref.run_eva=strcmp(one_work_info.ref.run_type, 'evaluate');
end

do_init_exp_info(opts, work_info_trn, work_info_eva_val, input_exp_info);
clear input_exp_info

work_info_trn.ref.task_info=imdb.ref.task_info_train;
work_info_eva_val.ref.task_info=imdb.ref.task_info_val;

start_epoch=work_info_trn.ref.current_epoch+1;

if opts.do_train

    my_clear_net_trn(net_config);

end

global_init_fn=net_run_opts.global_init_fn;
if ~isempty(global_init_fn)
    if opts.do_train
        global_init_fn(opts, imdb, work_info_trn, net_config);
    end

    if opts.do_eva_val
        global_init_fn(opts, imdb, work_info_eva_val, net_config);
    end
end

for epoch=start_epoch:opts.epoch_num

    if opts.do_train
        work_info_trn.ref.current_epoch=epoch;
        my_net_run_epoch(opts, imdb, work_info_trn, net_config);
    end

    if opts.do_eva_val
        work_info_eva_val.ref.current_epoch=epoch;
        if any(opts.eva_epoch_idxes==epoch)
            my_net_run_epoch(opts, imdb, work_info_eva_val, net_config);
        end
    end

  exp_info=[];
  exp_info.train=work_info_trn;
  exp_info.eva_val=work_info_eva_val;

  if ~isempty(net_run_opts.net_progress_disp_fn)
      net_run_opts.net_progress_disp_fn(opts, exp_info, imdb, net_config);
  end

  opts.epoch_model_cache_fn(opts, epoch, exp_info, net_config);

  my_diary_flush();

  if opts.do_train
      if work_info_trn.ref.epoch_valid_batch_count==0
          fprintf('\n\n===work_info_trn.ref.epoch_valid_batch_count==0, stop \n\n');
          break;
      end
  end

end

my_diary_flush();

end

function work_info=gen_work_info_basic()

work_info=[];
work_info.objective = [] ;
work_info.eva_results = cell(0, 1) ;
work_info.current_epoch=0;
work_info.valid_epoch_idxes=[];
work_info.total_iter_num=0;
work_info.tmp_cache=[];

make_ref_obj(work_info);

end

function do_init_exp_info(opts, work_info_trn, work_info_eva_val, input_exp_info)

    if ~isempty(input_exp_info)

          update_work_info_from_cache(work_info_trn, input_exp_info, 'train');
          update_work_info_from_cache(work_info_eva_val, input_exp_info, 'eva_val');

          current_epoch1=work_info_trn.ref.current_epoch;
          current_epoch2=work_info_eva_val.ref.current_epoch;
          current_epoch=max(current_epoch1, current_epoch2);

    else
          if isfield(opts, 'net_start_epoch') && ~isempty(opts.net_start_epoch)
            current_epoch=opts.work_info_missed_start_epoch-1;
          else
            current_epoch=0;
          end
    end

  work_info_trn.ref.current_epoch=current_epoch;
  work_info_eva_val.ref.current_epoch=current_epoch;

end

function update_work_info_from_cache(work_info, tmp_info, work_info_name)

if isfield(tmp_info, work_info_name)

      cached_work_info=tmp_info.(work_info_name);

      if isa(cached_work_info, 'ref_obj')
          cached_work_info=cached_work_info.ref;
      end

      work_info.ref=cached_work_info;
      work_info.ref.tmp_cache=[];

end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function C = mid(A)
%MID          Midpoint of long (only with error term)
%
%  C = mid(A)
%

% written  11/06/99     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_LONG_ERROR = getappdata(0,'INTLAB_LONG_ERROR');

  if ~INTLAB_LONG_ERROR
    error('mid called for long without error term')
  end

  C = A;
  C.error.mant = 0;
  C.error.exp = 0;

repo: OpenMAS


text:
function [pvec, pstruct] = tapas_ehgf_binary_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pvec    = NaN(1,length(ptrans));
pstruct = struct;

l = r.c_prc.n_levels;

pvec(1:l)         = ptrans(1:l);                           % mu_0
pstruct.mu_0      = pvec(1:l);
pvec(l+1:2*l)     = exp(ptrans(l+1:2*l));                  % sa_0
pstruct.sa_0      = pvec(l+1:2*l);
pvec(2*l+1:3*l)   = ptrans(2*l+1:3*l);                     % rho
pstruct.rho       = pvec(2*l+1:3*l);
pvec(3*l+1:4*l-1) = exp(ptrans(3*l+1:4*l-1));              % ka
pstruct.ka        = pvec(3*l+1:4*l-1);
pvec(4*l:5*l-1)   = ptrans(4*l:5*l-1);                     % om
pstruct.om        = pvec(4*l:5*l-1);

end

author: translationalneuromodeling

repo: tapas


text:
function varargout = aib(varargin)
% VL_AIB  Agglomerative Information Bottleneck
%   PARENTS = VL_AIB(PCX) runs Agglomerative Information Bottleneck (AIB)
%   on the class-feature co-occurrence matrix PCX and returns a vector
%   PARENTS representing the sequence of compressed VL_AIB alphabets.
%
%   PCX is the joint probability of the occurrence of the class label
%   C and the feature value X. PCX has one row for each class label
%   and one column for each feature value, non negative entires and
%   sums to one.  VL_AIB iteratively merges the pair of feature values
%   that decreases the mutual information I(X,C) the least. This
%   compresses the alphabet of the discrete random variable X in such
%   a way that the new variable is still informative about C.
%
%   Merge operations are represented by a binary tree. The nodes of
%   the tree correspond to the original feature values and any other
%   value obtained by merging.
%
%   The vector PARENTS represents the merge tree. The tree nodes are
%   numbered in breadth-first order, starting from the leaves. The
%   numbers associated to the tree leaves correspond to the original
%   feature values (so the first leaf has number one and correspond
%   to the first feature value).  In total there are 2*M-1 nodes,
%   where M is the number of feature values (the number of columns of
%   PCX). The internal nodes are numbered according to the order in
%   which VL_AIB generated them. It is therefore possible to recover from
%   the tree the state of the VL_AIB algorithm at each step (see VL_AIBCUT())
%   PARENTS is a UINT32 array with one element for each tree node storing the
%   index of the parent node. The root parent is conventionally set to 1.
%
%   Feature values with null probability (null columns of the PCX
%   matrix) are ignored by the VL_AIB algorithm and the corresponding
%   entries in the PARENTS vectors are set to zero.  Notice that this
%   causes the root of the tree to have index smaller of 2*M-1
%   (PARENTS has still 2*M-1 entries, but the last portion is
%   zero-padded).
%
%   Alternatively, the option ClusterNull can be used to assign the
%   null probability values to a special value. The result is similar
%   to pretending that the null probability nodes have indeed very
%   small probability, uniform across categories.
%
%   [PARENTS, COST] = VL_AIB(...) returns the values COST of the cost %
%   function being optimized by VL_AIB (i.e. the mutual information
%   I(X,C)). COST has M column. The first column is the initial value
%   of the cost function. The others correspond to the cost after each
%   of the M-1 merges. If less than M-1 merges are performed, the rest
%   of the vector is filled with NaNs.
%
%   VL_AIB accepts the following options:
%
%   Verbose::
%     If specified, increase verbosity level.
%
%   ClusterNull::
%     If specified, do not signal null nodes; instead cluster them.
%
%   See also: VL_AIBCUT(), VL_AIBHIST(), VL_AIBCUTHIST(), VL_AIBCUTPUSH(), VL_HELP().
[varargout{1:nargout}] = vl_aib(varargin{:});

author: jianxiongxiao

repo: ProfXkit


%% REW_boardingbridge
% %{
imfolder = 'images\REW_boardingbridge';
im_n = 2;
imfile = cell(im_n,1);
imfile{1} = [imfolder '\' 'boardingbridge_01.jpg'];
imfile{2} = [imfolder '\' 'boardingbridge_02.jpg'];

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = [1,2];

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 1, 'persp', 0, imfolder );
%}


text:
function tests = SE3Test
    tests = functiontests(localfunctions);
    clc
end

% we will assume that the primitives rotx,trotx, etc. all work

function constructor_test(tc)

    verifyClass(tc, SE3(), 'SE3');

    %% null
    tc.verifyEqual(SE3().double, eye(4,4));

    %% translation only
    t = [1 2 3];
    tc.verifyEqual(SE3(t).double, transl(t));
    tc.verifyEqual(SE3(t').double, transl(t));
    tc.verifyEqual(SE3(t(1), t(2), t(3)).double, transl(t));

    %% R,t
    R = roty(-pi/2);
    tc.verifyEqual(SE3(R,t).double, transl(t)*r2t(R));

    %% T
    T = rt2tr(R,t);
    tc.verifyEqual(SE3(T).double, T);

    %% copy constructor
    TT = SE3(T);
    tc.verifyEqual(SE3(TT).double, T);

    %% vectorised versions

    T1 = transl(1,2,3) * rpy2tr(0.2, 0.3, 0.4);
    T2 = transl(1,-2,3) * rpy2tr(0.2, -0.3, 0.4);

    TT = cat(3, T1, T2, T1, T2);

    tt = SE3(TT);
    tc.verifyEqual(length(tt), size(TT, 3) );
    tc.verifyEqual(tt.T, TT);

    t = rand(10,3);
    s = SE3(t);
    tc.verifyLength(s, 10);
    tc.verifyEqual(s(4).t, t(4,:)');

    %% symbolic version
    syms x y z

        verifyClass(tc, SE3(x, y, z), 'SE3');

end

function concat_test(tc)
    x = SE3();
    xx = [x x x x];

    tc.verifyClass(xx, 'SE3');
    tc.verifySize(xx, [1 4]);
end

function double_test(tc)

    a = SE3(1,2, 3);
    d = double(a);
    verifyInstanceOf(tc, d, 'double');
    verifySize(tc, d, [4 4]);
    tc.verifyEqual(d(1:3,4), [1 2 3]');

    aa = [a a a a a];
    d = double(aa);
    verifyInstanceOf(tc, d, 'double');
    verifySize(tc, d, [4 4 5]);
    tc.verifyEqual(d(1:3,4, 1), [1 2 3]');
    tc.verifyEqual(d(1:3,4, 5), [1 2 3]');
end

function staticconstructors_test(tc)
    %% rotation primitives
    for theta = [-pi/2 0 pi/2 pi]
        tc.verifyEqual(SE3.Rx(theta).T, trotx(theta), 'AbsTol', 1e-10  );
    end
    for theta = [-pi/2 0 pi/2 pi]
        tc.verifyEqual(SE3.Ry(theta).T, troty(theta), 'AbsTol', 1e-10  );
    end
    for theta = [-pi/2 0 pi/2 pi]
        tc.verifyEqual(SE3.Rz(theta).T, trotz(theta), 'AbsTol', 1e-10  );
    end

    for theta = [-pi/2 0 pi/2 pi]*180/pi
        tc.verifyEqual(SE3.Rx(theta, 'deg').T, trotx(theta, 'deg'), 'AbsTol', 1e-10  );
    end
    for theta = [-pi/2 0 pi/2 pi]
        tc.verifyEqual(SE3.Ry(theta, 'deg').T, troty(theta, 'deg'), 'AbsTol', 1e-10  );
    end
    for theta = [-pi/2 0 pi/2 pi]
        tc.verifyEqual(SE3.Rz(theta, 'deg').T, trotz(theta, 'deg'), 'AbsTol', 1e-10  );
    end

    %% 3 angle
    tc.verifyEqual(SE3.rpy( 0.1, 0.2, 0.3 ).T, rpy2tr( 0.1, 0.2, 0.3 ), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.rpy([ 0.1, 0.2, 0.3] ).T, rpy2tr( 0.1, 0.2, 0.3 ), 'AbsTol', 1e-10  );

    tc.verifyEqual(SE3.eul( 0.1, 0.2, 0.3 ).T, eul2tr( 0.1, 0.2, 0.3 ), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.eul([ 0.1, 0.2, 0.3] ).T, eul2tr( 0.1, 0.2, 0.3 ), 'AbsTol', 1e-10  );

    tc.verifyEqual(SE3.rpy( 10, 20, 30, 'deg' ).T, rpy2tr( 10, 20, 30, 'deg' ), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.rpy([ 10, 20, 30], 'deg' ).T, rpy2tr( 10, 20, 30, 'deg' ), 'AbsTol', 1e-10  );

    tc.verifyEqual(SE3.eul( 10, 20, 30, 'deg' ).T, eul2tr( 10, 20, 30, 'deg' ), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.eul([ 10, 20, 30], 'deg' ).T, eul2tr( 10, 20, 30, 'deg' ), 'AbsTol', 1e-10  );

    %% OA vectors
    tc.verifyEqual(SE3.oa([0 1 0], [0 0 1]).T, eye(4,4), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.oa([1 0 0], [0 1 0]).T, [0 0 1 0; 1 0 0 0; 0 1 0 0; 0 0 0 1]', 'AbsTol', 1e-10  );

    %% (theta, v)
    th = 0.2; v = unit([1 2 3]);
    tc.verifyEqual(SE3.angvec(th, v ).T, angvec2tr(th, v), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.angvec(-th, v ).T, angvec2tr(-th, v), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.angvec(-th, -v ).T, angvec2tr(-th, -v), 'AbsTol', 1e-10  );
    tc.verifyEqual(SE3.angvec(th, -v ).T, angvec2tr(th, -v), 'AbsTol', 1e-10  );

    %% exponential
    tc.verifyEqual(SE3.exp(zeros(4,4)).T, eye(4,4), 'AbsTol', 1e-10  );
    t= [1 2 3];
    tc.verifyEqual(SE3.exp([t 0 0 0]').T, transl(t), 'AbsTol', 1e-10  );
end

function delta_test(tc)
    T = SE3(1,2,3) * SE3.Ry(pi/2);

    d = [4 5 6 1 2 3]*1e-3;

    % delta
    Td = SE3.delta(d);
    tc.verifyEqual(Td.t, d(1:3)', 'AbsTol', 1e-10);
    tc.verifyEqual(Td.R, skew(d(4:6))+eye(3,3), 'AbsTol', 1e-10);

    % increment
    T2 = T.increment(d)
    tc.verifyEqual(double(T2), double(T.*Td), 'AbsTol', 1e-10 );
end

function isa_test(tc)

    verifyTrue(tc, SE3.isa(trotx(0)) );
    verifyTrue(tc, SE3.isa(trotx(0), 'valid') );

    verifyFalse(tc, SE3.isa(1) )
end

function resulttype_test(tc)

    t = SE3();
    verifyClass(tc, t, 'SE3');

    verifyClass(tc, t*t, 'SE3');

    verifyClass(tc, t.*t, 'SE3');
    % other combos all fail, test this?

    verifyClass(tc, t/t, 'SE3');

    verifyClass(tc, inv(t), 'SE3');
    end

function inverse_test(tc)

    R1 = rpy2r( randn(1,3) );   t1 = randn(3,1);
    T1 = rt2tr(R1, t1);

    TT1 = SE3(T1);

    % test inverse
    tc.verifyEqual(double(TT1.inv()), inv(T1), 'AbsTol', 1e-10  );

    tc.verifyEqual(double(TT1*TT1.inv()), eye(4,4), 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT1.inv()*TT1), eye(4,4), 'AbsTol', 1e-10  );

    % vector case
    tc.verifyEqual(double(TT1.inv()*TT1), eye(4,4), 'AbsTol', 1e-10  );
end

% function uminus_test(tc)
%     R1 = rpy2r( randn(1,3) );  t1 = randn(3,1); T1 = rt2tr(R1, t1);
%
%     TT1 = SE3(T1);
%
%
%     TM = - TT1;
%     tc.verifyEqual(TM.double, -T1, 'AbsTol', 1e-10  );
% end

function Rt_test(tc)

    R1 = rpy2r( randn(1,3) );  t1 = randn(3,1); T1 = rt2tr(R1, t1);

    TT1 = SE3(T1);

    tc.verifyEqual(TT1.T, T1, 'AbsTol', 1e-10  );
    tc.verifyEqual(TT1.R, R1, 'AbsTol', 1e-10  );
    tc.verifyEqual(TT1.t, t1, 'AbsTol', 1e-10  );

    tc.verifyEqual(TT1.transl, t1', 'AbsTol', 1e-10  );
    TT = [TT1 TT1 TT1];
    tc.verifyEqual(TT.transl, [t1 t1 t1]', 'AbsTol', 1e-10  );

    [a,b,c] = transl(TT1);
    tc.verifyEqual(a, t1(1));
    tc.verifyEqual(b, t1(2));
    tc.verifyEqual(c, t1(3));
end

function tv_test(tc)
    a = SE3(1,2,3);
    tc.verifyEqual(a.tv, [1 2 3]');

    aa = [a a a a];
    tc.verifyEqual(aa.tv, [1 2 3]'*[1 1 1 1]);

end

function arith_test(tc)

    R1 = rpy2r( randn(1,3) );  t1 = randn(3,1); T1 = rt2tr(R1, t1);
    R2 = rpy2r( randn(1,3) );  t2 = randn(3,1); T2 = rt2tr(R2, t2);

    TT1 = SE3(T1);
    TT2 = SE3(T2);

    I = SE3();

    %% SE3-SE3 product
    % scalar x scalar

    tc.verifyEqual(double(TT1*TT2), T1*T2, 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT2*TT1), T2*T1, 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT1*I), T1, 'AbsTol', 1e-10  );
    tc.verifyEqual(double(TT2*I), T2, 'AbsTol', 1e-10  );

    % vector x vector
    tc.verifyEqual([TT1 TT1 TT2] * [TT2 TT1 TT1], [TT1*TT2 TT1*TT1 TT2*TT1]);

    % scalar x vector
    tc.verifyEqual(TT1 * [TT2 TT1], [TT1*TT2 TT1*TT1]);

    % vector x scalar
    tc.verifyEqual([TT1 TT2]*TT2, [TT1*TT2 TT2*TT2]);

    %% SE3-vector product
    vx = [1 0 0]'; vy = [0 1 0]'; vz = [0 0 1]';

    % scalar x scalar

    tc.verifyEqual(TT1*vy, h2e( T1*e2h(vy) ), 'AbsTol', 1e-10);

    % vector x vector
    tc.verifyEqual([TT1 TT2 TT1] * [vz vx vy], [h2e(T1*e2h(vz)) h2e(T2*e2h(vx)) h2e(T1*e2h(vy))], 'AbsTol', 1e-10);

    % scalar x vector
    tc.verifyEqual(TT1 * [vx vy vz], h2e( T1*e2h([vx vy vz]) ), 'AbsTol', 1e-10);

    % vector x scalar
    tc.verifyEqual([TT1 TT2 TT1] * vy, [h2e(T1*e2h(vy)) h2e(T2*e2h(vy)) h2e(T1*e2h(vy))], 'AbsTol', 1e-10);

end

function adjoint_test(tc)
    R = rpy2r( randn(1,3) );  t = randn(3,1); T = rt2tr(R, t);
    TT = SE3(T);

    tc.verifyEqual(TT.Ad, [R skew(t)*R; zeros(3,3) R]);

    % velxform
    tc.verifyEqual(TT.velxform, [R zeros(3,3); zeros(3,3) R]);

end

function interp_test(tc)
    R = rpy2r( randn(1,3) );  t = randn(3,1); T = rt2tr(R, t);
    TT = SE3(T);
    I = SE3();

    tc.verifyEqual(double(interp(I, TT, 0)),   double(I), 'AbsTol', 1e-10 );
    tc.verifyEqual(double(interp(I, TT, 1)),   double(TT), 'AbsTol', 1e-4 );
    tc.verifyEqual(double(interp(I, TT, 0.5)), double(trinterp(T, 0.5)), 'AbsTol', 1e-10  );

end

function conversions_test(tc)

    %% 3 angle
    tc.verifyEqual(SE3.rpy( 0.1, 0.2, 0.3 ).torpy, [ 0.1, 0.2, 0.3], 'AbsTol', 1e-10  );

    tc.verifyEqual(SE3.eul( 0.1, 0.2, 0.3 ).toeul, [ 0.1, 0.2, 0.3 ], 'AbsTol', 1e-10  );

    tc.verifyEqual(SE3.rpy( 10, 20, 30, 'deg' ).R, rpy2r( 10, 20, 30, 'deg' ), 'AbsTol', 1e-10  );

    tc.verifyEqual(SE3.eul( 10, 20, 30, 'deg' ).R, eul2r( 10, 20, 30, 'deg' ), 'AbsTol', 1e-10  );

    %% (theta, v)
    th = 0.2; v = unit([1 2 3]);
    a = SE3.angvec(th, v ).toangvec;
    tc.verifyEqual(a, th, 'AbsTol', 1e-10  );

    [a,b] = SE3.angvec(th, v ).toangvec;
    tc.verifyEqual(a, th, 'AbsTol', 1e-10  );
    tc.verifyEqual(b, v, 'AbsTol', 1e-10  );

    %% quaternion

    R = rpy2r( 0.2, 0.3, 0.4);
    q = UnitQuaternion(R);
    T = SO3(R).SE3;
    verifyClass(tc, T.UnitQuaternion, 'UnitQuaternion');

    tc.verifyEqual(q, T.UnitQuaternion);

    %% Twist
    T = SE3.exp([1 2 3 0.2 0.3 0.4]);
    t = T.Twist();
    verifyInstanceOf(tc, t, 'Twist');
    tc.verifyEqual(t.v, [1 2 3]', 'AbsTol', 1e-10);
    tc.verifyEqual(t.w, [0.2 0.3 0.4]', 'AbsTol', 1e-10);

    %%  SE3                     convert to SE3 class

    T = rpy2tr(0.2, 0.3, 0.4);
    tc.verifyEqual(SO3(T).SE3, SE3(T), 'AbsTol', 1e-10 );
    tc.verifyEqual(SO3(T).T, T, 'AbsTol', 1e-10 );

end

function miscellany_test(tc)

    R = rpy2r( randn(1,3) );  t = randn(3,1); T = rt2tr(R, t);
    TT = SE3(T);

    tc.verifyEqual(dim(TT), 4);

    verifyEqual( tc, isSE(TT), true );

    verifyClass(tc, TT.new, 'SE3');

    verifyClass(tc, SE3.convert(TT), 'SE3');
    verifyClass(tc, SE3.convert(T), 'SE3');
    z = SE3.convert(TT);
    tc.verifyEqual(double(z), double(TT));

    z = SE3.convert(T);
    tc.verifyEqual(double(z), T);

end

function display_test(tc)

    T1 = SE3.rand;
    T2 = SE3.rand

    T1.print
    trprint(T1)   % old style syntax
end

function compatability_test(tc)

    T1 = SE3.rand();
    T2 = SE3.rand();

    T1.plot
    trplot(T1)   % old style syntax

    T1.animate
    T1.animate(T2)
    tranimate(T1)   % old style syntax
    tranimate(T1, T2)   % old style syntax
end

author: petercorke

repo: spatialmath-matlab


text:
function mbasis = basis_matrix_overhauser_uni ( )

%*****************************************************************************80
%
%% BASIS_MATRIX_OVERHAUSER_UNI sets up the uniform Overhauser spline basis matrix.
%
%  Discussion:
%
%    This basis matrix assumes that the data points P1, P2, P3 and
%    P4 are uniformly spaced in T, and that P2 corresponds to T = 0,
%    and P3 to T = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics: Principles and Practice,
%    page 505.
%
%  Parameters:
%
%    Output, real MBASIS(4,4), the basis matrix.
%
  mbasis(1,1) = - 1.0 / 2.0;
  mbasis(1,2) =   3.0 / 2.0;
  mbasis(1,3) = - 3.0 / 2.0;
  mbasis(1,4) =   1.0 / 2.0;

  mbasis(2,1) =   2.0 / 2.0;
  mbasis(2,2) = - 5.0 / 2.0;
  mbasis(2,3) =   4.0 / 2.0;
  mbasis(2,4) = - 1.0 / 2.0;

  mbasis(3,1) = - 1.0 / 2.0;
  mbasis(3,2) =   0.0;
  mbasis(3,3) =   1.0 / 2.0;
  mbasis(3,4) =   0.0;

  mbasis(4,1) =   0.0;
  mbasis(4,2) =   2.0 / 2.0;
  mbasis(4,3) =   0.0;
  mbasis(4,4) =   0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bwithde2(catalog)
    % BWITHDE2 plot b-values with depth
    % BWITHDE2(catalog);

    % turned into function by Celso G Reyes 2017

    ZG = ZmapGlobal.Data; % used by get_zmap_globals

    BV3 = [];
    Nmin = 50;
    zdlg = ZmapDialog();
    nEventsPerWindow = 150;
    overlap = 5;

    minMag = min(catalog.Magnitude);

    zdlg.AddEdit('nEventsPerWindow', 'Number of events in each window', nEventsPerWindow, 'tooltip');
    zdlg.AddEdit('overlap', 'Overlap Factor (advances 1/overlap)', overlap, 'tooltip');
    zdlg.AddPopup('mcCalculator', 'Determine Mc using:', {'Automatic','Fixed Mc=Mmin'},1, 'tooltip',...
        {@auto_mc, @(~)minMag });

    [res, okPressed] = zdlg.Create('Name','B-value with depth parameters');
    if ~okPressed
        return
    end
    watchon;

    [allDepths, I] = sort(catalog.Depth);
    allMagnitudes = catalog.Magnitude(I);

    nEventsPerWindow = res.nEventsPerWindow;
    overlap = res.overlap;

    stepSize = nEventsPerWindow / overlap;
    stepSize = max(stepSize, 1);
    startIdxs = 1 : stepSize : numel(allDepths) - nEventsPerWindow;
    endIdxs = startIdxs + nEventsPerWindow - 1;
    nSteps = numel(startIdxs);

    for t = 1 : nSteps
        % calculate b-value based an weighted LS
        startIdx = startIdxs(t);
        endIdx   = endIdxs(t);
        magnitudes = allMagnitudes(startIdx : endIdx);

        magco = res.mcCalculator(magnitudes);

        l = magnitudes >= magco - 0.05;
        if sum(l) >= Nmin
            [bv, stan, ~] = calc_bmemag(magnitudes(l), 0.1);
        else
            [bv, stan] = deal(nan);
        end
        minDepth = allDepths(startIdx);
        maxDepth = allDepths(endIdx);
        medianDepth = allDepths(startIdx + round(nEventsPerWindow/2));
        tripleidx = t.*3 - 2;
        bvalueTriplet(tripleidx : tripleidx+2, 1)     = [bv; bv; inf];
        mimaDepTriplet(tripleidx : tripleidx+2, 1) = [minDepth; maxDepth; inf];

        BV3(t,1:3) = [bv, medianDepth, stan];
    end

    watchoff

    plot_it(bvalueTriplet, mimaDepTriplet, BV3(:,1), BV3(:,2), BV3(:,3))
    return

    function plot_it(bvTriplet, minmaxDepthTriplet, bv, medianDepth, stan)
        myFigName='b-value with depth';
        errorColor = [0.5, 0.5, 0.5];
        % Find out if figure already exists
        %
        bdep = findobj('Type','Figure','-and','Name', myFigName);

        % Set up the Cumulative Number window

        if isempty(bdep)
            bdep = figure('Name'  , myFigName,...
                'NumberTitle'   , 'off', ...
                'NextPlot'      , 'add', ...
                'backingstore'  , 'on',...
                'Visible'       , 'on', ...
                'Position'      , position_in_current_monitor(ZG.map_len(1)-50, ZG.map_len(2)-20));
        else
            figure(bdep);
            clf(bdep);
        end
        orient tall
        rect    = [0.25, 0.15, 0.5, 0.75];
        ax      = axes(bdep, 'position',rect);
        ple     = errorbar(ax, medianDepth, bv, stan, stan, 'k');
        ple(1).Color = errorColor;

        ax.NextPlot = 'add';
        plot(ax, minmaxDepthTriplet, bvTriplet, 'color', errorColor);
        plot(ax, medianDepth, bv, 'sk', 'LineWidth', 1.0, 'MarkerSize', 4,...
            'MarkerFaceColor', 'w', 'MarkerEdgeColor', 'k', 'Marker','s');

        set(ax, 'box','on', 'SortMethod', 'childorder',...
            'TickDir', 'out', 'Ticklength', [ 0.02, 0.02],...
            'FontWeight', 'bold', 'FontSize', ZmapGlobal.Data.fontsz.m,...
            'Linewidth', 1)
        strib = sprintf('%s, ni = %g, Mmin = %g', catalog.Name, nEventsPerWindow, minMag);
        ylabel(ax, 'b-value')
        xlabel(ax, 'Depth [km]')
        title(ax, strib, 'FontWeight', 'bold', ...
            'FontSize', ZmapGlobal.Data.fontsz.m,...
            'Color', 'k',...
            'Interpreter', 'none')

        view(ax, [90 90])
    end

end

function [magco] = auto_mc(magnitudes)
    [Mc90, Mc95] = mcperc_ca3(magnitudes);
    if ~isnan(Mc95)
        magco = Mc95;
    elseif ~isnan(Mc90)
        magco = Mc90;
    else
        [magco] =  bvalca3(magnitudes, McAutoEstimate.auto);
    end
end

author: CelsoReyes


text:
%% An illustration of the use of MLFFIT2
%
% Igor Podlubny
%
% BERG Faculty, Technical University of Kosice, Slovakia

%% Introduction
%
% Several examples of the use of the
% two-parameter Mittag-Leffler function
%%
%
% $$y(x) = C t^{\beta -1} E_{\alpha,\beta}(a x^\alpha)$$
%
% for fitting data are provided below.
%

%% 1. Fitting back the noised Mittag-Leffler function:
%
% The first example is the "restoration" of the function
% $y(x) = 0.8 E_{1.5}(-0.2 x^{1.5})$. Three series of "measured" data
% are created by adding noise to the values of the Mittag-Leffler
% function at the same set of nodes x. Such moisy data
% are fitted using MLFFIT1.M function.
%

% Clear the workspace:
clear

% Define the set of nodes (x) and the parameters $\alpha$
% and $\beta$ of the % Mittag-Leffler function:

x = 0.1:0.35:16;
alfa = 1.5;
beta = 0.8;

% Now let us simulate measurements by adding noise to
% the exact values of the original function

y1 = 1 * x.^(beta-1) .* mlf(alfa, beta, -0.2*x.^alfa, 7) + (-.05 + .1*rand(size(x)));
y2 = 1 * x.^(beta-1) .* mlf(alfa, beta, -0.2*x.^alfa, 7) + (-.05 + .1*rand(size(x)));
y3 = 1 * x.^(beta-1) .* mlf(alfa, beta, -0.2*x.^alfa, 7) + (-.05 + .1*rand(size(x)));

% and fit these "measurements" by calling MLFFIT2:
[c, R2] = mlffit2([x x x], [y1 y2 y3], [2; 2; .5; -1], 6);

% Let us check if the coefficients of the fitting Mittag-Leffler
% function are close to the original coefficients:
alpha = c(1)
beta = c(2)
C = c(3)
a = c(4)

% Finally, we can plot the "measurements", the original function
% and the function that fits the "measurement":

xfine = x(1):0.01:x(end);
yfit = c(3)* xfine.^(c(2)-1) .* mlf(c(1), c(2), c(4)*xfine.^c(1), 6);
yorig = 1 * xfine.^(beta-1) .* mlf(alfa, beta, -0.2*xfine.^alfa, 6);

figure(1)
plot(xfine, yorig, 'r', x,y1, '.b', x, y2, '.g', x, y3, '.m', xfine, yfit, 'k')
grid on
legend('original function', 'noised data', 'noised data', 'noised data', 'fitting')

%% 2. Fitting the classics:
%
%

% Clear the workspace and define the set of nodes (x):
clear
x = 0.1:0.35:16;

% Now let us simulate measurements by adding noise to
% the exact values of the original function

y1 = exp(x).*erfc(sqrt(x))  + (-.02 + .04*rand(size(x)));
y2 = exp(x).*erfc(sqrt(x))  + (-.02 + .04*rand(size(x)));
y3 = exp(x).*erfc(sqrt(x))  + (-.02 + .04*rand(size(x)));

% and fit these "measurements" by calling MLFFIT2:
[c, R2] = mlffit2([x x x], [y1 y2 y3], [0.5; 0.5; .5; -0.5], 6);

% Let us check if the coefficients of the fitting Mittag-Leffler
% function:
alpha = c(1)
beta = c(2)
C = c(3)
a = c(4)

% Finally, we can plot the "measurements", the original function
% and the function that fits the "measurement":

xfine = x(1):0.01:x(end);
yfit = c(3)* xfine.^(c(2)-1) .* mlf(c(1), c(2), c(4)*xfine.^c(1), 6);
yorig = exp(xfine).*erfc(sqrt(xfine));

figure(2)
plot(xfine, yorig, 'r', x,y1, '.b', x, y2, '.g', x, y3, '.m', xfine, yfit, 'k')
grid on
legend('original function', 'noised data', 'noised data', 'noised data', 'fitting')

%% 3. Fitting the cosine function:

% Clear the workspace and define the set of nodes (x):
clear
x = 0.1:0.35:16;

% Now let us simulate measurements by adding noise to
% the exact values of the original function

y1 = cos(x)  + (-.1 + .2*rand(size(x)));
y2 = cos(x)  + (-.1 + .2*rand(size(x)));
y3 = cos(x)  + (-.1 + .2*rand(size(x)));

% and fit these "measurements" by calling MLFFIT2:
[c, R2] = mlffit2([x x x], [y1 y2 y3], [0.5; 0.5; .5; -0.5], 6);

% Let us check if the coefficients of the fitting Mittag-Leffler
% function:
alpha = c(1)
beta = c(2)
C = c(3)
a = c(4)

% Finally, we can plot the "measurements", the original function
% and the function that fits the "measurement":

xfine = x(1):0.01:x(end);
yfit = c(3)* xfine.^(c(2)-1) .* mlf(c(1), c(2), c(4)*xfine.^c(1), 6);
yorig = cos(xfine);

figure(3)
plot(xfine, yorig, 'r', x,y1, '.b', x, y2, '.g', x, y3, '.m', xfine, yfit, 'k')
grid on
legend('original function', 'noised data', 'noised data', 'noised data', 'fitting')

%% 4. Fitting the sine function:

% Clear the workspace and define the set of nodes (x):
clear
x = 0.1:0.35:16;

% Now let us simulate measurements by adding noise to
% the exact values of the original function

y1 = sin(x)  + (-.2 + .4*rand(size(x)));
y2 = sin(x)  + (-.2 + .4*rand(size(x)));
y3 = sin(x)  + (-.2 + .4*rand(size(x)));

% and fit these "measurements" by calling MLFFIT2:
[c, R2] = mlffit2([x x x], [y1 y2 y3], [1.5; 1.5; .5; -1], 6);

% Let us check if the coefficients of the fitting Mittag-Leffler
% function:
alpha = c(1)
beta = c(2)
C = c(3)
a = c(4)

% Finally, we can plot the "measurements", the original function
% and the function that fits the "measurement":

xfine = x(1):0.01:x(end);
yfit = c(3)* xfine.^(c(2)-1) .* mlf(c(1), c(2), c(4)*xfine.^c(1), 6);
yorig = sin(xfine);

figure(3)
plot(xfine, yorig, 'r', x,y1, '.b', x, y2, '.g', x, y3, '.m', xfine, yfit, 'k')
grid on
legend('original function', 'noised data', 'noised data', 'noised data', 'fitting')

%% 5. Dumped oscillation fitting:

% Clear the workspace and define the set of nodes (x):
clear
x = 0.1:0.35:16;

% Now let us simulate measurements by adding noise to
% the exact values of the original function

y1 = exp(-0.1*x).*cos(x) + (-.05 + .1*rand(size(x)));
y2 = exp(-0.1*x).*cos(x) + (-.05 + .1*rand(size(x)));
y3 = exp(-0.1*x).*cos(x) + (-.05 + .1*rand(size(x)));

% and fit these "measurements" by calling MLFFIT2:
[c, R2] = mlffit2([x x x], [y1 y2 y3], [2; 2; .5; -1], 6);

% Let us check if the coefficients of the fitting Mittag-Leffler
% function:
alpha = c(1)
beta = c(2)
C = c(3)
a = c(4)

% Finally, we can plot the "measurements", the original function
% and the function that fits the "measurement":

xfine = x(1):0.01:x(end);
yfit = c(3)* xfine.^(c(2)-1) .* mlf(c(1), c(2), c(4)*xfine.^c(1), 6);
yorig = exp(-0.1*xfine).*cos(xfine);

figure(3)
plot(xfine, yorig, 'r', x,y1, '.b', x, y2, '.g', x, y3, '.m', xfine, yfit, 'k')
grid on
legend('original function', 'noised data', 'noised data', 'noised data', 'fitting')

author: Sable

repo: mcbench-benchmarks


text:
% mu = weighted_mean(x,w)
% x is n x d, w is n x 1
function mu = weighted_mean(x,w)

%[n,d] = size(x);

z = sum(w);
x((isnan(x)|isinf(x)) & (w(:)==0)) = 0;
mu = sum(bsxfun(@times,x,w(:)),1)/z;

author: kristinbranson

repo: JAABA


text:
function [Pk,Phi,d,e] = intialize(p,q)
clc;
disp('Intaializing:');
Pk=p*q;
Phi=(p-1)*(q-1);
%Calculate the value of e
x=2;e=1;
while x > 1
    e=e+1;
    x=gcd(Phi,e);
end
%Calculate the value of d
i=1;
r=1;
while r > 0
    k=(Phi*i)+1;
    r=rem(k,e);
    i=i+1;
end
d=k/e;
clc;
disp(['The value of (N) is: ' num2str(Pk)]);
disp(['The public key (e) is: ' num2str(e)]);
disp(['The value of (Phi) is: ' num2str(Phi)]);
disp(['The private key (d)is: ' num2str(d)]);


text:
function [pvec, pstruct] = tapas_bayes_optimal_binary_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% This is a dummy function

pvec    = [];
pstruct = struct;

return;


text:
function [h_0,h_1] = daubcqf(N,TYPE)
%    [h_0,h_1] = daubcqf(N,TYPE);
%
%    Function computes the Daubechies' scaling and wavelet filters
%    (normalized to sqrt(2)).
%
%    Input:
%       N    : Length of filter (must be even)
%       TYPE : Optional parameter that distinguishes the minimum phase,
%              maximum phase and mid-phase solutions ('min', 'max', or
%              'mid'). If no argument is specified, the minimum phase
%              solution is used.
%
%    Output:
%       h_0 : Minimal phase Daubechies' scaling filter
%       h_1 : Minimal phase Daubechies' wavelet filter
%
%    Example:
%       N = 4;
%       TYPE = 'min';
%       [h_0,h_1] = daubcqf(N,TYPE)
%       h_0 = 0.4830 0.8365 0.2241 -0.1294
%       h_1 = 0.1294 0.2241 -0.8365 0.4830
%
%    Reference: "Orthonormal Bases of Compactly Supported Wavelets",
%                CPAM, Oct.89
%
%Author: Ramesh Gopinath  <ramesh@dsp.rice.edu>

if(nargin < 2),
  TYPE = 'min';
end;
if(rem(N,2) ~= 0),
  error('No Daubechies filter exists for ODD length');
end;

persistent cache;   % cache of precomputed wavelet solutions
key = [char('a'+N) TYPE];
try
    [h_0,h_1] = deal(cache.(key){:});
catch
    K = N/2;
    a = 1;
    p = 1;
    q = 1;
    h_0 = [1 1];
    for j  = 1:K-1,
      a = -a * 0.25 * (j + K - 1)/j;
      h_0 = [0 h_0] + [h_0 0];
      p = [0 -p] + [p 0];
      p = [0 -p] + [p 0];
      q = [0 q 0] + a*p;
    end;
    q = sort(roots(q));
    qt = q(1:K-1);
    if TYPE=='mid',
      if rem(K,2)==1,
        qt = q([1:4:N-2 2:4:N-2]);
      else
        qt = q([1 4:4:K-1 5:4:K-1 N-3:-4:K N-4:-4:K]);
      end;
    end;
    h_0 = conv(h_0,real(poly(qt)));
    h_0 = sqrt(2)*h_0/sum(h_0); 	%Normalize to sqrt(2);
    if(TYPE=='max'),
      h_0 = fliplr(h_0);
    end;
    if(abs(sum(h_0 .^ 2))-1 > 1e-4)
      error('Numerically unstable for this value of "N".');
    end;
    h_1 = rot90(h_0,2);
    h_1(1:2:N)=-h_1(1:2:N);
    cache.(key) = {h_0,h_1};
end

repo: REST