text:
function gpcf = gpcf_squared(varargin)
%GPCF_SQUARED  Create a squared (dot product) covariance function
%
%  Description
%    GPCF = GPCF_SQUARED('PARAM1',VALUE1,'PARAM2,VALUE2,...) creates
%    a squared (dot product) covariance function structure in which
%    the named parameters have the specified values. Any unspecified
%    parameters are set to default values. The squared covariance function
%    corresponds to x.^2 mean function and the respective covariance matrix
%    is given as C = x.^2*diag(gpcf.coeffSigma2)*(x'.^2);
%
%    GPCF = GPCF_SQUARED(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for squared (dot product) covariance function
%      interactions      - twoway interactions (default off)
%      coeffSigma2       - prior variance for regressor coefficients [10]
%                          This can be either scalar corresponding
%                          to a common prior variance or vector
%                          defining own prior variance for each
%                          coefficient.
%      coeffSigma2_prior - prior structure for coeffSigma2 [prior_logunif]
%      selectedVariables - vector defining which inputs are used [all]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*, MEAN_*
%
% Copyright (c) 2007-2016 Jarno Vanhatalo
% Copyright (c) 2008-2010 Jaakko Riihimäki
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2014 Arno Solin

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

  ip=inputParser;
  ip.FunctionName = 'GPCF_SQUARED';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('interactions', 'off', @(x) ismember(x,{'on' 'off'}))
  ip.addParamValue('coeffSigma2',10, @(x) isvector(x) && all(x>0));
  ip.addParamValue('coeffSigma2_prior',prior_logunif, @(x) isstruct(x) || isempty(x));
  ip.addParamValue('selectedVariables',[], @(x) isvector(x) && all(x>0));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

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

  if init || ~ismember('interactions',ip.UsingDefaults)
    gpcf.interactions=ip.Results.interactions;
  end

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

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

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_squared_pak;
    gpcf.fh.unpak = @gpcf_squared_unpak;
    gpcf.fh.lp = @gpcf_squared_lp;
    gpcf.fh.lpg = @gpcf_squared_lpg;
    gpcf.fh.cfg = @gpcf_squared_cfg;
    gpcf.fh.cfdg = @gpcf_squared_cfdg;
    gpcf.fh.cfdg2 = @gpcf_squared_cfdg2;
    gpcf.fh.ginput = @gpcf_squared_ginput;
    gpcf.fh.ginput2 = @gpcf_squared_ginput2;
    gpcf.fh.ginput3 = @gpcf_squared_ginput3;
    gpcf.fh.ginput4 = @gpcf_squared_ginput4;
    gpcf.fh.cov = @gpcf_squared_cov;
    gpcf.fh.trcov  = @gpcf_squared_trcov;
    gpcf.fh.trvar  = @gpcf_squared_trvar;
    gpcf.fh.recappend = @gpcf_squared_recappend;
    gpcf.fh.cf2ss = @gpcf_squared_cf2ss;
  end

end

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

  w = []; s = {}; h =[];
  if ~isempty(gpcf.p.coeffSigma2)
    w = log(gpcf.coeffSigma2);
    if numel(gpcf.coeffSigma2)>1
      s = [s; sprintf('log(squared.coeffSigma2 x %d)',numel(gpcf.coeffSigma2))];
    else
      s = [s; 'log(squared.coeffSigma2)'];
    end
    h = [h ones(1, numel(gpcf.coeffSigma2))];
    % Hyperparameters of coeffSigma2
    [wh, sh, hh] = gpcf.p.coeffSigma2.fh.pak(gpcf.p.coeffSigma2);
    sh=strcat(repmat('prior-', size(sh,1),1),sh);
    w = [w wh];
    s = [s; sh];
    h = [h 1+hh];
  end
end

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

  gpp=gpcf.p;

  if ~isempty(gpp.coeffSigma2)
    i2=length(gpcf.coeffSigma2);
    i1=1;
    gpcf.coeffSigma2 = exp(w(i1:i2));
    w = w(i2+1:end);

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

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

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

  if ~isempty(gpp.coeffSigma2)
    lp = gpp.coeffSigma2.fh.lp(gpcf.coeffSigma2, gpp.coeffSigma2) + sum(log(gpcf.coeffSigma2));
  end
end

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

  lpg = [];
  gpp=gpcf.p;

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

function DKff = gpcf_squared_cfg(gpcf, x, x2, mask, i1)
%GPCF_SQUARED_CFG  Evaluate gradient of covariance function
%                 with respect to the parameters
%
%  Description
%    DKff = GPCF_SQUARED_CFG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_SQUARED_CFG(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_SQUARED_CFG(GPCF, X, [], MASK) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the diagonal of gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_SQUARED_CFG(GPCF,X,X2,MASK,i) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradient of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    hyperparameter. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%   GPCF_SQUARED_PAK, GPCF_SQUARED_UNPAK, GPCF_SQUARED_LP, GP_G

  if nargin>2 && ~isempty(x2)
      if size(x,2) ~= size(x2,2)
          error('gpcf_squared -> _cfg: the number of columns of X1 and X2 has to be same')
      end
  end

  if isfield(gpcf, 'selectedVariables')
      x = x(:,gpcf.selectedVariables);
  end
  [n, m] =size(x);
  h = x.^2;
  if isequal(gpcf.interactions,'on')
      for xi1=1:m
          for xi2=xi1+1:m
              h = [h x(:,xi1).*x(:,xi2)];
          end
      end
  end

  DKff = {};

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

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

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

      if ~isempty(gpcf.p.coeffSigma2)
          if length(gpcf.coeffSigma2) == 1
              DKff{1}=gpcf.coeffSigma2*h*(h');
          else
              if isa(gpcf.coeffSigma2,'single')
                  epsi=eps('single');
              else
                  epsi=eps;
              end
              if ~savememory
                  i1=1:length(gpcf.coeffSigma2);
              end
              DKff=cell(1,length(i1));
              for ii1=i1
                  DD = gpcf.coeffSigma2(ii1)*h(:,ii1)*(h(:,ii1)');
                  DD(abs(DD)<=epsi) = 0;
                  DKff{ii1}= (DD+DD')./2;
              end
          end
      end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 || isempty(mask)

      error('this part of the subfunction has not been tested')

      if isfield(gpcf, 'selectedVariables')
          x2 = x2(:,gpcf.selectedVariables);
      end
      h2 = x2.^2;
      if isequal(gpcf.interactions,'on')
          for xi1=1:m
              for xi2=xi1+1:m
                  h2 = [h2 x2(:,xi1).*x2(:,xi2)];
              end
          end
      end

      if ~isempty(gpcf.p.coeffSigma2)
          if length(gpcf.coeffSigma2) == 1
              DKff{1}=gpcf.coeffSigma2*h*(h2');
          else
              if ~savememory
                  i1=1:m;
              end
              for ii1=i1
                  DKff{ii1}=gpcf.coeffSigma2(ii1)*h(:,ii1)*(h2(:,ii1)');
              end
          end
      end
      % Evaluate: DKff{1}    = d mask(Kff,I) / d coeffSigma2
    %           DKff{2...} = d mask(Kff,I) / d coeffSigma2
  elseif nargin == 4 || nargin == 5

      error('this part of the subfunction has not been tested')

      if ~isempty(gpcf.p.coeffSigma2)
        if length(gpcf.coeffSigma2) == 1
          DKff{1}=gpcf.coeffSigma2*sum(h.^2,2); % d mask(Kff,I) / d coeffSigma2
        else
          if ~savememory
            i1=1:m;
          end
          for ii1=i1
            DKff{ii1}=gpcf.coeffSigma2(ii1)*(h(:,ii1).^2); % d mask(Kff,I) / d coeffSigma2
          end
        end
      end
  end
  if savememory
    DKff=DKff{i1};
  end
end

function DKff = gpcf_squared_cfdg(gpcf, x, x2, dims)
%GPCF_SQUARED_CFDG  Evaluate gradient of covariance function, of
%                  which has been taken partial derivative with
%                  respect to x, with respect to parameters.
%
%  Description
%    DKff = GPCF_SQUARED_CFDG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of derivatived covariance matrix
%    dK(df,f)/dhyp = d(d k(X,X)/dx)/dhyp, with respect to the
%    parameters
%
%    Evaluate: DKff{1:m} = d Kff / d coeffSigma2
%    m is the dimension of inputs. This subfunction is needed when using
%    derivative observations.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_SQUARED_GINPUT

[~,m]=size(x);
if nargin<3
    x2=x;
end
if nargin < 4 || isempty(dims)
    dims = 1:m;
end
ii1=0;
DKff={};
if isfield(gpcf,'selectedVariables')
    selVars = gpcf.selectedVariables;
else
    selVars = 1:m;
end
c = gpcf.coeffSigma2;

h = x.^2;
h2 = x2.^2;

ii1=0;
DKff={};
if ~isempty(gpcf.p.coeffSigma2)
  if length(gpcf.coeffSigma2)==1
    % One coeffSigma2
    for i1=dims
      if isfield(gpcf, 'selectedVariables') && sum(gpcf.selectedVariables==i1)==0
        DK{i1}=zeros(size(x,1),size(x2,1));
      else
        DK{i1}=c(1).*2*x(:,i1)*h2(:,i1)';
        if isequal(gpcf.interactions,'on')
            for xi2=selVars
                if xi2~=i1
                    DK{i1} = DK{i1} + c(1)*x(:,xi2)*(x2(:,i1).*x2(:,xi2))';
                end
            end
        end
      end
    end
    ii1=ii1+1;
    DKff{ii1}=cat(1,DK{1:end});
  else
      % vector of coeffSigma2s
      for i1=1:length(selVars)
          for j=dims
              if selVars(i1)~=j %|| (isfield(gpcf, 'selectedVariables') && sum(gpcf.selectedVariables==i1)==0)
                  DK{j}=zeros(size(x,1),size(x2,1));
              else
                  DK{j}=c(i1).*2*x(:,j)*h2(:,j)';
              end
          end
          ii1=ii1+1;
          DKff{ii1}=cat(1,DK{1:end});
      end
      if isequal(gpcf.interactions,'on')
          for xi1=1:length(selVars)
              for xi2=xi1+1:length(selVars)
                  i1=i1+1;
                  for j=dims
                      if j==xi1
                          DK{j} = c(i1)*x(:,xi2)*(x2(:,xi1).*x2(:,xi2))';
                      elseif j==xi2
                          DK{j} = c(i1)*x(:,xi1)*(x2(:,xi1).*x2(:,xi2))';
                      else
                          DK{j}=zeros(size(x,1),size(x2,1));
                      end
                  end
                  ii1=ii1+1;
                  DKff{ii1}=cat(1,DK{1:end});
              end
          end
      end
  end
end

end

function DKff = gpcf_squared_cfdg2(gpcf, x, x2, dims1, dims2)
%GPCF_SQUARED_CFDG2  Evaluate gradient of covariance function, of which has
%                    been taken partial derivatives with respect to both
%                    input variables x, with respect to parameters.
%
%  Description
%    DKff = GPCF_SQUARED_CFDG2(GPCF, X) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of derivative covariance matrix
%    dK(df,df)/dhyp = d(d^2 k(X1,X2)/dX1dX2)/dhyp with respect to
%    the parameters
%
%    Evaluate: DKff{1:m} = d K(df,df) / d coeffSigma
%    m is the dimension of inputs.  This subfunction is needed when using
%    derivative observations.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%   GPCF_SQUARED_GINPUT, GPCF_SQUARED_GINPUT2

[~, m] =size(x);
if nargin <3 || isempty(x2)
    x2=x;
end
if nargin < 4 || isempty(dims1)
    %dims1 = 1:m;
    error('dims1 needs to be given')
end
if nargin < 5 || isempty(dims2)
    %dims2 = 1:m;
    error('dims2 needs to be given')
end
if isfield(gpcf,'selectedVariables')
    selVars = gpcf.selectedVariables;
else
    selVars = 1:m;
end
% NOTICE. AS OF NOW we assume that dims1 and dims2 are scalars

ii1=0;
if length(gpcf.coeffSigma2)==1
    c=repmat(gpcf.coeffSigma2,1,(1+length(selVars))*length(selVars)/2);
else
    c=gpcf.coeffSigma2;
end
if length(gpcf.coeffSigma2)==1
    % One coeffSigma2
    if dims1~=dims2 %|| (isfield(gpcf, 'selectedVariables') && (sum(gpcf.selectedVariables==j)==0 || sum(gpcf.selectedVariables==k)==0))
        DK=zeros(size(x,1),size(x2,1));
        if isequal(gpcf.interactions,'on') && (sum(selVars==dims1)>0 || sum(selVars==dims2)>0)
            DK = c(1)*x(:,dims1)*x2(:,dims2)';
        end
    else
        if sum(selVars==dims1)==0
            DK=zeros(size(x,1),size(x2,1));
        else
            DK=c(1).*4.*x(:,dims1)*x2(:,dims2)';
        end
        if isequal(gpcf.interactions,'on') && sum(selVars==dims1)>0
            for xi2=selVars
                if xi2~=dims1
                    DK = DK + c(1)*x(:,xi2)*x2(:,xi2)';
                end
            end
        end
    end
    ii1=ii1+1;
    DKff{ii1}=DK;
else
    % vector of coeffSigma2s
    for i1=1:length(selVars)
        if dims1~=dims2 || dims2~=selVars(i1)
            DK=zeros(size(x,1),size(x2,1));
        else
            DK=c(i1).*4.*x(:,dims1)*x2(:,dims2)';
        end
        ii1=ii1+1;
        DKff{ii1}=DK;
    end
    if isequal(gpcf.interactions,'on')
        for xi1=1:length(selVars)
            for xi2=xi1+1:length(selVars)
                i1=i1+1;
                %if k==xi1 && j==xi2
                if dims1==xi1 && dims2==xi2
                    DK = c(i1)*x(:,xi2)*x2(:,xi1)';
                %elseif k==xi2 && j==xi1
                elseif dims1==xi2 && dims2==xi1
                    DK = c(i1)*x(:,xi1)*x2(:,xi2)';
                %elseif k==j && k==xi1
                elseif dims1==dims2 && dims1==xi1
                    DK = c(i1)*x(:,xi2)*x2(:,xi2)';
                %elseif k==j && k==xi2
                elseif dims1==dims2 && dims1==xi2
                    DK = c(i1)*x(:,xi1)*x2(:,xi1)';
                else
                    DK=zeros(size(x,1),size(x2,1));
                end
                ii1=ii1+1;
                DKff{ii1}=DK;
            end
        end
    end
end

end

function DKff = gpcf_squared_ginput(gpcf, x, x2, i1)
%GPCF_SQUARED_GINPUT  Evaluate gradient of covariance function with
%                    respect to x.
%
%  Description
%    DKff = GPCF_SQUARED_GINPUT(GPCF, X, X2) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns DKff, the
%    gradients of covariance matrix Kff = k(X,X2) with respect to X (cell
%    array with matrix elements). If called with only two inputs
%    GPCF_SQUARED_GINPUT(GPCF, X), X2=X. This subfunction is needed when
%    computing gradients with respect to inducing inputs in sparse
%    approximations.
%
%    DKff = GPCF_SQUARED_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to ith covariate in X (matrix).
%    This subfunction is needed when using memory save option
%    in gp_set.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%   GPCF_SQUARED_PAK, GPCF_SQUARED_UNPAK, GPCF_SQUARED_LP, GP_G

if isfield(gpcf, 'selectedVariables')
    error('The selectedVariables option has not yet been implemented for gpcf_squared with derivobs=''on'' ')
    % notice, some parts of the code already take into account the
    % selectedVariables but the code has not been checked
end

[n, m] =size(x);

if nargin==4
    % Use memory save option
    savememory=1;
    if i1==0
        % Return number of covariates
        DKff=m;
        return
    end
else
    savememory=0;
end

if nargin == 2 || isempty(x2)

    xx = x.^2;
    if length(gpcf.coeffSigma2)==1
        s=gpcf.coeffSigma2.*ones(1,(1+m)*m/2);
    else
        s=gpcf.coeffSigma2;
    end
    ii1 = 0;
    if nargin<4
        i1=1:m;
    end
    for j = 1:n
        for i=i1
            DK = zeros(n);
            DK(j,:)=2*s(i)*x(j,i)*xx(:,i)';
            if isequal(gpcf.interactions,'on')
                for xi2=i1
                    if xi2~=i
                        DK(j,:) = DK(j,:) + s(i)*x(j,xi2).*(x(:,i).*x(:,xi2))';
                    end
                end
            end
            DK = DK + DK';
            ii1 = ii1 + 1;
            DKff{ii1} = DK;
        end
    end

elseif nargin == 3 || nargin == 4

    xx2 = x2.^2;
    if length(gpcf.coeffSigma2)==1
        s=gpcf.coeffSigma2.*ones(1,(1+m)*m/2);
    else
        s=gpcf.coeffSigma2;
    end
    ii1 = 0;
    if ~savememory
        i1=1:m;
    end
    for j = 1:n
        for i=i1
            DK = zeros(n, size(x2,1));
            DK(j,:)=2*s(i)*x(j,i)*xx2(:,i)';
            if isequal(gpcf.interactions,'on')
                for xi2=i1
                    if xi2~=i
                        DK(j,:) = DK(j,:) + s(i)*x(j,xi2).*(x2(:,i).*x2(:,xi2))';
                    end
                end
            end
            ii1 = ii1 + 1;
            DKff{ii1} = DK;
        end
    end
end

end

function DKff = gpcf_squared_ginput2(gpcf, x, x2, dims,takeOnlyDiag)
%GPCF_SQUARED_GINPUT2  Evaluate gradient of covariance function with
%                   respect to both input variables x and x2 in
%                   same dimension.
%
%  Description
%    DKff = GPCF_SQUARED_GINPUT2(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of twice derivatived covariance
%    matrix K(df,df) = dk(X1,X2)/dX1dX2 (cell array with matrix
%    elements). Input variable's dimensions are expected to be
%    same.  This subfunction is needed when using derivative
%    observations.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_SQUARED_GINPUT, GPCF_SQUARED_GINPUT2, GPCF_SQUARED_CFDG2

[~,m]=size(x);
ii1=0;
if nargin<4 || isempty(dims)
    dims=1:m;
end
if length(gpcf.coeffSigma2)==1
    c=repmat(gpcf.coeffSigma2,1,(1+m)*m/2);
else
    c=gpcf.coeffSigma2;
end
if isfield(gpcf, 'selectedVariables')
    sv = gpcf.selectedVariables;
else
    sv = 1:m;
end

if nargin==5 && isequal(takeOnlyDiag,'takeOnlyDiag')
    for i1=dims
        if ~any(sv==i1)%isfield(gpcf, 'selectedVariables') && sum(gpcf.selectedVariables==i1)==0
            DK=zeros(size(x,1),1);
        else
            DK=c(i1).*4.*(x(:,i1).*x2(:,i1));
            if isequal(gpcf.interactions,'on')
                i2=length(sv);
                for xi1=1:length(sv)
                    for xi2=xi1+1:length(sv)
                        if any(sv==xi1) && any(sv==xi2)
                            i2=i2+1;
                            if i1==xi1
                                DK = DK + c(i2)*x(:,xi2).*x2(:,xi2);
                            elseif i1==xi2
                                DK = DK + c(i2)*x(:,xi1).*x2(:,xi1);
                            end
                        end
                    end
                end
            end
        end
        ii1=ii1+1;
        DKff{ii1}=DK;
    end
else
    for i1=dims
        if ~any(sv==i1)%isfield(gpcf, 'selectedVariables') && sum(gpcf.selectedVariables==i1)==0
            DK=zeros(size(x,1),size(x2,1));
        else
            DK=c(i1).*4.*(x(:,i1)*x2(:,i1)');
            if isequal(gpcf.interactions,'on')
                i2=length(sv);
                for xi1=1:length(sv)
                    for xi2=xi1+1:length(sv)
                        if any(sv==xi1) && any(sv==xi2)
                            i2=i2+1;
                            if i1==xi1
                                DK = DK + c(i2)*x(:,xi2)*x2(:,xi2)';
                            elseif i1==xi2
                                DK = DK + c(i2)*x(:,xi1)*x2(:,xi1)';
                            end
                        end
                    end
                end
            end
        end
        ii1=ii1+1;
        DKff{ii1}=DK;
    end
end
end

function DKff = gpcf_squared_ginput3(gpcf, x, x2, dims1, dims2)
%GPCF_SQUARED_GINPUT3  Evaluate gradient of covariance function with
%                   respect to both input variables x and x2 (in
%                   different dimensions).
%
%  Description
%    DKff = GPCF_SQUARED_GINPUT3(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of twice derivatived covariance
%    matrix K(df,df) = dk(X1,X2)/dX1dX2 (cell array with matrix
%    elements). The derivative is calculated in multidimensional
%    problem between input's observation dimensions which are not
%    same. This subfunction is needed when using derivative
%    observations.
%
%    DKff is a cell array with the following elements:
%      DKff{1} = dk(X1,X2)/dX1_1dX2_2
%      DKff{2} = dk(X1,X2)/dX1_1dX2_3
%       ...
%      DKff{m-1} = dk(X1,X2)/dX1_1dX2_m
%      DKff{m} = dk(X1,X2)/dX1_2dX2_3
%       ...
%      DKff{m} = dk(X1,X2)/dX1_(m-1)dX2_m
%    where _m denotes the input dimension with respect to which the
%    gradient is calculated.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_SQUARED_GINPUT, GPCF_SQUARED_GINPUT2, GPCF_SQUARED_CFDG2

[~,m]=size(x);
if nargin < 3
    error('Needs at least 3 input arguments')
end
if nargin<4 || isempty(dims1)
    dims1=1:m;
end
if nargin<5 || isempty(dims2)
    dims2=1:m;
end
if isfield(gpcf, 'selectedVariables')
    sv = gpcf.selectedVariables;
else
    sv = 1:m;
end

ii1=0;
if length(gpcf.coeffSigma2)==1
  c=repmat(gpcf.coeffSigma2,1,(1+m)*m/2);
else
  c=gpcf.coeffSigma2;
end
DK=zeros(size(x,1),size(x2,1));
i2=m;
for i=dims1
    for j=dims2
        ii1=ii1+1;
        if isequal(gpcf.interactions,'on') && any(sv==i) && any(sv==j)
            i2 = min(i,j)*length(sv) - 0.5*min(i,j)*(min(i,j)-1) + max(i,j)-min(i,j);
            DKff{ii1} = c(i2)*x(:,j)*x2(:,i)';
        else
            DKff{ii1}=DK;
        end
    end
end

% i2=m;
% for i=1:m-1
%     for j=i+1:m
%         i2=i2+1;
%         if any(dims1==i) && any(dims2==j)
%             ii1=ii1+1;
%             if isequal(gpcf.interactions,'on') && any(sv==i) && any(sv==j)
% %             if isequal(gpcf.interactions,'on') &&...
% %                     ~(isfield(gpcf, 'selectedVariables') && (sum(gpcf.selectedVariables==i)==0 || sum(gpcf.selectedVariables==j)==0))
%                 DKff{ii1} = c(i2)*x(:,j)*x2(:,i)';
%             else
%                 DKff{ii1}=DK;
%             end
%         end
%     end
% end

end

function DKff = gpcf_squared_ginput4(gpcf, x, x2, dims)
%GPCF_SQUARED_GINPUT  Evaluate gradient of covariance function with respect
%                    to x. Simplified and faster version of squared_ginput,
%                    returns full matrices.
%
%  Description
%    DKff = GPCF_SQUARED_GINPUT4(GPCF, X, X2) takes a covariance function
%    structure GPCF, matrices X and X2 of input vectors and returns DKff,
%    the gradients of covariance matrix Kff = k(X,X2) with respect to X
%    (whole matrix); that is d k(X,X2)/dX. If called with only two inputs
%    GPCF_SQUARED_GINPUT4(GPCF, X), X2=X.
%
%    This subfunction is needed when using derivative observations.
%
%    Note! When coding the derivatives of the covariance function, remember
%    to double check them. See gp_cov for lines of code to check the
%    matrices
%
%  See also
%    GPCF_SQUARED_PAK, GPCF_SQUARED_UNPAK, GPCF_SQUARED_LP, GP_G

% if isfield(gpcf, 'selectedVariables')
%     error('The selectedVariables option has not yet been implemented for gpcf_squared with derivobs=''on'' ')
%     % notice, some parts of the code already take into account the
%     % selectedVariables but the code has not been checked
% end
%
% if isfield(gpcf, 'selectedVariables')
%     m = length(gpcf.selectedVariables);
%     sv = gpcf.selectedVariables;
% else
%     [~,m]=size(x);
%     sv = 1:m;
% end
% if nargin==2
%     x2=x;
% end
% if nargin<4
%     dims=1:size(x,2);
% end
% if length(gpcf.coeffSigma2)==1
%     c=repmat(gpcf.coeffSigma2,1,(1+m)*m/2);
% else
%     % If coeffSigma is vector, we trust it is of rigth length
%     c=gpcf.coeffSigma2;
% end
% h2=x2.^2;
% ii1=0;
% for i=dims
%     if ~any(sv==i) %isfield(gpcf, 'selectedVariables') && sum(gpcf.selectedVariables==i)==0
%         DK=zeros(size(x,1),size(x2,1));
%     else
%         DK=c(i).*2.*x(:,i)*h2(:,i)';
%         if isequal(gpcf.interactions,'on')
%             i2=m;
%             for xi1=sv
%                 sv2 = sv(sv>xi1);
%                 for xi2=sv2 %xi1+1:m
%                     i2 = i2+1;
%                     if i==xi1
%                         DK = DK + c(i2)*x(:,xi2)*(x2(:,i).*x2(:,xi2))';
%                     elseif i==xi2
%                         DK = DK + c(i2)*x(:,xi1)*(x2(:,i).*x2(:,xi1))';
%                     end
%                 end
%             end
%         end
%     end
%     ii1=ii1+1;
%     DKff{ii1}=DK;
% end

[n,m]=size(x);
if nargin<4
    dims=1:m;
end
if isfield(gpcf, 'selectedVariables')
    sv = gpcf.selectedVariables;
else
    sv = 1:m;
end
ii1=0;
if nargin==2
    x2=x;
end
h2=x2.^2;
if length(gpcf.coeffSigma2)==1
    c=repmat(gpcf.coeffSigma2,1,(1+m)*m/2);
else
    c=gpcf.coeffSigma2;
end
for i=dims
    if ~any(sv==i) % isfield(gpcf, 'selectedVariables') && sum(gpcf.selectedVariables==i)==0
        DK=zeros(size(x,1),size(x2,1));
    else
        DK=c(i).*2.*x(:,i)*h2(:,i)';
        if isequal(gpcf.interactions,'on')
            i2=length(sv);
            for xi1=1:length(sv)
                for xi2=xi1+1:length(sv)
                    if any(sv==xi1) && any(sv==xi2)
                        i2=i2+1;
%                        i2 = min(i,j)*length(sv) - 0.5*min(i,j)*(min(i,j)-1) + max(i,j)-min(i,j)
                        if i==xi1
                            DK = DK + c(i2)*x(:,xi2)*(x2(:,i).*x2(:,xi2))';
                        elseif i==xi2
                            DK = DK + c(i2)*x(:,xi1)*(x2(:,i).*x2(:,xi1))';
                        end
                    end
                end
            end
        end
    end
    ii1=ii1+1;
    DKff{ii1}=DK;
end
end

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

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);
  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end
  if isfield(gpcf, 'selectedVariables')
      x1 = x1(:,gpcf.selectedVariables);
      x2 = x2(:,gpcf.selectedVariables);
  end
  h1 = x1.^2;
  h2 = x2.^2;
  if isequal(gpcf.interactions,'on')
      m=size(x1,2);
      for xi1=1:m
          for xi2=xi1+1:m
              h1 = [h1 x1(:,xi1).*x1(:,xi2)];
              h2 = [h2 x2(:,xi1).*x2(:,xi2)];
          end
      end
  end
  C = h1*diag(gpcf.coeffSigma2)*(h2');
  C(abs(C)<=eps) = 0;

end

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

  if isfield(gpcf, 'selectedVariables')
      x = x(:,gpcf.selectedVariables);
  end
  h = x.^2;
  if isequal(gpcf.interactions,'on')
      m=size(x,2);
      for xi1=1:m
          for xi2=xi1+1:m
              h = [h x(:,xi1).*x(:,xi2)];
          end
      end
  end
  C = h*diag(gpcf.coeffSigma2)*(h');
  C(abs(C)<=eps) = 0;
  C = (C+C')./2;

end

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

  if isfield(gpcf, 'selectedVariables')
      x = x(:,gpcf.selectedVariables);
  end
  h = x.^2;
  if isequal(gpcf.interactions,'on')
      m=size(x,2);
      for xi1=1:m
          for xi2=xi1+1:m
              h = [h x(:,xi1).*x(:,xi2)];
          end
      end
  end
  if length(gpcf.coeffSigma2) == 1
      C=gpcf.coeffSigma2.*sum(h.^2,2);
  else
      C=sum(repmat(gpcf.coeffSigma2, size(x,1), 1).*h.^2,2);
  end
  C(abs(C)<eps)=0;

end

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

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

    % Initialize parameters
    reccf.coeffSigma2= [];

    % Set the function handles
    reccf.fh.pak = @gpcf_squared_pak;
    reccf.fh.unpak = @gpcf_squared_unpak;
    reccf.fh.lp = @gpcf_squared_lp;
    reccf.fh.lpg = @gpcf_squared_lpg;
    reccf.fh.cfg = @gpcf_squared_cfg;
    reccf.fh.cfdg = @gpcf_squared_cfdg;
    reccf.fh.cfdg2 = @gpcf_squared_cfdg2;
    reccf.fh.ginput = @gpcf_squared_ginput;
    reccf.fh.ginput2 = @gpcf_squared_ginput2;
    reccf.fh.ginput3 = @gpcf_squared_ginput3;
    reccf.fh.ginput4 = @gpcf_squared_ginput4;
    reccf.fh.cov = @gpcf_squared_cov;
    reccf.fh.trcov  = @gpcf_squared_trcov;
    reccf.fh.trvar  = @gpcf_squared_trvar;
    reccf.fh.recappend = @gpcf_squared_recappend;
    reccf.p=[];
    reccf.p.coeffSigma2=[];
    if ~isempty(ri.p.coeffSigma2)
      reccf.p.coeffSigma2 = ri.p.coeffSigma2;
    end

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

    reccf.interactions = gpcf.interactions;

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

    if isfield(gpcf, 'selectedVariables')
      reccf.selectedVariables = gpcf.selectedVariables;
    end
  end
end

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

  % Check arguments
  if nargin < 2 || isempty(x), x = 0; end

  % Scaling
  x0 = min(x);

  % Define the model
  F      = [0 1; 0 0];
  L      = [0; 1];
  Qc     = 0;
  H      = [1 0];
  Pinf   = [x0^2 x0; x0 1]*gpcf.coeffSigma2;
  dF     = zeros(2,2,1);
  dQc    = zeros(1,1,1);
  dPinf  = [x0^2 x0; x0 1];
  params = {};

  % Set params
  params.stationary = false;

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

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

end

author: gpstuff-dev

repo: gpstuff


text:
function [xs, ys, values]=centers2edges(xs, ys, values)
    % CENTERS2EDGES converts matrices of measurement centers to cells with measurements at center
    %
    % edges=CENTERS2EDGES(centers) returns the edges of each bin, for a vector CENTERS containing
    % the center of each bin.  for centers of length N, edges will be of length N+1. edges for the
    % first and last bin, are symmetrical about each bin's center.
    %
    % [xs, ys,values]=CENTERS2EDGES(xs, ys, values) 2-dimensional version. xs and ys are
    % matrices of the same size
    %
    % Pcolor typically uses the points as edges, and ignores the last values.
    %
    % warning, shading interp might look pretty, but it will shift items and is inaccurate.
    % instead, use 2interp
    %
    % ex.
    %  >> centers2edges([0 2 4 8])
    %  ans =
    %      -1     1     3     6    10
    %
    % see also EDGES2CENTERS, HISTCOUNTS, HISTOGRAM

    if nargout ~= nargin
        error('Each input should have a corresponding output');
    end

    if nargin==1 && isvector(xs)
        if isempty(xs)
            % nothing to do. return the empty value
        else
            if numel(xs) <2
                error('need at least 2 centers to determine bin size')
            end
            half_dx = diff(xs) ./ 2;
            if any(half_dx <= 0)
                error('values should be in ascending order');
            end
            half_dx = [-half_dx(1); half_dx(:); half_dx(end)];
            xs(2:end+1)=xs;
            xs = xs + reshape(half_dx,size(xs));
        end
    else
        if ~isequal(size(xs), size(ys))
            error('X and Y matrices should be same size');
        end
        if min(size(xs) ==1)
            error('X and Y must be matrices, not vectors');
        end
        % expand grid in each direction so that xs and ys are in the centers
        half_dx = diff(xs,[],2) ./ 2;
        if any(half_dx <= 0)
            error('X vector should be in ascending order')
        end
        half_dx = [-half_dx(:,1) , half_dx , half_dx(:,end)];
        xs=[xs(:,1), xs] + half_dx;
        xs(end+1,:)=xs(end,:);

        half_dy = diff(ys,[],1) ./ 2;

        if any(half_dx <= 0)
            error('Y vector should be in ascending order');
        end
        half_dy = [-half_dy(1,:) ; half_dy ; half_dy(end,:)];
        ys=[ys(1,:); ys] + half_dy;
        ys(:,end+1)=ys(:,end);

        if exist('values','var')
            assert(isequal(size(xs)-1,size(values)),'size of VALUES doesn''t match X and Y matrices');
            values(end+1,:)=nan;
            values(:,end+1)=nan;
        elseif nargout >=3
            error('Values will only be returned if provided');
        end
    end
end

author: CelsoReyes


text:
function CPD = learn_params(CPD, fam, data, ns, cnodes)
%function CPD = learn_params(CPD, local_data)
% LEARN_PARAMS Compute the ML/MAP estimate of the params of a tabular CPD given complete data
% CPD = learn_params(CPD, local_data)
%
% local_data(i,m) is the value of i'th family member in case m (can be cell array).

local_data = data(fam, :);
if iscell(local_data)
  local_data = cell2num(local_data);
end
counts = compute_counts(local_data, CPD.sizes);
switch CPD.prior_type
 case 'none', CPD.CPT = mk_stochastic(counts);
 case 'dirichlet', CPD.CPT = mk_stochastic(counts + CPD.dirichlet);
 otherwise, error(['unrecognized prior ' CPD.prior_type])
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function bessel_j0_int_values_test ( )

%*****************************************************************************80
%
%% BESSEL_J0_INT_VALUES_TEST demonstrates the use of BESSEL_JO_INT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_J0_INT_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_J0_INT_VALUES stores values of \n' );
  fprintf ( 1, '  the integral of the Bessel function J0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_j0_int_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function prtPlotUtilDataSetExploreGuiSimple(theObject,parent)
% Internal function,
% xxx Need Help xxx - see prtDataSetClass.explore

if nargin == 1
    windowSize = [754 600];
    pos = prtPlotUtilCenterFigure(windowSize);

    % Create the figure an UIControls
    figH = figure('Name','PRT Data Set Explorer','Menu','none','toolbar','figure','units','pixels','position',pos,'DockControls','off');
    try
        set(figH,'Number','Off');
    catch
        try
            set(figH,'NumberTitle','off');
        end
    end

    % Trim the toolbar down to just the zooming controls
    Toolbar.handle = findall(figH,'Type','uitoolbar');
    Toolbar.Children = findall(figH,'Parent',Toolbar.handle,'HandleVisibility','off');

    % Delete a bunch of things we dont need
    delete(findobj(Toolbar.Children,'TooltipString','New Figure',...
        '-or','TooltipString','Open File','-or','TooltipString','Save Figure',...
        '-or','TooltipString','Print Figure','-or','TooltipString','Edit Plot',...
        '-or','TooltipString','Data Cursor','-or','TooltipString','Brush/Select Data',...
        '-or','TooltipString','Link Plot','-or','TooltipString','Insert Colorbar',...
        '-or','TooltipString','Insert Legend','-or','TooltipString','Show Plot Tools and Dock Figure',...
        '-or','TooltipString','Hide Plot Tools'))
    bgc = get(figH,'Color');
else
    bgc = get(gcf,'Color');
    figH = parent;
end

popUpStrs = theObject.getFeatureNames;
popUpStrs = popUpStrs(:);

popX = uicontrol(figH,'Style','popup','units','normalized','FontUnits','Normalized','FontSize',0.6,'position',[0.15 0.01 0.19 0.04],'string',popUpStrs,'callback',{@plotSelectPopupCallback 1});
popXHead = uicontrol(figH,'Style','text','units','normalized','FontUnits','Normalized','FontSize',0.75,'position',[0.05 0.01 0.09 0.04],'string','X-Axis:','BackgroundColor',bgc,'HorizontalAlignment','Right'); %#ok

popY = uicontrol(figH,'Style','popup','units','normalized','FontUnits','Normalized','FontSize',0.6,'position',[0.45 0.01 0.19 0.04],'string',popUpStrs,'callback',{@plotSelectPopupCallback 2});
popYHead = uicontrol(figH,'Style','text','units','normalized','FontUnits','Normalized','FontSize',0.75,'position',[0.35 0.01 0.09 0.04],'string','Y-Axis:','BackgroundColor',bgc,'HorizontalAlignment','Right'); %#ok

popZ = uicontrol(figH,'Style','popup','units','normalized','FontUnits','Normalized','FontSize',0.6,'position',[0.75 0.01 0.19 0.04],'string',[{'None'}; popUpStrs],'callback',{@plotSelectPopupCallback 3});
popZHead = uicontrol(figH,'Style','text','units','normalized','FontUnits','Normalized','FontSize',0.75,'position',[0.65 0.01 0.09 0.04],'string','Z-Axis:','BackgroundColor',bgc,'HorizontalAlignment','Right'); %#ok

axisH = axes('Units','Normalized','outerPosition',[0.05 0.07 0.9 0.9]);

% Setup the PopOut Option
hcmenu = uicontextmenu;
hcmenuPopoutItem = uimenu(hcmenu, 'Label', 'Popout', 'Callback', @explorerPopOut); %#ok
set(axisH,'UIContextMenu',hcmenu);

if theObject.nFeatures > 1
    plotDims = [1 2 0];

    set(popX,'value',1); % Becase we have dont have a none;
    set(popY,'value',2); % Becase we have dont have a none;
    set(popZ,'value',1); % Becase we have a none;
else
    plotDims = [1 1 0];

    set(popX,'value',1); % Becase we have dont hvae a none;
    set(popY,'value',1); % Becase we have a none;
    set(popZ,'value',1); % Becase we have a none;
end
updatePlot;

    function plotSelectPopupCallback(myHandle, eventData, varargin) %#ok
        cVal = get(myHandle,'value');
        axisInd = varargin{1};
        if axisInd == 3
            % Z-axis we have a None option
            cVal = cVal - 1;
        end
        plotDims(axisInd) = cVal;
        updatePlot;
    end

    function updatePlot
        actualPlotDims = plotDims(plotDims>=1);
        axes(axisH); %#ok
        h = plot(theObject,actualPlotDims);
        set(h,'HitTest','off');
        set(axisH,'ButtonDownFcn',@axisOnClick);
    end
    function explorerPopOut(myHandle,eventData) %#ok
        figure;
        actualPlotDims = plotDims(plotDims>=1);
        plot(theObject,actualPlotDims);
    end
    function axisOnClick(myHandle,eventData)
        actualPlotDims = plotDims(plotDims>=1);
        data = theObject.getFeatures(actualPlotDims);

        [rP,rD] = rotateDataAndClick(data);

        dist = prtDistanceEuclidean(rP,rD);
        [~,i] = min(dist);
        obsName = theObject.getObservationNames(i);
        title(sprintf('Observation Closest To Last Click: %s',obsName{1}));

        debug = false;
        if debug
            hold on;
            d = theObject.getObservations(i);
            switch length(actualPlotDims)
                case 2
                    plot(d(1),d(2),'kx');
                case 3
                    plot3(d(1),d(2),d(3),'kx');
            end
            hold off;
        end
    end

    function [rotatedData,rotatedClick] = rotateDataAndClick(data)
        %[rotatedData,rotatedClick] = rotateDataAndClick(data)
        % Used internally; from Click3dPoint from matlab central; Copyright
        % follows

        %         Copyright (c) 2009, Babak Taati
        % All rights reserved.
        %
        % Redistribution and use in source and binary forms, with or without
        % modification, are permitted provided that the following conditions are
        % met:
        %
        %     * Redistributions of source code must retain the above copyright
        %       notice, this list of conditions and the following disclaimer.
        %     * 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 OWNER 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.
        point = get(gca, 'CurrentPoint'); % mouse click position
        camPos = get(gca, 'CameraPosition'); % camera position
        camTgt = get(gca, 'CameraTarget'); % where the camera is pointing to

        camDir = camPos - camTgt; % camera direction
        camUpVect = get(gca, 'CameraUpVector'); % camera 'up' vector

        % build an orthonormal frame based on the viewing direction and the
        % up vector (the "view frame")
        zAxis = camDir/norm(camDir);
        upAxis = camUpVect/norm(camUpVect);
        xAxis = cross(upAxis, zAxis);
        yAxis = cross(zAxis, xAxis);

        rot = [xAxis; yAxis; zAxis]; % view rotation

        if size(data,2) < 3
            data = cat(2,data,zeros(size(data,1),1));
        end

        % the point cloud represented in the view frame
        rotatedData = (rot * data')';
        rotatedData = rotatedData(:,1:2);
        % the clicked point represented in the view frame
        rotatedClick = rot * point' ;
        rotatedClick = rotatedClick(1:2,1)';
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function res = spm_eeg_artefact_zscore(S)
% Plugin for spm_eeg_artefact doing z-score thresholding
% S            - input structure
% fields of S:
%    S.D       - M/EEG object
%    S.chanind - vector of indices of channels that this plugin will look at
%
%    Additional parameters can be defined specific for each plugin.
%
% Output:
% res -
%    If no input is provided the plugin returns a cfg branch for itself.
%
%    If input is provided the plugin returns a matrix of size D.nchannels x D.ntrials
%    with zeros for clean channel/trials and ones for artefacts.
%__________________________________________________________________________
% Copyright (C) 2013-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_artefact_zscore.m 7132 2017-07-10 16:22:58Z guillaume $

%-This part if for creating a config branch that plugs into spm_cfg_eeg_artefact
% Any parameters can be specified and they are then passed to the plugin
% when it's called.
%--------------------------------------------------------------------------
if nargin == 0
    threshold         = cfg_entry;
    threshold.tag     = 'threshold';
    threshold.name    = 'Threshold';
    threshold.strtype = 'r';
    threshold.num     = [1 1];
    threshold.val     = {3};
    threshold.help    = {'Threshold value (in stdev)'};

    excwin         = cfg_entry;
    excwin.tag     = 'excwin';
    excwin.name    = 'Excision window';
    excwin.strtype = 'r';
    excwin.num     = [1 1];
    excwin.val     = {100};
    excwin.help    = {'Window (in ms) to mark as bad around each event (for mark mode only), 0 - do not mark data as bad.'};

    zscore      = cfg_branch;
    zscore.tag  = 'zscore';
    zscore.name = 'Threshold z-scored data';
    zscore.val  = {threshold, excwin};
    zscore.help = {''};

    res = zscore;

    return
end

SVNrev = '$Rev: 7132 $';

%-Startup
%--------------------------------------------------------------------------
spm('sFnBanner', mfilename, SVNrev);
spm('FigName','M/EEG zscore thresholding');

D = spm_eeg_load(S.D);

chanind = S.chanind;
threshold = S.threshold;

if isequal(S.mode, 'reject') && isequal(D.type, 'continuous')
    error('Rejection mode not for continuous data.');
end

res = zeros(D.nchannels, D.ntrials);

%-Artefact detection
%--------------------------------------------------------------------------

spm_progress_bar('Init', length(chanind), 'Channels checked');
if length(chanind) > 100, Ibar = floor(linspace(1, length(chanind),100));
else Ibar = [1:length(chanind)]; end

for j = 1:length(chanind)
    dat = squeeze(D(chanind(j), :, :));
    if size(dat, 1) == 1
        dat = dat';
    end

    dat = detrend(dat);

    zsdat = dat./repmat(std(dat), size(dat, 1), 1);

    if isequal(S.mode, 'reject')
        res(chanind(j), :) = any(abs(zsdat)>threshold);
    elseif isequal(S.mode, 'mark')
        bad  = abs(zsdat)>threshold;
        if ~any(bad(:))
            if isequal(D.type, 'continuous')
                if any(Ibar == j), spm_progress_bar('Set', j); end
            end
            continue;
        end

        if S.excwin>0
            excwin = ones(round(1e-3*S.excwin*D.fsample), 1);
            bad = ~~conv2(excwin, 1, double(bad), 'same');
        end

        for i = 1:D.ntrials
            res = [];

            onsets  = find(bad(:, i));
            offsets = find(~bad(:, i));
            onsets(find(diff(onsets)<2)+1) = [];

            if bad(end)
                offsets(end+1) = length(bad)+1;
            end

            for k = 1:length(onsets)
                res(end+1).type   = 'artefact_zscore';
                res(end).value    = char(D.chanlabels(chanind(j)));
                res(end).time     = D.time(onsets(k)+1) - D.time(1) + D.trialonset(i);
                res(end).duration = (min(offsets(offsets>onsets(k)))-onsets(k))./D.fsample;
            end

            if ~isempty(res)
                ev = D.events(i);
                if iscell(ev)
                    ev = ev{1};
                end

                if ~S.append
                    sel1 = strmatch('artefact_zscore', {ev.type});
                    sel2 = strmatch(char(D.chanlabels(chanind(j))), {ev(sel1).value});
                    ev(sel1(sel2)) = [];
                end

                D = events(D, i, spm_cat_struct(ev, res));
            end
        end
    end

    if any(Ibar == j), spm_progress_bar('Set', j); end
end

spm_progress_bar('Clear');

if isequal(S.mode, 'mark')
    res = D;
end

spm('FigName','M/EEG zscore thresholding: done');

author: spm


text:
function g = modelLogLikeGradients(model)

% MODELLOGLIKEGRADIENTS Compute a model's gradients wrt log likelihood.
% FORMAT
% DESC is a wrapper function to compute the gradients of the log
% likelihood of a given model.
% ARG model : the model for which likelihoods are computed.
% RETURN g : teh gradients of the likelihood with respect to the
% parameters.
%
% SEEALSO : modelCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2005

% MLTOOLS

fhandle = str2func([model.type 'LogLikeGradients']);
g = fhandle(model);

if isfield(model, 'paramGroups')
  g = g*model.paramGroups;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function inc = crtran ( a, c, c_size, m, k, n, critvl, i1, c1, c2, iswitch )

%*****************************************************************************80
%
%% CRTRAN determines the effect of moving an object to another class.
%
%  Discussion:
%
%    This computation is very inefficient.  It is only set up so that we
%    can compare algorithm ASA 113 to the K-means algorithms ASA 058 and
%    ASA 136.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Colin Banfield, LC Bassill,
%    Algorithm AS 113:
%    A transfer for non-hierarchichal classification,
%    Applied Statistics,
%    Volume 26, Number 2, 1977, pages 206-210.
%
%  Parameters:
%
%    Input, real A(M,N), the data values.  There are M objects,
%    each having spatial dimension N.
%
%    Input, integer C(M), the classification of each object.
%
%    Input, integer C_SIZE(K), the number of objects in each class.
%
%    Input, integer M, the number of objects.
%
%    Input, integer K, the number of classes.
%
%    Input, integer N, the number of spatial dimensions, or variates,
%    of the objects.
%
%    Input, real CRITVL, the current value of the criterion.
%
%    Input, integer I1, the object to be transferred.
%
%    Input, integer C1, C2, the current class of object I1, and the
%    class to which it may be transferred.
%
%    Input, integer ISWITCH:
%    1, indicates that I1 should be temporarily transferred, the
%       change in CRITVL should be computed, and then I1 restored.
%    2, indicates that I1 will be permanently transferred.
%
%    Output, real INC, the change to CRITVL that would occur if I1 were
%    transferred from class C1 to C2.  This is only computed for ISWITCH = 1.
%
  inc = 0.0;

  if ( iswitch == 2 )
    return
  end
%
%  Move object I from class C1 to class C2.
%
  c(i1) = c2;
  c_size(c1) = c_size(c1) - 1;
  c_size(c2) = c_size(c2) + 1;
%
%  Define the critical value as the sum of the squares of the distances
%  of the points to their cluster center.
%
  for i = 1 : k
    c_size(i) = 0;
    for j = 1 : n
      c_center(i,j) = 0.0;
    end
  end

  for i = 1 : m
    ci = c(i);
    c_size(ci) = c_size(ci) + 1;
    for j = 1 : n
      c_center(ci,j) = c_center(ci,j) + a(i,j);
    end
  end

  for i = 1 : k
    for j = 1 : n
      c_center(i,j) = c_center(i,j) / c_size(i);
    end
  end

  for i = 1 : k
    wss(i) = 0.0;
  end

  for i = 1 : m
    ci = c(i);
    for j = 1 : n
      wss(ci) = wss(ci) + ( a(i,j) - c_center(ci,j) )^2;
    end
  end

  critvl_new = 0.0;
  for i = 1 : k
    critvl_new = critvl_new + wss(i);
  end

  inc = critvl_new - critvl;
%
%  Move object I1 from class C2 to class C1.
%
  c(i1) = c1;
  c_size(c1) = c_size(c1) + 1;
  c_size(c2) = c_size(c2) - 1;

  return
end


text:
function outStrMat = mat2strDisplay(inNumMat,inNumFmt)
[nrow,ncol] = size(inNumMat);
dispFmt = ['%0.',num2str(inNumFmt),'f'];
outStrMat = '[';
for i=1:nrow
    for j=1:ncol
        if j~=ncol
            outStrMat = [outStrMat,sprintf(dispFmt,inNumMat(i,j)),','];
        else
            if nrow~=1
                outStrMat = [outStrMat,sprintf(dispFmt,inNumMat(i,j)),';'];
            else
                outStrMat = [outStrMat,sprintf(dispFmt,inNumMat(i,j))];
            end
        end
    end
end
outStrMat = [outStrMat,']'];
end

repo: paper-simulation


text:
 function ob = Gtomo2_wtfmex(file, chat, sp_group_type)
%function ob = Gtomo2_wtfmex(file, chat, sp_group_type)
%
% THIS IS OBSOLETE; USE Gtomo2_wtmex INSTEAD!
%
% Construct Gtomo2_wtfmex object, which does 2d forward and backprojection
% using the wtfmex cabilities for any system geometry.
% See Gtomo2_wtfmex_test.m for example usage.
% This object avoids the memory pigging property of matlab's sparse matrices
% using wtfmex's static system matrix capability.
% Basically, you create a wtfmex system matrix by calling:
%	G = Gtomo2_wtfmex(file)
% and then you can use it thereafter by typing commands like
%	y = G * x;
% which will auto-magically call wtfmex to do the multiplication
% in the compiled C mex program.
%
% Alternative usage: let "file" be geometry description from arg_pair().
%
% Copyright 2002-2-20, Jeff Fessler, The University of Michigan

if ~isvar('chat') | isempty(chat)
	chat = 0;
end

if ~isvar('sp_group_type') | isempty(sp_group_type)
	sp_group_type = 'row';
end

%
% create default object, as required by Mathworks
%
ob = Gtomo2(chat);
ob.file = '';
ob.arg = '';

if ~nargin % required by Mathworks
	ob = class(ob, 'Gtomo2_wtfmex');
return
end

warning 'Gtomo2_wtfmex is obsolete.  use Gtomo2_wtmex instead'

if exist(file, 'file')
	ob.file = file;

	[ob.nx, ob.ny, ob.nb, ob.na] = wtfmex('read', ob.file);
else
	ob.arg = file;
	wtfmex('gensys', ob.arg', sp_group_type);
	ob.nx = str2num(arg_get(ob.arg, 'nx'));
	ob.ny = str2num(arg_get(ob.arg, 'ny'));
	ob.nb = str2num(arg_get(ob.arg, 'nb'));
	ob.na = str2num(arg_get(ob.arg, 'na'));
end

ob.dims = [ob.nb * ob.na, ob.nx * ob.ny];

ob = class(ob, 'Gtomo2_wtfmex');

tmp = sum(ob) > 0;
ob.mask = reshape(tmp, ob.nx, ob.ny);

if ob.chat
	wtfmex('print'); % display header info
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [nu] = tapas_ti_init_u(u, ptheta, pars)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

nu = cell(1, numel(pars.T));
nu(:) = {u};

end

author: translationalneuromodeling

repo: tapas


text:
function [yr] = s2yr(s)
% Convert time from seconds to Julian years (365.25 days).
% Chad Greene 2012
yr = s*3.168808781403e-8;

author: Sable

repo: mcbench-benchmarks


text:
classdef(Abstract) AbstractLiftCoefficientModel < matlab.mixin.SetGet & matlab.mixin.Heterogeneous & matlab.mixin.Copyable
    %AbstractDragCoefficientModel Summary of this class goes here
    %   Detailed explanation goes here

    properties(Abstract, Constant)
        enum LiftCoefficientModelEnum
    end

    methods
        [ClS, liftUnitVectInertial] = getLiftCoeffAndDir(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEF, attState)

        useTf = openEditDialog(obj)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function m_demo(index)
% M_DEMO  Demonstration program showing various maps in M_Map package
%         Dig into this to look for examples of things you want to do.
%
%         M_DEMO runs all the demos.
%         M_DEMO(NUM) runs examples NUM (1<=NUM<=10), pausing between examples.
%
%         Some demos may require you to install GSHHS or TerrainBase datafiles
%         (see documentation)

% Rich Pawlowicz (rich@ocgy.ubc.ca) 7/May/1997
% (thanks to Art Newhall for putting these examples into an m-file,
% and the Chuck Denham for enhancing the interface)
%
% 27/July/98 - more examples.
% 17/Aug/98     "
% 15/Nov/98  - another example, better interface.
% 23/Dec/98  - another example.

%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

N_EXAMPLES=11;

if nargin==0
 index=1:N_EXAMPLES;
end

for i=index

clf;
switch i

  case 1

    m_proj('ortho','lat',48','long',-123');
    m_coast('patch','r');
    m_grid('linest','-','xticklabels',[],'yticklabels',[],'ytick',[-80:40:80]);
    xlabel('Orthographic Projection','visible','on');

  case 2

    m_proj('lambert','long',[-160 -40],'lat',[30 80]);
    m_coast('patch',[1 .85 .7]);
    m_elev('contourf',[500:500:6000]);
    m_grid('box','fancy','tickdir','in');
    colormap(flipud(copper));
    xlabel('Conic Projection of North America with elevations','visible','on');

  case 3

    m_proj('stereographic','lat',90,'long',30,'radius',25);
    m_elev('contour',[-3500:1000:-500],'edgecolor','b');
    m_grid('xtick',12,'tickdir','out','ytick',[70 80],'linest','-');
    m_coast('patch',[.7 .7 .7],'edgecolor','r');
    xlabel('Polar Stereographic Projection with bathymetry','visible','on');

  case 4

    subplot(211);
    Slongs=[-100 0;-75 25;-5 45; 25 145;45 100;145 295;100 290];
    Slats= [  8 80;-80  8; 8 80;-80   8; 8  80;-80   0;  0  80];
    for l=1:7
     m_proj('sinusoidal','long',Slongs(l,:),'lat',Slats(l,:));
     m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...
            'ytick',[-80:20:80],'yticklabels',[],'linest','-','color',[.9 .9 .9]);
     m_coast('patch','g');
    end
    xlabel('Interrupted Sinusoidal Projection of World Oceans');
    % In order to see all the maps we must undo the axis limits set by m_grid calls:
    set(gca,'xlimmode','auto','ylimmode','auto');

    subplot(212);
    Slongs=[-100 43;-75 20; 20 145;43 100;145 295;100 295];
    Slats= [  0  90;-90  0;-90   0; 0  90;-90   0;  0  90];
    for l=1:6
     m_proj('mollweide','long',Slongs(l,:),'lat',Slats(l,:));
     m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...
            'ytick',[-80:20:80],'yticklabels',[],'linest','-','color','k');
     m_coast('patch',[.6 .6 .6]);
    end
    xlabel('Interrupted Mollweide Projection of World Oceans');
    set(gca,'xlimmode','auto','ylimmode','auto');

  case 5

    %% Nice looking data
    [lon,lat]=meshgrid([-136:2:-114],[36:2:54]);
    u=sin(lat/6);
    v=sin(lon/6);

    m_proj('oblique','lat',[56 30],'lon',[-132 -120],'aspect',.8);

    subplot(121);
    m_coast('patch',[.9 .9 .9],'edgecolor','none');
    m_grid('tickdir','out','yaxislocation','right',...
	   'xaxislocation','top','xlabeldir','end','ticklen',.02);
    hold on;
    m_quiver(lon,lat,u,v);
    xlabel('Simulated surface winds');

    subplot(122);
    m_coast('patch',[.9 .9 .9],'edgecolor','none');
    m_grid('tickdir','out','yticklabels',[],...
	   'xticklabels',[],'linestyle','none','ticklen',.02);
    hold on;
    [cs,h]=m_contour(lon,lat,sqrt(u.*u+v.*v));
    clabel(cs,h,'fontsize',8);
    xlabel('Simulated something else');

  case 6

    % Plot a circular orbit
    lon=[-180:180];
    lat=atan(tan(60*pi/180)*cos((lon-30)*pi/180))*180/pi;

    m_proj('miller','lat',82);
    m_coast('color',[0 .6 0]);
    m_line(lon,lat,'linewi',3,'color','r');
    m_grid('linest','none','box','fancy','tickdir','out');

  case 7

    m_proj('lambert','lon',[-10 20],'lat',[33 48]);
    m_tbase('contourf');
    m_grid('linestyle','none','tickdir','out','linewidth',3);
    colormap(jet);

  case 8

    m_vec;

  case 9

    % Example showing the default coastline and all of the GSHHS coastlines.

    axes('position',[.35 .6 .37 .37]);
    m_proj('albers equal-area','lat',[40 60],'long',[-90 -50],'rect','on');
    m_coast('patch',[0 1 0]);
    m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right');
    m_text(-69,41,'Standard coastline','color','r','fontweight','bold');

    axes('position',[.09 .5 .37 .37]);
    m_proj('albers equal-area','lat',[40 54],'long',[-80 -55],'rect','on');
    m_gshhs_c('patch',[.2 .8 .2]);
    m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top');
    m_text(-80,52.5,'GSHHS\_C (crude)','color','m','fontweight','bold','fontsize',14);

    axes('position',[.13 .2 .37 .37]);
    m_proj('albers equal-area','lat',[43 48],'long',[-67 -59],'rect','on');
    m_gshhs_l('patch',[.4 .6 .4]);
    m_grid('linest','none','linewidth',2,'tickdir','out');
    m_text(-66.5,43.5,'GSHHS\_L (low)','color','m','fontweight','bold','fontsize',14);

    axes('position',[.35 .05 .37 .37]);
    m_proj('albers equal-area','lat',[45.8 47.2],'long',[-64.5 -62],'rect','on');
    m_gshhs_i('patch',[.5 .6 .5]);
    m_grid('linest','none','linewidth',2,'tickdir','out','yaxisloc','right');
    m_text(-64.4,45.9,'GSHHS\_I (intermediate)','color','m','fontweight','bold','fontsize',14);

    axes('position',[.55 .23 .37 .37]);
    m_proj('albers equal-area','lat',[46.375 46.6],'long',[-64.2 -63.7],'rect','on');
    m_gshhs_h('patch',[.6 .6 .6]);
    m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right');
    m_text(-64.18,46.58,'GSHHS\_H (high)','color','m','fontweight','bold','fontsize',14);

  case 10

    % Example showing a trackline plot

    clf
    m_proj('UTM','long',[-72 -68],'lat',[40 44]);
    m_gshhs_i('color','k');
    m_grid('box','fancy','tickdir','in');

    % fake up a trackline
    lons=[-71:.1:-67];
    lats=60*cos((lons+115)*pi/180);
    dates=datenum(1997,10,23,15,1:41,zeros(1,41));

    m_track(lons,lats,dates,'ticks',0,'times',4,'dates',8,...
           'clip','off','color','r','orient','upright');

  case 11

    % example showing range rings

    clf
    m_proj('hammer','clong',170);
    m_grid('xtick',[],'ytick',[],'linestyle','-');
    m_coast('patch','g');
    m_line(100.5,13.5,'marker','square','color','r');
    m_range_ring(100.5,13.5,[1000:1000:15000],'color','b','linewi',2);
    xlabel('1000km range rings from Bangkok');

  end

 if i<length(index)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   disp('  hit return to continue');
   pause
   disp('        ...drawing');
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 end
end


text:
clear
addpath('../PDM_helpers/');
addpath(genpath('../fitting/'));
addpath('../models/');
addpath(genpath('../face_detection'));
addpath('../CCNF/');

%% loading the patch experts

[clmParams, pdm] = Load_CLM_params_66();

% A CLM-Z model trained on Multi-PIE and BU-4DFE
[patches] = Load_Patch_Experts( '../models/clmz/', 'svr_patches_multi_pie_*.mat', '../models/clmz/', 'svr_depth_patches_*.mat', clmParams);

clmParams.multi_modal_types  = patches(1).multi_modal_types;

%%

images = {'sample_depth_imgs/1.jpg', 'sample_depth_imgs/2.jpg', 'sample_depth_imgs/3.jpg', 'sample_depth_imgs/4.jpg', 'sample_depth_imgs/5.jpg'};
images_depth = {'sample_depth_imgs/1d.png', 'sample_depth_imgs/2d.png', 'sample_depth_imgs/3d.png', 'sample_depth_imgs/4d.png', 'sample_depth_imgs/5d.png'};
verbose = true;

for img=1:numel(images)

    image_orig = imread(images{img});

    image_depth = imread(images_depth{img});

    % Need to convert from the disparity to depth values, and threshold
    image_depth = 10000./(image_depth);
    image_depth(image_depth > 300) = 0;

    % First attempt to use the Matlab one (fastest but not as accurate, if not present use yu et al.)
    [bboxs] = detect_faces(image_orig, {'cascade', 'zhu'});

    if(size(image_orig,3) == 3)
        image = rgb2gray(image_orig);
    end

    %%

    if(verbose)
        f = figure;
        if(max(image(:)) > 1)
            imshow(double(image_orig)/255, 'Border', 'tight');
        else
            imshow(double(image_orig), 'Border', 'tight');
        end
        axis equal;
        hold on;
    end

    for i=1:size(bboxs,2)

        % Convert from the initial detected shape to CLM model parameters
        bbox = bboxs(:,i);

        % Use the initial global and local params for clm fitting in the image
        [shape,~,~,lhood,lmark_lhood,view_used] = Fitting_from_bb(image, image_depth, bbox, pdm, patches, clmParams);

        % shape correction for matlab format
        shape = shape + 1;

        if(verbose)

            % valid points to draw (not to draw self-occluded ones)
            v_points = logical(patches(1).visibilities(view_used,:));

            try
                plot(shape(v_points,1), shape(v_points',2),'.r','MarkerSize',20);
                plot(shape(v_points,1), shape(v_points',2),'.b','MarkerSize',10);
            catch warn

            end
        end

    end
    hold off;

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function b = r8sp_mv ( m, n, nz_num, row, col, a, x )

%*****************************************************************************80
%
%% R8SP_MV multiplies an R8SP matrix by an R8VEC.
%
%  Discussion:
%
%    The R8SP storage format stores the row, column and value of each nonzero
%    entry of a sparse matrix.
%
%    It is possible that a pair of indices (I,J) may occur more than
%    once.  Presumably, in this case, the intent is that the actual value
%    of A(I,J) is the sum of all such entries.  This is not a good thing
%    to do, but I seem to have come across this in MATLAB.
%
%    The R8SP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP
%    ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of
%    the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero elements in
%    the matrix.
%
%    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and
%    column indices of the nonzero elements.
%
%    Input, real A(NZ_NUM), the nonzero elements of the matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(M), the product vector A*X.
%
  b = zeros(m,1);

  for k = 1 : nz_num

    i = row(k);
    j = col(k);
    b(i) = b(i) + a(k) * x(j);

  end

  return
end


text:
 function Xk = nufft_table_adj(st, X, order, flips, om)
%function Xk = nufft_table_adj(st, X, order, flips, om)
%| adjoint of table-based nufft interpolation.
%|
%| in
%|	st		structure from nufft_init
%|	X [M,nc]	DTFT values (usually nc=1)
%|	order	0|1	0th or 1st-order interpolation
%|			default is 0 for backward compatability
%|	flips	0|1	sign flips? (for real table with even N)
%|	om [M,1]	optional (default st.om)
%|
%| out
%|	Xk [*Kd,nc]	DFT coefficients
%|
%| Copyright 2004-3-30, Jeff Fessler and Yingying Zhang, University of Michigan

if nargin < 2, ir_usage, end

if ~isvar('order') || isempty(order)
	order = 0; % default 0th order for backward compatability
end

if ~isvar('flips') || isempty(flips)
	flips = zeros(size(st.Nd)); % default no flips for backward compatability
end

if nargin < 4
	om = st.om;
end

dd = length(st.Kd);

tm = zeros(size(om)); % must be double for mex file!
for id=1:dd
	gam = 2*pi / st.Kd(id);
	tm(:,id) = om(:,id) / gam; % t = omega / gamma
end

if size(X,1) ~= size(om,1)
	fail 'X size problem'
end

nc = size(X,2);

% adjoint of phase shift
class_X = class(X);
if isfield(st, 'phase_shift') && ~isempty(st.phase_shift)
	X = X .* repmat(conj(st.phase_shift), [1 nc]);
end

% convert X to complex double for mex file
if ~isa(X, 'double'), X = double(X); end

arg = {int32(st.Jd), int32(st.Ld), double(tm), int32(st.Kd(1:dd)), ...
	int32(order), int32(flips)};

X = complexify(X);

switch dd
case 1
	Xk = interp1_table_adj_mex(X, st.h{1}, arg{:});

case 2
	Xk = interp2_table_adj_mex(X, st.h{1}, st.h{2}, arg{:});

case 3
	Xk = interp3_table_adj_mex(X, st.h{1}, st.h{2}, st.h{3}, arg{:});

case 4
	Xk = interp4_table_adj_mex(X, st.h{1}, st.h{2}, st.h{3}, st.h{4}, arg{:});

otherwise
	fail '> 4d not done'
end
Xk = cast(Xk, class_X);


text:
%% Example
% An example script, it is used to show how to use the different
% functions of the KUKA iiwa matlab toolbox

% First run the following script in Matlab
% Then start the client on the KUKA iiwa controller

% Mohammad SAFEEA, 24th of April 2018

close all,clear all,clc
warning('off')
ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
else
    %% move to initial position
pinit={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % initial confuguration
relVel=0.15; % relative velocity
movePTPJointSpace( t_Kuka , pinit, relVel); % point to point motion in joint space

    %% Get the joints positions
      jPos  = getJointsPos( t_Kuka );
      fprintf('The joints positions of the robot are: \n');
      disp(jPos);
 % jPos =

 %   [-8.6886e-06]    [0.3491]    [9.9110e-05]    [-1.2218]    [-8.1109e-05]    [1.5708]    [7.7778e-05]

      %% Get position roientation of end effector
      fprintf('Cartesian position/orientation of end-effector:\n');
      pos=getEEFPos( t_Kuka );
      disp(pos);
% pos =

%    [536.8023]    [0.0255]    [563.8379]    [-3.1416]    [-1.1624e-04]    [-3.1415]

      %% Get position of end effector
      fprintf('Cartesian position of end-effecotr:\n')
      cpos=getEEFCartesianPosition( t_Kuka );
      disp(cpos);
%cpos =

%    [536.8023]    [0.0254]    [563.8380]
      %% Get orientation of end effector
      fprintf('Cartesian orientation of end-effecotr:\n')
      orie=getEEFCartesianOrientation( t_Kuka );
      disp(orie);
% orie =

%    [-3.1416]    [-1.1606e-04]    [-3.1415]
      %% Get force at end effector
      fprintf('Cartesian force acting at end effector:\n')
      f=getEEF_Force(t_Kuka);
      disp(f);
% f =

%    [-2.6678]    [-0.3954]    [10.5507]
      %% Get moment at end effector
      fprintf('moment at eef\n');
      m=getEEF_Moment(t_Kuka);
      disp(m);
% m =

%    [-0.1500]    [0.3462]    [-0.0949]
      %% turn off the server
       net_turnOffServer( t_Kuka );

       fclose(t_Kuka);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function check = reciprocal_check ( a, b )

%*****************************************************************************80
%
%% RECIPROCAL_CHECK checks the parameters of the Reciprocal CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A <= B.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( a <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'RECIPROCAL_CHECK - Fatal error!\n' );
    fprintf ( 1, '  A <= 0.0\n' );
    check = 0;
    return
  end

  if ( b < a ) then
    fprintf ( 1, '\n' );
    fprintf ( 1, 'RECIPROCAL_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B < A\n' );
    check = 0;
    return
  end

  check = 1;

  return
end


text:
classdef P2P < Algorithm

    methods (Access = public)

        function obj = P2P()
            obj.name = 'P2P';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        function result = compute(~,signal)
            result = peak2peak(signal);
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = 3 * n;
            memory = 1;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function display(SO3F,varargin)
% called by standard output

if check_option(varargin,'skipHeader')
  disp('  fibre component');
else
  displayClass(SO3F,inputname(1),[],'moreInfo',symChar(SO3F),varargin{:});
  if SO3F.antipodal, disp('  antipodal: true'); end
  disp(' ');
end

disp(['  kernel: ',char(SO3F.psi)]);
disp(['  fibre : ',char(round(SO3F.h)),' || ' char(round(SO3F.r))]);
disp(['  weight: ',xnum2str(mean(SO3F))]);
disp(' ');

end

author: mtex-toolbox

repo: mtex

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


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
%   - data                 MRI (head), Omega=(0,128)x(0,128), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             NGF
%   - pre-registration     rigid2D
%   - regularizer          mbElastic
%   - optimization         lBFGS
% version 2015/05/20
% ===============================================================================

close all, help(mfilename);
setup2DMRIData

imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e-1);
distance('reset','distance','NGF','edge',50);
trafo('reset','trafo','rigid2D');
regularizer('reset','regularizer','mbElastic','alpha',0.1,'mu',1,'lambda',0);

PIRpara = optPara('lBFGS','solver','backslash');
NPIRpara = optPara('lBFGS','solver',regularizer('get','solver'));

[yc,wc,his] = MLIR(ML,'PIRobj',@PIRBFGSobjFctn,'PIRpara',PIRpara,...
  'NPIRobj',@NPIRBFGSobjFctn,'NPIRpara',NPIRpara,...
  'minLevel',4,'maxLevel',7,'parametric',1,'plotMLiter',0);

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


text:
function [imBWGrid, imColorGridMasked, imBWBoardMask] = segColorGrid(whiteLight, colorGrid, camCorners, verbose)
%% Extracts the color grid on white board.
% See also: ImgProc.getMatchedNodes

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%

imLight = imread(whiteLight);
imColorGrid = imread(colorGrid);

% find seg mask
backSlashidx = strfind(whiteLight, '\');
whiteLightFileName = whiteLight(backSlashidx(end)+1:end);
setNumIdx = regexp(whiteLightFileName,'\d');

if(~isempty(camCorners))
    bwBoardName = fullfile(whiteLight(1:backSlashidx(end)), ['bwBoard', whiteLightFileName(setNumIdx),'.png']);
    if(verbose)
        figure('Name', 'segColorGrid', 'units','normalized','outerposition',[0 0 1 1]);
    end
else
    bwBoardName = fullfile(whiteLight(1:backSlashidx(end)), ['objectROI', whiteLightFileName(setNumIdx),'.png']);
end

%% Flood fill using the 1st checkerboard corner
% Convert RGB image into HSV color space.
% imHsv = rgb2hsv(imLight);
% imV = imHsv(:,:,3); % intensity channelim
% imLightEnhance = ImgProc.imadjust3(imLight);
% imGray = mat2gray(rgb2gray(imLightEnhance)); % normalized gray

% find color grid pixels
im1 = rgb2hsv(imColorGrid);
im2 = rgb2hsv(imLight);
s = mean(im1(:,:,3)) / mean(im2(:,:,3));
s = min(s, 1);
imGrid = imColorGrid-s*imLight;

%%
if(exist(bwBoardName,'file'))
    imBWBoard = imread(bwBoardName);
elseif(~isempty(camCorners))
    imGrayGrid = rgb2gray(imColorGrid-s*imLight);
    imGridROI = bwconvhull(bwareafilt(imGrayGrid > mean(imGrayGrid(:)), 1));
    imLightMasked = ImgProc.maskImage(imLight, imGridROI);
    imGrayMasked = rgb2gray(imLightMasked); % normalized gray

    imBWBoard = false(size(imGrayMasked));

    % get bounding box corners of the checkerboard
    maxXY = max(camCorners);
    minXY = min(camCorners);

    permsX = perms([minXY(:,1), maxXY(:,1)]);
    permsY = perms([minXY(:,2), maxXY(:,2)]);

    % check the best flood fill results when using the 4 corners
    winSize = 5;

    for i = 1:length(permsX)
        for j = 1:length(permsY)
            % init seed set to corner
            x = round(permsX(i));
            y = round(permsY(j));

            % find the brightest pixel index in a 5x5 local roi
            imRoi = imgaussfilt(imGrayMasked(y-winSize:y+winSize,x-winSize:x+winSize));
            [maxVal, maxIdx] = max(imRoi(:));

            % convert back to global row and col
            [row, col] = ind2sub(size(imRoi), maxIdx);
            col = col + x - winSize - 1;
            row = row + y - winSize - 1;

            % Flood fill using geodesic distance
            imCurSeg = ImgProc.segFloodFill(imLightMasked, row, col);
            imBWBoard = imBWBoard | imCurSeg;

            if(verbose)
                %             fs(imCurSeg);hold on; plot(col,row,'ro');
            end
        end
    end
else
    imBWBoard = true(size(imLight,1),size(imLight,2));
    % imBWBoard = bwconvhull(imBWBoard);

    % imWeight = graydiffweight(imLabNorm, col, row, 'GrayDifferenceCutoff', tol);
    % imBWBoard = imsegfmm(imWeight, col, row, 0.01);
end

if(verbose && ~isempty(camCorners))
    subplot(2,2,1);
    imshow(imBWBoard);
    title('Flood fill white board');
    drawnow
end

%% Crop checkerboard area
if(~isempty(camCorners))
    % find checkerboard area
    imBWCb = ~imBWBoard;
    winSize = 5;

    seeds = [];
    for i = 1:length(camCorners)
        x = round(camCorners(i,1));
        y = round(camCorners(i,2));
        imBWCb(y-winSize:y+winSize,x-winSize:x+winSize) = 1;

        % seeds for bwselect
        seeds = [seeds; [x,y]];
    end

    % select the checkerboard area using all corners as seeds
    imBWCb = bwselect(imBWCb, seeds(:,1), seeds(:,2), 8);

    imBWCb = bwconvhull(imBWCb);
    % imBWCb( round(minXY(:,2)):round(maxXY(:,2)), round(minXY(:,1)):round(maxXY(:,1))) = 1;
    imBWBoardMask = imfill(imBWBoard, 'holes');
    imBWBoardMask(imBWCb(:)) = 0;
else
    imBWCb = ~imBWBoard;
    imBWBoardMask = imBWBoard;
end

if(verbose && ~isempty(camCorners))
    subplot(2,2,2);
    imshow(imBWBoardMask);
    title('White board without checkerboard area');
    drawnow
end

%% Segment color grid
% im1 = rgb2hsv(imColorGrid);
% im2 = rgb2hsv(imLight);
% s = mean(im1(:,:,3)) / mean(im2(:,:,3));
% s = min(s, 1);
imColorGridMasked = ImgProc.maskImage(imGrid, imBWBoardMask);

% enhance color
imColorGridMasked = ImgProc.imadjust3(imColorGridMasked);

% to hsv
imHsv = rgb2hsv(imColorGridMasked);

% adaptive thresholding using Gaussian filter
sigma = 3;
imBWGrid = ImgProc.adaptiveThresh(imHsv(:,:,3), sigma);

% only keep the largest area
imBWGrid = bwareafilt(imBWGrid, 1);
% imBWGrid = bwmorph(imBWGrid,'clean', inf);

% close and open
se1 = strel('line',2,0);
se2 = strel('line',2,90);
imBWGrid = imclose(imBWGrid, se1);
imBWGrid = imclose(imBWGrid, se2);

% clean edges
imBWCb = imdilate(imBWCb, ones(8, 8));
imBWGrid(imBWCb(:)) = 0;

if(~isempty(camCorners))
    if(verbose)
        subplot(2,2,3);
        imshow(imColorGridMasked);
        title('Color grid masked');
        drawnow

        subplot(2,2,4);
        imshow(imBWGrid);
        title('Color grid binary mask');
        drawnow
    end
end

% imColorGridMasked = ImgProc.maskImage(imColorGridMasked, imBWGrid);
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
% Build and initialize the computational graph for regression based speech
% enhancement/dereverberation
%
function [layer, para] = Build_DereverbNet_Regression(Data_tr, para)
para.output = 'tmp';

if para.topology.useMasking
    layer = genNetworkDereverb_Masking(para.topology);     % generate the network graph
else
    layer = genNetworkDereverb_Regression(para.topology);
end
para.preprocessing{1} = {};                     % optional preprocessing for each data stream
para.preprocessing{2} = {};
para.cost_func.layer_idx = length(layer);       % specify which layers are cost function layers
para.cost_func.layer_weight = [1];              % set the weights of each cost function layer
para = ParseOptions2(para);

% generating the scaling factor for the input, as we will need to use a
% small constant in the logarithm. We need to make sure that the power of
% speech are larger than this constant most of the time.
scale = 1e4;        % we hard code the scale to be a constant so that all network will use the same number
scale = scale/2^16; % note that we are using int16 to store waveform samples, so need to scale down
if para.topology.useWav
    layer = InitWavScaleLayer(layer, scale);
end

if strcmpi(para.topology.RegressionNetType, 'DNN')    % if use DNN, splice the frames
    idx = ReturnLayerIdxByName(layer, 'splice');
else
    idx = ReturnLayerIdxByName(layer, 'delta'); % if use LSTM, use dynamic features
end
fft_net_length = idx(1);
fft_net = layer(1:fft_net_length);
paraTmp = para;
paraTmp.out_layer_idx = fft_net_length;
fprintf('Generate global MVN weights for mask subnet - %s\n', datestr(now));
[layer{fft_net_length+1}.W, layer{fft_net_length+1}.b] = computeGlobalCMVN(Data_tr, 100, paraTmp, fft_net);
VerifyPreprocessingTree(layer(1:fft_net_length+1), Data_tr, paraTmp, 100);

% set weight of static, velocity, and accelration features in the MSE cost
% function.
delta_idx = ReturnLayerIdxByName(layer, 'delta');
weight_idx = delta_idx+1;
for i=1:length(weight_idx)
    if isfield(layer{weight_idx(i)}, 'W') && numel(layer{weight_idx(i)}.W) == prod(layer{weight_idx(i)}.dim)
        continue;
    end
    layer{weight_idx(i)}.W = diag([ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(1); ...
        ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(2); ...
        ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(3)]);
    layer{weight_idx(i)}.b = zeros(para.topology.nFreqBin*3,1);
end

end

author: singaxiong

repo: SignalGraph


text:
function [fess,X] = fess_vbmc(vp,gp,X)
%FESS_VBMC Compute fractional effective sample size through importance sampling

if nargin < 3 || isempty(X); X = 100; end

% If a single number is passed, take it as the number of samples
if numel(X) == 1
    N = X;
    X = vbmc_rnd(vp,N,0);
else
    N = size(X,1);
end

% Can directly pass the estimated GP means instead of the full GP
if isstruct(gp)
    [~,~,fbar] = gplite_pred(gp,X,[],[],0,0);
else
    fbar = mean(gp,2);
end

if size(fbar,1) ~= size(X,1)
    error('Mismatch between number of samples from VP and GP.');
end

% Compute effective sample size (ESS) with importance sampling
vlnpdf = max(vbmc_pdf(vp,X,0,1),log(realmin));
logw = fbar - vlnpdf;
w = exp(logw - max(logw));
w = w/sum(w);
fess = 1/sum(w.^2) / N; % fractional ESS

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [fMedianMc, fMeanMc, fMc_org, v1Sigma, fStd_Mc, vMc] = calc_BstMc(mCatalog, fBinning, nSample, nMethod)
% function [fMedianMc, fMeanMc, fMc_org, v1Sigma, fStd_Mc, vMc] = calc_BstMc(mCatalog, fBinning, nSample, nMethod)
%---------------------------------------------------------------------------------------------------------
% Bootstrap EQ catalog and determine Mc choosing a specific method
%
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Magnitude binning interval
% nSample    : Number of bootstrap samples
% nMethod    : Method to determine the magnitude of completeness
%               1: Maximum curvature
%               2: Fixed Mc = minimum magnitude (Mmin)
%               3: Mc90 (90% probability)
%               4: Mc95 (95% probability)
%               5: Best combination (Mc95 - Mc90 - maximum curvature)
%               6: Mc using EMR-method
%               7: Mc due b using Shi & Bolt uncertainty
%               8: Mc due b using bootstrap uncertainty
%               9: Mc due b Cao-criterion
%
% Outgoing variables:
% fMedianMc  : Median (50 percentile) of vMc using nSample bootstraps
% fMeanMc    : Mean of Mc using nSample bootstraps
% fMc_org    : Mc estimate from original dataset
% v1Sigma    : [16 84] percentiles, according to standard deviation for normal distribution
% fStd_Mc    : Second moment of empirical Mc distribution (assuming normal distribution)
% vMc        : Mc estimates from bootstrap samples (empirical distribution)
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 28.05.03

% Check input
if nargin == 0, error('No catalog input'); end
if nargin == 1, fBinning = 0.1, nSample = 200, disp('Default Bin size 0.1, Default 200 bootstrap samples');end
if nargin == 2, nSample = 200, disp('Default 200 bootstrap samples'); end
if nargin == 3, nMethod = 6, disp('Default method: EMR-method'); end
if nargin > 4, error('Too many arguments!'); end

% Initialize
vMc = [];

% Get magnitudes
vMags = mCatalog(:,6);

% Create bootstrap samples using bootstrap matlab toolbox
mMag_bstsamp = bootrsp(vMags,nSample);
% Calculate Mc of original catalog
fMc_org = calc_Mc(mCatalog, nMethod);
% Determine Mc uncertainty
for nSamp=1:nSample
    mCatalog(:,6) = mMag_bstsamp(:,nSamp);
    fMc = calc_Mc(mCatalog, nMethod);
    vMc =  [vMc; fMc];
end

% Check for Nan and create output
vSel = isnan(vMc);
vNoNanMc = vMc(~vSel,:);
if ~isempty(vNoNanMc)
    fStd_Mc = calc_StdDev(vNoNanMc);
    fMeanMc = mean(vNoNanMc);
    fMedianMc = median(vNoNanMc);
else
    fStd_Mc = NaN;
    fMeanMc = NaN;
    fMedianMc = NaN;
end

if (~isempty(vNoNanMc) & length(vNoNanMc) > 1)
    v1Sigma = prctile(vNoNanMc,[16 84]);
elseif (~isempty(vNoNanMc)  &&  length(vNoNanMc) == 1)
    v1Sigma = prctile(vNoNanMc,[16 84]);
    v1Sigma = v1Sigma';
else
    v1Sigma = [NaN NaN];
end


text:
classdef log

methods(Static)

function display(message)
    % compute current time string
    curtime = datestr(now, 'HH:MM:SS');
    disp([curtime, ' ', message]);
end

end

end

author: indigits

repo: sparse-plex


text:
% Chapter 3 - Complex Iterative Maps.
% Program 3a - Julia Sets.
% Copyright Birkhauser 2013. Stephen Lynch.

% Plot the Julia set J(0,1.1) (Figure 3.1(d)).
clear
k=15;niter=2^k;
x=zeros(1,niter);y=zeros(1,niter);
x1=zeros(1,niter);y1=zeros(1,niter);
a=0;b=1.1;
x(1)=real(0.5+sqrt(0.25-(a+1i*b)));
y(1)=imag(0.5+sqrt(0.25-(a+1i*b)));
% Check that the point is unstable.
isunstable=2*abs(x(1)+1i*y(1))

hold on
for n=1:niter
    x1=x(n);y1=y(n);
    u=sqrt((x1-a)^2+(y1-b)^2)/2;v=(x1-a)/2;
    u1=sqrt(u+v);v1=sqrt(u-v);
    x(n+1)=u1;y(n+1)=v1;
    if y(n)<b
        y(n+1)=-y(n+1);
    end
    if rand < .5
       x(n+1)=-u1;y(n+1)=-y(n+1);
    end
end

fsize=15;
plot(x,y,'k.','MarkerSize',1)
set(gca,'XTick',-1.6:0.4:1.6,'FontSize',fsize)
set(gca,'YTick',-1.2:0.4:1.2,'FontSize',fsize)
xlabel('Re z','FontSize',fsize)
ylabel('Im z','FontSize',fsize)
hold off

% End of Program 3a.


text:
classdef S2AxisFieldHarmonic < S2AxisField
% a class represeneting a axis field on the sphere

properties
  sF;
end

properties(Dependent = true, Access = protected)
  xx;
  xy;
  yy;
  xz;
  yz;
  zz;
end

properties(Dependent = true)
  bandwidth
end

methods

  function sVF = S2AxisFieldHarmonic(sF, varargin)
    % initialize a spherical vector field
    if nargin == 0, return; end

    if length(sF) == 6
      sVF.sF = sF(:);
    end

  end

  function bw = get.bandwidth(sVF), bw = sVF.sF.bandwidth; end
  function xx = get.xx(sVF), xx = sVF.sF(1); end
  function xy = get.xy(sVF), xy = sVF.sF(2); end
  function yy = get.yy(sVF), yy = sVF.sF(3); end
  function xz = get.xz(sVF), xz = sVF.sF(4); end
  function yz = get.yz(sVF), yz = sVF.sF(5); end
  function zz = get.zz(sVF), zz = sVF.sF(6); end
  function sVF = set.xx(sVF, xx), sVF.sF(1) = xx; end
  function sVF = set.xy(sVF, xy), sVF.sF(2) = xy; end
  function sVF = set.yy(sVF, yy), sVF.sF(3) = yy; end
  function sVF = set.xz(sVF, xz), sVF.sF(4) = xz; end
  function sVF = set.yz(sVF, yz), sVF.sF(5) = yz; end
  function sVF = set.zz(sVF, zz), sVF.sF(6) = zz; end

end

methods(Static = true)
  sAF = quadrature(f, varargin)
  sAF = approximation(v, y, varargin)
  function sAF = normal
    sAF = S2AxisFieldHarmonic.quadrature(@(v) v(:),'bandwidth',2);
  end
end

end


text:
function PlotMap(cfig, map, path, scan, scanIdx)

world   = map.points;
scan = Transform(scan, path(:,end));

worldColor = [1 0 1];
scanColor = 'green';

% Plot
cfig(1); clf;
set(gca, 'color', [0,0,0]);
hold on; grid on; axis equal;

plot(world(:,1), world(:,2), '+', 'MarkerSize', 0.5, 'color', worldColor);
plot(scan(:,1),  scan(:,2),  '.', 'MarkerSize', 2, 'color', scanColor);
plot(path(1,:),  path(2,:),  '-.', 'LineWidth', 1, 'color', [0 1 1]);
for i = 1:20:length(scan)
    line([path(1,end), scan(i,1)], [path(2,end), scan(i,2)], 'color', [0.75 0.996 0.242]);
end
title(['Scan: ', num2str(scanIdx)]);
drawnow

% frame = getframe(cfig);
% cda = frame.cdata;
% frame.cdata = cda(30:1199, 270:2349, :);
% writeVideo(vip, frame);

author: meyiao

repo: LaserSLAM

save_path: github-repos/MATLAB/meyiao-LaserSLAM