text:
function [K] = kv1u1(x, y, xp, yp, hyp, i)

logsigma = hyp(1);
logthetax = hyp(2);
logthetay = hyp(3);

n_x = size(x,1);
n_y = size(y,1);
n_xp = size(xp,1);
n_yp = size(yp,1);

x = repmat(x,1,n_xp);
y = repmat(y,1,n_yp);
xp = repmat(xp',n_x,1);
yp = repmat(yp',n_y,1);

switch i

case 0

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(x+(-1).*xp).*(y+(-1).*yp);

case 1 % logsigma

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(x+(-1).*xp).*(y+(-1).*yp);

case 2 % logthetax

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*((-1)+(1/2).*exp(1).^((-1).*logthetax).*(x+(-1).*xp).^2).*(x+ ...
  (-1).*xp).*(y+(-1).*yp);

case 3 % logthetay

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(x+(-1).*xp).*((-1)+(1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(y+(-1).*yp);

otherwise

        K = zeros(n_x, n_xp);
end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
% v_rotro2eu_tab: Calculate tables needed for v_rotro2eu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 52 different rotation matrix patterns of -1,0,+1:                     %
%  1- 3: identity matrix rows in order: 123, 231, 312                   %
%  1- 3: negated identity matrix rows in order: 132, 213, 321           %
%  7-12: As 1-6 but with rows 2,3 negated                               %
% 13-18: As 1-6 but with rows 1,3 negated                               %
% 19-24: As 1-6 but with rows 1,2 negated                               %
% 25-33: +1 in position (i-24) and 0's in remainder of this row and col %
% 34-42: -1 in position (i-24) and 0's in remainder of this row and col %
% 43-51: 0 in position (i-42)                                           %
% 52: no special symmetry                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mes=[1:3 10:12 7:9 4:6]; % sign reversal look-up table
rtci=[2 3 5 6 8 9; 3 1 6 4 9 7; 1 2 4 5 7 8]';
rtsi=[3 2 6 5 9 8; 1 3 4 6 7 9; 2 1 5 4 8 7]';
rtr=[1 4 7 2 5 8 3 6 9]; % indices to transpose a vectorized 3x3 matrix
w6=ones(6,1); %
th6=3*w6;
x6=[2 1 2 1 2 1]'; % Index for sin components
scai=[0 0 0 1; 0 0 0 2; 0 0 0 3; 1 -1 0 1; 1 -1 0 2; 1 -1 0 3; 0 0 -1 1; 0 0 -1 2; 0 0 -1 3; -1 1 0 1; -1 1 0 2; -1 1 0 3]'; % [sin; -sin; cos; xyz] for fixed rotations
% create pattersn of non-zero entries
nzpatt=10*ones(3,3,52); % pattern of -1,0,+1
e3=eye(3);
nzpatt(:,:,1)=e3;
nzpatt(:,:,2)=e3([2 3 1],:);
nzpatt(:,:,3)=e3([3 1 2],:);
nzpatt(:,:,4)=-e3([1 3 2],:);
nzpatt(:,:,5)=-e3([2 1 3],:);
nzpatt(:,:,6)=-e3([3 2 1],:);
for j=1:3
    f3=-e3;
    f3(j,j)=1;
    for i=1:6
        nzpatt(:,:,i+6*j)=f3*nzpatt(:,:,i);
    end
end
for i=1:9
    ir=1+mod(i-1,3);
    ic=1+(i-ir)/3;
    nzpatt(:,ic,i+24)=0;
    nzpatt(ir,:,i+24)=0;
    nzpatt(ir,ic,i+24)=1;
        nzpatt(:,ic,i+33)=0;
    nzpatt(ir,:,i+33)=0;
    nzpatt(ir,ic,i+33)=-1;
    nzpatt(ir,ic,i+42)=0;
end
nzpattv=reshape(nzpatt,9,52); % vectorize the 3x3 matrices
nzpattc=reshape(sum(nzpatt~=0,1),3,52); % number of non-zero elements in each column
% now create transition map
trmap=zeros(52,12);  % result of applying transformation j to pattern i
zel=zeros(4,3,52); % elements to zero: [zero; non-zero; sine-sign; targ-sign],transformation,initial pattern
jm='xyz123456789'; % rotation patterns
for i=1:52
    for j=1:12
        rijv=reshape(v_rotqr2ro(v_roteu2qr(jm(j),pi/3))*nzpatt(:,:,i),9,1); % vectorized result of applying transformation
        rijv(abs(abs(abs(rijv)-0.5)-0.5)>1e-8)=10; % set entries to 10 unless close to -1,0,+1
        k=find(all(round(rijv)==nzpattv,1),1); % round to integers and find a match
        if isempty(k)
            error('cannot find match for (%d,%d)');
        else
            trmap(i,j)=k;
        end
        if j<=3
            icol=mod(find([nzpattc(:,i)==1 & nzpattc(:,k)==2;nzpattc(:,i)==2 & nzpattc(:,k)==3],1)-1,3)+1; % find the column to zero an element
            if ~isempty(icol)
                irow=(1:3)*(~nzpatt(:,icol,i) & nzpatt(:,icol,k)); % find zero that disappears
                jrow=6-j-irow; % find other row involved in rotation
                zel(:,j,i)=[3*icol-3+[irow; jrow]; mod(j-irow+1,3)-1; sign(nzpatt(jrow,icol,i))];
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print zel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('zel.txt','w');
fprintf(fid,'zel=reshape([');
for i=1:52
    if i>1
    fprintf(fid,';\n     ');
    end
    fprintf(fid,' %2d',zel(:,:,i));
end
fprintf(fid,']'',4,3,52);\n');
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print trmap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('trmap.txt','w');
fprintf(fid,'trmap=[');
for i=1:52
    if i>1
    fprintf(fid,';\n     ');
    end
    fprintf(fid,' %2d',trmap(i,:));
end
fprintf(fid,'];\n');
fclose(fid);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function gpcf = gpcf_periodic(varargin)
%GPCF_PERIODIC  Create a periodic covariance function for Gaussian Process
%
%  Description
%    GPCF = GPCF_PERIODIC('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates periodic covariance function structure in which the
%    named parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    GPCF = GPCF_PERIODIC(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Periodic covariance function with squared exponential decay
%    part as in Rasmussen & Williams (2006) Gaussian processes for
%    Machine Learning.
%
%    Parameters for periodic covariance function [default]
%      magnSigma2             - magnitude (squared) [0.1]
%      lengthScale            - length scale for each input [10]
%                               This can be either scalar
%                               corresponding isotropic or vector
%                               corresponding ARD
%      period                 - length of the periodic component(s) [1]
%      lengthScale_sexp       - length scale for the squared exponential
%                               component [10] This can be either scalar
%                               corresponding isotropic or vector
%                               corresponding ARD.
%      decay                  - determines whether the squared exponential
%                               decay term is used (1) or not (0).
%                               Not a hyperparameter for the function.
%      kalman_deg             - degree of approximation in type 'KALMAN' [6]
%      kalman_deg_sexp        - degree of sexp approximation in type
%                               'KALMAN' [6]
%      kalman_valid           - determines whether Bessel function is
%                               used (1) or not (0) in type 'KALMAN' [0]
%      magnSigma2_prior       - prior structure for magnSigma2 [prior_logunif]
%      lengthScale_prior      - prior structure for lengthScale [prior_t]
%      lengthScale_sexp_prior - prior structure for lengthScale_sexp
%                               [prior_fixed]
%      period_prior           - prior structure for period [prior_fixed]
%
%    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_*
%
% Copyright (c) 2009-2010 Heikki Peura
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2014 Arno Solin and Jukka Koskenranta

% 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_PERIODIC';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('magnSigma2',0.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('lengthScale',10, @(x) isvector(x) && all(x>0));
  ip.addParamValue('period',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('lengthScale_sexp',10, @(x) isvector(x) && all(x>0));
  ip.addParamValue('decay',0, @(x) isscalar(x) && (x==0||x==1));
  ip.addParamValue('kalman_deg',[], @(x) isscalar(x) && mod(x,1)==0);
  ip.addParamValue('kalman_deg_sexp',[], @(x) isscalar(x) && mod(x,1)==0);
  ip.addParamValue('kalman_valid',[], @(x) isscalar(x) && (x==0||x==1));
  ip.addParamValue('magnSigma2_prior',prior_logunif, @(x) isstruct(x) || isempty(x));
  ip.addParamValue('lengthScale_prior',prior_t, @(x) isstruct(x) || isempty(x));
  ip.addParamValue('lengthScale_sexp_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('period_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('selectedVariables',[], @(x) isempty(x) || ...
                   (isvector(x) && all(x>0)));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

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

  if init || ~ismember('magnSigma2',ip.UsingDefaults)
    gpcf.magnSigma2 = ip.Results.magnSigma2;
  end
  if init || ~ismember('lengthScale',ip.UsingDefaults)
    gpcf.lengthScale = ip.Results.lengthScale;
  end
  if init || ~ismember('period',ip.UsingDefaults)
    gpcf.period = ip.Results.period;
  end
  if init || ~ismember('lengthScale_sexp',ip.UsingDefaults)
    gpcf.lengthScale_sexp=ip . Results.lengthScale_sexp;
  end
  if init || ~ismember('decay',ip.UsingDefaults)
    gpcf.decay = ip.Results.decay;
  end
  if init || ~ismember('kalman_deg',ip.UsingDefaults)
      if ~isempty(ip.Results.kalman_deg)
          gpcf.kalman_deg = ip.Results.kalman_deg;
      end
  end
  if init || ~ismember('kalman_deg_sexp',ip.UsingDefaults)
      if ~isempty(ip.Results.kalman_deg_sexp)
          gpcf.kalman_deg_sexp = ip.Results.kalman_deg_sexp;
      end
  end
  if init || ~ismember('kalman_valid',ip.UsingDefaults)
      if ~isempty(ip.Results.kalman_valid)
          gpcf.kalman_valid = ip.Results.kalman_valid;
      end
  end
  if init || ~ismember('magnSigma2_prior',ip.UsingDefaults)
    gpcf.p.magnSigma2 = ip.Results.magnSigma2_prior;
  end
  if init || ~ismember('lengthScale_prior',ip.UsingDefaults)
    gpcf.p.lengthScale = ip.Results.lengthScale_prior;
  end
  if init || ~ismember('lengthScale_sexp_prior',ip.UsingDefaults)
    gpcf.p.lengthScale_sexp = ip.Results.lengthScale_sexp_prior;
  end
  if init || ~ismember('period_prior',ip.UsingDefaults)
    gpcf.p.period = ip.Results.period_prior;
  end

  if ~ismember('selectedVariables',ip.UsingDefaults)
    gpcf.selectedVariables = ip.Results.selectedVariables;
  end

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_periodic_pak;
    gpcf.fh.unpak = @gpcf_periodic_unpak;
    gpcf.fh.lp = @gpcf_periodic_lp;
    gpcf.fh.lpg = @gpcf_periodic_lpg;
    gpcf.fh.cfg = @gpcf_periodic_cfg;
    gpcf.fh.ginput = @gpcf_periodic_ginput;
    gpcf.fh.cov = @gpcf_periodic_cov;
    gpcf.fh.covvec = @gpcf_periodic_covvec;
    gpcf.fh.trcov  = @gpcf_periodic_trcov;
    gpcf.fh.trvar  = @gpcf_periodic_trvar;
    gpcf.fh.recappend = @gpcf_periodic_recappend;
    gpcf.fh.cf2ss = @gpcf_periodic_cf2ss;
  end

end

function [w, s, h] = gpcf_periodic_pak(gpcf)
%GPCF_PERIODIC_PAK  Combine GP covariance function parameters into
%                   one vector
%
%  Description
%    W = GPCF_PERIODIC_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.magnSigma2)
%             (hyperparameters of gpcf.magnSigma2)
%             log(gpcf.lengthScale(:))
%             (hyperparameters of gpcf.lengthScale)
%             log(gpcf.lengthScale_sexp)
%             (hyperparameters of gpcf.lengthScale_sexp)
%             log(gpcf.period)
%             (hyperparameters of gpcf.period)]'
%
%  See also
%    GPCF_PERIODIC_UNPAK

  if isfield(gpcf,'metric')
    error('Periodic covariance function not compatible with metrics.');
  else
    i1=0;i2=1;
    w = []; s = {}; h=[];

    if ~isempty(gpcf.p.magnSigma2)
      w = [w log(gpcf.magnSigma2)];
      s = [s; 'log(periodic.magnSigma2)'];
      h = [h 1];

      % Hyperparameters of magnSigma2
      [wh, sh, hh] = gpcf.p.magnSigma2.fh.pak(gpcf.p.magnSigma2);
      sh=strcat(repmat('prior-', size(sh,1),1),sh);
      w = [w wh];
      s = [s; sh];
      h = [h 1+hh];
    end

    if ~isempty(gpcf.p.lengthScale)
      w = [w log(gpcf.lengthScale)];
      s = [s; 'log(periodic.lengthScale)'];
      h = [h ones(1,numel(gpcf.lengthScale))];

      % Hyperparameters of lengthScale
      [wh, sh, hh] = gpcf.p.lengthScale.fh.pak(gpcf.p.lengthScale);
      sh=strcat(repmat('prior-', size(sh,1),1),sh);
      w = [w wh];
      s = [s; sh];
      h = [h 1+hh];
    end

    if gpcf.decay == 1 && ~isempty(gpcf.p.lengthScale_sexp)
      w = [w log(gpcf.lengthScale_sexp)];
      s = [s; 'log(periodic.lengthScale_sexp)'];
      h = [h ones(1,numel(gpcf.lengthScale_sexp))];

      % Hyperparameters of lengthScale_sexp
      [wh, sh, hh] = gpcf.p.lengthScale_sexp.fh.pak(gpcf.p.lengthScale_sexp);
      sh=strcat(repmat('prior-', size(sh,1),1),sh);
      w = [w wh];
      s = [s; sh];
      h = [h 1+hh];
    end

    if ~isempty(gpcf.p.period)
      w = [w log(gpcf.period)];
      s = [s; 'log(periodic.period)'];
      h = [h 1];

      % Hyperparameters of period
      [wh, sh, hh] = gpcf.p.period.fh.pak(gpcf.p.period);
      sh=strcat(repmat('prior-', size(sh,1),1),sh);
      w = [w wh];
      s = [s; sh];
      h = [h 1+hh];
    end
  end
end

function [gpcf, w] = gpcf_periodic_unpak(gpcf, w)
%GPCF_PERIODIC_UNPAK  Sets the covariance function parameters into
%                     the structure
%
%  Description
%    [GPCF, W] = GPCF_PERIODIC_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.magnSigma2)
%             (hyperparameters of gpcf.magnSigma2)
%             log(gpcf.lengthScale(:))
%             (hyperparameters of gpcf.lengthScale)
%             log(gpcf.lengthScale_sexp)
%             (hyperparameters of gpcf.lengthScale_sexp)
%             log(gpcf.period)
%             (hyperparameters of gpcf.period)]'
%
%  See also
%    GPCF_PERIODIC_PAK

  if isfield(gpcf,'metric')
    error('Covariance function not compatible with metrics');
  else
    gpp=gpcf.p;
    if ~isempty(gpp.magnSigma2)
      i1=1;
      gpcf.magnSigma2 = exp(w(i1));
      w = w(i1+1:end);
    end
    if ~isempty(gpp.lengthScale)
      i2=length(gpcf.lengthScale);
      i1=1;
      gpcf.lengthScale = exp(w(i1:i2));
      w = w(i2+1:end);
    end
    if ~isempty(gpp.lengthScale_sexp) && gpcf.decay == 1
      i2=length(gpcf.lengthScale_sexp);
      i1=1;
      gpcf.lengthScale_sexp = exp(w(i1:i2));
      w = w(i2+1:end);
    end
    if ~isempty(gpp.period)
      i2=length(gpcf.period);
      i1=1;
      gpcf.period = exp(w(i1:i2));
      w = w(i2+1:end);
    end
    % hyperparameters
    if ~isempty(gpp.magnSigma2)
      [p, w] = gpcf.p.magnSigma2.fh.unpak(gpcf.p.magnSigma2, w);
      gpcf.p.magnSigma2 = p;
    end
    if ~isempty(gpp.lengthScale)
      [p, w] = gpcf.p.lengthScale.fh.unpak(gpcf.p.lengthScale, w);
      gpcf.p.lengthScale = p;
    end
    if ~isempty(gpp.lengthScale_sexp)
      [p, w] = gpcf.p.lengthScale_sexp.fh.unpak(gpcf.p.lengthScale_sexp, w);
      gpcf.p.lengthScale_sexp = p;
    end
    if ~isempty(gpp.period)
      [p, w] = gpcf.p.period.fh.unpak(gpcf.p.period, w);
      gpcf.p.period = p;
    end

  end
end

function lp = gpcf_periodic_lp(gpcf)
%GPCF_PERIODIC_LP  Evaluate the log prior of covariance function parameters
%
%  Description
%    LP = GPCF_PERIODIC_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.
%
%    Also the log prior of the hyperparameters of the covariance
%    function parameters is added to E if hyperprior is
%    defined.
%
%  See also
%    GPCF_PERIODIC_PAK, GPCF_PERIODIC_UNPAK, GPCF_PERIODIC_LPG, GP_E

  lp = 0;
  gpp=gpcf.p;

  if isfield(gpcf,'metric')
    error('Covariance function not compatible with metrics');
  else
    % 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.

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

    if ~isempty(gpp.lengthScale_sexp) && gpcf.decay == 1
      lp = lp +gpp.lengthScale_sexp.fh.lp(gpcf.lengthScale_sexp, gpp.lengthScale_sexp) +sum(log(gpcf.lengthScale_sexp));
    end
    if ~isempty(gpcf.p.period)
      lp = gpp.period.fh.lp(gpcf.period, gpp.period) +sum(log(gpcf.period));
    end
  end

end

function lpg = gpcf_periodic_lpg(gpcf)
%GPCF_PERIODIC_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_PERIODIC_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_PERIODIC_PAK, GPCF_PERIODIC_UNPAK, GPCF_PERIODIC_LP, GP_G

  lpg = [];
  gpp=gpcf.p;
  if isfield(gpcf,'metric')
    error('Covariance function not compatible with metrics');
  end
  if ~isempty(gpcf.p.magnSigma2)
    lpgs = gpp.magnSigma2.fh.lpg(gpcf.magnSigma2, gpp.magnSigma2);
    lpg = [lpg lpgs(1).*gpcf.magnSigma2+1 lpgs(2:end)];
  end
  if ~isempty(gpcf.p.lengthScale)
    lll = length(gpcf.lengthScale);
    lpgs = gpp.lengthScale.fh.lpg(gpcf.lengthScale, gpp.lengthScale);
    lpg = [lpg lpgs(1:lll).*gpcf.lengthScale+1 lpgs(lll+1:end)];
  end
  if gpcf.decay == 1 && ~isempty(gpcf.p.lengthScale_sexp)
    lll = length(gpcf.lengthScale_sexp);
    lpgs = gpp.lengthScale_sexp.fh.lpg(gpcf.lengthScale_sexp, gpp.lengthScale_sexp);
    lpg = [lpg lpgs(1:lll).*gpcf.lengthScale_sexp+1 lpgs(lll+1:end)];
  end
  if ~isempty(gpcf.p.period)
    lpgs = gpp.period.fh.lpg(gpcf.period, gpp.period);
    lpg = [lpg lpgs(1).*gpcf.period+1 lpgs(2:end)];
  end
end

function DKff = gpcf_periodic_cfg(gpcf, x, x2, mask, i1)
%GPCF_PERIODIC_CFG  Evaluate gradient of covariance function
%                   with respect to the parameters
%
%  Description
%    DKff = GPCF_PERIODIC_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_PERIODIC_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_PERIODIC_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_PERIODIC_CFG(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), 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_PERIODIC_PAK, GPCF_PERIODIC_UNPAK, GPCF_PERIODIC_LP, GP_G

  gpp=gpcf.p;

  i2=1;
  gp_period=gpcf.period;
  DKff={};
  gprior=[];

  if nargin==5
    % Use memory save option
    savememory=1;
    if i1==0
      % Return number of hyperparameters
      i=0;
      if ~isempty(gpcf.p.magnSigma2)
        i=1;
      end
      if ~isempty(gpcf.p.lengthScale)
        i=i+length(gpcf.lengthScale);
      end
      if gpcf.decay==1 && ~isempty(gpcf.p.lengthScale_sexp)
        i=i+length(gpcf.lengthScale_sexp);
      end
      if ~isempty(gpcf.p.period)
        i=i+length(gpcf.lengthScale);
      end
      DKff=i;
      return
    end
  else
    savememory=0;
  end

  % Evaluate: DKff{1} = d Kff / d magnSigma2
  %           DKff{2} = d Kff / d lengthScale
  % 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))
    Cdm = gpcf_periodic_trcov(gpcf, x);

    ii1=0;
    if ~isempty(gpcf.p.magnSigma2) && (~savememory || all(i1==1));
        ii1=1;
        DKff{ii1} = Cdm;
    end

    if isfield(gpcf,'metric')
      error('Covariance function not compatible with metrics');
    else
      if isfield(gpcf,'selectedVariables')
        x = x(:,gpcf.selectedVariables);
      end
      [n, m] =size(x);
      if ~savememory
        i1=1:m;
      else
        if i1==1
          DKff=DKff{1};
          return
        else
          i1=i1-1;
        end
      end

      if ~isempty(gpcf.p.lengthScale) && (~savememory || all(i1 <= length(gpcf.lengthScale)))
        % loop over all the lengthScales
        if length(gpcf.lengthScale) == 1
          % In the case of isotropic PERIODIC
          s = 2./gpcf.lengthScale.^2;
          dist = 0;
          for i=1:m
            D = sin(pi.*bsxfun(@minus,x(:,i),x(:,i)')./gp_period);
            dist = dist + 2.*D.^2;
          end
          D = Cdm.*s.*dist;

          ii1 = ii1+1;
          DKff{ii1} = D;
        else
          % In the case ARD is used
          for i=i1
            s = 2./gpcf.lengthScale(i).^2;
            dist = sin(pi.*bsxfun(@minus,x(:,i),x(:,i)')./gp_period);
            D = Cdm.*s.*2.*dist.^2;

            ii1 = ii1+1;
            DKff{ii1} = D;
          end
        end
      end
      if savememory
        if length(DKff) == 1
          DKff=DKff{1};
          return
        end
        i1=i1-length(gpcf.lengthScale);
      end
      if gpcf.decay == 1
        if ~isempty(gpcf.p.lengthScale_sexp) && (~savememory || all(i1 <= length(gpcf.lengthScale_sexp)))
          if length(gpcf.lengthScale_sexp) == 1
            % In the case of isotropic PERIODIC
            s = 1./gpcf.lengthScale_sexp.^2;
            dist = 0;
            for i=1:m
              D = bsxfun(@minus,x(:,i),x(:,i)');
              dist = dist + D.^2;
            end
            D = Cdm.*s.*dist;

            ii1 = ii1+1;
            DKff{ii1} = D;
          else
            % In the case ARD is used
            for i=i1
              s = 1./gpcf.lengthScale_sexp(i).^2;
              dist = bsxfun(@minus,x(:,i),x(:,i)');
              D = Cdm.*s.*dist.^2;

              ii1 = ii1+1;
              DKff{ii1} = D;
            end
          end
        end
      end
      if savememory
        if length(DKff) == 1
          DKff=DKff{1};
          return
        end
        i1=i1-length(gpcf.lengthScale_sexp);
      end

      if ~isempty(gpcf.p.period)
        % Evaluate help matrix for calculations of derivatives
        % with respect to the period
        if length(gpcf.lengthScale) == 1
          % In the case of an isotropic PERIODIC
          s = repmat(1./gpcf.lengthScale.^2, 1, m);

          dist = 0;
          for i=1:m
            dist = dist + 2.*pi./gp_period.*sin(2.*pi.*bsxfun(@minus,x(:,i),x(:,i)')./gp_period).*bsxfun(@minus,x(:,i),x(:,i)').*s(i);
          end
          D = Cdm.*dist;
          ii1=ii1+1;
          DKff{ii1} = D;
        else
          % In the case ARD is used
          for i=i1
            s = 1./gpcf.lengthScale(i).^2;        % set the length
            dist = 2.*pi./gp_period.*sin(2.*pi.*bsxfun(@minus,x(:,i),x(:,i)')./gp_period).*bsxfun(@minus,x(:,i),x(:,i)');
            D = Cdm.*s.*dist;

            ii1=ii1+1;
            DKff{ii1} = D;
          end
        end
      end

    end
    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 || isempty(mask)
    if size(x,2) ~= size(x2,2)
      error('gpcf_periodic -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end

    K = gpcf.fh.cov(gpcf, x, x2);
    ii1=0;
    if ~isempty(gpcf.p.magnSigma2) && (~savememory || all(i1==1))
        ii1=1;
        DKff{ii1} = K;
    end

    if isfield(gpcf,'metric')
      error('Covariance function not compatible with metrics');
    else
      if isfield(gpcf,'selectedVariables')
        x = x(:,gpcf.selectedVariables);
        x2 = x2(:,gpcf.selectedVariables);
      end
      [n, m] =size(x);
      if ~savememory
        i1=1:m;
      else
        if i1==1
          DKff=DKff{1};
          return
        end
        i1=i1-1;
      end
      % Evaluate help matrix for calculations of derivatives with respect to the lengthScale
      if ~isempty(gpcf.p.lengthScale) && (~savememory || all(i1 <= length(gpcf.lengthScale)))
        if length(gpcf.lengthScale) == 1
          % In the case of an isotropic PERIODIC
          s = 1./gpcf.lengthScale.^2;
          dist = 0; dist2 = 0;
          for i=1:m
            dist = dist + 2.*sin(pi.*bsxfun(@minus,x(:,i),x2(:,i)')./gp_period).^2;
          end
          DK_l = 2.*s.*K.*dist;

          ii1=ii1+1;
          DKff{ii1} = DK_l;
        else
          % In the case ARD is used
          for i=i1
            s = 1./gpcf.lengthScale(i).^2;        % set the length
            dist = 2.*sin(pi.*bsxfun(@minus,x(:,i),x2(:,i)')./gp_period);
            DK_l = 2.*s.*K.*dist.^2;

            ii1=ii1+1;
            DKff{ii1} = DK_l;
          end
        end
      end
      if savememory
        if length(DKff) == 1
          DKff=DKff{1};
          return
        end
        i1=i1-length(gpcf.lengthScale);
      end

      if gpcf.decay == 1 && (~savememory || all(i1 <= length(gpcf.lengthScale_sexp)))
        % Evaluate help matrix for calculations of derivatives with
        % respect to the lengthScale_sexp
        if length(gpcf.lengthScale_sexp) == 1
          % In the case of an isotropic PERIODIC
          s = 1./gpcf.lengthScale_sexp.^2;
          dist = 0; dist2 = 0;
          for i=1:m
            dist = dist + bsxfun(@minus,x(:,i),x2(:,i)').^2;
          end
          DK_l = s.*K.*dist;

          ii1=ii1+1;
          DKff{ii1} = DK_l;
        else
          % In the case ARD is used
          for i=i1
            s = 1./gpcf.lengthScale_sexp(i).^2;        % set the length
            dist = bsxfun(@minus,x(:,i),x2(:,i)');
            DK_l = s.*K.*dist.^2;

            ii1=ii1+1;
            DKff{ii1} = DK_l;
          end
        end
      end
      if savememory
        if length(DKff) == 1
          DKff=DKff{1};
          return
        end
        i1=i1-length(gpcf.lengthScale_sexp);
      end

      if ~isempty(gpcf.p.period)
        % Evaluate help matrix for calculations of derivatives
        % with respect to the period
        if length(gpcf.lengthScale) == 1
          % In the case of an isotropic PERIODIC
          s = repmat(1./gpcf.lengthScale.^2, 1, m);
          dist = 0; dist2 = 0;
          for i=1:m
            dist = dist + 2.*pi./gp_period.*sin(2.*pi.*bsxfun(@minus,x(:,i),x2(:,i)')./gp_period).*bsxfun(@minus,x(:,i),x2(:,i)').*s(i);
          end
          DK_l = K.*dist;

          ii1=ii1+1;
          DKff{ii1} = DK_l;
        else
          % In the case ARD is used
          for i=i1
            s = 1./gpcf.lengthScale(i).^2;        % set the length
            dist = 2.*pi./gp_period.*sin(2.*pi.*bsxfun(@minus,x(:,i),x2(:,i)')./gp_period).*bsxfun(@minus,x(:,i),x2(:,i)');
            DK_l = s.*K.*dist;

            ii1=ii1+1;
            DKff{ii1} = DK_l;
          end
        end
      end
    end
    % Evaluate: DKff{1}    = d mask(Kff,I) / d magnSigma2
    %           DKff{2...} = d mask(Kff,I) / d lengthScale etc.
  elseif nargin == 4 || nargin == 5
    [n, m] =size(x);
    if isfield(gpcf,'metric')
      error('Covariance function not compatible with metrics');
    else
        ii1=0;
        if ~isempty(gpcf.p.magnSigma2) && (~savememory || all(i1==1))
            ii1=1;
            DKff{ii1} = gpcf.fh.trvar(gpcf, x);   % d mask(Kff,I) / d magnSigma2
        end
      for i2=1:length(gpcf.lengthScale)
        ii1 = ii1+1;
        DKff{ii1}  = 0;                          % d mask(Kff,I) / d lengthScale
      end
      if gpcf.decay == 1
        for i2=1:length(gpcf.lengthScale_sexp)
          ii1 = ii1+1;
          DKff{ii1}  = 0;                      % d mask(Kff,I) / d lengthScale_sexp
        end
      end
      if ~isempty(gpcf.p.period)
        ii1 = ii1+1;                             % d mask(Kff,I) / d period
        DKff{ii1}  = 0;
      end
    end
  end
  if savememory
    DKff=DKff{1};
  end
end

function DKff = gpcf_periodic_ginput(gpcf, x, x2, i1)
%GPCF_PERIODIC_GINPUT  Evaluate gradient of covariance function with
%                      respect to x
%
%  Description
%    DKff = GPCF_PERIODIC_GINPUT(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 X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_PERIODIC_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).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_PERIODIC_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), or k(X,X) if X2 is empty, with respect to ith
%    covariate in X. This subfunction is needed when using memory
%    save option in gp_set.
%
%  See also
%    GPCF_PERIODIC_PAK, GPCF_PERIODIC_UNPAK, GPCF_PERIODIC_LP, GP_G

  [n, m] =size(x);
  gp_period=gpcf.period;
  ii1 = 0;
  if length(gpcf.lengthScale) == 1
    % In the case of an isotropic PERIODIC
    s = repmat(1./gpcf.lengthScale.^2, 1, m);
    %gp_period = repmat(1./gp_period, 1, m);
  else
    s = 1./gpcf.lengthScale.^2;
  end
  if gpcf.decay == 1
    if length(gpcf.lengthScale_sexp) == 1
      % In the case of an isotropic PERIODIC
      s_sexp = repmat(1./gpcf.lengthScale_sexp.^2, 1, m);
    else
      s_sexp = 1./gpcf.lengthScale_sexp.^2;
    end
  end
  if nargin<4
    i1=1:m;
  else
    % Use memory save option
    if i1==0
      % Return number of covariates
      if isfield(gpcf,'selectedVariables')
        DKff=length(gpcf.selectedVariables);
      else
        DKff=m;
      end
      return
    end
  end

  if nargin == 2 || isempty(x2)
    K = gpcf.fh.trcov(gpcf, x);
    if isfield(gpcf,'metric')
      error('Covariance function not compatible with metrics');
    else
      for j = 1:n
        for i=i1
          DK = zeros(size(K));
          DK(j,:) = -s(i).*2.*pi./gp_period.*sin(2.*pi.*bsxfun(@minus,x(j,i),x(:,i)')./gp_period);
          if gpcf.decay == 1
            DK(j,:) = DK(j,:)-s_sexp(i).*bsxfun(@minus,x(j,i),x(:,i)');
          end
          DK = DK + DK';

          DK = DK.*K;      % dist2 = dist2 + dist2' - diag(diag(dist2));

          ii1 = ii1 + 1;
          DKff{ii1} = DK;
        end
      end
    end

  elseif nargin == 3
    K = gpcf.fh.cov(gpcf, x, x2);

    if isfield(gpcf,'metric')
      error('Covariance function not compatible with metrics');
    else
      ii1 = 0;
      for i=i1
        for j = 1:n
          DK= zeros(size(K));
          if gpcf.decay == 1
            DK(j,:) = -s(i).*2.*pi./gp_period.*sin(2.*pi.*bsxfun(@minus,x(j,i),x2(:,i)')./gp_period)-s_sexp(i).*bsxfun(@minus,x(j,i),x2(:,i)');
          else
            DK(j,:) = -s(i).*2.*pi./gp_period.*sin(2.*pi.*bsxfun(@minus,x(j,i),x2(:,i)')./gp_period);
          end
          DK = DK.*K;

          ii1 = ii1 + 1;
          DKff{ii1} = DK;
        end
      end
    end
  end
end

function C = gpcf_periodic_cov(gpcf, x1, x2)
%GP_PERIODIC_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_PERIODIC_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_PERIODIC_TRCOV, GPCF_PERIODIC_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
%   [n1,m1]=size(x1);
%   [n2,m2]=size(x2);
  gp_period=gpcf.period;

  if size(x1,2)~=size(x2,2)
    error('the number of columns of X1 and X2 has to be same')
  end

  if isfield(gpcf,'metric')
    error('Covariance function not compatible with metrics');
  else
    if isfield(gpcf,'selectedVariables')
      x1 = x1(:,gpcf.selectedVariables);
      x2 = x2(:,gpcf.selectedVariables);
    end
    [n1,m1]=size(x1);
    [n2,m2]=size(x2);

    C=zeros(n1,n2);
    ma2 = gpcf.magnSigma2;

    % Evaluate the covariance
    if ~isempty(gpcf.lengthScale)
      s = 1./gpcf.lengthScale.^2;
      if gpcf.decay == 1
        s_sexp = 1./gpcf.lengthScale_sexp.^2;
      end
      if m1==1 && m2==1
        dd = bsxfun(@minus,x1,x2');
        dist=2.*sin(pi.*dd./gp_period).^2.*s;
        if gpcf.decay == 1
          dist = dist + dd.^2.*s_sexp./2;
        end
      else
        % If ARD is not used make s a vector of
        % equal elements
        if size(s)==1
          s = repmat(s,1,m1);
        end
        if gpcf.decay == 1
          if size(s_sexp)==1
            s_sexp = repmat(s_sexp,1,m1);
          end
        end

        dist=zeros(n1,n2);
        for j=1:m1
          dd = bsxfun(@minus,x1(:,j),x2(:,j)');
          dist = dist + 2.*sin(pi.*dd./gp_period).^2.*s(:,j);
          if gpcf.decay == 1
            dist = dist +dd.^2.*s_sexp(:,j)./2;
          end
        end
      end
      dist(dist<eps) = 0;
      C = ma2.*exp(-dist);
    end

  end
end

function C = gpcf_periodic_trcov(gpcf, x)
%GP_PERIODIC_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_PERIODIC_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_PERIODIC_COV, GPCF_PERIODIC_TRVAR, GP_COV, GP_TRCOV

  if isfield(gpcf,'metric')
    error('Covariance function not compatible with metrics');
  else

    % Try to use the C-implementation
    C=trcov(gpcf, x);
    %     C = NaN;
    if isnan(C)
      % If there wasn't C-implementation do here
      if isfield(gpcf,'selectedVariables')
        x = x(:,gpcf.selectedVariables);
      end
      [n, m] =size(x);
      gp_period=gpcf.period;

      s = 1./(gpcf.lengthScale);
      s2 = s.^2;
      if size(s)==1
        s2 = repmat(s2,1,m);
        gp_period = repmat(gp_period,1,m);
      end
      if gpcf.decay == 1
        s_sexp = 1./(gpcf.lengthScale_sexp);
        s_sexp2 = s_sexp.^2;
        if size(s_sexp)==1
          s_sexp2 = repmat(s_sexp2,1,m);
        end
      end

      ma = gpcf.magnSigma2;

      C = zeros(n,n);
      for ii1=1:n-1
        d = zeros(n-ii1,1);
        col_ind = ii1+1:n;
        for ii2=1:m
          d = d+2.*s2(ii2).*sin(pi.*(x(col_ind,ii2)-x(ii1,ii2))./gp_period(ii2)).^2;
          if gpcf.decay == 1
            d=d+s_sexp2(ii2)./2.*(x(col_ind,ii2)-x(ii1,ii2)).^2;
          end
        end
        C(col_ind,ii1) = d;

      end
      C(C<eps) = 0;
      C = C+C';
      C = ma.*exp(-C);
    end
  end
end

function C = gpcf_periodic_trvar(gpcf, x)
%GP_PERIODIC_TRVAR  Evaluate training variance vector
%
%  Description
%    C = GP_PERIODIC_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_PERIODIC_COV, GP_COV, GP_TRCOV

  [n, m] =size(x);

  C = ones(n,1)*gpcf.magnSigma2;
  C(C<eps)=0;
end

function reccf = gpcf_periodic_recappend(reccf, ri, gpcf)
%RECAPPEND Record append
%
%  Description
%    RECCF = GPCF_PERIODIC_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_periodic';

    % Initialize parameters
    reccf.lengthScale= [];
    reccf.magnSigma2 = [];
    reccf.lengthScale_sexp = [];
    reccf.period = [];

    % Set the function handles
    reccf.fh.pak = @gpcf_periodic_pak;
    reccf.fh.unpak = @gpcf_periodic_unpak;
    reccf.fh.lp = @gpcf_periodic_lp;
    reccf.fh.lpg = @gpcf_periodic_lpg;
    reccf.fh.cfg = @gpcf_periodic_cfg;
    reccf.fh.cov = @gpcf_periodic_cov;
    reccf.fh.trcov  = @gpcf_periodic_trcov;
    reccf.fh.trvar  = @gpcf_periodic_trvar;
    reccf.fh.recappend = @gpcf_periodic_recappend;
    reccf.p=[];
    reccf.p.lengthScale=[];
    reccf.p.magnSigma2=[];
    if ri.decay == 1
      reccf.p.lengthScale_sexp=[];
      if ~isempty(ri.p.lengthScale_sexp)
        reccf.p.lengthScale_sexp = ri.p.lengthScale_sexp;
      end
    end

    reccf.p.period=[];
    if ~isempty(ri.p.period)
      reccf.p.period= ri.p.period;
    end
    if isfield(ri.p,'lengthScale') && ~isempty(ri.p.lengthScale)
      reccf.p.lengthScale = ri.p.lengthScale;
    end
    if ~isempty(ri.p.magnSigma2)
      reccf.p.magnSigma2 = ri.p.magnSigma2;
    end
  else
    % Append to the record

    gpp = gpcf.p;

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

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

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

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

    % record decay
    if ~isempty(gpcf.decay)
      reccf.decay(ri,:)=gpcf.decay;
    end

  end
end

function [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = gpcf_periodic_cf2ss(gpcf,x)
%GPCF_PERIODIC_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).
%
%  References:
%    Arno Solin and Simo Sarkka (2014). Explicit link between periodic
%    covariance functions and state space models. Accepted for
%    publication in Proceedings of the Seventeenth International
%    Conference on Artifcial Intelligence and Statistics (AISTATS 2014).
%

  % Check arguments
  if nargin < 2, x = []; end

  % Check if kalman_deg is given
  if isfield(gpcf,'kalman_deg')
      kalman_deg = gpcf.kalman_deg;
  else
      kalman_deg = [];
  end

  % Check if kalman_deg_sexp is given
  if isfield(gpcf,'kalman_deg_sexp')
      kalman_deg_sexp = gpcf.kalman_deg_sexp;
  else
      kalman_deg_sexp = [];
  end

  % Check if kalman_valid is given
  if isfield(gpcf,'kalman_valid')
      kalman_valid = gpcf.kalman_valid;
  else
      kalman_valid = [];
  end

  % Case squared exponential (i.e. quasi-periodic)
  if gpcf.decay

      % Return model matrices, derivatives and parameter information
      [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = ...
          cf_quasiperiodic_to_ss(gpcf.magnSigma2,gpcf.lengthScale, ...
          gpcf.period,gpcf.lengthScale_sexp,kalman_deg, ...
          inf,kalman_deg_sexp,kalman_valid);

      % Balance matrices for numerical stability
      [F,L,Qc,H,Pinf,dF,dQc,dPinf] = ...
          ss_balance(F,L,Qc,H,Pinf,dF,dQc,dPinf);

      % Check optimized parameters
      if isempty(gpcf.p.magnSigma2), ind(1) = false; else ind(1) = true; end
      if isempty(gpcf.p.lengthScale), ind(2) = false; else ind(2) = true; end
      if isempty(gpcf.p.lengthScale_sexp), ind(3) = false; else ind(3) = true; end
      if isempty(gpcf.p.period), ind(4) = false; else ind(4) = true; end

      % Change order: lengthScale_sexp <--> period
      % Use only optimized parameter gradients
      dF(:,:,[3,4])    = dF(:,:,[4,3]);
      dQc(:,:,[3,4])   = dQc(:,:,[4,3]);
      dPinf(:,:,[3,4]) = dPinf(:,:,[4,3]);

  % Case without squared exponential (i.e. purely periodic)
  else

      % Return model matrices, derivatives and parameter information
      [F,L,Qc,H,Pinf,dF,dQc,dPinf,params] = ...
          cf_periodic_to_ss(gpcf.magnSigma2,gpcf.lengthScale, ...
            gpcf.period,kalman_deg,kalman_valid);

      % Check optimization parameters
      if isempty(gpcf.p.magnSigma2), ind(1) = false; else ind(1) = true; end
      if isempty(gpcf.p.lengthScale), ind(2) = false; else ind(2) = true; end
      if isempty(gpcf.p.period), ind(3) = false; else ind(3) = true; end

  end

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

end

author: gpstuff-dev

repo: gpstuff

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


text:
function lf = leadfield_interpolate(pos, vol)

% LEADFIELD_INTERPOLATE interpolates the leadfield for a source at
% an arbitrary location given the pre-computed leadfields on a regular
% grid.
%
% Use as
%   lf = leadfield_interpolate(pos, vol)

% Copyright (C) 2013, Vladimir Litvak
%
% $Id$

% express the position in voxel coordinates
pos = ft_warp_apply(inv(vol.transform), pos);

lf = nan(length(vol.sens.label), 3*size(pos, 1));

use_splines = (vol.chan{1}.dat.dim(end) == 6);

for i = 1:length(vol.sens.label)
  for j = 1:3
    if use_splines
      lf(i, j:3:end) = spm_bsplins(squeeze(vol.chan{i}.dat(:, :, :, 3+j)), ...
        pos(:, 1), pos(:, 2), pos(:, 3), [4 4 4 0 0 0]);
    else
      lf(i, j:3:end) = spm_sample_vol(spm_vol([vol.filename{i} ',' num2str(j)]), ...
        pos(:, 1), pos(:, 2), pos(:, 3), 4);
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Make a linear dynamical system
%   X1 -> X2
%   |     |
%   v     v
%   Y1    Y2

intra = zeros(2);
intra(1,2) = 1;
inter = zeros(2);
inter(1,1) = 1;
n = 2;

X = 2; % size of hidden state
Y = 2; % size of observable state

ns = [X Y];
dnodes = [];
onodes = [2];
eclass1 = [1 2];
eclass2 = [3 2];
bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'eclass1', eclass1, 'eclass2', eclass2, ...
	      'observed', onodes);

x0 = rand(X,1);
V0 = eye(X);
C0 = rand(Y,X);
R0 = eye(Y);
A0 = rand(X,X);
Q0 = eye(X);

bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', x0, 'cov', V0, 'cov_prior_weight', 0);
bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(Y,1), 'cov', R0, 'weights', C0, ...
			   'clamp_mean', 1, 'cov_prior_weight', 0);
bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', zeros(X,1), 'cov', Q0, 'weights', A0, ...
			   'clamp_mean', 1, 'cov_prior_weight', 0);

T = 5; % fixed length sequences

clear engine;
engine{1} = kalman_inf_engine(bnet);
engine{2} = jtree_unrolled_dbn_inf_engine(bnet, T);
engine{3} = jtree_dbn_inf_engine(bnet);
N = length(engine);

% inference

ev = sample_dbn(bnet, T);
evidence = cell(n,T);
evidence(onodes,:) = ev(onodes, :);

t = 1;
query = [1 3];
m = cell(1, N);
ll = zeros(1, N);
for i=1:N
  [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence);
  m{i} = marginal_nodes(engine{i}, query, t);
end

% compare all engines to engine{1}
for i=2:N
  assert(approxeq(m{1}.mu, m{i}.mu));
  assert(approxeq(m{1}.Sigma, m{i}.Sigma));
  assert(approxeq(ll(1), ll(i)));
end

if 0
for i=2:N
  approxeq(m{1}.mu, m{i}.mu)
  approxeq(m{1}.Sigma, m{i}.Sigma)
  approxeq(ll(1), ll(i))
end
end

% learning

ncases = 5;
cases = cell(1, ncases);
for i=1:ncases
  ev = sample_dbn(bnet, T);
  cases{i} = cell(n,T);
  cases{i}(onodes,:) = ev(onodes, :);
end

max_iter = 2;
bnet2 = cell(1,N);
LLtrace = cell(1,N);
for i=1:N
  [bnet2{i}, LLtrace{i}] = learn_params_dbn_em(engine{i}, cases, 'max_iter', max_iter);
end

for i=1:N
  temp = bnet2{i};
  for e=1:3
    CPD{i,e} = struct(temp.CPD{e});
  end
end

for i=2:N
  assert(approxeq(LLtrace{i}, LLtrace{1}));
  for e=1:3
    assert(approxeq(CPD{i,e}.mean, CPD{1,e}.mean));
    assert(approxeq(CPD{i,e}.cov, CPD{1,e}.cov));
    assert(approxeq(CPD{i,e}.weights, CPD{1,e}.weights));
  end
end

% Compare to KF toolbox

data = zeros(Y, T, ncases);
for i=1:ncases
  data(:,:,i) = cell2num(cases{i}(onodes, :));
end
[A2, C2, Q2, R2, x2, V2, LL2trace] =  learn_kalman(data, A0, C0, Q0, R0, x0, V0, max_iter);

e = 1;
assert(approxeq(x2, CPD{e,1}.mean))
assert(approxeq(V2, CPD{e,1}.cov))
assert(approxeq(C2, CPD{e,2}.weights))
assert(approxeq(R2, CPD{e,2}.cov));
assert(approxeq(A2, CPD{e,3}.weights))
assert(approxeq(Q2, CPD{e,3}.cov));
assert(approxeq(LL2trace, LLtrace{1}))

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [outerproduct] = outer(tensor1, tensor2, squeezedimensions)
%C = OUTER(A, B, squeezedimensions) Computes the outer product
%C(i[1],...,i[m],j[1],...,j[n]) = A(i[1],...,i[m]) * B(j[1],...j[n]).
%Discards superfluous singleton dimensions if squeezedimensions ~= 0.
%
%Note: thanks to Emese Toth for suggesting the algorithm used here.
%
%Wynton Moore, January 2006

%store initial dimensions
dim1=size(tensor1);dim2=size(tensor2);

%evaluate
outerproduct=reshape(reshape(tensor1, [], 1)*reshape(tensor2, 1, []), [dim1 dim2]);

%discard superfluous singleton dimensions
if squeezedimensions
    outerproduct=squeeze(outerproduct);
end

author: Sable

repo: mcbench-benchmarks


text:
function [ zmat, seed ] = c8mat_uniform_01 ( m, n, seed )

%*****************************************************************************80
%
%% C8MAT_UNIFORM_01 returns a unit pseudorandom C8MAT.
%
%  Discussion:
%
%    The angles should be uniformly distributed between 0 and 2 * PI,
%    the square roots of the radius uniformly distributed between 0 and 1.
%
%    This results in a uniform distribution of values in the unit circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the matrix.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, double complex ZMAT(M,N), the pseudorandom complex matrix.
%
%    Output, integer SEED, a seed for the random number generator.
%
  for i2 = 1 : n
    for i1 = 1 : m

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      r = sqrt ( seed * 4.656612875E-10 );

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      theta = 2.0 * pi * seed * 4.656612875E-10;

      zmat(i1,i2) = r * ( cos ( theta ) + i * sin ( theta ) );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [vITRS,rotMat]=TIRS2ITRS(vTIRS,TT1,TT2,xpyp)
%%TIRS2ITRS Rotate vectors from the Terrestrial Intermediate Reference
%            System (TIRS) to the International Terrestrial Reference
%            System (ITRS). The ITRS is essentially the WGS-84 coordinate
%            system: it defines locations with respect to the crust of a
%            non-rotating Earth, where the z axis passes through a fixed
%            point on the surface. On the other hand, the TIRS is nearly
%            the same except the z axis is the axis of rotation of the
%            Earth, which slowly varies over time. Note that the velocity
%            conversion does not include the (small) centrifugal effect of
%            polar motion.
%
%INPUTS: x The NXnumVec collection of vectors in TIRS coordinates to
%          convert (units do not matter). N can be 3, or 6. If the vectors
%          are 3D, then they are position. 6D
%          vectors are assumed to be position and velocity. Since the TIRS
%          and ITRS co-rotate, there is no Coriolis effect to add. Also,
%          the accelerations due to the wobble of the rotation axis over
%          time are not considered. These accelerations are very small.
%          Thus, the function just rotates both halves of the vector.
% TT1, TT2 Two parts of a Julian date given in terrestrial time (TT).
%          The units of the date are days. The full date is the sum of
%          both terms. The date is broken into two parts to provide
%          more bits of precision. It does not matter how the date is
%          split.
%     xpyp xpyp=[xp;yp] are the polar motion coordinates in radians
%          including the effects of tides and librations. If this
%          parameter is omitted or if an empty matrix is passed, the
%          value from the function getEOP will be used.
%
%OUTPUTS: vITRS The NXnumVec vector of values of x rotated from the TIRS
%               into the ITRS.
%        rotMat The 3X3 rotation matrix used to rotate vectors from the
%               TIRS into the ITRS.
%
%The conversion functions from the International Astronomical Union's
%(IAU) Standard's of Fundamental Astronomy library are put together to get
%the necessary rotation matrix.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[vITRS,rotMat]=TIRS2ITRS(vTIRS,TT1,TT2)
%or if more parameters are known, using the format
%[vITRS,rotMat]=TIRS2ITRS(vTIRS,TT2,TT2,xpyp);
%
%Different celestial coordinate systems are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report,
%    Naval Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016,
%    173 pages.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function cbrt_test ( )

%*****************************************************************************80
%
%% CBRT_TEST tests R4_CBRT and R8_CBRT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CBRT_TEST:\n' );
  fprintf ( 1, '  Test CBRT_VALUES, R4_CBRT, R8_CBRT.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X        CBRT(X)\n' );
  fprintf ( 1, '                   R4_CBRT(X)         Diff\n' );
  fprintf ( 1, '                   R8_CBRT(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = cbrt_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_cbrt ( single ( x ) );
    fx3 = r8_cbrt ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end


text:
function value = r8_aid ( x )

%*****************************************************************************80
%
%% R8_AID evaluates the derivative of the Airy function Ai of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the derivative of the Airy function
%    Ai of X.
%
  persistent aifcs
  persistent aigcs
  persistent naif
  persistent naig
  persistent x2sml
  persistent x3sml

  if ( isempty ( naif ) )
    aifcs = [ ...
       0.105274612265314088088970057325134114, ...
       0.011836136281529978442889292583980840, ...
       0.000123281041732256643051689242469164, ...
       0.000000622612256381399016825658693579, ...
       0.000000001852988878441452950548140821, ...
       0.000000000003633288725904357915995625, ...
       0.000000000000005046217040440664768330, ...
       0.000000000000000005223816555471480985, ...
       0.000000000000000000004185745090748989, ...
       0.000000000000000000000002672887324883, ...
       0.000000000000000000000000001392128006, ...
       0.000000000000000000000000000000602653, ...
       0.000000000000000000000000000000000220 ]';
    aigcs = [ ...
       0.0212338781509186668523122276848937, ...
       0.0863159303352144067524942809461604, ...
       0.0017975947203832313578033963225230, ...
       0.0000142654998755506932526620687495, ...
       0.0000000594379952836832010488787064, ...
       0.0000000001524033664794478945214786, ...
       0.0000000000002645876603490435305100, ...
       0.0000000000000003315624296815020591, ...
       0.0000000000000000003139789757594792, ...
       0.0000000000000000000002325767379040, ...
       0.0000000000000000000000001384384231, ...
       0.0000000000000000000000000000676629, ...
       0.0000000000000000000000000000000276 ]';
    eta = 0.1 * r8_mach ( 3 );
    naif = r8_inits ( aifcs, 13, eta );
    naig = r8_inits ( aigcs, 13, eta );
    x3sml = r8_mach ( 3 )^0.3334;
    x2sml = sqrt ( r8_mach ( 3 ) );
  end

  if ( x < - 1.0 )
    [ xn, phi ] = r8_admp ( x );
    value = xn * cos ( phi );
  elseif ( abs ( x ) <= x2sml )
    x2 = 0.0;
    x3 = 0.0;
    value = ( x2 * ( 0.125 + r8_csevl ( x3, aifcs, naif ) ) ...
      - r8_csevl ( x3, aigcs, naig ) ) - 0.25;
  elseif ( abs ( x ) <= x3sml )
    x2 = x * x;
    x3 = 0.0;
    value = ( x2 * ( 0.125 + r8_csevl ( x3, aifcs, naif ) ) ...
      - r8_csevl ( x3, aigcs, naig ) ) - 0.25;
  elseif ( x <= 1.0 )
    x2 = x * x;
    x3 = x * x * x;
    value = ( x2 * ( 0.125 + r8_csevl ( x3, aifcs, naif ) ) ...
      - r8_csevl ( x3, aigcs, naig ) ) - 0.25;
  else
    value = r8_aide ( x ) ...
      * exp ( - 2.0 * x * sqrt ( x ) / 3.0 );
  end

  return
end


text:
function conv2conv()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    conv_layer_idx = get_conv_layer_idx_from_layer_idx(curr_layer_idx);
    P_reshape = reshape(mem.activations{curr_layer_idx-1}', config.feature_map_sizes{curr_layer_idx-1}(1), config.feature_map_sizes{curr_layer_idx-1}(2)*config.feature_map_sizes{curr_layer_idx-1}(3)*config.batch_size);
    P_expand = config.IM2COL(P_reshape, [config.kernel_size(conv_layer_idx, 1), config.kernel_size(conv_layer_idx, 2)]);
    P_expand(:, mem.conv2conv{curr_layer_idx}{2}) = 0;
    input_size = [0 0];
    input_size(1) = config.kernel_size(conv_layer_idx, 1)*config.kernel_size(conv_layer_idx, 2)*config.feature_map_sizes{curr_layer_idx-1}(3);
    input_size(2) = (config.feature_map_sizes{curr_layer_idx-1}(1)-config.kernel_size(conv_layer_idx, 1)+1)*(config.feature_map_sizes{curr_layer_idx-1}(2)-config.kernel_size(conv_layer_idx, 2)+1)*config.batch_size;
    mem.layer_inputs{curr_layer_idx} = reshape(accumarray(mem.conv2conv{curr_layer_idx}{1}, P_expand(:), [input_size(1)*input_size(2), 1]), input_size(1), input_size(2));
    % memory redundancy, need improve
    mem.orig_activation_size{curr_layer_idx-1} = size(mem.activations{curr_layer_idx-1});
    mem.activations{curr_layer_idx-1} = mem.layer_inputs{curr_layer_idx};
end

author: jimmy-ren


text:
function rtk=uddcb_ppp(rtk)

VAR_DCB=30^2;
idx=rtk.id+1;

if rtk.x(idx)==0
    rtk=initx(rtk,1e-6,VAR_DCB,idx);
end

return

repo: GINav


text:
%example_mat_dual

clear all;
clc
clear opt;
% add path
addpath(genpath('../../SLEP/'));

% load data and set regularization parameter
load('../../data/scene.mat');
lambda = 10^-4;

% center data
D = CenterRowData(D);
L = CenterRowData(L);

% call the main function
[W,fmin] = mat_dual(D,L,lambda);

author: zzstefan

repo: BrainNetClass


text:
%% Housekeeping
clear
close all
clc

%% read the models and their calibrations

m0=rise('frwz_nk');

%% solve the model

m=solve(m0,'steady_state_unique',true);

part_list={};

m2=solve(m0,'solve_perturbation_type',{'frwz',part_list});

%% print results
m.print_solution()

%% print solution for a subset of variables only
m.print_solution({'PAI','Y'})

%% Alternative calibrations

cal_2=struct('psi_a_2', 0.7);

m2=set(m,'parameters',cal_2);

%% construct a vector of models

bigm=[m,m2];

%% compute impulse responses for all models simultaneously

myirfs=irf(bigm,'irf_periods',20);

%% plot the impulse responses
close all
var_list=m.endogenous.name;
figure('name','Impulse responses to a monetary policy shock');
for ii=1:numel(var_list)
    subplot(3,1,ii)
    reg1=myirfs.EPS_R.regime_1.(var_list{ii});
    reg2=myirfs.EPS_R.regime_2.(var_list{ii});
    plot([reg1,reg2]);
    title(var_list{ii})
    if ii==1
        legend({'reg1-m1','reg1-m2','reg2-m1','reg2-m2'})
    end
    axis tight
end

author: jmaih


text:
function test_suite=test_moxunit_fieldtrip_test_suite_filtering
% #   For MOxUnit_fieldtrip's copyright information and license terms,   #
% #   see the COPYING file distributed with MOxUnit_fieldtrip.           #

    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_fieldtrip_suite_filter_basics()
% test that object can be instantiated and has expected properties

    % each element in the struct below containes these elements:
    % - the label for suite.filter ('maxwalltime','maxmem','loadfile')
    % - the value for suite.filter
    % - lines for which test cases should be kept    after the filter
    % - "                                  " removed "              "
    s=struct();
    s.walltime={'maxwalltime',...
                {20,'20'},...
                {'','% WALLTIME 10','% WALLTIME 00:00:10'},...
                {'% WALLTIME 30','% WALLTIME 00:05:00'}};

    s.mem={'maxmem',...
                {10000,'10000'},...
                {'','% MEM 10','% MEM 1kb'},...
                {'% MEM 20000','% MEM 1tb'}};
    s.dccnpath={'loadfile',...
                {false,'no'},...
                {'','% dccnpath'},...
                {'dccnpath','fullfile(dccnpath,''bar'')'}};

    keys=fieldnames(s);
    for k=1:numel(keys)
        key=keys{k};
        prop=s.(key);

        label=prop{1};
        values=prop{2};

        for skip_test=[false,true]
            prop_idx=3;
            if skip_test
                prop_idx=prop_idx+1;
            end
            line_cell=prop{prop_idx};

            for j=1:numel(line_cell)
                line=line_cell{j};
                for m=1:numel(values)
                    value=values{m};

                    test_dir=tempname();
                    fn=randtestfilename(test_dir);

                    helper_test_single_test_case_filter(fn,...
                                                        label,...
                                                        value,...
                                                        line,...
                                                        skip_test);
                end
            end
        end
    end

function test_fieldtrip_suite_filter_exclude_if_prefix_equals_failed()
    test_dir=tempname();

    prefix='failed_';

    arg='exclude_if_prefix_equals_failed';
    no_yes_values={{'no',false},...
                    {'yes',true}};
    for k=1:numel(no_yes_values)
        exclude=k>1;
        for j=1:numel(no_yes_values{k})
            value=no_yes_values{k}{j};

            fn=randtestfilename(test_dir,prefix);
            helper_test_single_test_case_filter(fn,arg,value,'',exclude);
        end
    end

function helper_test_single_test_case_filter(fn,arg,value,...
                                                line,...
                                                skip_test)

    % make test case
    test_dir=fileparts(fn);
    suite=make_basic_test_suite(test_dir);
    dir_cleaner=onCleanup(@()clean_dir(test_dir));

    tc=build_test_case(fn, {line});

    % suite is empty initially
    assertEqual(countTestCases(suite),0);

    % after adding the test it must be in the suite
    suite=addFromFile(suite,fn);
    assertEqual(countTestCases(suite),1);
    assertEqual(getLocation(getTestNode(suite,1)),...
                    getLocation(tc));

    % skip test
    suite=filter(suite,arg,value);

    % tests must still be there
    assertEqual(countTestCases(suite),1);

    report=MOxUnitTestReport(0,1,'');
    report=run(suite,report);

    assertEqual(countTestOutcomes(report),1);

    % get test outcome.
    outcome=getTestOutcome(report,1);

    % test should not have failed
    assertEqual(isNonFailure(outcome),true);

    % succesful if test was not skipped
    is_success=isSuccess(outcome);
    assertEqual(is_success,~skip_test);

    if skip_test
        assertTrue(isa(outcome,'MOxUnitSkippedTestOutcome'));
    end

function s=randstr
    s=char(ceil(rand(1,10)*26+64));

function suite=make_basic_test_suite(test_dir)
    suite=MOxUnitFieldTripTestSuite();
    mkdir(test_dir);

function fn=randtestfilename(parent_dir,prefix)
    if nargin<2
        prefix='test_';
    end
    test_name=randstr();
    func_name=sprintf('%s%s',prefix,test_name);
    fn=fullfile(parent_dir,sprintf('%s.m',func_name));

function clean_dir(dir_name)
    ds=dir(dir_name);
    n=numel(ds);
    for k=1:n
        d=ds(k);
        if ~d.isdir
            % delete file
            pth=fullfile(dir_name,d.name);
            delete(pth);
        end
    end

    rmdir(dir_name);

function tc=build_test_case(fn, body_lines)
    [unused,func_name]=fileparts(fn);
    lines=[sprintf('function %s',func_name),...
            body_lines];
    write_content(fn,lines);
    tc=MOxUnitFieldTripTestCase(fn);

function write_content(fn,lines)
    fid=fopen(fn,'w');
    cleaner=onCleanup(@()fclose(fid));
    fprintf(fid,'%s\n',lines{:});


text:
function [subs,vals,keepLocs] = nonMaxSuprWin( subs, vals, ...
  strLocs, endLocs, thresh, maxn )
% Nonmaximal suppression of values outside of a given window.
%
% Suppresses all location in subs that do not fall in given range (defined
% by strLocs and endLocs).  For example, if subs are 3D coordinates of
% maxes over an array of size siz,  "nonMaxSuprWin( subs, vals,
% [1,1,1]+10, siz-10 )" suppreses all locations within 10 pixels of the
% border of I.
%
% USAGE
%  [subs,vals,keepLocs] = nonMaxSuprWin( subs, vals, ...
%                                      strLocs, endLocs, thresh, maxn )
%
% INPUTS
%  subs        - subscripts of point locations (m x d)
%  vals        - values at point locations (m x 1)
%  strLocs     - locations at which to start cropping along each dim
%  endLocs     - locations at which to end cropping along each dim
%  thresh      - [] minimum value below which not to look fo
%  maxn        - [] return at most maxn of the largest vals
%
% OUTPUTS
%  subs        - subscripts of non-suppressed point locations (n x d)
%  vals        - values at non-suppressed point locations (n x 1)
%  keepLocs    - indicies of kept locations from subs (n x 1)
%
% EXAMPLE
%
% See also SUBSTOARRAY, NONMAXSUPR, NONMAXSUPRLIST
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

vals = vals(:);  nvals = length(vals);  nd=size(subs,2);
if( nargin<5 || isempty(thresh)); thresh=[]; end
if( nargin<6 || isempty(maxn)); maxn = []; end

[strLocs,er] = checkNumArgs( strLocs, [1 nd], 0, 0 ); error(er);
[endLocs,er] = checkNumArgs( endLocs, [1 nd], 0, 0 ); error(er);
if (any(strLocs>endLocs)); error('strLocs must be <= endLocs'); end

% discard vals below thresh
if (~isempty(thresh))
  keepLocs = vals > thresh;
else
  keepLocs = true( nvals, 1 );
end

% suppress all values outside of window defined by start and end locs
for d=1:nd
  if (strLocs(d)>0)
    keepLocsi = (subs(:,d)>=strLocs(d)) & (subs(:,d)<=endLocs(d));
    keepLocs = keepLocs & keepLocsi;
  end
end

% suppress all but the first maxn nonzero elts in keepLocs
if( ~isempty(maxn) && maxn>0 && maxn<sum(keepLocs) )
  [~,order] = sort( -vals ); [~,unorder]=sort(order);
  keepLocs = keepLocs(order);
  keepLocs( cumsum( keepLocs )>maxn ) = 0;
  keepLocs = keepLocs(unorder);
end

% discard locations where keepLocs==0
vals = vals( keepLocs ); subs = subs( keepLocs, : );

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
%islscalar - Test if parameter is a (pseudo) logical scalar.
%
%  USAGE
%
%    test = islscalar(x)
%
%    x              parameter to test
%
%  NOTE
%
%    To be considered logical, the scalar should be equal to 0 or 1, but it does
%    not need to actually be of class 'logical' (class(x) could be e.g. 'double').
%
%  SEE ALSO
%
%    See also islvector, islmatrix, isdmatrix, isdvector, isdscalar, isimatrix, isivector,
%    isstring_FMAT.
%

% Copyright (C) 2010-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function test = islscalar(x)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help islscalar'' for details).');
end

% Test: double, scalar
test = islogical(x) & isscalar(x);

function test = islogical(x)

test = builtin('islogical',x) | all(x(:)==0|x(:)==1);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function res = recenterTransform3d(transfo, center)
%RECENTERTRANSFORM3D Change the fixed point of an affine 3D transform.
%
%   TRANSFO2 = recenterTransform3d(TRANSFO, CENTER)
%   where TRANSFO is a 4x4 transformation matrix, and CENTER is a 1x3 row
%   vector, computes the new transformations that uses the same linear part
%   (defined by the upper-left 3x3 corner of the transformation matrix) as
%   the initial transform, and that will leave the point CENTER unchanged.
%
%
%
%   Example
%   % creating a re-centered rotation using:
%   rot1 = createRotationOx(pi/3);
%   rot2 = recenterTransform3d(rot1, [3 4 5]);
%   % will give the same result as:
%   rot3 = createRotationOx([3 4 5], pi/3);
%
%
%   See also
%   transforms3d, createRotationOx, createRotationOy, createRotationOz
%   createTranslation3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-27,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% remove former translation part
res = eye(4);
res(1:3, 1:3) = transfo(1:3, 1:3);

% create translations
t1 = createTranslation3d(-center);
t2 = createTranslation3d(center);

% compute translated transform
res = t2*res*t1;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function sc =  dbscale(q)

%tstoolbox/@unit/dbscale
%   returns scaling value when calculating decibel values from data of
%   this unit. dpscale returns either 10 (for power or energy units (e.g.
%   Watt)) or 20 (for all other units (e.g. Volt).
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

sc = q.dBScale;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%GOBJ_POINT Create point geometry object.
%
%   [ GOBJ ] = GOBJ_POINT( P, TAG ) Creates a point
%   geometry object. Accepts the following input parameters.
%
%       Parameter   Value/{Default}           Description
%       -----------------------------------------------------------------------------------
%       p           array  {[0,0]}            Coordinates of point
%       tag         string {P1}               Geometry object tag/name

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [X,Time,Z] = load_data(datafile,Spec,sample,loadExcel)
%loadData Load vintage of data from file and format as structure
%
% Description:
%
%   Load data from file
%
% Input Arguments:
%
%   datafile - filename of Microsoft Excel workbook file
%
% Output Arguments:
%
% Data - structure with the following fields:
%   .    X : T x N numeric array, transformed dataset
%   . Time : T x 1 numeric array, date number with observation dates
%   .    Z : T x N numeric array, raw (untransformed) dataset

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

if nargin < 4
    loadExcel = 0;
end

ext = datafile(find(datafile=='.',1,'last'):end); % file extension

idx = strfind(datafile,filesep); idx = idx(2);
datafile_mat = fullfile(datafile(1:idx),'mat',[strtok(datafile(idx+1:end),'.') '.mat']);
if exist(datafile_mat,'file')==2 && ~loadExcel
    % Load raw data from MATLAB formatted binary (.mat) file
    d = load(datafile_mat);
    Z    = d.Z;
    Time = d.Time;
    Mnem = d.Mnem;
elseif any(strcmpi(ext,{'.xlsx','.xls'}))
    % Read raw data from Excel file
    [Z,Time,Mnem] = readData(datafile);
    save(datafile_mat,'Z','Time','Mnem');
    %    Z : raw (untransformed) observed data
    % Time : observation periods for the time series data
    % Mnem : series ID for each variable
else
    error('Only Microsoft Excel workbook files supported.');
end

% Sort data based on model specification
Z = sortData(Z,Mnem,Spec);
clear Mnem; % since now Mnem == Spec.SeriesID

% Transform data based on model specification
[X,Time,Z] = transformData(Z,Time,Spec);

% Drop data not in estimation sample
if nargin == 3
    [X,Time,Z] = dropData(X,Time,Z,sample);
end

end

function [Z,Time,Mnem] = readData(datafile)
%readData Read data from Microsoft Excel workbook file

[DATA,TEXT] = xlsread(datafile,'data');
Mnem = TEXT(1,2:end);
if ispc
    Time = datenum(TEXT(2:end,1),'mm/dd/yyyy');
    Z    = DATA;
else
    Time = DATA(:,1) + datenum(1899,12,31) - 1;
    Z    = DATA(:,2:end);
end

end

function [Z,Mnem] = sortData(Z,Mnem,Spec)
%sortData Sort series by order of model specification

% Drop series not in Spec
inSpec = ismember(Mnem,Spec.SeriesID);
Mnem(~inSpec) = [];
Z(:,~inSpec)  = [];

% Sort series by ordering of Spec
N = length(Spec.SeriesID);
permutation = NaN(1,N);
for i = 1:N
    idxSeriesSpec = find(strcmp(Spec.SeriesID{i},Mnem));
    permutation(i) = idxSeriesSpec;
end

Mnem = Mnem(permutation);
Z    = Z(:,permutation);

end

function [X,Time,Z] = transformData(Z,Time,Spec)
%transformData Transforms each data series based on Spec.Transformation
%
% Input Arguments:
%
%      Z : T x N numeric array, raw (untransformed) observed data
%   Spec : structure          , model specification
%
% Output Arguments:
%
%      X : T x N numeric array, transformed data (stationary to enter DFM)

T = size(Z,1);
N = size(Z,2);

X = NaN(T,N);
for i = 1:N
    formula = Spec.Transformation{i};
    freq    = Spec.Frequency{i};
    switch freq  % time step for different frequencies based on monthly time
        case 'm' % Monthly
            step = 1;
        case 'q' % Quarterly
            step = 3;
    end
    t1 = step;    % assume monthly observations start at beginning of quarter
    n  = step/12; % number of years, needed to compute annual % changes
    series  = Spec.SeriesName{i};
    switch formula
        case 'lin' % Levels (No Transformation)
            X(:,i) = Z(:,i);
        case 'chg' % Change (Difference)
            X(t1:step:T,i) = [NaN; Z(t1+step:step:T,i) - Z(t1:step:T-t1,i)];
        case 'ch1' % Year over Year Change (Difference)
            X(12+t1:step:T,i) = Z(12+t1:step:T,i) - Z(t1:step:T-12,i);
        case 'pch' % Percent Change
            X(t1:step:T,i) = 100*[NaN; Z(t1+step:step:T,i) ./ Z(t1:step:T-t1,i) - 1];
        case 'pc1' % Year over Year Percent Change
            X(12+t1:step:T,i) = 100*(Z(12+t1:step:T,i) ./ Z(t1:step:T-12,i) - 1);
        case 'pca' % Percent Change (Annual Rate)
            X(t1:step:T,i) = 100*[NaN; (Z(t1+step:step:T,i) ./ Z(t1:step:T-step,i)).^(1/n) - 1];
        case 'log' % Natural Log
            X(:,i) = log(Z(:,i));
        otherwise
            warning(['Transformation ''' formula ''' not found for ' series '. ' 'Using untransformed data.']);
            X(:,i) = Z(:,i);
    end
end

% Drop first quarter of observations
% since transformations cause missing values
Time = Time(4:end);
Z    = Z(4:end,:);
X    = X(4:end,:);

end

function [X,Time,Z] = dropData(X,Time,Z,sample)
%dropData Remove data not in estimation sample

idxDrop = Time < sample;

Time(idxDrop) = [];
X(idxDrop,:)  = [];
Z(idxDrop,:)  = [];

end

author: FRBNY-TimeSeriesAnalysis

repo: Nowcasting


text:
function y = gprctile(x,p)
%GPRCTILE Percentile(s) of a grouped sample.
%  In some scientific works, once the data have been gathered from a
%  population of interest, it is often difficult to get a sense of what
%  the data indicate when they are presented in an unorganized fashion.
%  Assembling the raw data into a meaningful form, such as a frequency
%  distribution, makes the data easier to understand and interpret. It is
%  in the context of frequency distributions that the importance of
%  conveying in a succinct way numerical information contained in the data
%  is encountered.
%  So, grouped data is data that has been organized into groups known as
%  classes.  The raw dataset can be organized by constructing a table
%  showing the frequency distribution of the variable (whose values are
%  given in the raw dataset). Such a frequency table is often referred to
%  as grouped data.
%  Here, we developed a m-code to calculate the percentile(s) of a grouped
%  data. One can input the returns or modified vectors n and xout
%  containing the frequency counts and the bin locations of the hist
%  m-function, in a column form matrix.
%
%  Percentile calculation uses the straight forward formula,
%
%                  P = L + I*(N*P/100 - C)/F
%
%  where:
%  L = lower limit of the interval containing the percentile
%  I = width of the interval containing the percentile
%  N = total number of data
%  P = interested percentile
%  C = cumulative frequency corresponding to the previous percentile class
%  F = number of cases in the interval containing the percentile
%
%  Syntax: function y = gprctile(x,p)
%
%  Inputs:
%       x - data matrix (Size of matrix must be n-by-2; absolut frequency=
%           column 1, class mark=column 2)
%       p - is a scalar or a vector of percent values
%  Outputs:
%       y  - percentile(s) of the values in x
%
%  Example: From the example given at
%  http://www.emathzone.com/tutorials/basic-statistics/frequency-distributi
%  on-of-discrete-data.html
%  We are interested to get the percentile values 15,30,60,80.
%
%  Data: 2,4,6,1,3,5,3,7,8,6,4,7,4,4,2,1,3,6,4,2,5,7,9,1,2,10,1,8,9,2,3,1,
%  2,3,4,4,4,6,6,5,5,4,5,8,5,4,3,3,2,5,0,5,9,9,8,10,0,4,10,10,1,1,2,2,1,8,
%  6,9,10
%
%  Using [A,B] = hist(x,11)
%
%  x = [B' A']; p = [15,30,60,80];
%
%                     ----------------
%                        MC       F
%                     ----------------
%                      0.4545     2
%                      1.3636     8
%                      2.2727     9
%                      3.1818     7
%                      4.0909    11
%                      5.0000     8
%                      5.9091     6
%                      6.8182     3
%                      7.7273     5
%                      8.6364     5
%                      9.5455     5
%                     ----------------
%
%  Calling on Matlab the function:
%          y = gprctile(x,p)
%
%  Answer is:
%
%  y = 1.8535    2.9481    5.0455    7.4909
%
%  Created by A. Trujillo-Ortiz, R. Hernandez-Walls and R. Preciado-Perez
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.edu.mx
%
%  Copyright (C)  September 15, 2012.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A., R. Hernandez-Walls and R. Preciado-Perez. (2012).
%     gprctile:Percentile(s) of a grouped sample. [WWW document].
%     URL http://www.mathworks.com/matlabcentral/fileexchange/
%     38228-gprctile
%
%  Reference:
%  Cann, A. J. (2003), Maths from Scratch for Biologists. John Willey &
%              Sons Ltd. Chichester:England.
%

c = size(x,2);

if c ~= 2
    error('stats:gprctile:BadData','X must have two colums.');
end

if ~isvector(p) || numel(p) == 0
    error('stats:gprctile:BadProbs', ...
          'P must be a scalar or a non-empty vector.');
elseif any(p < 0 | p > 100) || ~isreal(p)
    error('stats:gprctile:BadPercents', ...
          'P must take real values between 0 and 100');
end

mc = x(:,1); %class mark
f = x(:,2); %absolut frequency
a = mc(2) - mc(1); %class interval
fa = cumsum(f); %cumulative absolut frequency
n = fa(end); %sample size

y = [];
for i = 1:length(p)
    r(i) = n*p(i)/100; %interested ordened data (iod)
    cl(i) = min(find(fa > r(i))); %class where the iod are
    l(i) = mc(cl(i));
    L(i) = l(i) - 0.5*a; %lower class boundary
    A(i) =(a/f(cl(i)));
    if (cl(i) - 1) == 0;
        B(i) = r(i);
    else
        B(i) = (r(i) - fa(cl(i) - 1));
    end

    yy = L(i) + A(i)*B(i); %calculation by a linear interpolation
    y = [y yy];
end

return,


text:
function [i,j] = argmin2(x)
%ARGMIN2  Index of minimum element of matrix.
% [i,j] = ARGMIN2(x) returns indices (i,j) such that x(i,j) == min(x(:)).
%
% See also ARGMIN, ARGMAX2.

[colmin,i] = min(x);
[ignore,j] = min(colmin);
i = i(j);

author: andrewssobral

repo: lrslibrary


text:
function randint_test ( )

%*****************************************************************************80
%
%% RANDINT_TEST shows how random integers are generated in MATLAB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDINT_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the MATLAB RANDINT function.\n' );

  randint_test01 ( );

  randint_test02 ( );

  seed = 123456789;
  randint_test03 ( seed );

  seed = 987654321;
  randint_test03 ( seed );

  seed = 123456789;
  randint_test03 ( seed );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDINT_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function randint_test01 ( )

%*****************************************************************************80
%
%% RANDINT_TEST01 simply calls the random integer generator a few times.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDINT_TEST01:\n' );
  fprintf ( 1, '  In MATLAB, random integers are generated by calling RANDINT:\n' );
  fprintf ( 1, '  If a range is not specified, the values are 0 or 1.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A = randint ( )      a random scalar value.\n' );
  fprintf ( 1, '  B = randint ( 5, 1 ) a random column vector of 5 entries.\n' );
  fprintf ( 1, '  C = randint ( 1, 5 ) a random row vector of 5 entries.\n' );
  fprintf ( 1, '  D = randint ( 3, 4 ) a 3 by 4 random matrix.\n' );
  fprintf ( 1, '  E = randint ( 5 )    a 5 by 5 random matrix.\n' );

  a = randint ( )
  b = randint ( 5, 1 )
  c = randint ( 1, 5 )
  d = randint ( 3, 4 )
  e = randint ( 5 )

  return
end
function randint_test02 ( )

%*****************************************************************************80
%
%% RANDINT_TEST02 specifies the range.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDINT_TEST02:\n' );
  fprintf ( 1, '  RANDINT allows the user to specify the numeric range.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A = randint (  1, 1, [ 5,    10 ] ) a random scalar value.\n' );
  fprintf ( 1, '  B = randint ( 10, 1, [ 7,     8 ] ) a random column vector of 5 entries.\n' );
  fprintf ( 1, '  C = randint (  1, 5, [ -1,   +1 ] ) a random row vector of 5 entries.\n' );
  fprintf ( 1, '  D = randint (  3, 4, [ -5,   +5 ] ) a 3 by 4 random matrix.\n' );
  fprintf ( 1, '  E = randint (  5, 5, [ 100, 200 ] ) a 5 by 5 random matrix.\n' );

  a = randint ( 1,  1, [   5,  10 ] )
  b = randint ( 10, 1, [   7,   8 ] )
  c = randint ( 1,  5, [  -1,  +1 ] )
  d = randint ( 3,  4, [  -5,  +5 ] )
  e = randint ( 5,  5, [ 100, 200 ] )

  return
end
function randint_test03 ( seed )

%*****************************************************************************80
%
%% RANDINT_TEST03 sets the seed before calling RANDINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RANDINT_TEST03:\n' );
  fprintf ( 1, '  By setting the random number seed, you can control\n' );
  fprintf ( 1, '  how the random number sequence begins.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The command "rng ( 123456789 )" sets the seed to 123456789.\n' );

  rng ( seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Seed has been set to %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now generate 5 random values.\n' );

  for i = 1 : 5
    a = randint ( 1, 1, [ 1, 100 ] );
    fprintf ( 1, '  RANDINT(1,1,[1,100]) = %g\n', a );
  end

  rng ( seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Seed has been reset to %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now generate 5 more random values.\n' );

  for i = 1 : 5
    a = randint ( 1, 1, [ 1, 100 ] );
    fprintf ( 1, '  RANDINT(1,1,[1,100]) = %g\n', a );
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
%--- help for format_transition_probabilities ---
%
%  INTERNAL FUNCTION
%
%


text:
function cy = zaxpy ( n, ca, cx, incx, cy, incy )

%*****************************************************************************80
%
%% ZAXPY computes a complex constant times a vector plus a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of elements in CX and CY.
%
%    Input, complex CA, the multiplier of CX.
%
%    Input, complex CX(*), the first vector.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Input, complex CY(*), the second vector.
%
%    Input, integer INCY, the increment between successive entries of CY.
%
%    Output, complex CY(*), the value of CY(*) + CA * CX(*).
%
  cy(1:incy:1+(n-1)*incy) = cy(1:incy:1+(n-1)*incy) ...
    + ca * cx(1:incx:1+(n-1)*incx);

  return
end


text:
% Demo for Structured Edge Detector (please see readme.txt first).
addpath('./toolbox/');
%% set opts for training (see edgesTrain.m)
opts=edgesTrain();                % default options (good settings)
opts.modelDir='models/';          % model will be in models/forest
opts.modelFnm='modelBsds';        % model name
opts.nPos=5e5; opts.nNeg=5e5;     % decrease to speedup training
opts.useParfor=0;                 % parallelize if sufficient memory

%% train edge detector (~20m/8Gb per tree, proportional to nPos/nNeg)
tic, model=edgesTrain(opts); toc; % will load model if already trained

%% set detection parameters (can set after training)
model.opts.multiscale=0;          % for top accuracy set multiscale=1
model.opts.sharpen=2;             % for top speed set sharpen=0
model.opts.nTreesEval=4;          % for top speed set nTreesEval=1
model.opts.nThreads=4;            % max number threads for evaluation
model.opts.nms=0;                 % set to true to enable nms

%% evaluate edge detector on BSDS500 (see edgesEval.m)
if(0), edgesEval(model, 'show',1, 'name','' ); end

%% detect edge and visualize results

path = './data/input/';
files = dir([path, '*.jpg']);
Imgnum = size(files, 1);
result_path = './data/edge_tone/';
mkdir(result_path);
for i=1:Imgnum
    I = imread(strcat(path, files(i,1).name));
    if size(I, 3) < 3
        I = repmat(I, [1,1,3]);
    end
    %%
    E=edgesDetect(I,model);
    imwrite(1-E, strcat(result_path, files(i,1).name(1:end-4),'_edge.jpg'));
    %%
    I=rgb2gray(I);
    nhoodSize = [11, 11];
    smoothValue=0.5*diff(getrangefromclass(I)).^2;
    I=imguidedfilter(I,I,'NeighborhoodSize', nhoodSize, 'DegreeOfSmoothing',smoothValue);
    imwrite(I,strcat(result_path,files(i,1).name(1:end-4),'_gf.jpg'));
end

author: Yijunmaverick


text:
function value = p27_autonomous ( )

%*****************************************************************************80
%
%% P27_AUTONOMOUS reports whether problem P27 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 0;

  return
end


text:
%% Copyright (C) 2014, 2016-2017, 2022-2023 Colin B. Macdonald
%% Copyright (C) 2016 Lagu
%% Copyright (C) 2016 Abhinav Tripathi
%% Copyright (C) 2017 NVS Abhilash
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@sym subsref {(@var{f}, @var{idx})}
%% @defopx Operator @@sym {@var{f}(@var{i})} {}
%% @defopx Operator @@sym {@var{f}(@var{i}, @var{j})} {}
%% @defopx Operator @@sym {@var{f}(@var{i}:@var{j})} {}
%% @defopx Operator @@sym {@var{f}.property} {}
%% Access entries of a symbolic array.
%%
%% Examples:
%% @example
%% @group
%% A = sym([10 11 12]);
%% A(2)
%%   @result{} (sym) 11
%%
%% A(2:3)
%%   @result{} (sym) [11  12]  (1×2 matrix)
%%
%% A(1, 1)
%%   @result{} (sym) 10
%% @end group
%%
%% @group
%% A.flat
%%   @result{} Matrix([[10, 11, 12]])
%% @end group
%% @end example
%%
%% @seealso{@@sym/subsasgn, @@sym/subsindex, @@sym/end}
%% @end defop

function varargout = subsref (f, idx)

  switch idx.type
    case '()'
      % sym(sym) indexing gets here
      for i = 1:length(idx.subs)
        if (isa(idx.subs{i}, 'sym'))
          idx.subs{i} = subsindex(idx.subs{i})+1;
        end
      end
      for i = 1:length(idx.subs)
        if (~ is_valid_index(idx.subs{i}))
          error('OctSymPy:subsref:invalidIndices', ...
                'invalid indices: should be integers or boolean');
        end
      end
      out = mat_access(f, idx.subs);

    case '.'
      fld = idx.subs;
      if (strcmp (fld, 'pickle'))
        out = f.pickle;
      elseif (strcmp (fld, 'flat'))
        out = f.flat;
      elseif (strcmp (fld, 'ascii'))
        out = f.ascii;
      elseif (strcmp (fld, 'unicode'))
        out = f.unicode;
      %elseif (strcmp (fld, 'extra'))
      %  out = f.extra;
      % not part of the interface
      %elseif (strcmp (fld, 'size'))
      %  out = f.size;
      else
        error ('@sym/subsref: invalid or nonpublic property ''%s''', fld);
      end

    otherwise
      error ('@sym/subsref: invalid subscript type ''%s''', idx.type);

  end

  % Octave 7.1, Issue #1112 and https://savannah.gnu.org/bugs/?61898
  varargout = {out};

end

%!shared a,b
%! b = [1:4];
%! a = sym(b);
%!assert(isequal( a(1), b(1) ))
%!assert(isequal( a(2), b(2) ))
%!assert(isequal( a(4), b(4) ))
%!assert(isempty( a([]) ))

%!shared a,b
%! b = [1:4];  b = [b; 3*b; 5*b];
%! a = sym(b);
%!assert(isequal( a(1), b(1) ))
%!assert(isequal( a(2), b(2) ))
%!assert(isequal( a(4), b(4) ))
%!assert(isequal( a(:,:), a ))
%!assert(isequal( a(1:2,1:3), a(1:2,1:3) ))
%!assert(isequal( a(1:2:3,[1 2 4]), b(1:2:3,[1 2 4]) ))
%!assert(isequal( a(1:2:3,[4 2 3 1]), b(1:2:3,[4 2 3 1]) ))
% repeats
%!assert(isequal( a(1:2:3,[4 1 1 1]), b(1:2:3,[4 1 1 1]) ))

%!assert(isequal( a([],:), b([],:) ))
%!assert(isequal( size(a([],:)), [0 4] ))
%!assert(isequal( a(1:2,[]), b(1:2,[]) ))
%!assert(isequal( size(a(1:2,[])), [2 0] ))
%!assert(isempty( a(1:2,[]) ))
%!assert(isempty( a([],[]) ))
%!assert(isequal( a([],[]), sym([]) ))
%!assert(~isequal( a(1:2,[]), sym([]) ))

%!shared e
%! e = sym([1 3 5; 2 4 6]);
%!assert(isequal( e(:), sym((1:6)') ))
%!assert(isequal( e([1 2 3]), sym([1 2 3]) ))
%!assert(isequal( e([1; 3; 4]), sym([1; 3; 4]) ))
%!assert(isempty( e([]) ))
%!assert(isempty( e('') ))
%!assert(isequal( e([]), sym([]) ))

%!shared a,b
%! b = 1:5; a = sym(b);
%!assert(isequal(  a([1 2 5]),  b([1 2 5])  ))
%!assert(isequal(  a([1; 2; 5]),  b([1; 2; 5])  ))

%!shared x
%! syms x

%!test
%! % logical with empty result
%! assert(isempty( x(false) ))
%! a = [x x];
%! assert(isempty( a([false false]) ))

%!test
%! % issue 18, scalar access
%! assert(isequal( x(1), x ))
%! assert(isequal( x(true), x ))

%!shared

%!test
%! % older access tests
%! syms x
%! f = [x 2; 3 4*x];
%! % element access
%! assert (logical(  f(1,1) == x  ))
%! assert (logical(  f(1,2) == 2  ))
%! % linear access of 2d array
%! assert (logical(  f(1) == x  ))
%! assert (logical(  f(2) == 3  ))  % column based
%! assert (logical(  f(3) == 2  ))

%!shared a,b
%! % effectively a random matrix
%! a = reshape( round(50*(sin(1:20)+1)),  5,4);
%! b = sym(a);

%!test
%! % older array refs test
%! assert (logical(b(1,1) == a(1,1)))
%! assert (logical(b(3,1) == a(3,1)))
%! assert (logical(b(1,3) == a(1,3)))
%! assert (logical(b(4,4) == a(4,4)))

%!test
%! % older array refs test: linear indices
%! assert (logical(b(1) == a(1)))
%! assert (logical(b(3) == a(3)))
%! assert (logical(b(13) == a(13)))

%!test
%! % older array refs test: end
%! assert (all(all(logical(  b(end,1) == a(end,1)  ))))
%! assert (all(all(logical(  b(2,end) == a(2,end)  ))))
%! assert (all(all(logical(  b(end,end) == a(end,end)  ))))
%! assert (all(all(logical(  b(end-1,1) == a(end-1,1)  ))))
%! assert (all(all(logical(  b(2,end-1) == a(2,end-1)  ))))
%! assert (all(all(logical(  b(end-1,end-1) == a(end-1,end-1)  ))))

%!shared

%!test
%! % older slicing tests
%! syms x
%! a = [1 2 3 4 5 6]; a = [a; 3*a; 5*a; 2*a; 4*a];
%! b = sym(a);
%! assert (isequal(  b(:,1), a(:,1)  ))
%! assert (isequal(  b(:,2), a(:,2)  ))
%! assert (isequal(  b(1,:), a(1,:)  ))
%! assert (isequal(  b(2,:), a(2,:)  ))
%! assert (isequal(  b(:,:), a(:,:)  ))
%! assert (isequal(  b(1:3,2), a(1:3,2)  ))
%! assert (isequal(  b(1:4,:), a(1:4,:)  ))
%! assert (isequal(  b(1:2:5,:), a(1:2:5,:)  ))
%! assert (isequal(  b(1:2:4,:), a(1:2:4,:)  ))
%! assert (isequal(  b(2:2:4,3), a(2:2:4,3)  ))
%! assert (isequal(  b(2:2:4,3), a(2:2:4,3)  ))

%!test
%! % 2D arrays
%! b = [1:4]; b = [b; 3*b; 5*b];
%! a = sym(b);
%! I = rand(size(b)) > 0.5;
%! assert (isequal (a(I), b(I)))
%! I = I(:);
%! assert (isequal (a(I), b(I)))
%! I = I';
%! assert (isequal (a(I), b(I)))
%! I = logical(zeros(size(b)));
%! assert (isequal (a(I), b(I)))

%!test
%! % 1D arrays, does right with despite warning
%! r = [1:6];
%! ar = sym(r);
%! c = r';
%! ac = sym(c);
%! Ir = rand(size(r)) > 0.5;
%! Ic = rand(size(c)) > 0.5;
%! assert (isequal (ar(Ir), r(Ir)))
%! assert (isequal (ac(Ic), c(Ic)))
%! assert (isequal (ar(Ic), r(Ic)))
%! assert (isequal (ac(Ir), c(Ir)))

%!test
%! % rccross tests
%! B = [1 2 3 4; 5 6 7 9; 10 11 12 13];
%! A = sym(B);
%! assert (isequal (A([1 3],[2 3]), B([1 3], [2 3])  ))
%! assert (isequal (A(1,[2 3]), B(1,[2 3])  ))
%! assert (isequal (A([1 2],4), B([1 2],4)  ))
%! assert (isequal (A([2 1],[4 2]), B([2 1],[4 2])  ))
%! assert (isequal (A([],[]), B([],[])  ))

%!error <integers>
%! % issue #445
%! A = sym([10 11]);
%! A(1.1)

%!error <integers>
%! % issue #445
%! A = sym([10 11]);
%! A(sym(4)/3)

%!error <integers>
%! % issue #445
%! A = sym([1 2; 3 4]);
%! A(1.1, 1)

%!error <integers>
%! % issue #445
%! A = sym([1 2; 3 4]);
%! A(1, sym(4)/3)

%!error <out of bound>
%! A = sym([1 2; 3 4]);
%! A(5)

%!shared a, b
%! a = [1 2 3 5; 4 5 6 9; 7 5 3 2];
%! b = sym (a);

%!test
%! c = true;
%! assert (isequal (a(c), b(c)))
%! c = false;
%! assert (isequal (a(c), b(c)))

%!test
%! c = [false true];
%! assert (isequal (a(c), b(c)))
%! d = c | true;
%! assert (isequal (a(d), b(d)))
%! d = c & false;
%! assert (isequal (a(d), b(d)))

%!test
%! c = [false true false true; true false true false; false true false true];
%! assert (isequal (a(c), b(c)))
%! d = c | true;
%! assert (isequal (a(d), b(d)))
%! d = c & false;
%! assert (isequal (a(d), b(d)))

%!test
%! c = [false true false true false];
%! assert (isequal (a(c), b(c)))
%! d = c | true;
%! assert (isequal (a(d), b(d)))
%! d = c & false;
%! assert (isequal (a(d), b(d)))

%!test
%! c = [false; true; false; true; false];
%! assert (isequal (a(c), b(c)))
%! d = c | true;
%! assert (isequal (a(d), b(d)))
%! d = c & false;
%! assert (isequal (a(d), b(d)))

%!test
%! c = [false true; false true; true false];
%! assert (isequal (a(c), b(c)))
%! d = c | true;
%! assert (isequal (a(d), b(d)))
%! d = c & false;
%! assert (isequal (a(d), b(d)))

%!shared

%!test
%! % Orientation of empty results of logical indexing on row or column vectors
%! r = [1:6];
%! c = r';
%! ar = sym(r);
%! ac = sym(c);
%! assert (isequal (ar(false), r(false)))
%! assert (isequal (ac(false), c(false)))
%! assert (isequal (ar(false (1, 6)), r(false (1, 6))))
%! assert (isequal (ac(false (1, 6)), c(false (1, 6))))
%! assert (isequal (ar(false (6, 1)), r(false (6, 1))))
%! assert (isequal (ac(false (6, 1)), c(false (6, 1))))

repo: octsympy