text:
function lik = lik_loglogistic(varargin)
%LIK_LOGLOGISTIC Create a right censored log-logistic likelihood structure
%
%  Description
%    LIK = LIK_LOGLOGISTIC('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a likelihood structure for right censored log-logistic
%    survival model in which the named parameters have the
%    specified values. Any unspecified parameters are set to
%    default values.
%
%    LIK = LIK_LOGLOGISTIC(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood structure with the named parameters
%    altered with the specified values.
%
%    Parameters for loggaussian likelihood [default]
%      shape       - shape parameter r [1]
%      shape_prior - prior for shape [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 [ (r/exp(f_i)*(y_i/exp(f_i))^(r-1)/
%                               (1+(y_i/exp(f_i))^r))^(1-z_i)
%                               *(1+(y_i/exp(f_i))^r)^(-z_i) ]
%
%
%    where r is the shape parameter of log-logistic distribution.
%    z is a vector of censoring indicators with z = 0 for uncensored event
%    and z = 1 for right censored event.
%
%    When using the loggaussian 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_LOGLOGISTIC';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('shape',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('shape_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

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

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

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_loglogistic_pak;
    lik.fh.unpak = @lik_loglogistic_unpak;
    lik.fh.lp = @lik_loglogistic_lp;
    lik.fh.lpg = @lik_loglogistic_lpg;
    lik.fh.ll = @lik_loglogistic_ll;
    lik.fh.llg = @lik_loglogistic_llg;
    lik.fh.llg2 = @lik_loglogistic_llg2;
    lik.fh.llg3 = @lik_loglogistic_llg3;
    lik.fh.tiltedMoments = @lik_loglogistic_tiltedMoments;
    lik.fh.siteDeriv = @lik_loglogistic_siteDeriv;
    lik.fh.invlink = @lik_loglogistic_invlink;
    lik.fh.predy = @lik_loglogistic_predy;
    lik.fh.recappend = @lik_loglogistic_recappend;
    lik.fh.predcdf=@lik_loglogistic_predcdf;
  end

end

function [w,s] = lik_loglogistic_pak(lik)
%LIK_LOGLOGISTIC_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_LOGLOGISTIC_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.shape)
%
%   See also
%   LIK_LOGLOGISTIC_UNPAK, GP_PAK

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

function [lik, w] = lik_loglogistic_unpak(lik, w)
%LIK_LOGLOGISTIC_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    [LIK, W] = LIK_LOGLOGISTIC_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.shape)
%
%   See also
%   LIK_LOGLOGISTIC_PAK, GP_UNPAK

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

function lp = lik_loglogistic_lp(lik, varargin)
%LIK_LOGLOGISTIC_LP  log(prior) of the likelihood parameters
%
%  Description
%    LP = LIK_LOGLOGISTIC_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_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, LIK_LOGLOGISTIC_LLG2, GPLA_E

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

end

function lpg = lik_loglogistic_lpg(lik)
%LIK_LOGLOGISTIC_LPG  d log(prior)/dth of the likelihood
%                parameters th
%
%  Description
%    E = LIK_LOGLOGISTIC_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_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, LIK_LOGLOGISTIC_LLG2, GPLA_G

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

function ll = lik_loglogistic_ll(lik, y, f, z)
%LIK_LOGLOGISTIC_LL  Log likelihood
%
%  Description
%    LL = LIK_LOGLOGISTIC_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_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, LIK_LOGLOGISTIC_LLG2, GPLA_E

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

  r = lik.shape;
  ll = sum((1-z).*(log(r)+(r-1).*log(y)-r.*f) +(z-2).*log(1+(y./exp(f)).^r));

end

function llg = lik_loglogistic_llg(lik, y, f, param, z)
%LIK_LOGLOGISTIC_LLG  Gradient of the log likelihood
%
%  Description
%    LLG = LIK_LOGLOGISTIC_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_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_LLG2, LIK_LOGLOGISTIC_LLG3, GPLA_E

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

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

function llg2 = lik_loglogistic_llg2(lik, y, f, param, z)
%LIK_LOGLOGISTIC_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_LOGLOGISTIC_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_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG3, GPLA_E

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

  r = lik.shape;
  switch param
    case 'param'

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

function llg3 = lik_loglogistic_llg3(lik, y, f, param, z)
%LIK_LOGLOGISTIC_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_LOGLOGISTIC_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_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_LLG, LIK_LOGLOGISTIC_LLG2, GPLA_E, GPLA_G

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

  r = lik.shape;
  switch param
    case 'param'

    case 'latent'
      llg3 = r.^3.*(z-2).*(y./exp(f)).^r.*(-1+(y./exp(f)).^r)./(1+(y./exp(f)).^r).^3;
    case 'latent2+param'
      llg3 = -(r.*(z-2).*(y./exp(f)).^r.*(-2-2.*(y./exp(f)).^r + ...
              r.*(-1+(y./exp(f)).^r).*log(y./exp(f))))./(1+(y./exp(f)).^r).^3;
      % correction due to the log transformation
      llg3 = llg3.*lik.shape;
  end
end

function [logM_0, m_1, sigm2hati1] = lik_loglogistic_tiltedMoments(lik, y, i1, sigm2_i, myy_i, z)
%LIK_LOGLOGISTIC_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_LOGLOGISTIC_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_loglogistic -> lik_loglogistic_tiltedMoments: missing z!'...
          'loglogistic likelihood needs the censoring            '...
          'indicators as an extra input z. See, for                 '...
          'example, lik_loglogistic and gpep_e.                       ']);
 end

  yy = y(i1);
  yc = 1-z(i1);
  r = lik.shape;
  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_loglogistic_norm(yy(i),myy_i(i),sigm2_i(i),yc(i),r);

    % 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_loglogistic_tilted_moments: sigm2hati1 >= sigm2_i');
      end
    end
    logM_0(i) = log(m_0);
  end
end

function [g_i] = lik_loglogistic_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_LOGLOGISTIC_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_LOGLOGISTIC_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_loglogistic -> lik_loglogistic_siteDeriv: missing z!'...
           'loglogistic likelihood needs the censoring        '...
           'indicators as an extra input z. See, for             '...
           'example, lik_loglogistic and gpla_e.                   ']);
  end

  yy = y(i1);
  yc = 1-z(i1);
  r = lik.shape;

  % get a function handle of an unnormalized tilted distribution
  % (likelihood * cavity = Log-Gaussian * Gaussian)
  % and useful integration limits
  [tf,minf,maxf]=init_loglogistic_norm(yy,myy_i,sigm2_i,yc,r);
  % 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.*r;

  function g = deriv(f)
    g = yc.*(1/r+log(yy)-f) + (-1-yc)./(1+(yy./exp(f)).^r).* ...
             (yy./exp(f)).^r.*(log(yy)-f);
  end
end

function [lpy, Ey, Vary] = lik_loglogistic_predy(lik, Ef, Varf, yt, zt)
%LIK_LOGLOGISTIC_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_LOGLOGISTIC_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_LOGLOGISTIC_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_loglogistic -> lik_loglogistic_predy: missing zt!'...
           'loglogistic likelihood needs the censoring    '...
           'indicators as an extra input zt. See, for         '...
           'example, lik_loglogistic and gpla_e.               ']);
  end

  yc = 1-zt;
  r = lik.shape;

  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_loglogistic_norm(...
        yt(i1),Ef(i1),Varf(i1),yc(i1),r);
      % integrate over the f to get posterior predictive distribution
      lpy(i1) = log(quadgk(pdf, minf, maxf));
    end
  end
end

function [df,minf,maxf] = init_loglogistic_norm(yy,myy_i,sigm2_i,yc,r)
%INIT_LOGLOGISTIC_NORM
%
%  Description
%    Return function handle to a function evaluating
%    loglogistic * Gaussian which is used for evaluating
%    (likelihood * cavity) or (likelihood * posterior) Return
%    also useful limits for integration. This is private function
%    for lik_loglogistic. This subfunction is needed by subfunctions
%    tiltedMoments, siteDeriv and predy.
%
%  See also
%    LIK_LOGLOGISTIC_TILTEDMOMENTS, LIK_LOGLOGISTIC_SITEDERIV,
%    LIK_LOGLOGISTIC_PREDY

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

  % Create function handle for the function to be integrated
  df = @loglogistic_norm;
  % use log to avoid underflow, and derivates for faster search
  ld = @log_loglogistic_norm;
  ldg = @log_loglogistic_norm_g;
  ldg2 = @log_loglogistic_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 loglogistic likelihood and Gaussian
    mu=-log(yc.^(1/r)./yy);
    s2=-r^2.*yc.*(-1-yc)./(1+yc).^2;
