text:
function lik = lik_binomial(varargin)
%LIK_BINOMIAL  Create a Binomial likelihood structure
%
%  Description
%    LIK = LIK_BINOMIAL creates Binomial likelihood structure.
%
%    The likelihood is defined as follows:
%                  __ n
%      p(y|f, z) = || i=1 [ p_i^(y_i)*(1-p_i)^(z_i-y_i)) *
%                           gamma(z_i+1)/(gamma(y_i+1)*gamma(z_i-y_i+1))]
%    where p_i = exp(f_i)/ (1+exp(f_i)) is the succes probability,
%    which is a function of the latent variable f_i and z is a
%    vector of numbers of trials.
%
%    When using Binomial likelihood you need to give the vector z
%    as an extra parameter to each function that requires y also.
%    For example, you should call gpla_e as follows
%      gpla_e(w, gp, x, y, 'z', z)
%
%  See also
%    GP_SET, LIK_*
%

% Copyright (c) 2009-2010 Jaakko Riihim�ki & Jarno Vanhatalo
% Copyright (c) 2010-2011 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'LIK_BINOMIAL';
  ip.addOptional('lik', [], @isstruct);
  ip.parse(varargin{:});
  lik=ip.Results.lik;

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

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_binomial_pak;
    lik.fh.unpak = @lik_binomial_unpak;
    lik.fh.ll = @lik_binomial_ll;
    lik.fh.llg = @lik_binomial_llg;
    lik.fh.llg2 = @lik_binomial_llg2;
    lik.fh.llg3 = @lik_binomial_llg3;
    lik.fh.tiltedMoments = @lik_binomial_tiltedMoments;
    lik.fh.predy = @lik_binomial_predy;
    lik.fh.predprcty = @lik_binomial_predprcty;
    lik.fh.invlink = @lik_binomial_invlink;
    lik.fh.recappend = @lik_binomial_recappend;
  end

end

function [w,s] = lik_binomial_pak(lik)
%LIK_BINOMIAL_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_BINOMIAL_PAK(LIK) takes a likelihood structure LIK
%    and returns an empty verctor W. If Binomial likelihood had
%    parameters this would combine them into a single row vector
%    W (see e.g. likelih_negbin). This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%  See also
%    LIK_NEGBIN_UNPAK, GP_PAK

  w = []; s = {};
end

function [lik, w] = lik_binomial_unpak(lik, w)
%LIK_BINOMIAL_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_BINOMIAL_UNPAK(W, LIK) Doesn't do anything.
%
%    If Binomial likelihood had parameters this would extracts
%    them parameters from the vector W to the LIK structure.
%    This is a mandatory subfunction used for example in energy
%    and gradient computations.
%
%  See also
%    LIK_BINOMIAL_PAK, GP_UNPAK

  lik=lik;
  w=w;

end

function ll = lik_binomial_ll(lik, y, f, z)
%LIK_BINOMIAL_LL  Log likelihood
%
%  Description
%    LL = LIK_BINOMIAL_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, succes counts Y, numbers of trials 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_BINOMIAL_LLG, LIK_BINOMIAL_LLG3, LIK_BINOMIAL_LLG2, GPLA_E

  if isempty(z)
    error(['lik_binomial -> lik_binomial_ll: missing z!'...
           'Binomial likelihood needs the expected number of   '...
           'occurrences as an extra input z. See, for         '...
           'example, lik_binomial and gpla_e.             ']);
  end

  expf = exp(f);
  p = expf ./ (1+expf);
  N = z;
  ll =  sum(gammaln(N+1)-gammaln(y+1)-gammaln(N-y+1)+y.*log(p)+(N-y).*log(1-p));
end

function llg = lik_binomial_llg(lik, y, f, param, z)
%LIK_BINOMIAL_LLG    Gradient of the log likelihood
%
%  Description
%    LLG = LIK_BINOMIAL_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, succes counts Y, numbers of trials 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_BINOMIAL_LL, LIK_BINOMIAL_LLG2, LIK_BINOMIAL_LLG3, GPLA_E

  if isempty(z)
    error(['lik_binomial -> lik_binomial_llg: missing z!'...
           'Binomial likelihood needs the expected number of   '...
           'occurrences as an extra input z. See, for         '...
           'example, lik_binomial and gpla_e.             ']);
  end

  switch param
    case 'latent'
      expf = exp(f);
      N = z;

      llg = y./(1+expf) - (N-y).*expf./(1+expf);
  end
end

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

  if isempty(z)
    error(['lik_binomial -> lik_binomial_llg2: missing z!'...
           'Binomial likelihood needs the expected number of    '...
           'occurrences as an extra input z. See, for          '...
           'example, lik_binomial and gpla_e.              ']);
  end

  switch param
    case 'latent'
      expf = exp(f);
      N = z;

      llg2 = -N.*expf./(1+expf).^2;
  end
end

function llg3 = lik_binomial_llg3(lik, y, f, param, z)
%LIK_BINOMIAL_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_BINOMIAL_LLG3(LIK, Y, F, PARAM) takes a
%    likelihood structure LIK, succes counts Y, numbers of trials
%    Z and latent values F and returns the third gradients of the
%    log likelihood with respect to PARAM. At the moment PARAM
%    can be only 'latent'. G3 is a vector with third gradients.
%    This subfunction is needed when using Laplace appoximation
%    for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_BINOMIAL_LL, LIK_BINOMIAL_LLG, LIK_BINOMIAL_LLG2, GPLA_E, GPLA_G

  if isempty(z)
    error(['lik_binomial -> lik_binomial_llg3: missing z!'...
           'Binomial likelihood needs the expected number of    '...
           'occurrences as an extra input z. See, for          '...
           'example, lik_binomial and gpla_e.              ']);
  end

  switch param
    case 'latent'
      expf = exp(f);
      N = z;
      llg3 = N.*(expf.*(expf-1))./(1+expf).^3;
  end
end

function [logM_0, m_1, sigm2hati1] = lik_binomial_tiltedMoments(lik, y, i1, sigm2_i, myy_i, z)
%LIK_BINOMIAL_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_BINOMIAL_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, succes counts Y,
%    numbers of trials 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_binomial -> lik_binomial_tiltedMoments: missing z!'...
%           'Binomial likelihood needs the expected number of               '...
%           'occurrences as an extra input z. See, for                     '...
%           'example, lik_binomial and gpla_e.                         ']);
%  end

  yy = y(i1);
  N = z(i1);
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    % Create function handle for the function to be integrated
    % (likelihood * cavity) and useful integration limits
    [tf,minf,maxf]=init_binomial_norm(yy(i),myy_i(i),sigm2_i(i),N(i));

    % 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 = m_2 - m_1(i).^2;
      %    if sigm2hati1 >= sigm2_i
      %      error('lik_binomial_tilted_moments: sigm2hati1 >= sigm2_i');
      %    end
    end
    logM_0(i) = log(m_0);
  end
end

function [lpy, Ey, Vary] = lik_binomial_predy(lik, Ef, Varf, yt, zt)
%LIK_BINOMIAL_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    [LPY] = LIK_BINOMIAL_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | y, zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires also the succes counts YT, numbers of trials ZT.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_BINOMIAL_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
%    GPEP_PRED, GPLA_PRED, GPMC_PRED

  if isempty(zt)
    error(['lik_binomial -> lik_binomial_predy: missing z!'...
           'Binomial likelihood needs the expected number of       '...
           'occurrences as an extra input z. See, for             '...
           'example, lik_binomial and gpla_e.                 ']);
  end

  if nargout > 1
    nt=length(Ef);
    Ey=zeros(nt,1);
    EVary = zeros(nt,1);
    VarEy = zeros(nt,1);
    for i1=1:nt
      ci = sqrt(Varf(i1));
      F  = @(x)zt(i1)./(1+exp(-x)).*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
      Ey(i1) = quadgk(F,Ef(i1)-6*ci,Ef(i1)+6*ci);

      F2  = @(x)zt(i1)./(1+exp(-x)).*(1-1./(1+exp(-x))).*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
      EVary(i1) = quadgk(F2,Ef(i1)-6*ci,Ef(i1)+6*ci);

      F3  = @(x)(zt(i1)./(1+exp(-x))).^2.*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
      VarEy(i1) = quadgk(F3,Ef(i1)-6*ci,Ef(i1)+6*ci) - Ey(i1).^2;
    end
    Vary = EVary+VarEy;
  end

  nt=length(yt);
  lpy=zeros(nt,1);
  for i1=1:nt
    ci = sqrt(Varf(i1));
    F  = @(x)exp(gammaln(zt(i1)+1)-gammaln(yt(i1)+1)-gammaln(zt(i1)-yt(i1)+1) + yt(i1).*log(1./(1+exp(-x))) + (zt(i1)-yt(i1)).*log(1-(1./(1+exp(-x))))).*norm_pdf(x,Ef(i1),sqrt(Varf(i1)));
    lpy(i1) = log(quadgk(F,Ef(i1)-6*ci,Ef(i1)+6*ci));
  end

end

function prctys = lik_binomial_predprcty(lik, Ef, Varf, zt, prcty)
%LIK_BINOMIAL_PREDPRCTY  Returns the percentiled of predictive density of y
%
%  Description
%    PRCTY = LIK_BINOMIAL_PREDPRCTY(LIK, EF, VARF YT, ZT)
%    Returns percentiles of the predictive density PY of YT, that is
%    This requires also the succes counts YT, numbers of trials ZT. This
%    subfunction is needed when using function gp_predprcty.
%
%  See also
%    GP_PREDPCTY

  if isempty(zt)
    error(['lik_binomial -> lik_binomial_predprcty: missing z!'...
           'Binomial likelihood needs the expected number of       '...
           'occurrences as an extra input z. See, for             '...
           'example, lik_binomial and gpla_e.                 ']);
  end

  opt=optimset('TolX',.5,'Display','off');
  nt=size(Ef,1);
  prctys = zeros(nt,numel(prcty));
  prcty=prcty/100;
  for i1=1:nt
    ci = sqrt(Varf(i1));
    for i2=1:numel(prcty)
      a=floor(fminbnd(@(a) (quadgk(@(f) binocdf(a,zt(i1),logitinv(f)).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)-prcty(i2)).^2,binoinv(prcty(i2),zt(i1),logitinv(Ef(i1)-1.96*ci)),binoinv(prcty(i2),zt(i1),logitinv(Ef(i1)+1.96*ci)),opt));
      if quadgk(@(f) binocdf(a,zt(i1),logitinv(f)).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)<prcty(i2)
        a=a+1;
      end
      prctys(i1,i2)=a;
    end
  end
end

function [df,minf,maxf] = init_binomial_norm(yy,myy_i,sigm2_i,N)
%INIT_LOGIT_NORM
%
%  Description
%    Return function handle to a function evaluating Binomial *
%    Gaussian which is used for evaluating (likelihood * cavity)
%    or (likelihood * posterior) Return also useful limits for
%    integration. This is private function for lik_binomial. This
%    subfunction is needed by subfunctions tiltedMoments and predy.
%
% See also
%   LIK_BINOMIAL_TILTEDMOMENTS, LIK_BINOMIAL_PREDY

% avoid repetitive evaluation of constant part
  ldconst = gammaln(N+1)-gammaln(yy+1)-gammaln(N-yy+1) - log(sigm2_i)/2 - log(2*pi)/2;