%     s2=1;
    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_loglogistic -> init_loglogistic_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_loglogistic -> init_loglogistic_norm: ' ...
             'integration interval maximun not found ' ...
             'even after looking hard!'])
    end
  end

  function integrand = loglogistic_norm(f)
  % loglogistic * Gaussian
    integrand = exp(ldconst ...
                     - yc.*r.*f +(-1-yc).*log(1+(yy./exp(f)).^r) ...
                    -0.5*(f-myy_i).^2./sigm2_i);
  end

  function log_int = log_loglogistic_norm(f)
  % log(loglogistic * Gaussian)
  % log_loglogistic_norm is used to avoid underflow when searching
  % integration interval
    log_int = ldconst ...
               -yc.*r.*f +(-1-yc).*log(1+(yy./exp(f)).^r) ...
              -0.5*(f-myy_i).^2./sigm2_i;
  end

  function g = log_loglogistic_norm_g(f)
  % d/df log(loglogistic * Gaussian)
  % derivative of log_loglogistic_norm
    g = -r.*yc - (-1-yc).*r.*(yy./exp(f)).^r./(1+(yy./exp(f)).^r) ...
        + (myy_i - f)./sigm2_i;
  end

  function g2 = log_loglogistic_norm_g2(f)
  % d^2/df^2 log(loglogistic * Gaussian)
  % second derivate of log_loglogistic_norm
    g2 =  r.^2.*(-1-yc).*(yy./exp(f)).^r./(1+(yy./exp(f)).^r).^2 ...
              -1/sigm2_i;
  end

end

function cdf = lik_loglogistic_predcdf(lik, Ef, Varf, yt)
%LIK_LOGLOGISTIC_PREDCDF  Returns the predictive cdf evaluated at yt
%
%  Description
%    CDF = LIK_LOGLOGISTIC_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

  r = lik.shape;

  % 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-logistic * Gaussian)
    % and useful integration limits.
    % yc=0 when evaluating predictive cdf
    [pdf,minf,maxf]=init_loglogistic_norm(...
      yt(i1),Ef(i1),Varf(i1),0,r);
    % integrate over the f to get posterior predictive distribution
    cdf(i1) = 1-quadgk(pdf, minf, maxf);
  end
end

function p = lik_loglogistic_invlink(lik, f)
%LIK_loglogistic Returns values of inverse link function
%
%  Description
%    P = LIK_LOGLOGISTIC_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_LOGLOGISTIC_LL, LIK_LOGLOGISTIC_PREDY

p = exp(f);
end

function reclik = lik_loglogistic_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = GPCF_LOGLOGISTIC_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-Logistic';

    % Initialize parameter
    reclik.shape = [];

    % Set the function handles
    reclik.fh.pak = @lik_loglogistic_pak;
    reclik.fh.unpak = @lik_loglogistic_unpak;
    reclik.fh.lp = @lik_loglogistic_lp;
    reclik.fh.lpg = @lik_loglogistic_lpg;
    reclik.fh.ll = @lik_loglogistic_ll;
    reclik.fh.llg = @lik_loglogistic_llg;
    reclik.fh.llg2 = @lik_loglogistic_llg2;
    reclik.fh.llg3 = @lik_loglogistic_llg3;
    reclik.fh.tiltedMoments = @lik_loglogistic_tiltedMoments;
    reclik.fh.invlink = @lik_loglogistic_invlink;
    reclik.fh.predy = @lik_loglogistic_predy;
    reclik.fh.predcdf=@lik_loglogistic_predcdf;
    reclik.fh.recappend = @lik_loglogistic_recappend;
    reclik.p=[];
    reclik.p.shape=[];
    if ~isempty(ri.p.shape)
      reclik.p.shape = ri.p.shape;
    end
  else
    % Append to the record
    reclik.shape(ri,:)=lik.shape;
    if ~isempty(lik.p)
      reclik.p.shape = lik.p.shape.fh.recappend(reclik.p.shape, ri, lik.p.shape);
    end
  end
end


text:
% NUMEL   Number of elements in an array or subscripted array expression.
%    N = NUMEL(A) returns the number of elements, N, in array A, equivalent
%    to PROD(SIZE(A)).
%
%    See also SIZE, PROD, SUBSREF.
%
%    Reference page in Doc Center
%       doc numel
%
%    Other functions named numel
%
%       categorical/numel      distributed/numel    tabular/numel
%       codistributed/numel    gpuArray/numel       tall/numel
%       Composite/numel        sym/numel            ts/numel
%       dataset/numel          symfun/numel
%

author: jmaih


text:
function [eventID] = ma_GetEventIDByEventNum(eventNum, script)
%getEventIDByEventNum Summary of this function goes here
%   Detailed explanation goes here
    try
        eventID = script{eventNum}.id;
    catch
        eventID = -1;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
close all
clear
clc

%% overall score
load('./AADBinfo.mat');
score = [testScore, trainScore];
figure;
hist(score, 10);
xlabel('score');
ylabel('#images');

%% attributes
attNameList = {'BalacingElements', 'ColorHarmony', 'Content', 'DoF', 'Light', 'MotionBlur', 'Object', 'Repetition', 'RuleOfThirds', 'Symmetry', 'VividColor' };
folderpath = './AADB_prepare';
attScore = cell(11,1);
for attID = 1:length(attNameList)
    fprintf('attribute-%s\n', attNameList{attID});
    attScore{attID} = [];
    filename = fullfile(folderpath, ['imgListTrainRegression_' attNameList{attID} '.txt']);
    numImage = numel(textread(filename,'%1c%*[^\n]'));
    scoreTMP = zeros(1,numImage);
    fn = fopen( filename, 'r' );
    for i = 1:numImage
        tline = fgets(fn);
        C = strsplit(tline, ' ');
%         imgName = C{1};
        imgLabel = str2double(C{2});
%         nameList{i} = imgName;
        scoreTMP(i) = imgLabel;
    end
    fclose(fn);
    attScore{attID} = [attScore{attID}, scoreTMP];

    filename = fullfile(folderpath, ['imgListTestRegression_' attNameList{attID} '.txt']);
    numImage = numel(textread(filename,'%1c%*[^\n]'));
    scoreTMP = zeros(1,numImage);
    fn = fopen( filename, 'r' );
    for i = 1:numImage
        tline = fgets(fn);
        C = strsplit(tline, ' ');
%         imgName = C{1};
        imgLabel = str2double(C{2});
%         nameList{i} = imgName;
        scoreTMP(i) = imgLabel;
    end
    fclose(fn);
    attScore{attID} = [attScore{attID}, scoreTMP];
end

%% show attribute distribution

A = cell2mat(attScore);
countA = zeros(11,3);
for attID = 1:length(attNameList)
    a = find(A(attID,:)==1);
    countA(attID, 1) = length(a);
    a = find(A(attID,:)==2);
    countA(attID, 2) = length(a);
    a = find(A(attID,:)==3);
    countA(attID, 3) = length(a);
end

figure;
h = bar(countA);
set(gca,'XTickLabel',{'BalacingElements'; 'ColorHarmony'; 'Content'; 'DoF'; 'Light'; 'MotionBlur'; 'Object'; 'Repetition'; 'RuleOfThirds'; 'Symmetry'; 'VividColor' })
l = {'negative', 'null', 'positive'};
legend(h,l);

author: aimerykong

repo: deepImageAestheticsAnalysis


text:
%Q2TR	Convert unit-quaternion to homogeneous transform
%
%	T = q2tr(Q)
%
%	Return the rotational homogeneous transform corresponding to the unit
%	quaternion Q.
%
%	See also TR2Q

%	Copyright (C) 1993 Peter Corke
function t = q2tr(q)
	s = q(1);
	x = q(2);
	y = q(3);
	z = q(4);

	r = [	1-2*(y^2+z^2)	2*(x*y-s*z)	2*(x*z+s*y)
		2*(x*y+s*z)	1-2*(x^2+z^2)	2*(y*z-s*x)
		2*(x*z-s*y)	2*(y*z+s*x)	1-2*(x^2+y^2)	];
	t = eye(4,4);
	t(1:3,1:3) = r;
	t(4,4) = 1;

repo: ekfmonoslam


text:
function f = dubna(varargin)
% Construct the odf from the dubna data set as example for an SO3Fun.

f = mtexdata('dubnaODF');

if nargin > 0 && isa(varargin{1},'rotation')
  ori = varargin{1};
  f = f.eval(ori);
end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function gui_adv(x, y)
%
% This file is needed by FASTICAG

% This is the advanced options -dialog