%   ldconst = log(factorial(N)/(factorial(yy)*factorial(N-yy))-log(sigm2_i)/2 -log(2*pi)/2;

 % Create function handle for the function to be integrated
  df = @binomial_norm;
 % use log to avoid underflow, and derivates for faster search
  ld = @log_binomial_norm;
  ldg = @log_binomial_norm_g;
  ldg2 = @log_binomial_norm_g2;

  % Set the limits for integration
  % Binomial likelihood is log-concave so the binomial_norm
  % function is unimodal, which makes things easier
  if yy==0 || yy==N
    % with yy==0 or yy==N 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
    % binomial likelihood and Gaussian
    mean_app = log(yy./(N-yy));
    ld0=1/(1+exp(-mean_app));
    ld1=(1-ld0)*ld0;
    ld2=ld0-3*ld0^2+2*ld0^3;
    var_app=inv(-( yy*(ld2*ld0-ld1^2)/ld0^2 + (N-yy)*(ld2*(ld0-1)-ld1^2)/(ld0-1)^2 ));

    modef = (myy_i/sigm2_i + mean_app/var_app)/(1/sigm2_i + 1/var_app);
%     sigm_app = sqrt((1/sigm2_i + 1/var_app)^-1);
  end
  % find the mode of the integrand using Newton iterations
  % few iterations is enough, since the first guess in the right direction
  niter=3;       % 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-4*modes;
  maxf=modef+4*modes;
  modeld=ld(modef);
  iter=0;
  % check that density at end points is low enough
  lddiff=12; % 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_negbin -> init_negbin_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_negbin -> init_negbin_norm: ' ...
             'integration interval maximun not found ' ...
             'even after looking hard!'])
    end
  end

  function integrand = binomial_norm(f)
  % Logit * Gaussian
    integrand = exp(ldconst + yy*log(1./(1.+exp(-f)))+(N-yy)*log(1-1./(1.+exp(-f)))...
                   - 0.5 * (f-myy_i).^2./sigm2_i);
%     integrand = exp(ldconst ...
%                     +yy*log(x)+(N-yy)*log(1-x) ...
%                     -0.5*(f-myy_i).^2./sigm2_i);
    integrand(isnan(integrand))=0;
  end

  function log_int = log_binomial_norm(f)
  % log(Binomial * Gaussian)
  % log_binomial_norm is used to avoid underflow when searching
  % integration interval

    log_int = ldconst + yy*log(1./(1.+exp(-f)))+(N-yy)*log(1-1./(1.+exp(-f)))...
                   - 0.5 * (f-myy_i).^2./sigm2_i;
%     log_int = ldconst ...
%               -log(1+exp(-yy.*f)) ...
%               -0.5*(f-myy_i).^2./sigm2_i;
  end

  function g = log_binomial_norm_g(f)
  % d/df log(Binomial * Gaussian)
  % derivative of log_logit_norm
    g = -(f-myy_i)./sigm2_i - exp(-f).*(N-yy)./((1+exp(-f)).^2.*(1-1./(1+exp(-f)))) ...
        + exp(-f).*yy./(1+exp(-f));
%     g = yy./(exp(f*yy)+1)...
%         + (myy_i - f)./sigm2_i;
  end

  function g2 = log_binomial_norm_g2(f)
  % d^2/df^2 log(Binomial * Gaussian)
  % second derivate of log_logit_norm
    g2 = - (1+exp(2.*f)+exp(f).*(2+N*sigm2_i)./((1+exp(f))^2*sigm2_i));
%     a=exp(f*yy);
%     g2 = -a*(yy./(a+1)).^2 ...
%          -1/sigm2_i;
  end

end

function p = lik_binomial_invlink(lik, f, z)
%LIK_BINOMIAL_INVLINK  Returns values of inverse link function
%
%  Description
%    P = LIK_BINOMIAL_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 gp_predprctmu.
%
%     See also
%     LIK_BINOMIAL_LL, LIK_BINOMIAL_PREDY

  p = logitinv(f);
end

function reclik = lik_binomial_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = GPCF_BINOMIAL_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
    reclik.type = 'Binomial';

    % Set the function handles
    reclik.fh.pak = @lik_binomial_pak;
    reclik.fh.unpak = @lik_binomial_unpak;
    reclik.fh.ll = @lik_binomial_ll;
    reclik.fh.llg = @lik_binomial_llg;
    reclik.fh.llg2 = @lik_binomial_llg2;
    reclik.fh.llg3 = @lik_binomial_llg3;
    reclik.fh.tiltedMoments = @lik_binomial_tiltedMoments;
    reclik.fh.invlink = @lik_binomial_invlink;
    reclik.fh.predprcty = @lik_binomial_predprcty;
    reclik.fh.predy = @lik_binomial_predy;
    reclik.fh.recappend = @likelih_binomial_recappend;
    return
  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function toggle_aspectratio(src, ~, ax, options)
    % TOGGLE_ASPECTRATIO toggles axes aspect ratio and sets checked property
    % uses the checked state to determine currentstate.
    % aspect ratio is based on current latitude (Y)
    %
    % TOGGLE_ASPECTRATIO(source, event, axesHandle, option)
    % options can be:
    %    "SetGlobal",in which case it affects ZmapGlobal
    %    "UpdateOnly",in which case it does not toggle values, but fixes daspect.
    %
    % to set multiple options, use either cell of chars, or array of strings

    if ~exist('options','var')
        options = "";
    end
    assert(all(ismember(options,["", "UpdateOnly","SetGlobal","on","off"]) ),...
        'options must be member of ["UpdateOnly", "SetGlobal","on","off"]');

    if ~any(options == "UpdateOnly")
        if any(options=="on")
            src.Checked="on";
        elseif any(options =="off")
            src.Checked ="off";
        else
            src.Checked=toggleOnOff(src.Checked);
        end
    end

    switch src.Checked
        case 'on'
            for i=1:numel(ax)
                daspect(ax(i), [1 cosd(mean(ax(i).YLim)) 10]);
            end
        case 'off'
            for i=1:numel(ax)
                daspect(ax(i),'auto');
            end
    end

    if any(options == "SetGlobal")
        ZG = ZmapGlobal.Data;
        ZG.lock_aspect = matlab.lang.OnOffSwitchState(src.Checked);
    end
end

author: CelsoReyes


text:
function plotPendulumCart(t,x,param)

Cart_Width = 0.2;
Cart_Height = 0.1;
Limits = param.axis;

Cart = x(1:2);
Bob = x(3:4);

if param.clearFigure    %THEN RUNNING ANIMATION
    clf;
    hold on;

    %Plot Cart
    x = Cart(1) - 0.5*Cart_Width;
    y = -0.5*Cart_Height;
    w = Cart_Width;
    h = Cart_Height;
    h = rectangle('Position',[x,y,w,h],'LineWidth',4,'Curvature',[0.3,0.3]);
    set(h,'EdgeColor',[0.1,0.8,0.1])

    %Plot Pendulum
    Rod_X = [Cart(1), Bob(1)];
    Rod_Y = [Cart(2), Bob(2)];
    plot(Rod_X,Rod_Y,'k-','LineWidth',4)

    %Plot Bob
    plot(Bob(1),Bob(2),'k.','MarkerSize',50)

    %Plot Rails
    plot([Limits(1) Limits(2)],-0.5*Cart_Height*[1,1],'k-','LineWidth',2)

    %Title
    title(['Simulation Time: ' num2str(t) ' s'])

    %These commands keep the window from automatically rescaling in funny ways.
    axis(Limits);
    axis('equal');
    axis manual;
    axis off;

else    %THEN RUNNING STOP ACTION

    hold on;

    %Plot Cart
    x = Cart(1) - 0.5*Cart_Width;
    y = -0.5*Cart_Height;
    w = Cart_Width;
    h = Cart_Height;
    rectangle('Position',[x,y,w,h],'LineWidth',2);

    %Plot Pendulum
    Rod_X = [Cart(1), Bob(1)];
    Rod_Y = [Cart(2), Bob(2)];
    plot(Rod_X,Rod_Y,'k-','LineWidth',2)

    %Plot Bob and hinge
    plot(Bob(1),Bob(2),'k.','MarkerSize',30)
    plot(Cart(1),Cart(2),'k.','MarkerSize',30)

    %Plot Rails
    plot([Limits(1) Limits(2)],-0.5*Cart_Height*[1,1],'k-','LineWidth',2)

    %Plot trace of the Bob's path:
    plot(param.Bob(1,:),param.Bob(2,:),'k:','LineWidth',1)

    %These commands keep the window from automatically rescaling in funny ways.
    axis(Limits);
    axis('equal');
    axis manual;
    axis off;

end
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

% You need to return the following variables correctly
g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the sigmoid of each value of z (z can be a matrix,
%               vector or scalar).

g = 1 ./ (1 + e.^(-z));

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function node_xyz = sphere_grid_q9_node_xyz ( nelemx, nelemy )

%*****************************************************************************80
%
%% SPHERE_GRID_Q9_NODE_XYZ produces node coordinates for a Q9 sphere grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.
%
%    Output, real NODE_XYZ(3,4*NELEMX*NELEMY-2*NELEMX+2),
%    the node coordinates.
%
  node = 0;

  node = node + 1;
  node_xyz(1,node) =  0.0;
  node_xyz(2,node) =  0.0;
  node_xyz(3,node) = -1.0;

  for j = 2 * nelemy : -1 : 2

    phi = ( j - 1 ) * pi / ( 2 * nelemy );

    for i = 1 : 2 * nelemx

      theta = ( i - 1 ) * 2.0 * pi / ( 2 * nelemx );

      node = node + 1;
      node_xyz(1,node) = cos ( theta ) * sin ( phi );
      node_xyz(2,node) = sin ( theta ) * sin ( phi );
      node_xyz(3,node) =                 cos ( phi );

    end
  end

  node = node + 1;
  node_xyz(1,node) =  0.0;
  node_xyz(2,node) =  0.0;
  node_xyz(3,node) =  1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 % Copyright 2001, Brown University, Providence, Rhode Island.
 %
 % All Rights Reserved
 %
 % Permission to use this software for noncommercial research and
 % educational purposes is hereby granted without fee.
 % Redistribution, sale, or incorporation of this software into a
 % commercial product is prohibited.
 %
 % BROWN UNIVERSITY DISCLAIMS ANY AND ALL WARRANTIES WITH REGARD TO
 % THIS SOFTWARE,INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
 % AND FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL BROWN
 % UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 % DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE,
 % DATA OR PROFITS.

% -----------------------------------------------------------------
%   jacobi() - jacobi polynomials
%
%   Get a vector 'poly' of values of the n_th order Jacobi polynomial
%   P^(alpha,beta)_n(z) alpha > -1, beta > -1 at the np points in z
%   -----------------------------------------------------------------

function jf = JACOBI1D(z, n, alpha, beta)

  dims = size(z);
  jf   = zeros(dims);
  one = 1.0;
  two = 2.0;

  if(n == 0)
   jf(:) = one;
  elseif (n == 1)
   jf(:) = 0.5*(alpha - beta + (alpha + beta + two)*z);
  else
   two = 2.0;
   apb = alpha + beta;

   poly   = zeros(dims);
   polyn2 = ones(dims);
   polyn1 = 0.5*(alpha - beta + (alpha + beta + two)*z);

   for k = 2:n
    a1 =  two*k*(k + apb)*(two*k + apb - two);
    a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
    a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
    a4 =  two*(k + alpha - one)*(k + beta - one)*(two*k + apb);

    a2 = a2/a1;
    a3 = a3/a1;
    a4 = a4/a1;

    poly   = (a2 + a3*z).*polyn1 - a4*polyn2;
    polyn2 = polyn1;
    polyn1 = poly;
   end
   jf = poly;
  end