% @(#)$Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables

% Handle to the window
global hf_FastICA_AdvOpt;

% Handles to some of the controls in window
global hpm_FastICA_finetune;
global he_FastICA_a1;
global he_FastICA_a2;
global he_FastICA_myy;
global he_FastICA_epsilon;
global he_FastICA_maxIterations;
global he_FastICA_maxFinetune;
global he_FastICA_sampleSize
global hpm_FastICA_initState;
global hb_FastICA_initGuess;
global ht_FastICA_initGuess;
global hpm_FastICA_displayMode;
global he_FastICA_displayInterval;
global hpm_FastICA_verbose;

% Some of the main variables needed
global g_FastICA_initGuess;
global g_FastICA_finetune;
global g_FastICA_a1;
global g_FastICA_a2;
global g_FastICA_myy;
global g_FastICA_epsilon;
global g_FastICA_maxNumIte;
global g_FastICA_maxFinetune;
global g_FastICA_initState;
global g_FastICA_sampleSize;
global g_FastICA_displayMo;
global g_FastICA_displayIn;
global g_FastICA_verbose;

global c_FastICA_appr_strD;
global c_FastICA_appr_strV;
global c_FastICA_finetune_strD;
global c_FastICA_finetune_strV;
global c_FastICA_iSta_strD;
global c_FastICA_iSta_strV;
global c_FastICA_dMod_strD;
global c_FastICA_dMod_strV;
global c_FastICA_verb_strD;
global c_FastICA_verb_strV;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Configuration options
FIGURENAME = 'FastICA: advanced options';
FIGURETAG = 'f_FastICAAdvOpt';
FIGURESIZE = [x y 450 280];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check to see if this figure is already open - it should not!
% Can't have more than one copy - otherwise the global
% variables and handles can get mixed up.
if ~isempty(findobj('Tag',FIGURETAG))
  error('Error: advanced options dialog already open!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize some of the controls' values

% Did we already load some initial guess
pm_initState_Value = g_FastICA_initState;
if isempty(g_FastICA_initGuess) | (g_FastICA_initGuess == 1)
  t_initGuess_String = 'Not loaded';
else
  t_initGuess_String = 'Loaded';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the figure
a = figure('Color',[0.8 0.8 0.8], ...
	   'PaperType','a4letter', ...
	   'Name', FIGURENAME, ...
	   'NumberTitle', 'off', ...
	   'Tag', FIGURETAG, ...
	   'Position', FIGURESIZE, ...
	   'MenuBar', 'none');
set (a, 'Resize', 'off');

hf_FastICA_AdvOpt = a;

set(hf_FastICA_AdvOpt, 'HandleVisibility', 'callback');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From here on it get's ugly as I have not had time to clean it up

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the frames
pos_l=2;
pos_w=FIGURESIZE(3)-4;
pos_h=FIGURESIZE(4)-4;
pos_t=2;
h_f_adv_background = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_adv_background');

pos_w=120;
pos_l=FIGURESIZE(3)-(pos_w+2+2);
pos_h=FIGURESIZE(4)-2*4;
pos_t=4;
h_f_adv_side = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_adv_side');

pos_l=4;
pos_w=FIGURESIZE(3)-8-pos_w-2;
pos_h=FIGURESIZE(4)-8;
pos_t=4;
h_f_advopt = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_advopt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_advopt
bgc = get(h_f_advopt, 'BackgroundColor');

pos_w1=230;
pos_w2=70;

pos_frame=get(h_f_advopt, 'Position');
pos_h = 20;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - 6;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Fine-tune (g)', ...
  'Style','text', ...
  'Tag','t_727');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_finetune = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_finetune_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_finetune', ...
  'Value',g_FastICA_finetune);

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Parameter a1 (g = ''tanh'')', ...
  'Style','text', ...
  'Tag','t_22');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_a1 = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc Checka1', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_a1), ...
  'Style','edit', ...
  'Tag','e_a1');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Parameter a2 (g = ''gauss'')', ...
  'Style','text', ...
  'Tag','t_222');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_a2 = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc Checka2', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_a2), ...
  'Style','edit', ...
  'Tag','e_a2');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','mu (step size)', ...
  'Style','text', ...
  'Tag','t_223');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_myy = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc CheckMyy', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_myy), ...
  'Style','edit', ...
  'Tag','e_myy');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','epsilon (stopping criterion)', ...
  'Style','text', ...
  'Tag','t_23');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_epsilon = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_epsilon), ...
  'Style','edit', ...
  'Tag','e_epsilon');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Maximum number of iterations', ...
  'Style','text', ...
  'Tag','t_24');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_maxIterations = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_maxNumIte), ...
  'Style','edit', ...
  'Tag','e_maxIterations');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Maximum iterations in fine-tuning', ...
  'Style','text', ...
  'Tag','t_2412');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_maxFinetune = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_maxFinetune), ...
  'Style','edit', ...
  'Tag','e_maxFinetune');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Sample size (proportion)', ...
  'Style','text', ...
  'Tag','t_224');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_sampleSize = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Callback','gui_advc CheckSampleSize', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_sampleSize), ...
  'Style','edit', ...
  'Tag','e_sampleSize');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Initial state', ...
  'Style','text', ...
  'Tag','t_25');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_initState = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_iSta_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_initState', ...
  'Value',pm_initState_Value);

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
hb_FastICA_initGuess = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc loadGuess', ...
  'Position',[pos_l pos_t (pos_w1-60) pos_h], ...
  'String','Load initial guess', ...
  'UserData', g_FastICA_initGuess, ...
  'Tag','b_LoadGuess');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
ht_FastICA_initGuess = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[(pos_l) pos_t (pos_w2) pos_h], ...
  'String',t_initGuess_String, ...
  'Style','text', ...
  'Tag','t_initGuess');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Display mode', ...
  'Style','text', ...
  'Tag','t_27');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_displayMode = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_dMod_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_displayMode', ...
  'Value',g_FastICA_displayMo);

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Iterations between displays', ...
  'Style','text', ...
  'Tag','t_28');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
he_FastICA_displayInterval = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','right', ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',num2str(g_FastICA_displayIn), ...
  'Style','edit', ...
  'Tag','e_displayInterval');

pos_t = pos_t - pos_h;
pos_l = pos_frame(1) + 6;
b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w1 pos_h], ...
  'String','Verbose', ...
  'Style','text', ...
  'Tag','t_29');

pos_l = pos_frame(1) + pos_frame(3) - 6 - pos_w2;
hpm_FastICA_verbose = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'Position',[pos_l pos_t pos_w2 pos_h], ...
  'String',c_FastICA_verb_strD, ...
  'Style','popupmenu', ...
  'Tag','pm_verbose', ...
  'Value',g_FastICA_verbose);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_adv_side
pos_vspace = 6;
pos_hspace = 10;
pos_frame = get(h_f_adv_side, 'Position');
pos_w = 100;
pos_h = 30;
pos_l = pos_frame(1) + pos_hspace;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc OK', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','OK', ...
  'Tag','b_advOK');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Cancel', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Cancel', ...
  'Tag','b_advCancel');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Default', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Default', ...
  'Tag','b_advDefault');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Apply', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Apply', ...
  'Tag','b_advApply');

pos_t = pos_frame(2) + pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_advc Help', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Help', ...
  'Tag','b_advHelp');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  FILE:   cnn_widerface.m
%
%    This function serves as the main function for training a model.
%
%  INPUT:  configuration     (see code for details)
%
%  OUTPUT: net               (trained network)
%          info              (training stats)

function [net, info] = cnn_widerface(varargin)
startup;

opts.keepDilatedZeros = false;
opts.inputSize = [500, 500];
opts.learningRate = 1e-4;