repo: Aero-matlab


text:
function data_num = p02_data_num ( )

%*****************************************************************************80
%
%% P02_DATA_NUM returns the dimension of the data vector for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the dimension of the data vector.
%
  data_num = 11;

  return
end


text:
function [R, movement, verbose] = tapas_physio_create_movement_regressors(...
    movement, verbose)
% Reads realignment parameters, creates derivative/squared & outlier regressor
%
% [R, movement, verbose] = ...
%       tapas_physio_create_movement_regressors(movement, verbose)
%
% The 6 realignment parameters can be augmented by their derivatives (in
% total 12 parameters) + the squares of parameters and derivatives (in
% total 24 parameters), leading to a Volterra expansion as in
%
%               Friston KJ, Williams S, Howard R, Frackowiak
%               RS, Turner R. Movement-related effects in fMRI
%               time-series. Magn Reson Med. 1996;35:346?355.)
%
% IN
%   movement    physio.model.movement
%   verbose     physio.verbose
% OUT
%   R                   [nScans, (6|12|24)+nOutliers] regressor matrix
%                       from movement modeling
%   movement            structure (as defined in tapas_physio_new) with the
%                       following fields
%   rp                  [nScans, 6] realignment parameters
%   dRp                 [nScans, 6] temporal difference of realignment
%                       parameters, prepending one line of zeros (for 1st
%                       scan)
%   quality_measures.   structure holding the suggested quality control
%                       measures for subject motion of Power et al., 2014,
%                       Fig. 2; all rotational parameters are transformed
%                       into translations by multiplying with rHead (arc
%                       length)
%       FD              framewise displacement (FD)
%       absTransDisplacement
%                       sum of absolute values of translational (x,y,z)
%                       realignment estimates, reflecting absolute
%                       displacement of the head
%       absRotDisplacement
%                       sum of absolute values of rotational (rotation
%                       around x,y,z-axis = pitch/roll/yaw)
%                       realignment estimates, reflecting absolute
%                       displacement of the head on its surface (arc
%                       length)
%       rmsdTrans       root mean squared framewise translational displacement (over x,y,z)
%                       equals to euclidean distance of center of
%                       mass of head between consecutive volumes
%       rmsdRot         root mean squared framewise rotational displacement
%                       (rotation around x,y,z axis)
%                       equals to euclidean distance of between angles
%                       between consecutive volumes
%       meanFD          mean (over volumes) of framewise displacement
%                       summary measure for subject (or session)
%       rmsMovement     root mean square (over scans) of detrended
%                       realignment estimates
%                       summary measure for subject (or session)
%
%   censoring.          structure with censoring information (detected outliers,
%                       computed quality values for censoring method maxval/fd:
%       nOutliers       number of detected outliers
%       R_outlier       stick/spike regressors indicating location of
%                       detected outliers
%       iOutlierTrans   volume indices of translation-related outliers
%       iOutlierRot     volume indices of rotation-related outliers
%       iOutlierArray   volume indices of all detected outliers
%
%
% EXAMPLE
%   tapas_physio_create_movement_regressors
%
%   See also tapas_physio_new tapas_physio_main_create_regressors

% Author: Lars Kasper
% Created: 2015-07-10
% Copyright (C) 2015 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% License (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/>.

rHead = 50; % head radius in mm for FD computation (Power et al., 2012)

[fp, fn, ext] = fileparts(movement.file_realignment_parameters);

if ~exist(movement.file_realignment_parameters, 'file')
    verbose = tapas_physio_log('No input multiple regressors found', verbose, 1);
    R = [];
else
    tmp = load(movement.file_realignment_parameters);
    if strcmp('.txt', ext) % text-file
        rp = tmp;
    else % mat file
        rp = tmp.R;
    end

    %% Motion Quality control

     [quality_measures, dRp] = ...
         tapas_physio_get_movement_quality_measures(rp, rHead);

    %% Motion 6/12/24 model
    R = rp;

    % Include derivatives
    if movement.order > 6
        R = [rp, dRp];
    end

    % Include squared regressors/derivatives
    if movement.order > 12
        R = [R, R.^2];
    end

    % Sanity check if order misspecified
    movement.order = size(R, 2);

    %% Motion Censoring
    % Include outlier movements exceeding thresholds as stick regressors
    % Euclidean distance used!
    iOutlierArray = [];
    iOutlierTrans = [];
    iOutlierRot = [];

    switch lower(movement.censoring_method)
        case 'none' % done
        case 'fd' % framewise displacement
            iOutlierArray = find(quality_measures.FD > ...
                movement.censoring_threshold);
        case 'maxval'   % tresholds for max abrupt translation/rotation,
                        % even per axis
            switch numel(movement.censoring_threshold)
                case {1,2}
                    outlier_translation_mm = movement.censoring_threshold(1);
                    outlier_rotation_deg = movement.censoring_threshold(end);
                    iOutlierTrans = find(quality_measures.rmsdTrans > outlier_translation_mm);
                    iOutlierRot = find(quality_measures.rmsdRot > outlier_rotation_deg/180*pi);
                case 6
                    ct = movement.censoring_threshold;
                    ct(4:6) = ct/180*pi;
                    iOutlierTrans = dRp(:,1) > ct(1) || dRp(:,2) > ct(2) || dRp(:,3) > ct(3);
                    iOutlierRot = dRp(:,4) > ct(4) || dRp(:,5) > ct(5) || dRp(:,6) > ct(6);
                otherwise
                    error('censoring threshold has to be 1,2 or 6 element vector. See tapas_physio_new');
            end

            iOutlierArray = unique([iOutlierTrans; iOutlierRot]);
        case 'dvars' %DVARS, as in Power et al, 2012
            % TODO
    end

    %% Construct censoring regressor matrix
    nOutliers = numel(iOutlierArray);

    nScans = size(rp,1);
    R_outlier = zeros(nScans, nOutliers);
    for iOutlier = 1:nOutliers
        R_outlier(iOutlierArray(iOutlier), iOutlier) = 1;
    end

    censoring = struct('nOutliers', nOutliers, 'R_outlier', R_outlier, ...
        'iOutlierTrans', iOutlierTrans, 'iOutlierRot', iOutlierRot, ...
        'iOutlierArray', iOutlierArray);

    if verbose.level >= 2
        switch lower(movement.censoring_method)
            case 'fd'
                verbose.fig_handles(end+1) = tapas_physio_plot_movement_outliers_fd(rp, ...
                    quality_measures, censoring, movement.censoring_threshold);
            case 'maxval'
                verbose.fig_handles(end+1) = tapas_physio_plot_movement_outliers_maxval(rp, ...
                    quality_measures, censoring, movement.censoring_threshold);
        end
    end

    %% Gather return values
    movement.censoring = censoring;
    movement.rp = rp;
    movement.dRp = dRp;
    movement.quality_measures = quality_measures;

    R = [R, R_outlier];

end

author: translationalneuromodeling

repo: tapas


text:
function v = rowNorm(A)
%ROWNORM  norm of each row
%   V = ROWNORM(A) returns the Euclidean norm of each row of A. When A is
%   an M*N matrix, the output, V, will be a M*1 vector.
%
%   Babak taati May 4, 2005
%   revised May 19, 2009

if nargin ~= 1
    error('Requires one input arguments.')
end

v = sqrt(sum(A.*A, 2));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function Pout = scale_imgs_by_csf(hP)
% Takes a string matrix of image file names
% finds the mean and std of the CSF space
% specified in a mask (hard-coded)
% and standardizes images by these values
%
% :Usage:
% ::
%
%     Pout = scale_imgs_by_csf(hP)
%
% Writes SC* images (SCaled)
%
% assumes images are spatially normalized.
% uses a canonical CSF mask!
%
% ..
%    tor wager
% ..

mP = which('canonical_ventricles.img');
[tmp,mP] = reslice_imgs(deblank(hP(1,:)),mP,0);

M = spm_general_hist(hP,mP,'canonical_ventricles');
disp(' ')
disp(['Means are: ' num2str(M(:,1)')])
disp(['Var    is: ' num2str(M(:,2)')])

for i = 1:size(hP,1)

	V = spm_vol(deblank(hP(i,:)));
	v = spm_read_vols(V);

	v = v - M(i,1);
	v = v ./ M(i,2);

	[d,f,e] = fileparts(V.fname);
	V.fname = [d filesep 'SC' f e];

	if i == 1, Pout = V.fname;, else, Pout = str2mat(Pout,V.fname);, end
	spm_write_vol(V,v);
	disp(['Written ' V.fname])

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [w,x]=fclencurt(N1,a,b)
% Fast Clenshaw Curtis Quadrature

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% fclencurt.m - Fast Clenshaw Curtis Quadrature
%
% Compute the N nodes and weights for Clenshaw-Curtis
% Quadrature on the interval [a,b]. Unlike Gauss
% quadratures, Clenshaw-Curtis is only exact for
% polynomials up to order N, however, using the FFT
% algorithm, the weights and nodes are computed in linear
% time. This script will calculate for N=2^20+1 (1048577
% points) in about 5 seconds on a normal laptop computer.
%
% Written by: Greg von Winckel - 02/12/2005
% Contact: gregvw(at)chtm(dot)unm(dot)edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N=N1-1; bma=b-a;
c=zeros(N1,2);
c(1:2:N1,1)=(2./[1 1-(2:2:N).^2 ])'; c(2,2)=1;
f=real(ifft([c(1:N1,:);c(N:-1:2,:)]));
w=bma*([f(1,1); 2*f(2:N,1); f(N1,1)])/2;
x=0.5*((b+a)+N*bma*f(1:N1,2));

author: mtex-toolbox

repo: mtex

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


text:
%% Machine Learning Online Class
%  Exercise 1: Linear regression with multiple variables
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear regression exercise.
%
%  You will need to complete the following functions in this
%  exericse:
%
%     warmUpExercise.m
%     plotData.m
%     gradientDescent.m
%     computeCost.m
%     gradientDescentMulti.m
%     computeCostMulti.m
%     featureNormalize.m
%     normalEqn.m
%
%  For this part of the exercise, you will need to change some
%  parts of the code below for various experiments (e.g., changing
%  learning rates).
%

%% Initialization

%% ================ Part 1: Feature Normalization ================

%% Clear and Close Figures
clear ; close all; clc

fprintf('Loading data ...\n');

%% Load Data
data = load('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

% Print out some data points
fprintf('First 10 examples from the dataset: \n');
fprintf(' x = [%.0f %.0f], y = %.0f \n', [X(1:10,:) y(1:10,:)]');

fprintf('Program paused. Press enter to continue.\n');
pause;

% Scale features and set them to zero mean
fprintf('Normalizing Features ...\n');

[X mu sigma] = featureNormalize(X);

% Add intercept term to X
X = [ones(m, 1) X];

%% ================ Part 2: Gradient Descent ================

% ====================== YOUR CODE HERE ======================
% Instructions: We have provided you with the following starter
%               code that runs gradient descent with a particular
%               learning rate (alpha).
%
%               Your task is to first make sure that your functions -
%               computeCost and gradientDescent already work with
%               this starter code and support multiple variables.
%
%               After that, try running gradient descent with
%               different values of alpha and see which one gives
%               you the best result.
%
%               Finally, you should complete the code at the end
%               to predict the price of a 1650 sq-ft, 3 br house.
%
% Hint: By using the 'hold on' command, you can plot multiple
%       graphs on the same figure.
%
% Hint: At prediction, make sure you do the same feature normalization.
%

fprintf('Running gradient descent ...\n');

% Choose some alpha value
alpha = 0.01;
num_iters = 400;

% Init Theta and Run Gradient Descent
theta = zeros(3, 1);
[theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters);

% Plot the convergence graph
figure;
plot(1:numel(J_history), J_history, '-b', 'LineWidth', 2);
xlabel('Number of iterations');
ylabel('Cost J');

% Display gradient descent's result
fprintf('Theta computed from gradient descent: \n');
fprintf(' %f \n', theta);
fprintf('\n');

% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================
% Recall that the first column of X is all-ones. Thus, it does
% not need to be normalized.

x = [1 1650 3]';
price = theta' * x;

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

fprintf(['Predicted price of a 1650 sq-ft, 3 br house ' ...
         '(using gradient descent):\n $%f\n'], price);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 3: Normal Equations ================

fprintf('Solving with normal equations...\n');

% ====================== YOUR CODE HERE ======================
% Instructions: The following code computes the closed form
%               solution for linear regression using the normal
%               equations. You should complete the code in
%               normalEqn.m
%
%               After doing so, you should complete this code
%               to predict the price of a 1650 sq-ft, 3 br house.
%

%% Load Data
data = csvread('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

% Add intercept term to X
X = [ones(m, 1) X];

% Calculate the parameters from the normal equation
theta = normalEqn(X, y);

% Display normal equation's result
fprintf('Theta computed from the normal equations: \n');
fprintf(' %f \n', theta);
fprintf('\n');

% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================

x = [1 1650 3]';
price = theta' * x;

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

fprintf(['Predicted price of a 1650 sq-ft, 3 br house ' ...
         '(using normal equations):\n $%f\n'], price);

author: zsiciarz

repo: ml-coursera


text:
%Polypeptide_Test
% Polypeptide test class.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 1/4/2011
classdef Polypeptide_Test < edu.stanford.covert.cell.sim.CellStateTestCase
    methods
        function this = Polypeptide_Test(name)
            this = this@edu.stanford.covert.cell.sim.CellStateTestCase(name);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [metabolite_structure,IDsAdded] = getCas2Echa(metabolite_structure)
%
% The input file was downloaded from https://echa.europa.eu/documents/10162/13629/ec_inventory_en.xlsx
%
%
% INPUT
% metabolite_structure  metabolite structure
%
% OUTPUT
% metabolite_structure  Updated metabolite structure
%
%
% Ines Thiele, 2020-2021
% first column contains echa_id,4th col is cas registry
[NUM,TXT,RAW]=xlsread('ec_inventory_en.xlsx');

annotationSource = 'Echa_id mapping from CasRegistry';
annotationType = 'automatic';

Mets = fieldnames(metabolite_structure);
for i = 1 : length(Mets)
    Cas{i,1} = Mets{i};
    Cas{i,2} = num2str(metabolite_structure.(Mets{i}).casRegistry);
end

for i = 4 : size(RAW,1)
    casRAW{i-3,1} = RAW{i,4};
    casRAWList{i-3,1} = RAW{i,1};
end
a = 1;
IDsAdded = '';

for i = 1 : size(Cas,1)
    % no keggId exists
    if ~isempty(metabolite_structure.(Cas{i,1}).casRegistry) && isempty(find(isnan(metabolite_structure.(Cas{i,1}).casRegistry),1))
        if isempty(metabolite_structure.(Cas{i,1}).echa_id) || ~isempty(find(isnan(metabolite_structure.(Cas{i,1}).echa_id),1))
            match = strmatch(Cas{i,2},casRAW,'exact');
            if ~isempty(match)
                metabolite_structure.(Cas{i,1}).echa_id = casRAWList{match};
                metabolite_structure.(Cas{i,1}).echa_id_source =  [annotationSource,':',annotationType,':',datestr(now)];
                IDsAdded{a,1} = Cas{i,1};
                IDsAdded{a,2} = 'echa_id';
                IDsAdded{a,3} =  metabolite_structure.(Cas{i}).echa_id;
                a = a + 1;
            end
        end
    end
end
[metabolite_structure] = cleanUpMetabolite_structure(metabolite_structure);

author: opencobra

repo: cobratoolbox


text:
function inside = p15_inside ( m, n, point )

%*****************************************************************************80
%
%% P15_INSIDE reports if a point is inside the region in problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, logical INSIDE(N), is TRUE if the point is in the region.
%
  inside(1:n) =                                          ...
       ( -8.0 <= point(1,1:n) & point(1,1:n) <= 2.0 &    ...
         -1.0 <= point(2,1:n) & point(2,1:n) <= 0.0    ) ...
    |                                                    ...
       ( -2.0 <= point(1,1:n) & point(1,1:n) <= 8.0 &    ...
          0.0 <= point(2,1:n) & point(2,1:n) <= 1.0 );

  return
end


text:
function [bar] = dynpcm22bar(dynpcm2)
% Convert pressure from dynes per square centimeter to bar
% Chad Greene 2012
bar = dynpcm2*0.00000100000;

author: Sable

repo: mcbench-benchmarks


text:
function [AGC_atten_dB] = extract_AGC( file_in, Nlines )
%
% -------------------------------------------------------------------------
% This program reads the RADARSAT aux data and calculates the AGC settings
% -------------------------------------------------------------------------
%  file_in      : input file name including it's path
%  Nlines       : number of azimuth lines
%  AGC_atten_dB : Receiver attenuation in dB for each range line
% -------------------------------------------------------------------------
% Created  : May 05, 2004 by Kaan Ersahin & Millie Sikdar
% Modified : Nov 22, 2004 by Kaan Ersahin
% -------------------------------------------------------------------------

% The parameters encoded in the auxilary data bits are defined in RSI-D6
% also known as RSCSA-IC0009 (X-band ICD)
% -------------------------------------------------------------------------
% PARAMETER NAME                          LOCATION         LENGTH   ID
% -------------------------------------------------------------------------
% aux_sync_marker         = aux_bits (:,   1:  32);     % 32 bit -  1
% image_ref_id            = aux_bits (:,  33:  64);     % 32 bit -  2
% payload_status          = aux_bits (:,  65:  80);     % 16 bit -  3
% replica_AGC             = aux_bits (:,  81:  86);     %  6 bit -  4
% CALN_atten_LPT_pow_set  = aux_bits (:,  89:  96);     %  8 bit -  6
% pulse_waveform_number   = aux_bits (:,  97: 100);     %  4 bit -  7
% temperature             = aux_bits (:, 113: 144);     % 32 bit -  9
% beam_sequence           = aux_bits (:, 145: 160);     % 16 bit - 10
% ephemeris               = aux_bits (:, 161: 176);     % 16 bit - 11
% number_of_beams         = aux_bits (:, 177: 178);     %  2 bit - 12
% ADC_rate                = aux_bits (:, 179: 180);     %  2 bit - 13
% pulse_count_1           = aux_bits (:, 185: 192);     %  8 bit - 15
% pulse_count_2           = aux_bits (:, 193: 200);     %  8 bit - 16
% PRF_beam                = aux_bits (:, 201: 213);     % 13 bit - 17
% beam_select             = aux_bits (:, 214: 215);     %  2 bit - 18
% Rx_window_start_time    = aux_bits (:, 217: 228);     % 12 bit - 20
% Rx_window_duration      = aux_bits (:, 233: 244);     % 12 bit - 22
% altitude                = aux_bits (:, 249: 344);     % 96 bit - 24
% time                    = aux_bits (:, 345: 392);     % 48 bit - 25
% SC_T02_defaults         = aux_bits (:, 393: 393);     %  1 bit - 26
% first_replica           = aux_bits (:, 394: 394);     %  1 bit - 27
% Rx_AGC_setting          = aux_bits (:, 395: 400);     %  6 bit - 28
% -------------------------------------------------------------------------
%                                    Total  => 50 bytes (400 bits)
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% read the binary aux data file
% -------------------------------------------------------------------------
file_in
fid_aux  = fopen(file_in,'r');
aux_bytes = fread(fid_aux,[50, Nlines],'uint8');
aux_bytes = aux_bytes';
fclose(fid_aux);

aux_bits = char(zeros(Nlines,400));
% -------------------------------------------------------------------------
% convert 50 bytes to 400 bits in general (it is only done for 50th byte)
% -------------------------------------------------------------------------
for byte = 50:50 %1:50 %% in general
   aux_bits(:,8*(byte-1)+1:8*byte) = num2str(dec2bin(aux_bytes(:,byte),8));
end
% -------------------------------------------------------------------------
%  Convert last 6 bits from binary to decimal
% -------------------------------------------------------------------------
d_Rx_AGC = bin2dec( aux_bits (:, 395: 400) );

% -------------------------------------------------------------------------
% For values greater than 31, the binary representation defined in the
% documentation is different and there is the need to substract 24 from
% decimal values to get the AGC setting in dB
% -------------------------------------------------------------------------
AGC_atten_dB = d_Rx_AGC - ( 24 * (d_Rx_AGC > 31));

author: denkywu

repo: SAR-Synthetic-Aperture-Radar


text:
function S = pad_strn (strg,padchar,n,loc)
%	S = PAD_STRN (strg,padchar,n,loc)
%This function pad the string 'strg' with 'padchar' characters
%  such that 'S' is 'n' characters long.
% if loc = 1 the 'padchar' are added to the left of 'strg'
% if loc = -1 the 'padchar' are added to the right of 'strg'

if (~ischar(strg) | ~ischar(padchar)),
	err = ErrEval(mfilename, 'Err_First two arguments must be characters');
	S = '';
	return;
end

S = strg;
while (length(S) < n)

	if (loc == 1),
		S = sprintf ('%s%s',padchar,S);
	elseif (loc == -1),
		S = sprintf ('%s%s',S,padchar);
	end

end;

return;


text:
% DEMSVARGPLVMGENERIC A generic demo for the MRD method.
%
% COPYRIGHT: Andreas C. Damianou, 2012
%
% SEEALO: MRDEmbed.m
%
% VARGPLVM

%------
% This is a generic demo for svargplvm. It only requires a cell-array
% Yall to be a priori created. Each cell should contain one dataset.
% This demo will use default options and configurations for the model.
% These can be changed by investigating the fields of globalOpt structure,
% which are set during 'svargplvm_init'.

%{
%eg:
clear
dataSetNames = 'toy'; toyDataCreate = 'fols'; svargplvmPrepareData
baseKern = {{'linard2','white'},{'linard2','white'}};
latentDimPerModel = 4;
init_X = 'separately';
initVardistIters = 200;
itNo = 200;
demSvargplvmGeneric
%}

% % Fix seeds
% randn('seed', 1e5);
% rand('seed', 1e5);

if ~exist('Yall', 'var')
    error('# This demo requires that a cell-array Yall is a priori created. This array should contain one dataset per cell.')
end
if ~exist('trainModel', 'var'), trainModel = true; end

M = length(Yall);

% This script initialises the options structure 'globalOpt'.
svargplvm_init;

% Don't use more inducing points than data
globalOpt.indPoints = min(globalOpt.indPoints, size(Yall{1},1));

%-- Load datasets
for i=1:M
    Y = Yall{i};
    dims{i} = size(Y,2);
    if i == 1
        Nall = size(Y,1);
    else
        if size(Y,1) ~= Nall
            error('The number of observations in each dataset must be the same!');
        end
    end
    indTr = globalOpt.indTr;
    if indTr == -1
        indTr = 1:Nall;
    end
    if ~exist('Yts')
        indTs = setdiff(1:size(Y,1), indTr);
        Yts{i} = Y(indTs,:);
    end
    Ytr{i} = Y(indTr,:);
end
clear('Y','Yall')

%--- Create model
if ~exist('options','var')
    options = svargplvmOptions(Ytr, globalOpt);
else
    warning('options field exists already!');
end

if ~isempty(globalOpt.dynamicsConstrainType)
    if exist('optionsDyn', 'var')
        warning('optionsDyn field exists already!');
    else
        if exist('timeStampsTraining', 'var')
            t = timeStampsTraining;
        else
            t = linspace(0, 2*pi, Nall+1)'; t = t(1:end-1, 1);
        end
        timeStampsTraining = t(indTr,1); %timeStampsTest = t(indTs,1);
        optionsDyn.type = 'vargpTime';
        optionsDyn.inverseWidth=30;
        optionsDyn.initX = globalOpt.initX;
        optionsDyn.constrainType = globalOpt.dynamicsConstrainType;
        if exist('timeStampsTraining', 'var')
            optionsDyn.t = timeStampsTraining;
        end
        if exist('labelsTrain', 'var') && ~isempty(labelsTrain)
            optionsDyn.labels = labelsTrain;
        end
    end
else
    optionsDyn= [];
end

model = svargplvmModelCreate(Ytr, globalOpt, options, optionsDyn);
if exist('diaryFile'), model.diaryFile = diaryFile; end

if ~isfield(globalOpt, 'saveName') || isempty(globalOpt.saveName)
    model.saveName = vargplvmWriteResult([], model.type, globalOpt.dataSetName, globalOpt.experimentNo, [], globalOpt.saveModelDir);
end
model.globalOpt = globalOpt;
model.options = options;

%-- Define what level of parallelism to use (w.r.t submodels or/and w.r.t
% datapoints).
%{
fprintf('# Parallel computations w.r.t the submodels!\n');
model.parallel = 1;
model = svargplvmPropagateField(model,'parallel', 1);
%
fprintf('# Parallel computations w.r.t the datapoints!\n');
model.vardist.parallel = 1;
for i=1:model.numModels
    model.comp{i}.vardist.parallel = 1;
end
%}

% Force kernel computations
params = svargplvmExtractParam(model);
model = svargplvmExpandParam(model, params);

if ~isempty(globalOpt.dynamicsConstrainType)
    fprintf('# Median of vardist. covars: %d \n',median(median(model.vardist.covars)));
    fprintf('# Min of vardist. covars: %d \n',min(min(model.vardist.covars)));
    fprintf('# Max of vardist. covars: %d \n',max(max(model.vardist.covars)));
end

%%

if trainModel
    model = svargplvmOptimiseModel(model);
    svargplvmShowScales(model)
end

return

%% PREDICTIONS

% obsMod = 1 means that the 1st modality is considered to be observed and
% the other modality is considered to be the unobserved one, during test
% time.
obsMod = 1; % one of the involved sub-models (the one for which we have the data)
infMod = setdiff(1:2, obsMod);

[sharedDims, privateDims] = svargplvmFindSharedDims(model);

% We can either use a totally new test set (a matrix Yts must be created
% where e.g. Yts{1} and Yts{2} are the two test modalities) or using the
% training set. In the later case, we solve the correspondance problem, ie
% for a given y we find the K most similar z's in the other modality (in
% the code, K = numberOfNN). See the Yale faces example.
fprintf('\n# PREDICTIONS: \n\n');
if ~exist('testOnTraining')
    testOnTraining=0;
end

numberOfNN = 5;

%--------------------%
svargplvmPredictions %---- Script returning: ZpredMuAll and mini(the indices for NN)
%--------------------%

author: SheffieldML

repo: vargplvm


text:
% Run file to produce the Poisson experiment from the paper
%   S. Dolgov, D. Savostyanov,
%   "Alternating minimal energy methods for linear systems in higher
%   dimensions"

try
    maxNumCompThreads(1);
catch
    % Just skip. Sometimes if you specify -singleCompThread in the command
    % line, MATLAB will fail at maxNumCompThreads with scary, so tell him
    % it's okay.
end;

d = 16;
L = 6;

%!!! Last run: on dx360-11, E5-2670 0 @ 2.60GHz  !!!!!!!

tol = 1e-6;
nswp = 5;
kickrank = 4;

tol_gmres = 1e-3; % I would not set anything lower

% Matrix, RHS
A = tt_qlaplace_dd(L*ones(1,d))*(2^L+1)^2;
A = tt_reshape(A, 2^L*ones(d,2));
f = tt_ones(2^L, d);

% Reference solution
u_ex = amen_solve2(A, f, 1e-10, 'nswp', 30, 'max_full_size', 1500);

% Compare the solvers
% testdata format: {d,swp}
[u_svd,sd_svd] = amen_solve2(A, f, tol, 'x0', f, 'verb', 4, 'kicktype', 'svd', 'kickrank', kickrank, 'nswp', nswp, 'max_full_size', 50, 'trunc_norm', 'fro', 'ismex', true);
[u_als,sd_als] = amen_solve2(A, f, tol, 'x0', f, 'verb', 4, 'kicktype', 'als', 'kickrank', kickrank, 'nswp', nswp, 'max_full_size', 50, 'trunc_norm', 'fro', 'ismex', true);
[u_alstpz,sd_alstpz] = alstpz_solve(A,f,tol, 'x0',f, 'verb', 4, 'max_full_size', 50, 'kickrank', kickrank, 'nswp', nswp, 'kicktype', 'svd', 'symm', false, 'ismex', true);
[u_dmrg,sd_dmrg] = dmrg_solve3(A, f, tol, 'x0', f, 'kickrank', 0, 'min_dpow', 0.5, 'step_dpow', 0, 'step_drank', 0, 'max_full_size', 50, 'nswp', nswp, 'verb', 4, 'ismex', true, 'trunc_norm', 'fro', 'dirfilter', 1);
[U_gmres,td_gmres] = tt_gmres(core(A), core(f), tol_gmres, 4, 15, tol_gmres, tol_gmres, [], [], [], [], 3);

% Measure the A-errors
err0 = sqrt(dot(u_ex, A*u_ex));
errs_svd = zeros(d,nswp);
errs_als = zeros(d,nswp);
errs_alstpz = zeros(d,nswp);
errs_dmrg = zeros(d-1,nswp);
for i=1:nswp
    for j=1:d
        errs_svd(j,i) = sqrt(dot(sd_svd{2}{j,i}-u_ex, A*(sd_svd{2}{j,i}-u_ex)))/err0;
        errs_als(j,i) = sqrt(dot(sd_als{2}{j,i}-u_ex, A*(sd_als{2}{j,i}-u_ex)))/err0;
        errs_alstpz(j,i) = sqrt(dot(sd_alstpz{2}{j,i}-u_ex, A*(sd_alstpz{2}{j,i}-u_ex)))/err0;
        if (j<d)
            errs_dmrg(j,i) = sqrt(dot(sd_dmrg{2}{j,2*i-1}-u_ex, A*(sd_dmrg{2}{j,2*i-1}-u_ex)))/err0;
        end;
    end;
end;

% Prepare the data in the TikZ-readable form
%   1       2       3       4       5       6       7
% iter, t_alstpz, e_alstpz, t_svd, e_svd, t_als, e_als
dat_amens = [(1/d:1/d:nswp)', sd_alstpz{1}(:), errs_alstpz(:), sd_svd{1}(:), errs_svd(:), sd_als{1}(:), errs_als(:)];
dat_dmrg = [(1/(d-1):1/(d-1):nswp)', reshape(sd_dmrg{1}(1:d-1,1:2:nswp*2), [],1), errs_dmrg(:)];

% Process the GMRES output. It is different...
if (~isempty(whos('td_gmres')))
    iter_gmres = min([find(td_gmres{1}==0, 1), numel(td_gmres{1})+1])-1;
    err_gmres = zeros(iter_gmres, 1);
    for i=1:iter_gmres
        err_gmres(i) = sqrt(dot(tt_tensor(td_gmres{2}{i})-u_ex, A*(tt_tensor(td_gmres{2}{i})-u_ex)))/err0;
    end;
    dat_gmres = [(1:iter_gmres)', td_gmres{1}(1:iter_gmres)', err_gmres];
end;

% Draw 'em for humans
% iter
figure(1);
semilogy(dat_amens(:,1), dat_amens(:,[3,5,7]), dat_dmrg(:,1), dat_dmrg(:,3), dat_gmres(1:min(iter_gmres,nswp),1), dat_gmres(1:min(iter_gmres,nswp),3));
legend('alstpz', 'amen-svd', 'amen-als', 'dmrg', 'gmres');
% time
figure(2);
loglog(dat_amens(:,2), dat_amens(:,3), dat_amens(:,4), dat_amens(:,5), dat_amens(:,6), dat_amens(:,7), dat_dmrg(:,2), dat_dmrg(:,3), dat_gmres(:,2), dat_gmres(:,3));

% % Uncomment this if you want to draw the data elsewhere
% save('conv_lp16_amens.dat', '-ascii', 'dat_amens');
% save('conv_lp16_dmrg.dat', '-ascii', 'dat_dmrg');
% save('conv_lp16_gmres.dat', '-ascii', 'dat_gmres');

author: oseledets

repo: TT-Toolbox


text:
% interp1 举例
clear; clc;

xi=0:pi/5:2*pi; % 将插值区间分成若干等距小区间
yi=sin(xi); % 插值节点处的函数值

xh=0:pi/50:2*pi; % 需要插值的点

% 分段线性插值
subplot(2,2,1)
yh=interp1(xi,yi,xh); % 根据插值函数求出的近似值
plot(xi,yi,'+r', xh,yh,'o-','LineWidth',1.5);
title('分段线性插值','FontSize',18)

% 分段零次插值
subplot(2,2,2)
yh=interp1(xi,yi,xh,'nearst'); % 用邻近的值近似
plot(xi,yi,'+r', xh,yh,'o-','LineWidth',1.5);
title('分段零次插值','FontSize',18)

% 分段三次 Hermite
subplot(2,2,3)
yh=interp1(xi,yi,xh,'pchip');
plot(xi,yi,'+r', xh,yh,'o-','LineWidth',1.5);
title('分段三次Hermite插值','FontSize',18)

% 分段三次样条插值
subplot(2,2,4)
yh=interp1(xi,yi,xh,'spline');
plot(xi,yi,'+r', xh,yh,'o-','LineWidth',1.5);
title('分段三次样条插值','FontSize',18)

repo: NumericalAnalysis


text:
function tests = test_ft_connectivity_mutualinformation

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_connectivity_mutualinformation

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  func = localfunctions;
  for i=1:numel(func)
    fprintf('evaluating %s\n', func2str(func{i}));
    feval(func{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_chan_time(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan  = 6;
ntime  = 100000;
dat    = randn(nchan, ntime);

%   histmethod = The way that histograms are generated from the data. Possible values
%                are 'eqpop' (default), 'eqspace', 'ceqspace', 'gseqspace'.
%                See the help of the 'binr' function in the ibtb toolbox for more information.
%   numbin     = scalar value. The number of bins used to create the histograms needed for
%                the entropy computations
%   opts       = structure that is passed on to the 'information' function in the ibtb
%                toolbox. See the help of that function for more information.
%   refindx    = scalar value or 'all'. The channel that is used as 'reference channel'.

result = {};
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'eqpop'    , 'numbin', 5);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'eqspace'  , 'numbin', 5);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'ceqspace' , 'numbin', 5);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'gseqspace', 'numbin', 5);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'eqpop'    , 'numbin', 50);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'eqspace'  , 'numbin', 50);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'ceqspace' , 'numbin', 50);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'ibtb', 'histmethod', 'gseqspace', 'numbin', 50);
result{end+1} = ft_connectivity_mutualinformation(dat, 'method', 'gcmi');

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequaln(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

repo: fieldtrip


text:
%% Demo for "Fitit"
% Copyright (c) 2011, The MathWorks, Inc.

%% Generate a set of random data

clear all
clc

s = RandStream('mt19937ar','seed',1971);
RandStream.setDefaultStream(s);

X = linspace(1,10,100);
X = X';

% Specify the parameters for a second order Fourier series

w = .6067;
a0 = 1.6345;
a1 = -.6235;
b1 = -1.3501;
a2 = -1.1622;
b2 = -.9443;

% Fourier2 is the true (unknown) relationship between X and Y
Y = a0 + a1*cos(X*w) + b1*sin(X*w) + a2*cos(2*X*w) + b2*sin(2*X*w);

% Add in a noise vector

K = max(Y) - min(Y);
noisy = Y +  .2*K*randn(100,1);

%% Generate a set of random data

s = RandStream('mt19937ar','seed',1973);
RandStream.setDefaultStream(s);

hold off

X1 = linspace(1,10,100);
X1 = X1';

% Specify the parameters for a second order Fourier series

w = rand;
a0 = randn;
a1 = randn;
b1 = randn;
a2 = randn;
b2 = randn;

% Fourier2 is the true (unknown) relationship between X and Y
Y1 = a0 + a1*cos(X*w) + b1*sin(X*w) + a2*cos(2*X*w) + b2*sin(2*X*w);

% Add in a noise vector

K = max(Y1) - min(Y1);
noisy1 = Y1 +  .2*K*randn(100,1);

foo = fitit(X1,noisy1)

% Compare the output of "fitit" with the original model

hold on
plot(X1,Y1, 'r', 'linewidth', 2)

%%  Working with the fit object

figure
plot(My_Fit)
My_Fit(noisy(1:10))

figure;
subplot(2,1,1); plot(X,My_Fit(X)); title('My Fit');
subplot(2,1,2); plot(X, differentiate(My_Fit, X)); refline(0,0),title('First derivative of My Fit');

methods(My_Fit)

%%  Add confidence intervals into the mix

[avgfit,fitstd1,fitstd2] = fitit(X,noisy,1000);

%%  Explain Localized Regression

figure
scatter(X, noisy)

fit1 = smooth(X,noisy,0.02,'loess')
hold on
plot(X, fit1, 'r', 'linewidth', 2)

fit2 = smooth(X,noisy,0.95,'loess')
plot(X, fit2, 'k', 'linewidth', 2)

fit3 = smooth(X,noisy,0.2,'loess')
plot(X, fit3, 'b', 'linewidth', 2)

%% Explain Cross Validation

% Divide data set into a test set and a training set

cp = cvpartition(size(noisy,1),'k',10);

trainingset1 = X(training(cp,1));
trainingset2 = noisy(training(cp,1));

testset1 = X(test(cp,1));
testset2 = noisy(test(cp,1));

figure
trainingscatter = scatter(trainingset1,trainingset2,'+', 'r');
hold on
testscatter = scatter(testset1, testset2,'FaceColor','b','Filled');
legend('Training Set','Test Set', 'location', 'NorthWest');

%% Use Training Set to create a LOESS smooth

delete(testscatter);
Z = smooth(trainingset1,trainingset2,.5,'loess');
ZPlot = plot(trainingset1, Z, 'color','r','linestyle','-','linewidth',2);

%% Use the test set to evaluate goodness of fit

foobar = fit(trainingset1, Z, 'cubicinterp')
scatter(testset1, testset2,'FaceColor','b','Filled')

for i = 1: length(testset1)

    plot([testset1(i) testset1(i)], [foobar(testset1(i)) testset2(i) ])

end

%%  Explain Bootstrap

figure
scatter(X, noisy, 'b')

% Create a new data set
index = randsample(length(X),length(X), 'true');
hold on
scatter(X(index), noisy(index), 'filled', 'r')
bar = fit(X(index), noisy(index), 'fourier2');
plot(bar)
legend('Original Data Set', 'Random Sample', 'Fit', 'location', 'NorthWest')

%%  Repeat

index = randsample(length(X),length(X), 'true');
scatter(X(index), noisy(index), 'filled', 'k')
bar = fit(X(index), noisy(index), 'fourier2');
plot(bar, 'k')
legend('off')

%% Repeat 20 more times

for i = 1:20

    index = randsample(length(X),length(X), 'true');
    bar = fit(X(index), noisy(index), 'fourier2');
    plot(bar, 'b')

end

legend off


text:
% TRANSFORMCOORDS - Select nazion and inion in anatomical MRI images.
%
% Usage:
%   mewcoords = transformcoords(coords, rotate, scale, center, reverse);
%
% Inputs:
%   coords    - array of 3-D coordinates (3 by N or N by 3)
%   rotate    - [pitch roll yaw] rotate in 3-D using pitch (x plane),
%               roll (y plane) and yaw (z plane). An empty array does
%               not perform any rotation.
%   scale     - [scalex scaley scalez] scale axis. A single numeric
%               input scale all the dimensions the same. Default 1
%               does not scale.
%   shifts    - [x y z] shift coordinates (after rotation and scaling).
%               Default [0 0 0] does not move the center.
%   reverse   - [0|1] when set to 1 perform the reverse transformation,
%               first moving to the old center, unscaling, and unrotating.
%               Default is 0.
%
% Output:
%   newcoords - coordinates after rotating, scaling and recentering
%
% Author: Arnaud Delorme, Salk, SCCN, UCSD, CA, March 23, 2004

% Copyright (C) 2004 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function coords = transformcoords(coords, rotate, scale, center, reverse);

    if nargin < 2
        help transformcoords;
        return;
    end
    if nargin < 3
        scale =1;
    end
    if nargin < 4
        center = [0 0 0];
    end
    if nargin < 5
        reverse = 0;
    end
    if size(coords, 1) ~= 3
        coords = coords';
        trp    = 1;
    else
        trp = 0;
    end
    if size(coords, 1) ~= 3
        error('Number of columns must be 3 for the coordinate input');
    end
    if length(rotate) > 0 && length(rotate) ~= 3
        error('rotate parameter must have 3 values');
    end

    % decode parameters
    % -----------------
    centx = -center(1);
    centy = -center(2);
    centz = -center(3);
    if length(scale) == 1
        scale = [scale scale scale];
    end
    scalex = scale(1);
    scaley = scale(2);
    scalez = scale(3);
    if length(rotate) < 3
        rotate = [0 0 0]
    end
    pitch = rotate(1);
    roll  = rotate(2);
    yaw   = rotate(3);

    if ~reverse
        % pitch roll yaw rotation
        % -----------------------
        % pitch (x-axis); roll = y axis rotation; yaw = z axis
        % see http://bishopw.loni.ucla.edu/AIR5/homogenous.html
        cp = cos(pitch); sp = sin(pitch);
        cr = cos(roll);  sr = sin(roll);
        cy = cos(yaw);   sy = sin(yaw);
        rot3d = [ cy*cr+sy*sp*sr    sy*cr-cy*sp*sr     cp*sr  ;
                  -sy*cp            cy*cp               sp     ;
                  sy*sp*cr-cy*sr    -cy*sp*cr-sy*sr     cp*cr  ];
        coords = rot3d*coords;

        % scaling and centering
        % ---------------------
        coords(1,:) = coords(1,:)*scalex-centx;
        coords(2,:) = coords(2,:)*scaley-centy;
        coords(3,:) = coords(3,:)*scalez-centz;
    else
        % unscaling and uncentering
        % -------------------------
        coords(1,:) = (coords(1,:)+centx)/scalex;
        coords(2,:) = (coords(2,:)+centy)/scaley;
        coords(3,:) = (coords(3,:)+centz)/scalez;

        % pitch roll yaw rotation
        % -----------------------
        cp = cos(-pitch); sp = sin(-pitch);
        cr = cos(-roll);  sr = sin(-roll);
        cy = cos(-yaw);   sy = sin(-yaw);
        rot3d = [ cy*cr+sy*sp*sr    sy*cr-cy*sp*sr     cp*sr  ;
                  -sy*cp            cy*cp               sp     ;
                  sy*sp*cr-cy*sr    -cy*sp*cr-sy*sr     cp*cr  ];
        coords = rot3d*coords;
    end

    if trp
        coords = coords';
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function TestPhoneRecognizerLSTM_TIMIT()

modeldir = my_dir('nnet');
modelfiles = findFiles(['nnet/' modeldir{1}], 'mat');
modelfiles = sort_nnet_by_itr(modelfiles);      % get the last iteration model

dnn = load(modelfiles{1});
para = dnn.para;
para.IO = RemoveIOStream(para.IO, 2);   % during evaluation, we don't have the label usually. So remove that stream from the configuration
layer = dnn.layer(1:end-2);     % we discard the last two layers that is not useful for ASR

clean_cond = 1;
[Data_cv, ~, para] = LoadData_TIMIT(para, 'test');

para.out_layer_idx = length(layer) + [0];   % you can specify which layers' activation will be outputed

for i=1:length(Data_cv(1).data)
    tmpData(1).data{1} = Data_cv(1).data{i};
    output = FeatureTree2(tmpData, para, layer);
    posterior = output{1}{1};

    subplot(2,1,1); imagesc(Data_cv(1).data{i});    title('Feature');
    subplot(2,1,2); imagesc(posterior);       title('Posteriorgram and true label'); hold on
    plot(Data_cv(2).data{i}, 'r'); hold off;
    pause
end

end

author: singaxiong

repo: SignalGraph


text:
function tri_surface_write ( node_file_name, triangle_file_name, dim_num, ...
  node_num, order_num, triangle_num, node_xyz, triangle_node )

%*****************************************************************************80
%
%% TRI_SURFACE_WRITE writes graphics information to a pair of TRI_SURFACE files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILE_NAME, the name of the node file.
%
%    Input, string TRIANGLE_FILE_NAME, the name of the triangle file.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, integer ORDER_NUM, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, real NODE_XYZ(DIM_NUM,NODE_NUM), the node coordinates.
%
%    Input, integer TRIANGLE_NODE(ORDER_NUM,TRIANGLE_NUM),
%    the nodes that form the triangles.
%
  header = 0;

  r8mat_write ( node_file_name, dim_num, node_num, node_xyz );

  i4mat_write ( triangle_file_name, order_num, triangle_num, ...
    triangle_node, header );

  return
end

repo: jburkardt-m


text:
function F = conj( F )
%CONJ Complex conjugate of a SPHEREFUNV.
%   CONJ(F) returns the complex conjugate of F. For a complex F, CONJ(F) =
%   REAL(F) - i*IMAG(F).
%
%   Since only real SPHEREFUNV objects are supported, the result of this
%   function is always just the original function F.
%
%   See also SPHEREFUNV/IMAG, SPHEREFUNV/REAL.

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

% Empty check:
if ( isempty( F ) )
    return
end

% Take conj part of each component:
F.components = cellfun( @conj, F.components, 'UniformOutput', false );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y = reflect(x,minx,maxx)

% function y = reflect(x,minx,maxx)
% Reflect the values in matrix x about the scalar values minx and maxx.
% Hence a vector x containing a long linearly increasing series
% is converted into a waveform which ramps linearly up and down
% between minx and maxx.
% If x contains integers and minx and maxx are (integers + 0.5),
% the ramps will have repeated max and min samples.
%
% Nick Kingsbury, Cambridge University, January 1999.

y = x;

% Reflect y in maxx.
t = find(y > maxx);
y(t) = 2*maxx - y(t);

% Reflect y in minx.
t = find(y < minx);
while ~isempty(t)
   y(t) = 2*minx - y(t);

   % Repeat until no more values out of range.
   t = find(y > maxx);
   if ~isempty(t), y(t) = 2*maxx - y(t); end
   t = find(y < minx);
end

return;

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%A case of Buck DC-DC converter
C=440;
L=20;
RC=0.002;
RL=0.001;
Vor=5;
Vin=24;
Vc0=5;
ILmin=0;


text:
function h = gaussianKernel3D( r, sigma)
%Create pre-defined 3-D Gaussian kernel.
%
%   Description :
%
%   h = gaussianKernel3D(r,sigma) returns a rotationally symmetric Gaussian
%   kernel h of size 2*r+1 with standard deviation sigma (positive). The
%   default size of r is 1 and the default sigma is 0.5.
%
%   Author : Mohammad Mustafa
%   By courtesy of University of Nottingham and Mirada Medical Limited,
%   Oxford, UK
%
%   Published under a Creative Commons Attribution-Non-Commercial-Share Alike
%   3.0 Unported Licence http://creativecommons.org/licenses/by-nc-sa/3.0/
%
%   June 2012

    if nargin<1
        r=1; sigma=0.5;
    elseif nargin==1
        sigma=0.5;
    end
    [x,y,z] = ndgrid(-r:r,-r:r,-r:r);
    arg   = -(x.*x + y.*y + z.*z)/(2*sigma*sigma);

    h = exp(arg);

    if sum(h(:)) ~= 0,
      h  = h/sum(h(:));
    end;

end


text:
function [ v_weight, predecessor ] = bellman_ford ( v_num, e_num, source, ...
  e, e_weight )

%*****************************************************************************80
%
%  Purpose:
%
%    BELLMAN_FORD finds shortest paths from a given vertex of a weighted directed graph.
%
%  Discussion:
%
%    The Bellman-Ford algorithm is used.
%
%    Each edge of the graph has a weight, which may be negative.  However,
%    it should not be the case that there is any negative loop, that is,
%    a circuit whose total weight is negative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer V_NUM, the number of vertices.
%
%    Input, integer E_NUM, the number of edges.
%
%    Input, integer SOURCE, the vertex from which distances will
%    be calculated.
%
%    Input, integer E(2,E_NUM), the edges, given as pairs of
%    vertex indices.
%
%    Input, real E_WEIGHT(E_NUM), the weight of each edge.
%
%    Output, real V_WEIGHT(V_NUM), the weight of each node,
%    that is, its minimum distance from SOURCE.
%
%    Output, integer PREDECESSOR(V_NUM), a list of predecessors,
%    which can be used to recover the shortest path from any node back to SOURCE.
%
  r8_big = 1.0E+30;
%
%  Step 1: initialize the graph.
%
  v_weight(1:v_num) = r8_big;
  v_weight(source) = 0.0;

  predecessor(1:v_num) = -1;
%
%  Step 2: Relax edges repeatedly.
%
  for i = 1 : v_num - 1
    for j = 1 : e_num
      u = e(2,j);
      v = e(1,j);
      t = v_weight(u) + e_weight(j);
      if ( t < v_weight(v) )
        v_weight(v) = t;
        predecessor(v) = u;
      end
    end
  end
%
%  Step 3: check for negative-weight cycles
%
  for j = 1 : e_num
    u = e(2,j);
    v = e(1,j);
    if ( v_weight(u) + e_weight(j) < v_weight(v) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'BELLMAN_FORD - Fatal error!\n' );
      fprintf ( 1, '  Graph contains a cycle with negative weight.\n' );
      error ( 'BELLMAN_FORD - Fatal error!' );
    end
  end

  return
end


text:
function initialize
% add paths for Thomas' subdirectories
addpath('~/zmap/src/thomas:', ...
     '~/zmap/src/thomas/synthetic:', ...
     '~/zmap/src/thomas/decluster:', ...
     '~/zmap/src/thomas/decluster/gardnerknopoff:', ...
     '~/zmap/src/thomas/decluster/marsan:', ...
     '~/zmap/src/thomas/decluster/reasen:', ...
     '~/zmap/src/thomas/etas:', ...
     '~/zmap/src/thomas/gui:', ...
     '~/zmap/src/thomas/plot:', ...
     '~/zmap/src/thomas/pov:', ...
     '~/zmap/src/thomas/seismicrates:', ...
     '~/zmap/src/thomas/seismicrates/ex:', ...
     '~/zmap/src/thomas/synthetic:', ...
     '~/zmap/src/thomas/ex:', ...
     '~/zmap/results:', ...
     '~/zmap/eq_data:', ...
     '~/zmap/eq_data/Alaska:', ...
     '~/zmap/eq_data/Adak:', ...
     '~/zmap/eq_data/SCEC:', ...
     '~/zmap/eq_data/Shearer:', ...
     '~/zmap/eq_data/Taiwan:')

author: CelsoReyes


text:
function nlp = Step1Constraints(nlp, bounds, varargin)
    domain = nlp.Plant;

    %% virtual constraints
    nlp = CstrFcns.ImposeVirtualConstraint(nlp, bounds);

    %% initial position
    q0 = bounds.initialPos;
    dq0 = bounds.initialVel;

    updateVariableProp(nlp,'x','first','lb',q0,'ub',q0,'x0',q0);
    updateVariableProp(nlp,'dx','first','lb',dq0,'ub',dq0,'x0',dq0);

    %% torso

    torso = CoordinateFrame('Name','torso',...
        'Reference',nlp.Plant.Joints(end),...
        'Offset',[0,0,0.3],...
        'R',[0,0,0]);
    rpy = getEulerAngles(domain, torso);
    constr_fun = SymFunction(['torso_',domain.Name],rpy(2),{domain.States.x});
    lb = bounds.torso.lb;
    ub = bounds.torso.ub;
    addNodeConstraint(nlp,constr_fun,{'x'},'all', lb, ub,'Linear');
    %% initial position
    %     q0 = bounds.restPos;
    %     dq_minus = zeros(6,1);
    %     %     dq0 = bounds.restVel;
    %     %     ddq0 = bounds.restAcc;
    %     robot = nlp.Plant;
    %     D = robot.calcMassMatrix(q0);
    %     Jfoot = Jh_foot_acrobot(q0);
    %     Jimp = jac_push_acrobot(q0);
    %
    %     E = [D -Jfoot';Jfoot, zeros(3,3)];
    %     F = [D*dq_minus + Jimp'*bounds.push_force;
    %         zeros(3,1)];
    %
    %     y = E\F;
    %
    %     dq0 = y(1:6);
    %
    %
    %     updateVariableProp(nlp,'x','first','lb',q0,'ub',q0,'x0',q0);
    %     updateVariableProp(nlp,'dx','first','lb',dq0,'ub',dq0,'x0',dq0);
    %     updateVariableProp(nlp,'ddx','first','lb',ddq0,'ub',ddq0,'x0',ddq0);

    %     updateVariableProp(nlp,'x','last','lb',q0,'ub',q0,'x0',q0);
    %     updateVariableProp(nlp,'dx','last','lb',dq0,'ub',dq0,'x0',dq0);
    %     updateVariableProp(nlp,'ddx','last','lb',ddq0,'ub',ddq0,'x0',ddq0);

    %     D = nlp.Plant.Mmat;
    %     q = nlp.Plant.States.x;
    %     dq = nlp.Plant.States.dx;
    %
    %     f_push = CoordinateFrame('Name','push',...
    %         'Reference',nlp.Plant.Joints(end),...
    %         'Offset',[0,0,0.3],...
    %         'R',[0,0,0]);
    %
    %     %     p_push_sym = getCartesianPosition(robot,f_push);
    %     J_push_sym = getSpatialJacobian(nlp.Plant,f_push);
    %     J = J_push_sym([1,3],:);
    %
    %     F = SymVariable('f',[2,1]);
    %
    %     nlp.addVariable('fg','first','Name','fg','Dimension',3,...
    %         'lb',[-inf,0,-inf],'ub',[inf,inf,inf],'x0',[0,0,0]);
    %     Fg = SymVariable('fg',[3,1]);
    %     Jg = nlp.Plant.Gmap.ConstraintWrench.ffoot;
    %     constr = D*dq - J.'*F + Jg*Fg;
    %     constr_fun = SymFunction(['impact_',nlp.Plant.Name],constr,{q,dq,Fg},{F});
    %     addNodeConstraint(nlp,constr_fun,{'x','dx','fg'},'first',0,0,'Nonlinear',bounds.push_force);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [x,P,K] = UnscentedKalmanFilterX_Update(xPred,PPred,y,yPred,S,Pxy)
% UNSCENTEDKALMANFILTERX_UPDATE Perform the discrete-time UKF update step,
% under the assumption of additive process noisem for a single measurement.
% This is essentially identical to a KF Update step.
%
% Parameters
% ----------
% xPred: column vector
%   The (xDim x 1) predicted state estimate.
% PPred: matrix
%   The (xDim x xDim) predicted state covariance matrix.
% y: column vector
%   The (yDim x 1) measurement vector.
% yPred: column vector
%   The (yDim x 1) predicted measurement estimate.
% S: matrix
%   The (yDim x yDim) innovation covariance matrix.
% Pxy: matrix
%   The (xDim x yDim) cross-covariance matrix.
%
% Returns
% -------
% x: column vector
%   The (xDim x 1) state estimate at the current time-step.
% P: matrix
%   The (xDim x xDim) state covariance matrix at the current
%   time-step.
% K: matrix
%   The (xDim x yDim) Kalman gain matrix at the current
%   time-step.
%
% October 2017 Lyudmil Vladimirov, University of Liverpool.

    [x,P,K] = KalmanFilterX_Update(xPred,PPred,y,yPred,S,Pxy);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [xUpdate,SUpdate,innov,Szz,W]=sqrtCubKalUpdate(xPred,SPred,z,SR,h,xi,w,innovTrans,measAvgFun,stateDiffTrans,stateTrans)
%SQRTCUBKALUPDATE Perform the measurement update step in the square root
%                 cubature Kalman filter with additive measurement noise.
%                 Unlike the non-square root version, the covariance update
%                 is such that finite precision problems cannot lead to a
%                 non-positive (semi-)definite covariance. However, unlike
%                 the non-square root version, all of the cubature weights
%                 must be positive.
%
%INPUTS: xPred The xDim X 1 predicted target state.
%        SPred The xDim X xDim lower-triangular square root predicted state
%              covariance matrix.
%            z The zDim X 1  measurement vector.
%           SR The zDim X zDim lower-triangular square root of the
%              measurement covariance matrix in the native coordinate
%              system of the measurement.
%            h A function handle for the measurement function that the
%              state as its argument.
%           xi An xDim X numCubPoints matrix of cubature points. If this
%              and the next parameter are omitted or empty matrices are
%              passed, then fifthOrderCubPoints(xDim). It is suggested that
%              xi and w be provided to avoid needless recomputation of the
%              cubature points.
%            w A numCubPoints X 1 vector of the weights associated with the
%              cubature points.
%   innovTrans An optional function handle that computes and optionally
%              transforms the value of the difference between the
%              observation and any predicted points. This is called as
%              innovTrans(a,b) and the default if omitted or an empty
%              matrix is passed is @(a,b)bsxfun(@minus,a,b). This must be
%              able to handle sets of values. For a zDimX1 measurement,
%              either of the inputs could be zDimXN in size while one of
%              the inputs could be zDimX1 in size.  This only needs to be
%              supplied when a measurement difference must be restricted
%              to a certain range. For example, the innovation between two
%              angles will be 2*pi if one angle is zero and the other
%              2*pi, even though they are the same direction. In such an
%              instance, a function handle to the
%              wrapRange(bsxfun(@minus,a,b),-pi,pi) function with the
%              appropriate parameters should be passed for innovTrans.
%   measAvgFun An optional function handle that, when given N measurement
%              values with weights, produces the weighted average. This
%              function only has to be provided if the domain of the
%              measurement is not linear. For example, when averaging
%              angular values, then the function meanAng should be used.
% stateDiffTrans An optional function handle that takes an xDimXN matrix of
%              N differences between states and transforms them however
%              might be necessary. If not transformation is necessary, this
%              parameter can be omitted or an empty matrix passed.
%   stateTrans An optional function that takes a state estimate and
%              transforms it. This is useful if one wishes the elements of
%              the state to be bound to a certain domain. For example, if
%              an element of the state is an angle, one might generally
%              want to bind it to the region +/-pi.
%
%OUTPUTS: xUpdate The xDim X 1 updated state vector.
%         SUpdate The updated xDim X xDim lower-triangular square root
%                 state covariance matrix.
%      innov, Szz The zDimX1 innovation and the zDimXzDim square root
%                 innovation covariance matrix are returned in case one
%                 wishes to analyze the consistency of the estimator or use
%                 those values in gating or likelihood evaluation.
%               W The xDimXzDim gain used in the update. This can be
%                 useful when gating and using the function
%                 calcMissedGateCov.
%
%If the function h needs additional parameters beyond the state, then the
%parameters can be passed by using an anonymous function as the function
%handle. For example, suppose that the measurement function is measFunc and
%it needs the additional parameters param1 and param2. In this instance,
%rather than using
%h=@measFunc
%one should use
%h=@(x)measFunc(x,param1,param2)
%This way, every time sqrtCubKalUpdate calls measFunc (via h) with a
%different x, those two parameters are always passed.
%
%The mathematics behind the function sqrtCubKalUpdate are described in more
%detail in Section IX of [1] and in [2]. Note that this is essentially one
%type of square root "unscented Kalman filter" with additive noise. One
%simply has to provide the filter with the appropriate cubature points and
%weights.
%
%The optional parameters innovTrans and measAvgFun are not described in
%references [1] and [2], but allow for possible modifications to the filter
%as described in [3] The parameters have been added to allow the filter to
%be used with angular quantities. For example, if the measurement consisted
%of range and angle, z=[r;theta], then
%innovTrans=@(a,b)[bsxfun(@minus,a(1,:),b(1,:));
%                  wrapRange(bsxfun(@minus,a(2,:),b(2,:)),-pi,pi)];
%measAvgFun=@(z,w)[calcMixtureMoments(z(1,:),w);
%                  meanAng(z(2,:),w')];
%should be used to approximately deal with the circular nature of the
%measurements.
%
%REFERENCES:
%[1] D. F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems Magazine,
%    vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%[2] I. Arasaratnam and S. Haykin, "Cubature Kalman filters," IEEE
%    Transactions on Automatic Control, vol. 54, no. 6, pp. 1254-1269,
%    Jun. 2009.
%[3] D. F. Crouse, "Cubature/ unscented/ sigma point Kalman filtering with
%    angular measurement models," in Proceedings of the 18th International
%    Conference on Information Fusion, Washington, D.C., 6-9 Jul. 2015.
%
%July 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    xDim=size(xPred,1);

    if(nargin<6||isempty(xi))
        [xi,w]=fifthOrderCubPoints(xDim);
    end

    if(nargin<8||isempty(innovTrans))
        %The function just returns the input.
        innovTrans=@(a,b)bsxfun(@minus,a,b);
    end

    if(nargin<9||isempty(measAvgFun))
        measAvgFun=@(zPoints,w)calcMixtureMoments(zPoints,w);
    end

    if(nargin<10||isempty(stateDiffTrans))
        stateDiffTrans=@(x)x;
    end

    if(nargin<11||isempty(stateTrans))
        stateTrans=@(x)x;
    end

    zDim=size(z,1);
    numCubPoints=size(xi,2);
    sqrtW=sqrt(w);
    %Calculate the updated estimates

    %Predicted cubature state points
    xPredPoints=stateTrans(transformCubPoints(xi,xPred,SPred));

    %Predicted, centered cubature state points
    xPredCenPoints=bsxfun(@times,stateDiffTrans(bsxfun(@minus,xPredPoints,xPred)),sqrtW');

    %Predicted cubature measurement points
    zPredPoints=zeros(zDim,numCubPoints);
    for curP=1:numCubPoints
        zPredPoints(:,curP)=h(xPredPoints(:,curP));
    end

    %Measurement prediction.
    zPred=measAvgFun(zPredPoints,w);

    %Centered, predicted cubature measurement points, transformed as
    %necessary to keep the values within a desired range.
    zPredCenPoints=bsxfun(@times,innovTrans(zPredPoints,zPred),sqrtW');

    %Root innovation covariance
    Szz=tria([zPredCenPoints,SR]);

    %The cross covariance.
    Pxz=xPredCenPoints*zPredCenPoints';

    %The filter gain
    W=(Pxz/Szz')/Szz;

    %The innovation, transformed as necessary to keep values in a desired
    %range.
    innov=innovTrans(z,zPred);

    %Updated state estimate
    xUpdate=stateTrans(xPred+W*innov);

    %Updated state root covariance
    SUpdate=tria([stateDiffTrans(xPredCenPoints-W*zPredCenPoints),W*SR]);
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:
function SO3VF = rdivide(SO3VF1, SO3VF2)
% overloads |SO3VF1 ./ SO3VF2|
%
% Syntax
%   SO3VF = SO3VF1 ./ SO3VF2
%   SO3VF = SO3VF1 ./ a
%   SO3VF = a ./ SO3VF2
%   SO3VF = SO3F .* SO3VF1
%   SO3VF = SO3VF1 .* SO3F
%
% Input
%  SO3VF1, SO3VF2 - @SO3VectorFieldHarmonic
%  a - double
%  SO3F - @SO3Fun
%
% Output
%  SO3VF - @SO3VectorFieldHarmonic
%

if isa(SO3VF2,'vector3d')
  SO3VF2 = SO3VF2.xyz.';
end

if isnumeric(SO3VF1)
  SO3VF = SO3VectorFieldHandle(@(rot) SO3VF1 ./ SO3VF2.eval(rot),SO3VF2.SRight,SO3VF2.SLeft);
  SO3VF = SO3VectorFieldHarmonic(SO3VF, 'bandwidth', min(getMTEXpref('maxSO3Bandwidth'),2*SO3VF2.bandwidth));
  return
end
if isnumeric(SO3VF2)
  SO3VF = times(SO3VF1,1./SO3VF2);
  return
end

ensureCompatibleSymmetries(SO3VF1,SO3VF2);
SO3VF = SO3VectorFieldHandle(@(rot) SO3VF1.eval(rot)./ SO3VF2.eval(rot),SO3VF1.SRight,SO3VF1.SLeft);
SO3VF = SO3VectorFieldHarmonic(SO3VF, 'bandwidth', min(getMTEXpref('maxSO3Bandwidth'),2*max(SO3VF1.bandwidth, SO3VF2.bandwidth)));

end

repo: mtex


text:
function fiff_write_int(fid,kind,data)
%
% fiff_write_int(fid,kind,data)
%
% Writes a 32-bit integer tag to a fif file
%
%     fid           An open fif file descriptor
%     kind          Tag kind
%     data          The integers to use as data
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.3  2006/04/27 22:38:37  msh
%   Splitting an empty list now results in an empty output.
%   Added fiff_write_double and fiff_write_short
%   Write an array of floats, ints, and doubles instead of just one value.
%   Fixed help text of fiff_pick_channels.
%
%   Revision 1.2  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%

me='MNE:fiff_write_int';

if nargin ~= 3
        error(me,'Incorrect number of arguments');
end

FIFFT_INT=3;
FIFFV_NEXT_SEQ=0;
nel=numel(data);
datasize=nel*4;
count = fwrite(fid,int32(kind),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFT_INT),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(datasize),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFV_NEXT_SEQ),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(data),'int32');
if count ~= nel
    error(me,'write failed');
end

return;

author: fieldtrip


text:
function [files, status] = urldir(page, tag)
%
% Lists the files within a URL
%
% Collecion of folders:
% files = urldir(page, 'DIR')
% Collection of files:
% files = urldir(page, 'TXT')

% Created: A. Torralba, 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LabelMe, the open annotation tool
% Contribute to the database by labeling objects using the annotation tool.
%
% CSAIL, MIT
% 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin == 1
    tag = '[dir]';
    tag = '/'
else
    tag = lower(tag);
    if strcmp(tag, 'dir')
        tag = '/';
    end
    if strcmp(tag, 'txt')
        tag = 'xml';
    end
    if strcmp(tag, 'img')
        tag = 'jpg';
    end
end

nl = length(tag);
nfiles = 0;
files = [];

% Read page
page = strrep(page, '\', '/');
status =0;
trials = 0;
while status == 0 && trials<10
    [webpage, status] = urlread(page, 'Timeout', 5);
    if status == 0
        disp(sprintf('Warning: Failure urldir %s. Trying again...',page))
        drawnow
        trials = trials+1;
    end
end

if trials == 10
    disp('warning: skipping this folder...')
    status = 1;
    return
end

if status
    % Parse page
    j1 = findstr(lower(webpage), '<a href="');
    j2 = findstr(lower(webpage), '');
    Nelements = length(j1);
    if Nelements>0
        for f = 1:Nelements
            % get HREF element
            chain = webpage(j1(f):j2(f));
            jc = findstr(lower(chain), '">');
            if ~isempty(jc)
                chain = deblank(chain(10:jc(1)-1));

                % check if it is the right type
                if length(chain)>length(tag)-1
                    if strcmp(chain(end-nl+1:end), tag)
                        nfiles = nfiles+1;
                        chain = strrep(chain, '%20', ' '); % replace space character
                        files(nfiles).name = chain;
                        files(nfiles).bytes = 1;
                    end
                end
            end
        end
    end
end

author: CSAILVision

repo: LabelMeToolbox