%% use customized training function ie. adam
opts.trainFn = '@cnn_train_dag';
opts.batchGetterFn = '@cnn_get_batch_logistic_zoom';
opts.freezeResNet = false;
opts.tag = '';
opts.clusterNum = 25;
opts.clusterName = '';
opts.bboxReg = true;
opts.skipLRMult = [1, 0.001, 0.0001, 0.00001];
opts.sampleSize = 256;
opts.posFraction = 0.5;
opts.posThresh = 0.7;
opts.negThresh = 0.3;
opts.border = [0, 0];
opts.pretrainModelPath = 'matconvnet/pascal-fcn8s-tvg-dag.mat';
opts.dataDir = fullfile('data','widerface') ;
opts.modelType = 'pascal-fcn8s-tvg-dag' ;
opts.networkType = 'dagnn' ;
opts.batchNormalization = true ;
opts.weightInitMethod = 'gaussian' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.minClusterSize = [2, 2];
opts.maxClusterSize = opts.inputSize;
[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.modelType ;
if opts.freezeResNet, sfx = ['freezed-' sfx] ; end
sfx = [sfx '-' 'sample' num2str(opts.sampleSize)] ;
sfx = [sfx '-' 'posfrac' num2str(opts.posFraction)] ;
sfx = [sfx '-' 'N' num2str(opts.clusterNum)];
if opts.bboxReg, sfx = [sfx '-' 'bboxreg']; end

if any(opts.inputSize~=500),
    sz = opts.inputSize;
    sfx = [sfx '-input' num2str(sz(1)) 'x' num2str(sz(2))];
end
if ~isempty(opts.clusterName)
    sfx = [sfx '-' 'cluster-' opts.clusterName];
end

opts.expDir = fullfile('models', ['widerface-' sfx]) ;
if ~isempty(opts.tag)
    opts.expDir = [opts.expDir '-' opts.tag];
end
[opts, varargin] = vl_argparse(opts, varargin) ;
opts.expDir

opts.batchSize = 10;
opts.numSubBatches = 1;
opts.numEpochs = 50;
opts.gpus = [1];
opts.numFetchThreads = 8;
opts.lite = false ;
opts.imdbPath = fullfile(opts.dataDir, 'imdb.mat');
opts = vl_argparse(opts, varargin) ;

opts.train = struct() ;
opts = vl_argparse(opts, varargin) ;

opts.train.gpus = opts.gpus;
opts.train.batchSize = opts.batchSize;
opts.train.numSubBatches = opts.numSubBatches;
opts.train.numEpochs = opts.numEpochs;
opts.train.learningRate = opts.learningRate;

opts.train.keepDilatedZeros = opts.keepDilatedZeros;

%% initialize model structure
fprintf('Trying to initialize the structure of %s\n', opts.modelType);
net = cnn_init('model', opts.modelType, ...
               'batchNormalization', opts.batchNormalization, ...
               'weightInitMethod', opts.weightInitMethod, ...
               'networkType', opts.networkType) ;

%% load pretrained weights
if ~isempty(opts.pretrainModelPath)
  fprintf('Loading pretrained weights from %s\n', opts.pretrainModelPath);
  if ~exist(opts.pretrainModelPath)
    [~,model_name,~] = fileparts(opts.pretrainModelPath);
    download_path = sprintf('./trained_models/%s.mat', model_name);
    download_url = sprintf('http://www.vlfeat.org/matconvnet/models/%s.mat', model_name);
    fprintf('Downloading model %s from %s to %s ... \n', model_name, download_url, download_path);
    cmd = sprintf('wget -O %s %s', download_path, download_url);
    system(cmd);
  end
  net = cnn_load_pretrain(net, opts.pretrainModelPath);
end

net.meta.inputSize = opts.inputSize;
net.meta.normalization.inputSize = opts.inputSize;
net.meta.normalization.border = opts.border;
net.meta.augmentation.transformation = 'none';
net.meta.augmentation.rgbVariance = [];

if ~exist(opts.expDir), mkdir(opts.expDir); end;

%% prepare data
if exist(opts.imdbPath)
    imdb = load(opts.imdbPath) ;
    fprintf('Loaded imdb from %s\n', opts.imdbPath);
else
    imdb = cnn_setup_imdb('dataDir', opts.dataDir);
    save(opts.imdbPath, '-struct', 'imdb') ;
    fprintf('Saved imdb to %s\n', opts.imdbPath);
end

%% save model options
optpath = fullfile(opts.expDir, 'opts.mat');
if ~exist(optpath), save(optpath, 'opts'); end

%% define batch getter function
%batchGetter = @cnn_get_batch_logistic_zoom;
if ~isempty(opts.batchGetterFn)
    batchGetter = str2func(opts.batchGetterFn);
end

%% compute image stats
imageStatsPath = fullfile(opts.dataDir, 'imageStats.mat') ;
if exist(imageStatsPath)
    load(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
else
    [averageImage, rgbMean, rgbCovariance] = getImageStats(batchGetter, ...
                                                      opts, net.meta, imdb) ;
    save(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
end
net.meta.augmentation.transformation = 'f5';
net.meta.normalization.averageImage = rgbMean ;
[v,d] = eig(rgbCovariance) ;
net.meta.augmentation.rgbVariance = 0.1*sqrt(d)*v' ;
clear v d ;

%% clustering
minh = opts.minClusterSize(1); minw = opts.minClusterSize(2);
maxh = opts.maxClusterSize(1); maxw = opts.maxClusterSize(2);
if isempty(opts.clusterName)
    tmp_str = ['%s_N%d_' sprintf('min%dx%d_max%dx%d.mat',minh,minw,maxh,maxw)];
else
    tmp_str = ['%s_N%d_' opts.clusterName '.mat'];
end
clusterName = @(name,nd)(sprintf(tmp_str, name, nd));

%clusterName = @(name,nd,minh,minw,maxh,maxw)(sprintf( ...
%    '%s_N%d_min%dx%d_max%dx%d.mat',name,nd,minh,minw,maxh,maxw));
%clusterPath = fullfile(opts.dataDir, clusterName(...
%    'RefBox',opts.clusterNum,minh,minw,maxh,maxw));
boxName = 'RefBox';
clusterPath = fullfile(opts.dataDir, clusterName(boxName,opts.clusterNum));

fprintf('cluster path: %s\n', clusterPath);

if ~exist(clusterPath)
  clusters = cluster_rects(imdb, opts.clusterNum, [minh minw], [maxh maxw]);
  save(clusterPath, 'clusters');
else
  load(clusterPath);
end
net.meta.clusters = clusters;

%% add predictors/losses
switch opts.modelType
  case 'resnet-101-simple'
    net = cnn_add_loss_fcn8s_resnet101_simple(opts, net);
  otherwise
    error(sprintf('Not Implemented: model type %s', opts.modelType));
end

%% compute receptive fields and canonical variable sizes
var2idx = containers.Map;
for i = 1:numel(net.vars)
    var2idx(net.vars(i).name) = i;
end
net.meta.var2idx = var2idx;

sz_ = opts.inputSize;
net.meta.varsizes = net.getVarSizes({'data',[sz_,3,1]});
net.meta.recfields = net.getVarReceptiveFields('data');

%% configure sampling hyperparameters
net.meta.sampleSize = opts.sampleSize;
net.meta.posFraction = opts.posFraction;
net.meta.posThresh = opts.posThresh;
net.meta.negThresh = opts.negThresh;

%% assign training function
trainFn = str2func(opts.trainFn);

%% print out training options
opts
opts.train

%% start training (no validation)
derOutputs = {'loss_cls', 1, 'loss_reg', 1};
[net, info] = trainFn(net, imdb, getBatchFn(batchGetter, opts, net.meta), ...
                      'expDir', opts.expDir,...
                      'derOutputs', derOutputs, ...
                      'val', nan,...
                      opts.train) ;

%% wrapper for batch getter function
function fn = getBatchFn(batchGetter, opts, meta)
useGpu = numel(opts.train.gpus) > 0 ;

bopts.numThreads = opts.numFetchThreads ;
bopts.inputSize = meta.normalization.inputSize ;
bopts.border = meta.normalization.border ;
bopts.averageImage = meta.normalization.averageImage ;
bopts.rgbVariance = meta.augmentation.rgbVariance ;
bopts.transformation = meta.augmentation.transformation ;

if isfield(meta, 'sampleSize'), bopts.sampleSize = meta.sampleSize; end;
if isfield(meta, 'posFraction'), bopts.posFraction = meta.posFraction; end;
if isfield(meta, 'posThresh'), bopts.posThresh = meta.posThresh; end;
if isfield(meta, 'negThresh'), bopts.negThresh = meta.negThresh; end;

if isfield(meta, 'clusters'), bopts.clusters = meta.clusters; end;
if isfield(meta, 'var2idx'), bopts.var2idx = meta.var2idx; end;
if isfield(meta, 'varsizes'), bopts.varsizes = meta.varsizes; end;
if isfield(meta, 'recfields'), bopts.recfields = meta.recfields; end;
switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(batchGetter, bopts,x,y) ;
  case 'dagnn'
    fn = @(x,y) getDagNNBatch(batchGetter, bopts,useGpu,x,y) ;
end

%% interface to batch getter function
function inputs = getDagNNBatch(batchGetter, opts, useGpu, imdb, batch)
imagePaths = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
imageSizes = imdb.images.size(batch, :);
%isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;
isVal = 0;
if isfield(opts, 'clusters') && ~isVal
    labelRects = imdb.labels.rects(batch);
else
    labelRects = [];
end

[images, clsmaps, regmaps] = batchGetter(imagePaths, imageSizes, labelRects, ...
                                       opts, 'prefetch', nargout == 0) ;

if nargout > 0
    % if we are training
    if useGpu && ~isempty(clsmaps) && ~isempty(regmaps)
        images = gpuArray(images) ;
    end
    inputs = {'data', images};

    if ~isempty(clsmaps)
        inputs(end+1:end+2) = {'label_cls', clsmaps};
    end
    if ~isempty(regmaps)
        inputs(end+1:end+2) = {'label_reg', regmaps};
    end
end

%% function to collect image statistics
function [averageImage, rgbMean, rgbCovariance] = ...
    getImageStats(batchGetter, opts, meta, imdb)
train = find(imdb.images.set == 1) ;
train = train(1: 101: end);
bs = 256 ;
opts.networkType = 'dagnn' ;
fn = getBatchFn(batchGetter, opts, meta) ;
avg = {}; rgbm1 = {}; rgbm2 = {};

for t=1:bs:numel(train)
    batch_time = tic ;
    batch = train(t:min(t+bs-1, numel(train))) ;
    fprintf('collecting image stats: batch starting with image %d ...', batch(1)) ;
    temp = fn(imdb, batch) ;
    temp = gather(temp{2});
    z = reshape(permute(temp,[3 1 2 4]),3,[]) ;
    n = size(z,2) ;
    avg{end+1} = mean(temp, 4) ;
    rgbm1{end+1} = sum(z,2)/n ;
    rgbm2{end+1} = z*z'/n ;
    batch_time = toc(batch_time) ;
    fprintf(' %.2f s (%.1f images/s)\n', batch_time, numel(batch)/ batch_time) ;
end
averageImage = mean(cat(4,avg{:}),4) ;
rgbm1 = mean(cat(2,rgbm1{:}),2) ;
rgbm2 = mean(cat(3,rgbm2{:}),3) ;
rgbMean = rgbm1 ;
rgbCovariance = rgbm2 - rgbm1*rgbm1' ;

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function sphere_grid_test10 ( )

%*****************************************************************************80
%
%% SPHERE_GRID_TEST10 tests SPHERE_CUBED_POINTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_GRID_TEST10\n' );
  fprintf ( 1, '  SPHERE_CUBED_POINTS computes points on a cubed sphere grid.\n' );

  n = 10;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of divisions on each face = %d\n', n );

  ns = sphere_cubed_point_num ( n );
  fprintf ( 1, '  Total number of points = %d\n', ns );

  xyz = sphere_cubed_points ( n, ns );

  r8mat_transpose_print_some ( 3, ns, xyz, 1, 1, 3, 20, ...
    '  Initial part of XYZ array:' );
%
%  Display the nodes on a sphere.
%
  figure ( )
  clf
  hold on
  [ x, y, z ] = sphere ( 20 );
  c = ones ( size ( z ) );
  surf ( x, y, z, c );
  plot3 ( xyz(1,:), xyz(2,:), xyz(3,:), 'b.', 'Markersize', 20 );
  axis equal
  grid on
  view ( 3 )
  xlabel ( '<--X-->' )
  ylabel ( '<--Y-->' )
  zlabel ( '<--Z-->' )
  title ( 'Cubed sphere points' );
  hold off
  filename = sprintf ( 'sphere_cubed_points_f%d.png', n );
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file saved to "%s".\n', filename );
%
%  Write the nodes to a file.
%
  if ( 1 )

    filename = sprintf ( 'sphere_cubed_grid_f%d.xyz', n );

    file_unit = fopen ( filename, 'wt' );
    for node = 1 : ns
      for i = 1 : 3
        fprintf ( file_unit, '  %f', xyz(i,node ) );
      end
      fprintf ( file_unit, '\n' );
    end
    fclose ( file_unit );

    fprintf ( 1, '  Wrote grid nodes to "%s".\n', filename );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Example 7.4: Binary hypothesis testing
% Figure 7.4
% Boyd & Vandenberghe "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Michael Grant, 2005-12-19

% Generate the data
P = [0.70  0.10
     0.20  0.10
     0.05  0.70
     0.05  0.10];
[n,m] = size(P);

% Construct the tradeoff curve by finding the
% the Pareto optimal deterministic detectors,
% which are the curve's vertices

nopts   = 1000;
weights = logspace(-5,5,nopts);
obj     = [0;1];
inds    = ones(n,1);

% minimize  -t1'*q1 - w*t2'*q2
% s.t.      t1+t2 = 1,  t1,t2 \geq 0

next = 2;
for i = 1 : nopts,
   PW = P * diag( [ 1 ; weights(i) ] );
   [ maxvals, maxinds ] = max( PW' );  % max elt in each row
   if (~isequal(maxinds', inds(:,next-1)))
       inds(:,next) = maxinds';
       T = zeros(m,n);
       for j=1:n
          T(maxinds(1,j),j) = 1;
       end;
       obj(:,next) = 1-diag(T*P);
       next = next+1;
   end;
end;
plot(obj(1,:), obj(2,:),[0 1], [0 1],'--');
grid on
for i=2:size(obj,2)-1
   text(obj(1,i),obj(2,i),['a', num2str(i-1)]);
end;

% Minimax detector: not deterministic

cvx_begin
    variables T( m, n ) D( m, m )
    minimize max( D(1,2), D(2,1) )
    subject to
        D == T * P;
        sum( T, 1 ) == 1;
        T >= 0;
cvx_end

objmp = 1 - diag( D );
text( objmp(1), objmp(2), 'b' );
xlabel('P_{fp}'); ylabel('P_{fn}');

%print -deps roc.eps

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function z=orthoAlt4Pres(P,P0,T0)
%%ORTHOALT4PRES Get the orthometric height (pressure altitude) for a
%               point up to 86km altitude using the U.S. Standard
%               Atmosphere 1976 model given a reading of the air pressure
%               at the point. The model is built into many altimeters.
%
%INPUTS: P A vector or matrix of the pressures in Pascals that are to be
%          converted into orthometric heights in meters.
%       P0 An optional parameter specifying the pressure at sea level at
%          the location in question in Pascals. This pressure might have
%          to be interpolated down below terrain given a measurement on
%          the terrain. If omitted, the value of
%          Constants.standardAtmosphericPressure is used for the pressure.
%          Except in certain instances, aircraft over 18,000ft altitude in
%          the United States have to have their altimeter calibrated to
%          the standard atmospheric pressure and not the true pressure on
%          the ground/ at sea level.
%       T0 An optional parameter specifying the temperature at sea level
%          in degrees Kelvin. The standard atmosphere only specifies one
%          value, which is in Constants.standardTemp and is the default
%          value used if this parameter is omitted.
%
%OUTPUTS: z The orthometric height of the point implied by the specified air
%           pressure.The standard does not specify a reference geoid/
%           permanent tide system.
%
%The atmospheric model is defined in [1]. and the full model goes up to
%1000km altitude. Up to an altitude of 32km, the U.S. standard atmosphere
%is identical to the International Civil Aviation Organization's (ICAO)
%standard atmosphere. Regulations related to how altimeters should be set
%with respect to a standard pressure are at
%https://www.faa.gov/air_traffic/publications/atpubs/aim/aim0702.html
%
%REFERENCES:
%[1] U.S. Standard Atmosphere, 1976, National Oceanic and Atmospheric
%   Administration Std., Oct. 1976.
%
%June 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    P0=Constants.standardAtmosphericPressure;%N/m
end

if(nargin<3)
    T0=Constants.standardTemp;
end

%%%Constants defined in the standard atmosphere.
%The radius of the spherical Earth used in the standard atmosphere. THis is
%not the same as the semi-major axis of the WGS-84 reference ellipsoid.
rE=6356766;%Meters
%The 1962 estimate of the UNiversal Gas Constant that is used in the
%standard.
R=8.31432;%J/(mole*K)
%The mean moleacular mass of air at sea level.
M0=28.9644e-3;%kg/mole
%Standard acceleration due to gravity.
g0=9.80665;%m/s

%Hb and LMb are from table 4 in the U.S. Standard Atmosphere 1976.
%Hb is the base geopotential height and LMb is the molecular scale
%temperature gradient from that base.
HbTable=[0;
    11000;
    20000;
    32000;
    47000;
    51000;
    71000;
    84852];
LMbTable=[-6.5e-3;
    0;
    1e-3;
    2.8e-3;
    0;
    -2.8e-3;
    -2e-3;
    0];

%TMb is computed using Equation 23 starting from T0. It isn't just
%tabulated in case the user wants to specify a non-standard T0.
TMbTable=zeros(8,1);
TMbTable(1)=T0;
for n=2:8
    TMbTable(n)=TMbTable(n-1)+LMbTable(n-1)*(HbTable(n)-HbTable(n-1));
end

numP=size(P(:));
z=zeros(size(P));
for curP=1:numP
    PCur=P(curP);

    %To determine the correct formula to use, we must first determine which
    %band of altitudes the pressure falls. Thus, we will first evaluate
    %pres4OrthoAlt for the pressure at all of the bands that separate, knowing
    %that the pressure monotonically decreases with altitude in the model.
    H=11000;
    zt=rE*(H/(rE-H));
    [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
    if(PCur<Ph)
        H=20000;
        zt=rE*(H/(rE-H));
        [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
        if(PCur<Ph)
            H=32000;
            zt=rE*(H/(rE-H));
            [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
            if(PCur<Ph)
                H=47000;
                zt=rE*(H/(rE-H));
                [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                if(PCur<Ph)
                    H=51000;
                    zt=rE*(H/(rE-H));
                    [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                    if(PCur<Ph)
                        H=71000;
                        zt=rE*(H/(rE-H));
                        [Ph,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                        if(PCur<Ph)
                            H=72000;
                            zt=rE*(H/(rE-H));
                            [~,~,Hb,Pb,TMb,LMb]=presTemp4OrthoAlt(zt,P0,T0);
                        end
                    end
                end
            end
        end
    end

    %Solve the inverse of the appropriate air pressure equation to determine
    %the geopotential height.
    if(LMb==0)
        %The inverse of the exponential expression for air pressure as a
        %function of altitude in Equation 33b
        H=Hb-R*TMb*log(PCur/Pb)/(g0*M0);
    else
        %The inverse of air pressure according to equation 33a, which is the
        %linear formulation.
        H=Hb+(TMb/LMb)*((PCur/Pb)^(-LMb*R/(g0*M0))-1);
    end

    %Convert the geopotential height into an orthometric height using the
    %approximation in the standard.
    z(curP)=rE*(H/(rE-H));
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef nnbnorm < nntest
  properties (TestParameter)
    rows = {2 8 13}
    cols = {2 8 17}
    numDims = {1 3 4}
    batchSize = {2 7}
  end
  methods (Test)
    function regression(test, rows,cols, numDims, batchSize)
      r = rows ;
      c = cols ;
      nd = numDims ;
      bs = batchSize ;
      x = test.randn(r, c, nd, bs) ;
      g = test.randn(nd, 1) / test.range ;
      b = test.randn(nd, 1) / test.range ;
      epsilon = 0.001 ;

      [y,m] = vl_nnbnorm(x,g,b,'epsilon',epsilon) ;
      n = numel(x) / size(x,3) ;
      mu = sum(sum(sum(x,1),2),4) / n ;
      sigma2 = sum(sum(sum(bsxfun(@minus,x,mu).^2,1),2),4) / n + epsilon ;
      sigma = sqrt(sigma2) ;
      m_ = [mu(:),sigma(:)] ;
      test.eq(m,m_) ;

      g = reshape(g,1,1,[]) ;
      b = reshape(b,1,1,[]) ;
      y_ = bsxfun(@plus,b,bsxfun(@times,g,bsxfun(@rdivide,bsxfun(@minus,x,mu),sigma))) ;
      test.eq(y,y_) ;
    end

    function basic(test, rows, cols, numDims, batchSize)
      r = rows ;
      c = cols ;
      nd = numDims ;
      bs = batchSize ;
      x = test.randn(r, c, nd, bs) ;
      %g = test.randn(1, 1, nd, 1) ;
      %b = test.randn(1, 1, nd, 1) ;
      g = test.randn(nd, 1) / test.range ;
      b = test.randn(nd, 1) / test.range ;

      y = vl_nnbnorm(x,g,b) ;
      dzdy = test.randn(size(y)) ;
      [dzdx,dzdg,dzdb] = vl_nnbnorm(x,g,b,dzdy) ;

      test.der(@(x) vl_nnbnorm(x,g,b), x, dzdy, dzdx, test.range * 1e-3) ;
      test.der(@(g) vl_nnbnorm(x,g,b), g, dzdy, dzdg, 1e-2) ;
      test.der(@(b) vl_nnbnorm(x,g,b), b, dzdy, dzdb, 1e-3) ;
    end

    function givenMoments(test, rows, cols, numDims, batchSize)
      r = rows ;
      c = cols ;
      nd = numDims ;
      bs = batchSize ;
      x = test.randn(r, c, nd, bs) ;
      %g = test.randn(1, 1, nd, 1) ;
      %b = test.randn(1, 1, nd, 1) ;
      g = test.randn(nd, 1) / test.range ;
      b = test.randn(nd, 1) / test.range ;

      [y,m] = vl_nnbnorm(x,g,b) ;
      [y_,m_] = vl_nnbnorm(x,g,b,'moments',m) ;

      test.eq(y,y_) ;
      test.eq(m,m_) ;

      dzdy = test.randn(size(y)) ;
      [dzdx,dzdg,dzdb,m__] = vl_nnbnorm(x,g,b,dzdy) ;
      [dzdx_,dzdg_,dzdb_,m___] = vl_nnbnorm(x,g,b,dzdy,'moments',m) ;

      test.eq(m,m__)
      test.eq(m,m___)
      test.eq(dzdx,dzdx_) ;
      test.eq(dzdg,dzdg_) ;
      test.eq(dzdb,dzdb_) ;
    end

  end
end

repo: UDT


text:
function t = uminus(t)
%UMINUS Unary minus for ttensor.
%
%   See also TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

t.core = -t.core;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function test085 ( dim_num, n, z )

%*****************************************************************************80
%
%% TEST085 tests R0_MEASURE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST085\n' );
  fprintf ( 1, '  R0_MEASURE computes the R0 measure of quality.\n' );
  fprintf ( 1, '  The Riesz S = 0 energy measure  R0 = %14f\n', ...
    r0_measure ( dim_num, n, z ) );

  return
end


text:
% read_erplabevent() - import ERPLAB dataset events
%
% Usage:
%     >> event = read_erplabevent(filename, ...);
%
% Inputs:
%     filename - [string] file name
%
% Optional inputs:
%     'header' - FILEIO structure header
%
% Outputs:
%     event - FILEIO toolbox event structure
%
% Modified from read_eeglabevent

%123456789012345678901234567890123456789012345678901234567890123456789012
%
% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

function event = read_erplabevent(filename, varargin)

if nargin < 1
  help read_erplabheader;
  return;
end;

hdr = ft_getopt(varargin, 'header');

if isempty(hdr)
  hdr = read_erplabheader(filename);
end

event    = [];                % these will be the output in FieldTrip format
oldevent = hdr.orig.bindescr;    % these are in ERPLAB format

for index = 1:length(oldevent)
    event(end+1).type     = 'trial';
    event(end  ).sample   = (index-1)*hdr.nSamples + 1;
    event(end  ).value    = oldevent{index};
    event(end  ).offset   = -hdr.nSamplesPre;
    event(end  ).duration =  hdr.nSamples;
end;

author: spm


text:
function msm_to_mm_test16 ( )

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

  output_filename = 'msm_to_mm_test16.mm';

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

  return
end


text:
function ellipsoid_monte_carlo_test03 ( )

%*****************************************************************************80
%
%% ELLIPSOID_MONTE_CARLO_TEST03 uses ELLIPSOID_SAMPLE on a 3D ellipse centered at (1,2,3).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2014
%
%  Author:
%
%    John Burkardt
%
  m = 3;

  a = [ ...
    9.0, 6.0, 3.0; ...
    6.0, 5.0, 4.0; ...
    3.0, 4.0, 9.0 ]';
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    0, 0, 3 ]';
  r = 0.5;
  v = [ 1.0, 2.0, 3.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSOID_MONTE_CARLO_TEST03\n' );
  fprintf ( 1, '  Use ELLIPSOID_SAMPLE to estimate integrals\n' );
  fprintf ( 1, '  in a 3D ellipse (x-v)'' * A * (x-v) <= r^2.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ellipsoid radius R = %g\n', r );
  r8vec_print ( m, v, '  Ellipsoid center V:' )
  r8mat_print ( m, m, a, '  Ellipsoid matrix A:' )

  volume = ellipsoid_volume ( m, a, v, r );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ellipsoid volume = %g\n', volume );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N        1              X               Y  ' );
  fprintf ( 1, '              Z                X^2            YZ              Z^3\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = ellipsoid_sample ( m, n, a, v, r, seed );

    fprintf ( 1, '  %8d', n );

    for j = 1 : 7

      e(1:m) = e_test(1:m,j);

      value = monomial_value ( m, n, e, x );

      result = volume * sum ( value(1:n) ) / n;

      fprintf ( 1, '  %14.6g', result );

    end

    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  return
end


text:
function Clg=cct2clg(Cct,lat,lon)
% CCT2CLG  Convert CT covariance matrix to local geodetic.
%   Vectorized. Requires rotct2lg. See also CLG2CCT.
% Version: 2011-04-05
% Useage:  Clg=cct2clg(Cct,lat,lon)
% Input:   Cct - CT covariance matrix
%          lat - vector of station latitudes (rad)
%          lon - vector of station longitudes (rad)
% Output:  Clg - LG covariance matrix

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

n=length(lat);
if (n*3 ~= max(size(Cct)) )
  error('Size of lat,lon does not match size of Cct');
end

for i=1:n
  sinlat=sin(lat(i));
  coslat=cos(lat(i));
  sinlon=sin(lon(i));
  coslon=cos(lon(i));
  Ji=rotct2lg(lat(i),lon(i));
  indi=(i-1)*3+[1:3];
  for j=1:n
    sinlat=sin(lat(j));
    coslat=cos(lat(j));
    sinlon=sin(lon(j));
    coslon=cos(lon(j));
    Jj=rotct2lg(lat(j),lon(j));
    indj=(j-1)*3+[1:3];
    Clg(indi,indj)=Ji*Cct(indi,indj)*Jj';
  end
end

author: Sable

repo: mcbench-benchmarks


text:
%

% @
% Copyright (C) 2017 Jonas Ruesch
%
% 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/>.
% @

function [ac, ... % the aerodynamic center of the aircraft in meters measured from the Tornado frame of reference origin.
    staticMarginMAC, ... % the static margin of the aircraft with the given center of gravity as a factor of the mean chord.
    staticMarginInMeters]...
    = computeStaticMargin(...
nameTornadoAircraftGeometry, ...
estimatedVelocity, ...
airDensity, ...
tornadoDirectory, ...
centerOfGravity)

% Define an initial Tornado state with given estimated velocity.
state.AS=      estimatedVelocity; %Airspeed m/s
state.alpha=   0.0;       %Angle of attack, radians
state.betha=   0;         %Angle of sideslip, radians
state.P=       0;         %Rollrate, rad/s
state.Q=       0;         %pitchrate, rad/s
state.R=       0;         %yawrate, rad/s
state.adot=    0;         %Alpha time derivative rad/s
state.bdot=    0;         %Betha time derivative rad/s
state.ALT=     0;         %Altitude, m
state.rho=     airDensity;%Desity, kg/m^3
state.pgcorr=  0;         %Apply prandtl glauert compressibility correction

% Get some properties from the aircraft geometry
cd(tornadoDirectory);
settings=config('startup');
% Load the aircraft
cd(settings.acdir);
load(nameTornadoAircraftGeometry); % aircraft geometry is now in variable 'geo'.
% Overwrite center of gravity with given value (!)
geo.CG = [centerOfGravity 0 0];
% Overwrite reference point for moment computation with center of gravity.
geo.ref_point = [centerOfGravity 0 0];

% Find static margin
% Since aerodynamic center for entire aircraft should correspond with
% neutral point, the static margin should add up (cg + margin = neutral
% point = aerodynamic center entire plane).
cd(tornadoDirectory);
out=fFindstaticmargin(geo,state);
staticMarginMAC = out.h(1);
mean_chord = geo.c(1);
staticMarginInMeters = out.h(1)*mean_chord;
ac = out.ac(1);

cd(fileparts(which(mfilename))); % Back to the directory of this script.

end

author: jrgenerative

repo: fixed-wing-sim


text:
function sendmail = spm_cfg_sendmail
% SPM Configuration file for sendmail
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_sendmail.m 6929 2016-11-14 13:07:31Z guillaume $

% ---------------------------------------------------------------------
% Recipient
% ---------------------------------------------------------------------
recipient         = cfg_entry;
recipient.tag     = 'recipient';
recipient.name    = 'Recipient';
recipient.help    = {'User to receive mail.'};
recipient.strtype = 's';
recipient.num     = [1 Inf];
% ---------------------------------------------------------------------
% Subject
% ---------------------------------------------------------------------
subject         = cfg_entry;
subject.tag     = 'subject';
subject.name    = 'Subject';
subject.val     = {['[SPM] [%DATE%] On behalf of ' spm('Ver')]};
subject.help    = {'The subject line of the message. %DATE% will be replaced by a string containing the time and date when the email is sent.'};
subject.strtype = 's';
subject.num     = [1 Inf];
% ---------------------------------------------------------------------
% Message
% ---------------------------------------------------------------------
message         = cfg_entry;
message.tag     = 'message';
message.name    = 'Message';
message.val     = {'Hello from SPM!'};
message.help    = {'A string containing the message to send.'};
message.strtype = 's';
message.num     = [1 Inf];
% ---------------------------------------------------------------------
% Attachments
% ---------------------------------------------------------------------
attachments         = cfg_files;
attachments.tag     = 'attachments';
attachments.name    = 'Attachments';
attachments.val     = {{}};
attachments.help    = {'List of files to attach and send along with the message.'};
attachments.filter  = '.*';
attachments.ufilter = '.*';
attachments.num     = [0 Inf];
% ---------------------------------------------------------------------
% SMTP Server
% ---------------------------------------------------------------------
smtp         = cfg_entry;
smtp.tag     = 'smtp';
smtp.name    = 'SMTP Server';
smtp.help    = {'Your SMTP server. If not specified, look for sendmail help.'};
smtp.strtype = 's';
try
    smtp.val = {getpref('Internet','SMTP_Server')};
end
smtp.num     = [1 Inf];
% ---------------------------------------------------------------------
% E-mail
% ---------------------------------------------------------------------
email         = cfg_entry;
email.tag     = 'email';
email.name    = 'E-mail';
email.help    = {'Your e-mail address. Look in sendmail help how to store it.'};
email.strtype = 's';
try
    email.val = {getpref('Internet','E_mail')};
end
email.num     = [1 Inf];
% ---------------------------------------------------------------------
% Zip attachments
% ---------------------------------------------------------------------
zip         = cfg_menu;
zip.tag     = 'zip';
zip.name    = 'Zip attachments';
zip.val     = {'No'};
zip.help    = {'Zip attachments before being sent along with the message.'};
zip.labels  = {'Yes' 'No'}';
zip.values  = {'Yes' 'No'}';
% ---------------------------------------------------------------------
% Parameters
% ---------------------------------------------------------------------
params         = cfg_branch;
params.tag     = 'params';
params.name    = 'Parameters';
params.val     = { smtp email zip};
params.help    = {'Preferences for your e-mail server (Internet SMTP server) and your e-mail address. MATLAB tries to read the SMTP mail server from your system registry. This should work flawlessly. If you encounter any error, identify the outgoing mail server for your electronic mail application, which is usually listed in the application''s preferences, or, consult your e-mail system administrator, and update the parameters. Note that this function does not support e-mail servers that require authentication.'};
% ---------------------------------------------------------------------
% Sendmail
% ---------------------------------------------------------------------
sendmail       = cfg_exbranch;
sendmail.tag   = 'sendmail';
sendmail.name  = 'Sendmail';
sendmail.val   = { recipient subject message attachments params};
sendmail.help  = {'Send a mail message (attachments optionals) to an address.'};
sendmail.prog  = @spm_sendmail;
%_______________________________________________________________________

%_______________________________________________________________________
function spm_sendmail(job)

try
    setpref('Internet','SMTP_Server',job.params.smtp);
    setpref('Internet','E_mail',job.params.email);

    subj = strrep(job.subject,'%DATE%',datestr(now));
    mesg = strrep(job.message,'%DATE%',datestr(now));
    mesg = [mesg 10 10 '-- ' 10 10 'Statistical Parametric Mapping'];

    if ~isempty(job.attachments)
        if strcmpi(job.params.zip,'Yes')
            zipfile = fullfile(tempdir,'spm_sendmail.zip');
            zip(zipfile,job.attachments);
            job.attachments = {zipfile};
        end
        sendmail(job.recipient,subj,mesg,job.attachments);
    else
        sendmail(job.recipient,subj,mesg);
    end
catch
    %- not an error to prevent an analysis to crash because of just that...
    fprintf('Sendmail failed...\n');
end


text:
%Script to compile mex files stored in "private" folder
%Written by: Prateek Jain (pjain@cs.utexas.edu) and Raghu Meka (raghu@cs.utexas.edu)
%Last updated: October, 2009

cd private;
%mex -O -lmwlapack -lmwblas bdsqr_mex.c dbdqr.f -UDBDQR_IN_C -output bdsqr
%mex -O -lmwlapack -lmwblas reorth_mex.c reorth.f -output reorth
%mex -O -lmwlapack -lmwblas tqlb_mex.c tqlb.f -output tqlb
mex compute_X_Omega.c
mex compOmegaYx.c
mex compOmegaYtx.c
cd ..

author: andrewssobral

repo: lrslibrary


text:
function f2 = p02_f2 ( x )

%*****************************************************************************80
%
%% P02_F2 evaluates the second derivative for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    LE Scales,
%    Introduction to Non-Linear Optimization,
%    Springer, 1985.
%
%  Parameters:
%
%    Input, real X, the values of the variables.
%
%    Output, real F2, the second derivative.
%
  f2 = 2.0 + exp ( -x );

  return
end


text:
function pdf = cosine_pdf ( x, a, b )

%*****************************************************************************80
%
%% COSINE_PDF evaluates the Cosine PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) = ( 1 / ( 2 * PI * B ) ) * COS ( ( X - A ) / B )
%    for A - PI * B <= X <= A + PI * B
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < a - pi * b )

    pdf = 0.0;

  elseif ( x <= a + pi * b )

    y = ( x - a ) / b;

    pdf = 1.0 / ( 2.0 * pi * b ) * cos ( y );

  elseif ( a + pi * b < x )

    pdf = 0.0;

  end

  return
end


text:
function Dij = TS_PairwiseDist(tsOrOps,whatData,distanceMetric,doSave)
% TS_PairwiseDist   Pairwise distances between either TimeSeries or Operations
%
% (If doSave), stores pairwise distances back the data back in the local .mat
% file. Warning: this can be huge for large datasets.
%
%---INPUTS:
%
% tsOrOps: Compute pairwise distances between all pairs of time series
%               ('ts', default), or operations ('ops')
%
% whatData: what hctsa data source to use (cf. TS_LoadData)
%
% distanceMetric: what distance metric to use, e.g., 'euclidean' (default for
%                 time series) or 'corr_fast' (default for operations)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
% Check Inputs, set defaults:
% ------------------------------------------------------------------------------

if nargin < 1 || isempty(tsOrOps)
    tsOrOps = 'ts';
    fprintf(1,'Time series by default.\n');
end
if ~ismember(tsOrOps,{'ts','ops'})
    error('Invalid specifier ''%s'' (should be ''ts'' or ''ops'')',tsOrOps);
end

if nargin < 2
    whatData = 'norm';
end

if nargin < 3
    switch tsOrOps
    case 'ts'
        distanceMetric = 'euclidean';
    case 'ops'
        distanceMetric = 'corr_fast';
    end
end

% doSave
if nargin < 4
    % By default, save back to the file specified as whatData
    doSave = 1;
end

% ------------------------------------------------------------------------------
% Load data
% ------------------------------------------------------------------------------

if ischar(whatData)
    [TS_DataMat,~,~,theFile] = TS_LoadData(whatData);
else
    % If provided a matrix, assume it's from HCTSA_N.mat by default:
    TS_DataMat = whatData;
    theFile = 'HCTSA_N.mat';
end

% ------------------------------------------------------------------------------
% Compute pairwise distances between all objects
% ------------------------------------------------------------------------------
if strcmp(tsOrOps,'ops')
    % Transpose if we want columns ('ops')
    TS_DataMat = TS_DataMat';
end

fprintf(1,'Computing %s distances between %u objects...',distanceMetric,size(TS_DataMat,1));
Dij = BF_pdist(TS_DataMat,distanceMetric,1,[],1);
fprintf(1,' Done.\n');

% ------------------------------------------------------------------------------
% Save back to the local hctsa .mat file
% ------------------------------------------------------------------------------
if doSave
    switch tsOrOps
    case 'ts'
        theWhat = 'time series';
    case 'ops'
        theWhat = 'operations';
    end
    fprintf(1,'Saving pairwise distance information for %s back to %s...',...
                        theWhat,theFile);

    switch tsOrOps
    case 'ts'
        ts_clust = struct('Dij',Dij,'distanceMetric',distanceMetric);
        save(theFile,'ts_clust','-append')
    case 'ops'
        op_clust = struct('Dij',Dij,'distanceMetric',distanceMetric);
        save(theFile,'op_clust','-append')
    end

    fprintf(1,' Done.\n');
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

    Update = 0;
    for i=1:m,
        Update = Update + alpha/m * (theta' * X(i,:)' - y(i)) *  X(i,:)';
    end

    theta = theta - Update;

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

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

repo: machine-learning-coursera


text:
function [AccountVal, PayoutVal, FeeVal] = calcValuePayoutAndFees_orig( ...
    Holdings_0, assetPrice, annualWithdrawalRate, annualFee)
% CALCVALUEPAYOUTANDFEES_ORIG estimates the account value, the payout
% amounts, and the fees collected for a Monte Carlo simulation of the GMWB.
%
% Inputs:
% Holdings_0: An nAssets-by-1 vector of initial portfolio holdings, in
%   shares
% assetPrice: An nYears-by-nAssets-by-nTrials matrix of simulated asset
%   prices
% annualWithdrawalRate: A scalar indicating the fraction of the initial
%   portfolio value that the customer is guaranteed to withdraw each year
% annualFee: A scalar indicating the fraction of the current portfolio
%   value that the insurer may assess each year
%
% Outputs: All have size nYears+1-by-nTrials.  The first row is for the
%   initial year (t = 0)
% AccountVal: The total portofolio value after withdrawals.
% PayoutVal: The amount paid by the insurer to the customer if there is a
%   shortfall
% FeeVal: The amount collected by the insurer if funds remain in the
%   account

[nYears, ~, nTrials] = size(assetPrice);

AccountVal  = nan(nYears, nTrials);
PayoutVal   = nan(nYears, nTrials);
FeeVal      = nan(nYears, nTrials);

for iTrial = 1:nTrials
    % Initialize values at time 0.
    currentHoldings      = Holdings_0;
    currentPrice         = assetPrice(1, :, iTrial);
    AccountVal(1,iTrial) = currentPrice * currentHoldings;
    PayoutVal(1, iTrial) = 0;
    FeeVal(   1, iTrial) = 0;
    guaranteedWithdrawal = currentPrice * currentHoldings * annualWithdrawalRate;

    for iYear = 2:nYears
        currentPrice    = assetPrice(iYear, :, iTrial);
        curPortVal      = currentPrice * currentHoldings;
        totalWithdrawal = curPortVal * annualFee + guaranteedWithdrawal;

        FeeVal(iYear, iTrial) = curPortVal * annualFee;

        [currentHoldings, shortfall] = updateHoldings(currentPrice, ...
            currentHoldings, totalWithdrawal);

        PayoutVal(iYear, iTrial) = shortfall;
        AccountVal(iYear, iTrial) = currentPrice * currentHoldings;
    end
end
end

function [holdings, shortfall] = updateHoldings(currentPrice, ...
    currentHoldings, totalWithdrawal)
% Adjusts the current holdings and calculates the shortfall given the
% pre-withdrawal holdings, the withdrawal, and the current asset prices.

holdings = currentHoldings;
currentPrice = currentPrice';
curPortVal = currentHoldings .* currentPrice;

if any(curPortVal > 0)
    % We still have some funds in the portfolio: apply withdrawals in the
    % best way possible by distributing withdrawals across the holdings
    % proportionately to their values:
    newPortVal = curPortVal - curPortVal.*totalWithdrawal/sum(curPortVal);

    if sum(newPortVal) > 0      % Still have fund left in the account
        shortfall = 0;
        holdings = newPortVal ./ currentPrice;
        holdings(~isfinite(holdings)) = 0;
    else
        % We've just gone bankrupt.  Set all holdings to 0 and pay out the
        % shortfall.
        shortfall = -sum(newPortVal);
        holdings(:) = 0;
    end

else
    % Every asset has zero value-- the calculation is straightforward.
    shortfall   = totalWithdrawal;
    holdings(:) = 0;
end

end


text:
function cmap = CorrPlot(coeffs,isPlotText,ylabels)
im = coeffs;
im(isnan(im)) = 0;

% red-white-blue colormap
cmap = zeros(64,3);
cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)];
cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)];
cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)];
minlim = -1; %min(min(im));
maxlim = 1; %max(max(im));

RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap

image(RGB); axis equal; axis tight;
set(gca,'XTick',1:length(im),'YTick',1:length(im));
set(gca,'TickLength',[0,0]);

if exist('isPlotText','var'),
    if isPlotText,
        for i = 1:size(im,1),
            for j = 1:size(im,2),
                if im(j,i)~=0
                    text(i-0.3, j, num2str(round(im(j,i)*100)/100));%, 'Units', 'data') % flipped i/j necessary
                end
            end
        end
    end
end
if exist('ylabels','var'),
    set(gca,'YTickLabel',ylabels);
end
end

function RGB = ImageToRGB(im,cmap,minlim,maxlim)
L = size(cmap,1);
ix = round(interp1(linspace(minlim,maxlim,L),1:L,im,'linear','extrap'));
RGB = reshape(cmap(ix,:),[size(ix) 3]); % Make RGB image from scaled.
end

author: xiuyechen

repo: FishExplorer


text:
function [ node_num, tet_num ] = tet_mesh_order4_example_size ( )

%*****************************************************************************80
%
%% TET_MESH_ORDER4_EXAMPLE_SIZE sizes an example linear tet mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer TET_NUM, the number of tetrahedrons.
%
  node_num = 63;
  tet_num = 144;

  return
end