text:
function gpcf = gpcf_covar(varargin)
%   GPCF_COVAR creates a gpstuff covariance structure for the multivariate
%   Gaussian Process model based on the coregionalization models (LMC).
%   See Gelfand et al. (2004) in nonstationary multivariate process modelling
%   through spatially varying coregionalization. Sociedad de Estadistica e
%   Investigacion Operativa, Test.
%
%   The coregionalization matrix is a covariance matrix and its
%   parametrization is given via variances and the correlation matrix. The
%   correlation matrix is parametrizated with the transformation proposed by
%   Lewandowski et al. (2009) in generating correlation matrices based on
%   vines and extended onion method. Journal of multivariate analysis.
%
%   Description :
%
%   ─ GPCF = GPCF_COVAR('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%     creates correlation matrix structure in which
%     the named parameters have the specified values. Any
%     unspecified parameters are set to default values.
%
%   ─ GPCF = GPCF_COVAR(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%   ─ Parameters (of the gpstuff structure) for the covariance structure [default]
%     R_prior             - prior for correlation matrix  [prior_corrunif()]
%     V_prior             - prior for the variances. Must be a structure
%                           with each component being also a structure (the
%                           priors). Otherwise, just use prior_fixed.
%     corrFun             - structure containing other covariance
%                           functions. More specifically, each element of
%                           the structure must be a correlation function
%                           (σ²_f = 1). Each element of the structure can
%                           also be a combination of covariance functions,
%                           which will result to a correlation function.
%                           You always need to follow this structure when
%                           using that specific type of covariance function
%                           in the gpstructure (gp_set).
%
%     numberClass         - number of classes (GPs, categories, species, ...)
%                           being modelled (this is required to initialize
%                           other parameters)
%
%     degreeFreedom_prior - prior for degree of freedoms 'nu' [prior_corrunif()]
%     classVariables      - value defining which column of x is used to
%                           identify the class variables. If this is not
%                           given returns an error. They have to be given
%                           increasing downwards, otherwise you get an
%                           error.
%
%   »» 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_*
%
% ─────────────── 2018 Marcelo Hartmann.

% If you use this file. Learn to recognize its author. Copyright is good but
% misleading.

% 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.

  if nargin > 0 && ischar(varargin{1}) && ismember(varargin{1}, {'init' 'set'})
    % remove init and set if isfield(gpcf,'metric')
    varargin(1) = [];
  end

  ip = inputParser;
  ip.FunctionName = 'GPCF_COVAR';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('R', 0.01, @(x) isvector(x) && ~any(abs(x) > 1));
  ip.addParamValue('V', 1, @(x) isvector(x) && ~any(x < 0));
  ip.addParamValue('R_prior', prior_corrunif(), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('V_prior', prior_t('s2', 2), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('corrFun', {}, @(x) ~isempty(x) && iscell(x));
  ip.addParamValue('classVariables', [], @(x) ~isempty(x) && mod(x, 1) == 0 && x > 0);
  ip.addParamValue('numberClass', [], @(x) mod(x, 1) == 0 && x > 1);
  ip.addParamValue('aValue', 1, @(x) isreal(x) &&  x > 0);

  ip.parse(varargin{:});

  gpcf = ip.Results.gpcf;

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

  if init
      % Set the function handles to the subfunctions
      gpcf.fh.realtoz = @gpcf_covar_realtoz;
      gpcf.fh.ztoreal = @gpcf_covar_ztoreal;
      gpcf.fh.rhotoreal = @gpcf_covar_rhotoreal;
      gpcf.fh.realtorho = @gpcf_covar_realtorho;
      gpcf.fh.pak = @gpcf_covar_pak;
      gpcf.fh.unpak = @gpcf_covar_unpak;
      gpcf.fh.lp = @gpcf_covar_lp;
      gpcf.fh.lpg = @gpcf_covar_lpg;
      gpcf.fh.cfg = @gpcf_covar_cfg;
      gpcf.fh.cov = @gpcf_covar_cov;
      gpcf.fh.sigma = @gpcf_covar_sigma;
      gpcf.fh.trcov = @gpcf_covar_trcov;
      gpcf.fh.trvar = @gpcf_covar_trvar;
      gpcf.fh.recappend = @gpcf_covar_recappend;
  end

  % initialize class variables and selected variables
  if isempty(ip.Results.classVariables) || isempty(ip.Results.numberClass)
      error('class variables and number of classes must be given')
  end
  gpcf.classVariables = ip.Results.classVariables;
  gpcf.numberClass = ip.Results.numberClass;
  gpcf.vectSize = (gpcf.numberClass^2 - gpcf.numberClass)/2;
  gpcf.aValue = ip.Results.aValue;

  % create indexes which will easy our lives from now on
  seq = 1:gpcf.vectSize;
  i = ceil(0.5 + 0.5 * sqrt(1 + 8 * seq)); % ith column
  j = seq - (i - 2).*(i - 1)/2;            % jth column
  ind1 = (j - 1) * gpcf.numberClass + i;   % linear indexes for lower triangular
  ind2 = (i - 1) * gpcf.numberClass + j;   % linear indexes for upper triangular
  gpcf.index = [i; j; ind1; ind2; seq]';

  % Initialize correlations functions
  if init || ~ismember('corrFun', ip.UsingDefaults)
      gpcf.corrFun = ip.Results.corrFun;
  end

  if length(gpcf.corrFun) ~= gpcf.numberClass
      error('the number of correlation functions must match the number of classes')
  end

  % Initialize correlations parameters
  if init && ~ismember('R', ip.UsingDefaults)
      gpcf.R = ip.Results.R;
      if length(gpcf.R) ~= gpcf.vectSize
          error('length of the vector does not agree with the number of classes')
      end
  elseif init
      gpcf.R = (-0.99 + 1.98.*rand(1, gpcf.vectSize)) * 0.001;
  end

  % Initialize variance parameters
  if init && ~ismember('V', ip.UsingDefaults)
      gpcf.V = ip.Results.V;
      if length(gpcf.V) ~= gpcf.numberClass
          error('length of the variance vector does not agree with the number of classes')
      end
  elseif init
      % takes one if nothing is given.
      gpcf.V = ones(1, gpcf.numberClass);
  end

  % initialize prior structure
  if init
      gpcf.p = [];
  end

  % initialize correlations
  if init || ~ismember('R_prior', ip.UsingDefaults)
    gpcf.p.R = ip.Results.R_prior;
    if ~isempty(gpcf.p.R)
        % corrects the prior
        gpcf.p.R.numberClass = gpcf.numberClass;
        gpcf.p.R.vectSize = gpcf.vectSize;
        gpcf.p.R.aValue = gpcf.aValue;
        if gpcf.p.R.numberClass > gpcf.p.R.nu
            error('degrees of freedom must be greater than number of classes')
        end
    end
  end

  % initialize variances
  if init && ~ismember('V_prior', ip.UsingDefaults)
      gpcf.p.V = ip.Results.V_prior;
      if ~isempty(ip.Results.V_prior)
          gpcf.varfields = (fields(gpcf.p.V)');
          if length(structfun(@numel, gpcf.p.V)) ~= gpcf.numberClass
              error('prior missing for some variance')
          end
      end
  else
      var = [repmat('var', gpcf.numberClass , 1), num2str((1:gpcf.numberClass)')];
      for ind = 1:size(var, 1)
          var_aux = var(ind, :);
          var_aux = var_aux(~isspace(var_aux));
          Vpriors.(var_aux).p = ip.Results.V_prior;
          %Vpriors.(var(ind, :)).p = ip.Results.V_prior;
      end
      gpcf.p.V = Vpriors;
      gpcf.varfields = (fields(gpcf.p.V)');
  end

end

function y = gpcf_covar_realtoz(x, a, d)
% Description :
%   · Transforms the real line to the interval (-1, 1) using the shifted
%   logistic function
%
%   · rho = 2 ./ (1 + exp(-a*x)) - 1;
%   · x is the vector on the real line
%   · d equals to [], 1 or 2, and indicates whether the output corresponds
%     no derivatives, first or second derivative.

if isempty(d)
    % modified logistic distribution function
    y = 2 ./ (1 + exp(-a*x)) - 1;
elseif d == 1
    % repeated terms
    eminax = exp(-a*x);
    % 1st derivative, drho/dx
    y = 2 * a * (1 + eminax).^(-2) .* eminax;
elseif d == 2
    % repeated terms
    a2 = a ^ 2;
    eminax = exp(-a*x);
    % 2nd derivative, d2rho/dx2
    y = 4 * a2 * (1 + eminax).^(-3) .* eminax.^2 - 2 * a2 * (1 + eminax).^(-2) .* eminax;
else
    error('Something is wrong. Maybe d should be [], 1 or 2')
end
end

function y = gpcf_covar_ztoreal(x, a)
%  Description :
%   · Transforms the interval (-1, 1) to the real line, using the inverse
%   hyperbolic tangent function
%
%   · y = - 1/a * log(2/(x+1) - 1)
%   · x in (-1, 1)^(dim(x))

if any(abs(x) >= 1)
    error('domain error');
else
    y = -1/a * log(2./(x+1) - 1);
end
end

function y = gpcf_covar_rhotoreal(gpcf)
%  Description :
%   · Transforms the matrix R to the matrix Z (composed by the elements
%   gamma on the real line)

ind1 = gpcf.index(:, 3);
ind2 = gpcf.index(:, 4);

Z = zeros(gpcf.numberClass, gpcf.numberClass);
R = eye(gpcf.numberClass);
R([ind1, ind2]) = [gpcf.R'; gpcf.R'];

W = chol(R);
Z(1, 2:gpcf.numberClass) = W(1, 2:gpcf.numberClass);
for i1 = 2:(gpcf.numberClass-1)
    for j1 = (i1+1):gpcf.numberClass
        Z(i1, j1) = W(i1, j1) .* exp(-0.5*sum(log(1-Z(1:(i1-1), j1).^2)));
        % Z(i1, j1) = W(i1, j1) .* 1/sqrt(prod(1-Z(1:(i1-1), j1).^2));
    end
end

y = gpcf.fh.ztoreal(Z(ind2), gpcf.aValue);
end

function y = gpcf_covar_realtorho(gpcf, w, d)
%  Description :
%   · Transforms the vector gamma on the real line to the matrix R (in a vector);

%   · z = 2 ./ (1 + exp(-a*x)) - 1;
%   · x is the vector on the real line
%   · d equals to [], 1 or 2, and indicates whether the output corresponds
%     no derivatives, first or second derivative.

i = gpcf.index(:, 1);
j = gpcf.index(:, 2);
ind2 = gpcf.index(:, 4);
seq = gpcf.index(:, end);
i2 = gpcf.vectSize;

Z = zeros(gpcf.numberClass, gpcf.numberClass);
W = eye(gpcf.numberClass, gpcf.numberClass);

Z(ind2) = gpcf.fh.realtoz(w(1:i2), gpcf.aValue, []);
W(1, 2:gpcf.numberClass) = Z(1, 2:gpcf.numberClass);

for i1 = 2:gpcf.numberClass
    for j1 = i1:gpcf.numberClass
        ztmp = 0.5.*sum(log(1 - Z(1:(i1-1), j1).^2));
        W(i1, j1) = ((i1 == j1) + (i1 ~= j1)*Z(i1, j1)) * exp(ztmp);
        % ztmp = sqrt(1-Z(1:(i1-1), j1).^2);
        % W(i1, j1) = ((i1 == j1) + ~(i1 == j1)*Z(i1, j1)) * prod(ztmp);
    end
end

if isempty(d)
    % correlation matrix
    R = W'*W;
    y = R(ind2);

elseif d == 1
    % indexes
    ind3 = [j i seq];

    % derivatives of the cholesk decomposition w.r.t values on the real line
    dW = zeros(gpcf.numberClass, gpcf.numberClass, i2);
    dL = zeros(gpcf.numberClass, gpcf.numberClass, i2);
    z = Z(ind2);
    dz = gpcf.fh.realtoz(w, gpcf.aValue, 1);

    fz = sqrt(1 - z.^2);
    dfz = - 1./sqrt(1 - z.^2) .*  z .* dz;

    for k = 1:i2
        % take the column
        j1 = ind3(ind3(:, 3) == k, 2);

        % take the indexes for that column
        ind4 = ind3(ind3(:, 2) == j1, end);

        fzaux = fz(ind4);
        dfzaux = dfz(ind4);

        for i1 = 1:j1
            if i1 == 1
                dW(i1, j1, k) = (k == ind4(1))*dz(k);

            elseif i1 ~= j1
                kaux = ind4(i1);

                if kaux < k
                    continue
                else
                    ind5 = ind4(1:(i1-1));
                    fzaux2 = fz(ind5);
                    dfzaux2 = dfz(ind5);
                    e = (k == ind5);

                    if  kaux == k
                        %dW(i1, j1, k) = dz(k).*prod(fzaux2(~e));
                        dW(i1, j1, k) = dz(k).*exp(sum(log(fzaux2(~e))));
                    end

                    if kaux > k
                        %dW(i1, j1, k) = z(kaux).*dfzaux2(e)*prod(fzaux2(~e));
                        dW(i1, j1, k) = z(kaux).*dfzaux2(e)*exp(sum(log(fzaux2(~e))));
                    end
                end

            elseif i1 == j1
                e = (k == ind4);
                %dW(i1, j1, k) = dfzaux(e)*prod(fzaux(~e));
                dW(i1, j1, k) = dfzaux(e)*exp(sum(log(fzaux(~e))));

            end
        end
        dL(:, :, k) = dW(:, :, k)';

    end
    y = dL;

elseif d == 2
    % lower cholesk decompostion of R
    y = W';

end
end

function [w, s, h] = gpcf_covar_pak(gpcf, w)
%  GPCF_COVAR_PAK Combine GP correlation matrix values into one vector
%   keeping the order of input below.
%
%  Description :
%   · W = GPCF_COVAR_PAK(GPCF) takes the non-diagonal elements of
%     the correlation matrix and put them into a single row vector w
%     This is a mandatory subfunction used for example
%     in energy and gradient computations.
%
%   · w = [rho_(1, 2), ..., rho_(1, J), rho_(2, 3), ..., rho_(2, J), ..., rho_((J-1), J)]
%
%  * See also
%    GPCF_covar_UNPAK

  w = []; s = {}; h = [];

  % pak parameters of correlation functions
  for i = 1:gpcf.numberClass
      [wh, sh, hh] = gpcf.corrFun{i}.fh.pak(gpcf.corrFun{i});
      w = [w wh];
      s = [s; sh];
      h = [h 1+hh];
  end

  % pak correlations
  if ~isempty(gpcf.p.R)
      i = gpcf.index(:, 1);
      j = gpcf.index(:, 2);
      waux = gpcf.fh.rhotoreal(gpcf);
      w = [w waux'];
      S = [repmat('realcorr.', gpcf.vectSize, 1), num2str(j), num2str(i)];
      s = [s; cellstr(S)];
      h = [h 1];

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

  % pak variances
  if ~isempty(gpcf.p.V)
    w = [w log(gpcf.V)];
    S = [repmat('log var.', gpcf.numberClass, 1), num2str((1:gpcf.numberClass)')];
    s = [s; cellstr(S)];
    h = [h 1];

    % pak hyperparameters of V
    for i = 1:gpcf.numberClass
        cfvf = gpcf.varfields{i};
        [wh, sh, hh] = gpcf.p.V.(cfvf).p.fh.pak(gpcf.p.V.(cfvf).p);
        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_covar_unpak(gpcf, w)
%  GPCF_covar_UNPAK Sets the covariance function parameters into the
%   structure
%
%  Description :
%   · [GPCF, W] = GPCF_COVAR_UNPAK(GPCF, W) takes a covariance
%     function structure GPCF and a hyperparameter vector W,
%     and returns a covariance function structure identical to
%     the input, except that the covariance hyperparameters 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.
%
%  * See also
%    GPCF_covar_PAK

  % unpak parameters of correlations functions
  for i = 1:gpcf.numberClass
      cf = gpcf.corrFun{i};
      [gpcf.corrFun{i}, w] = cf.fh.unpak(cf, w);
  end

  gpp = gpcf.p;

  % unpak for Σ
  if ~ isempty(gpp.R)
      gpcf.R = gpcf.fh.realtorho(gpcf, w(1:gpcf.vectSize), [])';
      w = w((gpcf.vectSize + 1):end);

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

  if ~ isempty(gpp.V)
      i3 = gpcf.numberClass;
      gpcf.V = exp(w(1:i3));
      w = w((i3 + 1):end);

      % unpak hyperparameters of V
      % loop for the hyperpriors
      for i = 1:gpcf.numberClass
          cfvf = gpcf.varfields{i};
          [p, w] = gpcf.p.V.(cfvf).p.fh.unpak(gpcf.p.V.(cfvf).p, w);
          gpcf.p.V.(cfvf).p = p;
      end
  end
end

function lp = gpcf_covar_lp(gpcf)
%  GPCF_COVAR_LP Evaluate the log prior of covariance function parameters
%
%  Description :
%   · LP = GPCF_COVAR_LP(GPCF, X, T) takes a correlation function
%     structure GPCF and evaluates log-prior(R)
%
%  * See also
%    GPCF_COVAR_PAK, GPCF_covar_UNPAK, GPCF_covar_LPG, GP_E

  lp = 0;

  % add log-prior for correlations functions parameters
  for i = 1:gpcf.numberClass
      lp = lp + gpcf.corrFun{i}.fh.lp(gpcf.corrFun{i});
  end

  gpp = gpcf.p;

  % log-prior for correlations
  if ~isempty(gpp.R)
     % mapping to the real line
     x  = gpcf.fh.rhotoreal(gpcf);
     z  = gpcf.fh.realtoz(x, gpcf.aValue, []);
     dz = gpcf.fh.realtoz(x, gpcf.aValue, 1);

     % add log-prior on rho parameterization
     lp = lp + gpp.R.fh.lp(gpcf.R, gpp.R);

     % add log(|det J|) on the real line parameterization. Equation 11 in
     % Lewandowski et al (2009).
     pwr = gpcf.numberClass - gpcf.index(:, 2) - 1;
     lp = lp + 0.5 .* sum(pwr .* log(1 - z.^2)) + sum(log(dz));

  end

  % log-prior for variances
  if ~isempty(gpp.V)
      for i = 1:gpcf.numberClass
          cfvf = gpcf.varfields{i};

          lp = lp + gpp.V.(cfvf).p.fh.lp(gpcf.V(i), ...
              gpp.V.(cfvf).p) + log(gpcf.V(i));
      end
  end
end

function lpg = gpcf_covar_lpg(gpcf)
%  GPCF_COVAR_LPG Evaluate gradient of the log prior with respect
%   to the parameters.
%
%  Description :
%   · LPG = GPCF_COVAR_LPG(GPCF) takes a covariance function
%     structure GPCF and returns LPG = d log (p(th))/dth, where th
%     is the parametric vector of correlations and variances.
%     This is a mandatory subfunction used for example in gradient
%     computations.
%
%  * See also
%    GPCF_covar_PAK, GPCF_covar_UNPAK, GPCF_covar_LP, GP_G

  lpg = [];

  % stack grad vector of the correlation functions parameters
  for i = 1:gpcf.numberClass
      lpg = [lpg gpcf.corrFun{i}.fh.lpg(gpcf.corrFun{i})];
  end

  gpp = gpcf.p;

  % for correlations
  if ~isempty(gpcf.p.R)
      % mapping to the real line
      real = gpcf.fh.rhotoreal(gpcf);
      z   = gpcf.fh.realtoz(real, gpcf.aValue, [])';
      dz  = gpcf.fh.realtoz(real, gpcf.aValue, 1)';
      d2z = gpcf.fh.realtoz(real, gpcf.aValue, 2)';
      pwr = (gpcf.numberClass - gpcf.index(:, 2) - 1)';

      % cholesk decomposition
      L = gpcf.fh.realtorho(gpcf, real, 2);

      % transforming to first derivative of L w.r.t gammas
      dL = gpcf.fh.realtorho(gpcf, real, 1);

      b = zeros(gpcf.vectSize, gpcf.vectSize);
      for k = 1:size(dL, 3)
          dRaux = dL(:, :, k)*L' + L*dL(:, :, k)';
          b(:, k) = dRaux(gpcf.index(:, 4));
      end

      % dlogp(R(gamma))/dgamma
      dlpdgm = gpp.R.fh.lpg(gpcf.R, gpp.R) * b;

      lpgs = dlpdgm - pwr.*(z.*dz./(1 - z.^2)) + d2z./dz;
      lpg = [lpg lpgs];
  end

  % for variances
  if ~isempty(gpp.V)
      % building the grad vector on the real line parameterization
      for i = 1:gpcf.numberClass
          cfvf = gpcf.varfields{i};
          lpgs = gpp.V.(cfvf).p.fh.lpg(gpcf.V(i), gpp.V.(cfvf).p);
          lpg = [lpg lpgs(1).*gpcf.V(i)+1 lpgs(2:end)];
      end
  end
end

function C = gpcf_covar_sigma(gpcf, type)
%  GP_COVAR_COV  Evaluate the covariance matrix Σ of the linear model of
%  coregionalization and its Cholesky decomposition.
%
%  Description :
%   · C = GP_COVAR_COV(GP, i1) takes the gpcf structure and
%     returns the covariance matrix Σ and its Cholesky decomposition.
%     Every element (i, j) of Σ contains the covariance (parameterised by
%     correlation and variances) between class i and class j. This is a
%     mandatory subfunction.
%
%  * See also
%    GPCF_COVAR_TRCOV, GPCF_COVAR_TRVAR, GP_COV, GP_TRCOV

% linear indexes
ind1 = gpcf.index(:, 3);
ind2 = gpcf.index(:, 4);

% corr matrix
R = eye(gpcf.numberClass);

% filling elements in lower and upper part
R([ind1, ind2]) = [gpcf.R, gpcf.R];

% matrix Σ
Sig = bsxfun(@times, bsxfun(@times, sqrt(gpcf.V)', R), sqrt(gpcf.V));
% Sig = S * R * S;

% other options
if nargin == 2
    switch type
        case 'cov'
            % empty matrix
            C{1} = Sig;

            % calcule cholesk decomposition of Sigma
            [L, notpositivedefinite] = chol(Sig, 'lower');
            % test whether is positive-definite or not
            if notpositivedefinite
                L = NaN(gpcf.numberClass, gpcf.numberClass);
                L(ind2) = 0;
            end

            % take cholesky decomposition
            C{2} = L;

        case 'corr'
            x = gpcf.fh.rhotoreal(gpcf);
            L = gpcf.fh.realtorho(gpcf, x, 2);

            C{1} = L*L';
            C{2} = L;
    end
else
    C = Sig;

end
end

function C = gpcf_covar_cov(gpcf, x1, x2)
%  GP_COVAR_COV  Evaluate the covariance matrix between two input vectors
%
%  Description
%   · C = GP_COVAR_COV(GP, TX, X) takes in correlation structure
%     and two matrixes TX and X that contain input vectors to GP.
%     Returns covariance matrix C.
%     Every element ij of C contains correlation 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_COVAR_TRCOV, GPCF_COVAR_TRVAR, GP_COV, GP_TRCOV

 if isempty(x2)
     x2 = x1;
 end

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

 xClass1 = x1(:, gpcf.classVariables);
 xClass2 = x2(:, gpcf.classVariables);

 if ~issorted(xClass1) || ~issorted(xClass2)
     error('you need to give the class variable increasing downwards');
 end

 if max(xClass1) > gpcf.numberClass || max(xClass2) > gpcf.numberClass
     error('more classes than the given number of classes');
 end

 n1 = size(xClass1, 1);
 n2 = size(xClass2, 1);

 % information of the classes in the data
 ind1 = unique(xClass1)'; nind1 = size(ind1, 2);
 ind2 = unique(xClass2)'; nind2 = size(ind2, 2);

 % positions of observations in each vector
 nb1 = find(diff([-inf xClass1' inf]));
 nb2 = find(diff([-inf xClass2' inf]));

 % chol(Σ)
 M = gpcf.fh.sigma(gpcf, 'cov');
 L = M{2};

 % M = gpcf.fh.sigma(gpcf, 'corr');
 % L = bsxfun(@times, sqrt(gpcf.V)', M{2});

 C = zeros(n1, n2);
 for k = 1:gpcf.numberClass

     % kth column of the cholesk decomposition of Σ = SRS
     T = L(:, k) * L(:, k)';

     for i = 1:nind1
         for j = 1:nind2
             inb1 = nb1(i) : nb1(i + 1) - 1;
             inb2 = nb2(j) : nb2(j + 1) - 1;
             K(inb1, inb2) = T(ind1(i), ind2(j));
         end
     end

     cf = gpcf.corrFun{k};
     C = C + cf.fh.cov(cf, x1, x2) .* K;
end

 C(abs(C) < eps) = 0;
end

function dCRV = gpcf_covar_cfg(gpcf, x, x2, mask, i1)
%  GPCF_COVAR_CFG  Evaluate the derivarive of the covar matrix
%
%  Description :
%   · dRff = GPCF_COVAR_CFG(GPCF, X) takes a
%     covariance structure GPCF, a matrix of inputs
%     vectors and returns dSRS/drho and/or dSRS/dV, the gradients of
%     SRS matrix
%
%   · dRff = GPCF_COVAR_CFG(GPCF, X, X2) takes a
%     covariance structure GPCF, a matrix X of input
%     vectors and returns  dSRS/drho and/or dSRS/dV, the gradients of
%     SRS matrix with respect to rho (cell array with matrix
%     elements). This subfunction is needed when using sparse
%     approximations (e.g. FIC).
%
%  * See also
%    GPCF_COVAR_PAK, GPCF_COVAR_UNPAK, GPCF_COVAR_LP, GP_G

dCRV = {};

% covariance matrix Σ and chol(Σ)
M = gpcf.fh.sigma(gpcf, 'corr'); % R = M{1};
L = M{2};

% covariance matrix
%Sig = sVL * sVL';

% \sqrt(diag(var)) times L
sVL = bsxfun(@times, sqrt(gpcf.V)', L);

% take class variables
xC = x(:, gpcf.classVariables);

% positions of the observations
nb = diff([0 xC' xC(end) + 1]) .* (1:(length(xC) + 1));
nb = nb(nb ~= 0);

% number of observation in each class
diffnb = diff(nb);

% number of classes
nc = gpcf.numberClass;

% number of observations
n = length(xC);

% use savememory option
savememory = nargin == 5;

if nargin == 2 || (isempty(x2) && isempty(mask))

    % Evaluate the values (dSRS/dgamma_12, dSRS/dgamma_13, ..., dSRS/dgamma_(J-1)J
    % and/or dSRS/dV_1 ...  dSRS/dV_J

    % build the big Tj's
    T = cell(1, gpcf.numberClass);

    for  k = 1:gpcf.numberClass
        Taux = sVL(:, k) * sVL(:, k)';

        for i = 1:gpcf.numberClass
            inb1 = nb(i) : nb(i + 1) - 1;
            T{k}(inb1, inb1) = Taux(i, i);

            for j = (i + 1):gpcf.numberClass
                inb2 = nb(j):nb(j + 1) - 1;

                T{k}(inb1, inb2) = Taux(i, j);
                T{k}(inb2, inb1) = Taux(j, i);
            end
        end
    end

    % correlations matrices
    dC = {};

    for i1 = 1:nc
        dKK = gpcf.corrFun{i1}.fh.cfg(gpcf.corrFun{i1}, x);

        if ~isempty(dKK)
            % number of derivatives of the correlation function
            i2 = length(dKK);

            % derivatives of the specific correlation function
            for i3 = 1:i2
                dC{end + 1} = dKK{i3} .* T{i1};
            end
        end
    end

    % correlations and variances in Σ = sqrt(S) L L sqrt(S)
    sVdL = [];
    dsVL = [];

    if ~isempty(gpcf.p.R) || ~isempty(gpcf.p.V)

        % for gammas in Σ
        if ~isempty(gpcf.p.R)
            % transforming to the real line
            real = gpcf.fh.rhotoreal(gpcf);

            % transforming to first derivative of L w.r.t gammas
            dL = gpcf.fh.realtorho(gpcf, real, 1);

            % variances times first derivatives w.r.t correlations
            % \sqrt(diag(var)) (dLj/dΘ)
            sVdL = bsxfun(@times, sqrt(gpcf.V)', dL);
        end

        % for variances in Σ
        if ~isempty(gpcf.p.V)
            dV = 0.5 .* sqrt(gpcf.V);
            dsVLm = bsxfun(@times, dV', L); % take the lines ...
            dsVL = zeros(nc, nc, nc);

            for k = 1:nc
                dsVL(k, :, k) = dsVLm(k, :);
            end
        end

        % derivatives w.r.t to correlations or/and variances
        dSig = cat(3, sVdL, dsVL);

        if ~isempty(dSig)

            nSig = length(dSig);
            dL = cell(1, nSig);
            % dCorr = zeros(n, n);

            for k = 1:nSig
                % sum_j R_j kron (daj/dΘ * aj' + aj * daj/dΘ')
                dL{k} = zeros(n, n);

                for j = 1:nc
                    aj  = repelem(sVL(:, j).', diffnb);
                    daj = repelem(dSig(:, j, k), diffnb);
                    dd = daj * aj;

                    cf = gpcf.corrFun{j};
                    dCorr = cf.fh.trcov(cf, x) .* (dd + dd.');
                    dL{k} = dL{k} + dCorr;
                end
            end
        end
    end

    % return matrices
    dCRV = [dC dL];

% related to sparse approximations
elseif nargin == 3 || isempty(mask)
    error('nargin == 3 || isempty(mask) not implemented yet');

elseif nargin == 4 || nargin == 5
    error('nargin == 4 || nargin == 5 not implemented yet');

end
end

function C = gpcf_covar_trcov(gpcf, x)
%  GP_COVAR_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description :
%   · C = GP_COVAR_TRCOV(GP, TX) takes in correlation structure
%     matrix TX that contains training input vectors.
%     Returns covariance matrix C. Every element ij of C contains
%     the correlation between inputs i and j in TX.
%     This is a mandatory subfunction used for example in
%     prediction and energy computations.
%
%  * See also
%    GPCF_COVAR_COV, GPCF_COVAR_TRVAR, GP_COV, GP_TRCOV
%
% HERE IS TO TAKE ADVANTAGE OF THE SYMMETRIC MATRIX

% takes the class variables
xC = x(:, gpcf.classVariables);

% checking
if ~issorted(xC)
    error('class variable should increasing downwards');
end

% number of observations
n = size(xC, 1);

% getting the information of the classes in the data
a = unique(xC); % na = size(a, 1);

% checking
if max(a) > gpcf.numberClass
% if na ~= gpcf.numberClass || max(a) ~= gpcf.numberClass
    error('more or less classes than given to the model');
end

% getting the information of the classes in the data
ind = unique(xC)'; nind = size(ind, 2);

% number of observations in each class
nb = find(diff([-inf xC' inf]));

% Try to use the C implementation
% C = trcov_corr(gpcf, x1);

C = NaN;

if isnan(C)
    % take chol(Σ)
    M = gpcf.fh.sigma(gpcf, 'cov');
    L = M{2};

    % full matrices
    C = zeros(n, n);  K = zeros(n, n);

    for k = 1:gpcf.numberClass
        % T_k = a_k * a_k' matrices,
        T =  L(:, k) * L(:, k)';

        for i = 1:nind
            inb1 = nb(i) : nb(i + 1) - 1;
            K(inb1, inb1) = T(ind(i), ind(i));

            for j = (i + 1):nind
                inb2 = nb(j) : nb(j + 1) - 1;

                K(inb1, inb2) = T(ind(i), ind(j));
                K(inb2, inb1) = T(ind(j), ind(i));
            end
        end

        % corrFun kronecker T_k
        cf = gpcf.corrFun{k};
        C = C + cf.fh.trcov(cf, x) .* K;
    end

    C(abs(C) < eps) = 0;
end
end

function C = gpcf_covar_trvar(gpcf, x)
%  GP_COVAR_VECTOR  Evaluate training variance vector
%
%  Description:
%   · C = GP_COVAR_TRVAR(GPCF, TX) takes in correlation structure
%     of and matrix TX that contains training inputs.
%     Returns correlation vector C, which are ones.
%     Every element i of C contains the correlation of the input in TX.
%     This is a mandatory subfunction used for example in prediction and
%     energy computations.
%
%  * See also:
%    GPCF_COVAR_COV, GP_COV, GP_TRCOV

 % take indicator class
 x = x(:, gpcf.classVariables);

 % information of the classes in the data
 a = unique(x);

 % locate each class
 nb = find(diff([-inf x' inf]));

 % number of observation in each class
 diffnb = diff(nb);

 C = repelem(gpcf.V(a), diffnb)';
 end

function reccf = gpcf_covar_recappend(reccf, ri, gpcf)
%  RECAPPEND  Record append
%
%  Description:
%   · RECCF = GPCF_COVAR_RECAPPEND(RECCF, RI, GPCF) takes a
%     correlation structure record structure RECCF, record index RI
%     and correlation 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_covar';

      % initialize parameters for correlations functions
      for i = 1:ri.numberClass
          reccf.corrFun{i} = ri.corrFun{i}.fh.recappend([], ri.corrFun{i});
      end

      % Initialize parameters
      reccf.R = [];
      reccf.V = [];

      % Set the function handles
      reccf.fh.realtoz = @gpcf_covar_realtoz;
      reccf.fh.ztoreal = @gpcf_covar_ztoreal;
      reccf.fh.rhotoreal = @gpcf_covar_rhotoreal;
      reccf.fh.realtorho = @gpcf_covar_realtorho;
      reccf.fh.pak = @gpcf_covar_pak;
      reccf.fh.unpak = @gpcf_covar_unpak;
      reccf.fh.lp = @gpcf_covar_lp;
      reccf.fh.lpg = @gpcf_covar_lpg;
      reccf.fh.cfg = @gpcf_covar_cfg;
      reccf.fh.cov = @gpcf_covar_cov;
      reccf.fh.sigma = @gpcf_covar_sigma;
      reccf.fh.trcov = @gpcf_covar_trcov;
      reccf.fh.trvar = @gpcf_covar_trvar;
      reccf.fh.recappend = @gpcf_covar_recappend;
      reccf.p = [];
      reccf.p.R = [];
      reccf.p.V = [];

      if isfield(ri.p, 'R') && ~isempty(ri.p.R)
          reccf.p.R = ri.p.R;
      end

      if isfield(ri.p, 'V') && ~isempty(ri.p.V)
          for i = 1:ri.numberClass
              reccf.p.V.(ri.varfields{i}).p = ri.p.V.(ri.varfields{i}).p;
          end
      end

      if isfield(ri, 'classVariables')
          reccf.classVariables = ri.classVariables;
          reccf.numberClass = ri.numberClass;
          reccf.vectSize = ri.vectSize;
          reccf.aValue = ri.aValue;
          reccf.varfields = ri.varfields;
          reccf.index = ri.index;
      end

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

      if ~ isfield(gpcf,'metric')
          % record correlation function elements
          for i = 1:gpcf.numberClass
              cf = gpcf.corrFun{i};
              reccf.corrFun{i} = cf.fh.recappend(reccf.corrFun{i}, ri, cf);
          end

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

          % record V elements
          reccf.V(ri, :) = gpcf.V;
          if isfield(gpp, 'V') && ~isempty(gpp.V)
              for i = 1:gpcf.numberClass
                  cfvf = gpcf.varfields{i};
                  reccf.p.V.(cfvf).p = ...
                      gpp.V.(cfvf).p.fh.recappend(reccf.p.V.(cfvf).p, ...
                      ri, gpcf.p.V.(cfvf).p);
              end
          end
      end

  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
classdef dme_line3p < mp.dm_element
%MP.DME_LINE3P  MATPOWER data model class for 3-phase line data

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        fbus    %% bus index vector for "from" port (port 1) (all lines)
        tbus    %% bus index vector for "to" port (port 2) (all lines)
        freq    %% system frequency, in Hz
        lc      %% index into lc_tab for lines that are on
        len     %% length for lines that are on
        lc_tab  %% line construction table
        ys      %% cell array of 3x3 series admittance matrices for lc rows
        yc      %% cell array of 3x3 shunt admittance matrices for lc rows
    end     %% properties

    methods
        function name = name(obj)
            name = 'line3p';
        end

        function label = label(obj)
            label = '3-ph Line';
        end

        function label = labels(obj)
            label = '3-ph Lines';
        end

        function name = cxn_type(obj)
            name = 'bus3p';
        end

        function name = cxn_idx_prop(obj)
            name = {'fbus', 'tbus'};
        end

        function names = main_table_var_names(obj)
            names = horzcat( main_table_var_names@mp.dm_element(obj), ...
                {'bus_fr', 'bus_to', 'lc', 'len', ...
                 'pl1_fr', 'ql1_fr', 'pl2_fr', 'ql2_fr', 'pl3_fr', 'ql3_fr', ...
                 'pl1_to', 'ql1_to', 'pl2_to', 'ql2_to', 'pl3_to', 'ql3_to' ...
                 });
        end

%         function vars = export_vars(obj)
%             vars = {'pl1_fr', 'ql1_fr', ...
%                     'pl2_fr', 'ql2_fr', ...
%                     'pl3_fr', 'ql3_fr', ...
%                     'pl1_to', 'ql1_to', ...
%                     'pl2_to', 'ql2_to', ...
%                     'pl3_to', 'ql3_to' };
%         end

        function obj = initialize(obj, dm)
            initialize@mp.dm_element(obj, dm); %% call parent

            %% get bus mapping info
            b2i = dm.elements.bus3p.ID2i;   %% bus num to idx mapping

            %% set bus index vectors for port connectivity
            obj.fbus = b2i(obj.tab.bus_fr);
            obj.tbus = b2i(obj.tab.bus_to);
        end

        function obj = update_status(obj, dm)
            %% get bus status info
            bs = dm.elements.bus3p.tab.status;  %% bus status

            %% update status of branches connected to isolated/offline buses
            obj.tab.status = obj.tab.status & bs(obj.fbus) & ...
                                              bs(obj.tbus);

            %% call parent to fill in on/off
            update_status@mp.dm_element(obj, dm);
        end

        function obj = build_params(obj, dm)
            nlc = size(obj.lc, 1);
            obj.ys = zeros(nlc, 6);
            obj.yc = zeros(nlc, 6);
            obj.lc = obj.tab.lc(obj.on);
            obj.len = obj.tab.len(obj.on);

            %% build Ys and Yc for relevant lines
            idx = unique(obj.lc);
            rr = obj.lc_tab.r(idx, :);
            xx = obj.lc_tab.x(idx, :);
            cc = obj.lc_tab.c(idx, :);
%             rr = obj.lc_tab{idx, 2:7};
%             xx = obj.lc_tab{idx, 8:13};
%             cc = obj.lc_tab{idx, 14:19};
            for k = 1:length(idx)
                R = obj.vec2symmat(rr(k, :));
                X = obj.vec2symmat(xx(k, :));
                C = obj.vec2symmat(cc(k, :));
                Ys = inv(R + 1j * X);
                Yc = 1j * 2*pi * obj.freq * 1e-9 * C;
                obj.ys(k, :) = obj.symmat2vec(Ys);
                obj.yc(k, :) = obj.symmat2vec(Yc);
            end
        end

        function M = vec2symmat(obj, v)
            % makes a symmetric matrix from a vector of 6 values
            M = [v(1) v(2) v(3);
                 v(2) v(4) v(5);
                 v(3) v(5) v(6) ];
        end

        function v = symmat2vec(obj, M)
            % extracts a vector of 6 values from a matrix assumed to be symmetric
            v = [M(1, :) M(2,2:3) M(3,3)];
        end

        function obj = pretty_print(obj, dm, section, out_e, mpopt, fd, pp_args)
            switch section
                case 'det'
                    %% compute currents/powers for pp_args
                    s_fr = [  obj.tab.pl1_fr + 1j * obj.tab.ql1_fr ...
                            obj.tab.pl2_fr + 1j * obj.tab.ql2_fr ...
                            obj.tab.pl3_fr + 1j * obj.tab.ql3_fr    ];
                    s_to = [  obj.tab.pl1_to + 1j * obj.tab.ql1_to ...
                            obj.tab.pl2_to + 1j * obj.tab.ql2_to ...
                            obj.tab.pl3_to + 1j * obj.tab.ql3_to    ];
                    t = dm.elements.bus3p.tab;    %% bus3p table
                    vm = [ t.vm1 t.vm2 t.vm3 ] .* (t.base_kv/sqrt(3) * [1 1 1]);
                    va = [ t.va1 t.va2 t.va3 ];
                    v_ = vm .* exp(1j * va * pi/180);
                    i_fr = conj( s_fr ./ v_(obj.fbus, :));
                    i_to = conj( s_to ./ v_(obj.tbus, :));
                    c_fr = struct('cm', abs(i_fr), 'ca', angle(i_fr) * 180/pi);
                    c_to = struct('cm', abs(i_to), 'ca', angle(i_to) * 180/pi);

                    pp_args.line3p = {'c', 'f', abs(i_fr), angle(i_fr) * 180/pi};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);

                    pp_args.line3p = {'c', 't', abs(i_to), angle(i_to) * 180/pi};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);

                    pp_args.line3p = {'s', 'f', real(s_fr), imag(s_fr)};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);

                    pp_args.line3p = {'s', 't', real(s_to), imag(s_to)};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);
                otherwise
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);
            end
        end

        function TorF = pp_have_section_sum(obj, mpopt, pp_args)
            TorF = true;
        end

        function obj = pp_data_sum(obj, dm, rows, out_e, mpopt, fd, pp_args)
            %% call parent
            pp_data_sum@mp.dm_element(obj, dm, rows, out_e, mpopt, fd, pp_args);

            %% print generation summary
            t = obj.tab;
            ploss = [ t.pl1_fr t.pl2_fr t.pl3_fr ] + ...
                    [ t.pl1_to t.pl2_to t.pl3_to ];
            qloss = [ t.ql1_fr t.ql2_fr t.ql3_fr ] + ...
                    [ t.ql1_to t.ql2_to t.ql3_to ];
            fprintf(fd, '  %-29s %12.1f kW %12.1f kVAr\n', 'Total 3-ph line loss', ...
                sum(sum(ploss(obj.on, :))), sum(sum(qloss(obj.on, :))) );
        end

        function TorF = pp_have_section_det(obj, mpopt, pp_args)
            TorF = true;
        end

        function h = pp_get_headers_det(obj, dm, out_e, mpopt, pp_args)
            cs = pp_args.line3p{1};
            ft = pp_args.line3p{2};
            if cs == 'c' && ft == 'f'
                h1 = pp_get_headers_det@mp.dm_element(obj, dm, out_e, mpopt, pp_args);
            else
                h1 = {};
            end
            if cs == 'c'
                if ft == 'f'
                    h2 = {'-->  Current Injections at "From" Bus' };
                else    %% ft == 't'
                    h2 = {'', ...
                        '<--  Current Injections at "To" Bus' };
                end
            else    %% cs == 's'
                if ft == 'f'
                    h2 = {'', ...
                        '-->  Power Injections at "From" Bus' };
                else    %% ft == 't'
                    h2 = {'', ...
                        '<--  Power Injections at "To" Bus' };
                end
            end
            switch cs
                case 'c'
                    h = [ h1 h2 ...
                        {   '  3-ph    3-ph Bus  3-ph Bus          Phase A Current  Phase B Current  Phase C Current', ...
                            'Line ID    From ID   To ID    Status   (A)    (deg)     (A)    (deg)     (A)    (deg)', ...
                            '--------  --------  --------  ------  ------  ------   ------  ------   ------  ------' } ];
                    %%       1234567 123456789 123456789 -----1 123456.89 12345.7 12345.78 12345.7 12345.78 12345.7
                case 's'
                    h = [ h1 h2 ...
                        {   '  3-ph    3-ph Bus  3-ph Bus          Phase A Power    Phase B Power    Phase C Power', ...
                            'Line ID    From ID   To ID    Status   (kW)   (kVAr)    (kW)   (kVAr)    (kW)   (kVAr)', ...
                            '--------  --------  --------  ------  ------  ------   ------  ------   ------  ------' } ];
                    %%       1234567 123456789 123456789 -----1 1234567.9 12345.7 123456.8 12345.7 123456.8 12345.7
            end
        end

        function str = pp_data_row_det(obj, dm, k, out_e, mpopt, fd, pp_args)
            switch pp_args.line3p{1}    %% cs
                case 'c'
                    cm = pp_args.line3p{3}(k, :);
                    ca = pp_args.line3p{4}(k, :);
                    str = sprintf('%7d %9d %9d %6d %9.2f %7.1f %8.2f %7.1f %8.2f %7.1f', ...
                        obj.tab.uid(k), obj.tab.bus_fr(k), obj.tab.bus_to(k), ...
                        obj.tab.status(k), ...
                        cm(:, 1), ca(:, 1), ...
                        cm(:, 2), ca(:, 2), ...
                        cm(:, 3), ca(:, 3) );
                case 's'
                    p = pp_args.line3p{3}(k, :);
                    q = pp_args.line3p{4}(k, :);
                    str = sprintf('%7d %9d %9d %6d %9.1f %7.1f %8.1f %7.1f %8.1f %7.1f', ...
                        obj.tab.uid(k), obj.tab.bus_fr(k), obj.tab.bus_to(k), ...
                        obj.tab.status(k), ...
                        p(:, 1), q(:, 1), ...
                        p(:, 2), q(:, 2), ...
                        p(:, 3), q(:, 3) );
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [data, w] = normalize_data(data, K, scale, w)
[m,n] = size(data.A);

MIN_SCALE = 1e-3;
MAX_SCALE = 1e3;
minRowScale = MIN_SCALE * sqrt(n);
maxRowScale = MAX_SCALE * sqrt(n);
minColScale = MIN_SCALE * sqrt(m);
maxColScale = MAX_SCALE * sqrt(m);

D = ones(m,1);
E = ones(n,1);
NN = 1; % NN = 1, other choices bad
for j=1:NN
    %% D scale:
    Dt = twonorms(data.A(1:K.f,:)')';
    idx = K.f;
    Dt = [Dt;twonorms(data.A(idx+1:idx+K.l,:)')'];
    idx = idx + K.l;
    for i=1:length(K.q)
        if (K.q(i) > 0)
            nmA = mean(twonorms(data.A(idx+1:idx+K.q(i),:)'));
            Dt = [Dt;nmA*ones(K.q(i),1)];
            idx = idx + K.q(i);
        end
    end
    for i=1:length(K.s)
        if (K.s(i) > 0)
            nmA = mean(twonorms(data.A(idx+1:idx+get_sd_cone_size(K.s(i)),:)'));
            Dt = [Dt;nmA*ones(get_sd_cone_size(K.s(i)),1)];
            idx = idx + get_sd_cone_size(K.s(i));
        end
    end
    for i=1:K.ep
        nmA = mean(twonorms(data.A(idx+1:idx+3,:)'));
        Dt = [Dt;nmA*ones(3,1)];
        idx = idx + 3;
    end
    for i=1:K.ed
        nmA = mean(twonorms(data.A(idx+1:idx+3,:)'));
        Dt = [Dt;nmA*ones(3,1)];
        idx = idx + 3;
    end
    for i=1:length(K.p)
        nmA = mean(twonorms(data.A(idx+1:idx+3,:)'));
        Dt = [Dt;nmA*ones(3,1)];
        idx = idx + 3;
    end

    Dt(Dt < minRowScale) = 1;
    Dt(Dt > maxRowScale) = maxRowScale;
    data.A = sparse(diag(1./Dt))*data.A;

    %% E Scale
    Et = twonorms(data.A)';
    Et(Et < minColScale) = 1;
    Et(Et > maxColScale) = maxColScale;
    data.A = data.A*sparse(diag(1./Et));

    %%
    D = D.*Dt;
    E = E.*Et;
end

nmrowA = mean(twonorms(data.A'));
nmcolA = mean(twonorms(data.A));

data.A = data.A*scale;

data.b = data.b./D;
sc_b = nmcolA/ max(norm(data.b), MIN_SCALE);
data.b = data.b * sc_b * scale;

data.c = data.c./E;
sc_c = nmrowA/max(norm(data.c), MIN_SCALE);
data.c = data.c * sc_c * scale;

w.D = D;
w.E = E;
w.sc_b = sc_b;
w.sc_c = sc_c;

    function twoNorms = twonorms(A)
        twoNorms = sqrt(sum(A.^2,1));
    end

end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% POS conversion

posfile = dir('*pos');
pos = importdata(posfile.name);

% determine if timestamps are first or last column
[a b] = min(nanstd(diff(pos))); % find the column with smallest variability..

behav.timestamps = pos(:,b);
pos(:,b) = []; % remove timestamps from pos mat

if size(pos,2) > 5  % if optitrack
    columns = [7 9 8 3 5 4 6 10 1 2];
    pos = pos(:,columns);
    behav.position.x = pos(:,2);
    behav.position.y = pos(:,3);
    behav.position.z = pos(:,4);
    behav.orientation.rx = pos(:,5);
    behav.orientation.ry = pos(:,6);
    behav.orientation.rz = pos(:,7);
    behav.orientation.rw = pos(:,8);
    behav.timestamps = pos(:,1);
    behav.errorPerMarker = pos(:,9);
    behav.frameCount = pos(:,10);

elseif size(pos,2) < 5  % if LED tracking
    behav.position.x = nanmean(pos(:,[2 4]));
    behav.position.y = nanmean(pos(:,[1 3]));
    dx = pos(:,3) - pos(:,5);
    dy = pos(:,2) - pos(:,4);
	ang = atan2(dy,dx)-angOffset;
	ang = mod(ang,2*pi);
    behav.orientation.z = ang;
    warning('come up with a better head dir calculation...')
end

%% LFP conversion

%% SPIKE conversion

% times

spktimes.timestamps = ;

save('spktimes.cellinfo.mat')

% waveforms

save('waveforms.cellinfo.mat')

% features

save('features.cellinfo.mat')

% metadata

save('metadata.cellinfo.mat')

%% METADATA conversion

%% EVENT conversion

%% move old files to FMAT_format folder

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function output = strjoin(input, separator)
%STRJOIN Concatenate an array into a single string.
%
%     S = strjoin(C)
%     S = strjoin(C, separator)
%
% Description
%
% S = strjoin(C) takes an array C and returns a string S which concatenates
% array elements with comma. C can be a cell array of strings, a character
% array, a numeric array, or a logical array. If C is a matrix, it is first
% flattened to get an array and concateneted. S = strjoin(C, separator) also
% specifies separator for string concatenation. The default separator is comma.
%
% Examples
%
%     >> str = strjoin({'this','is','a','cell','array'})
%     str =
%     this,is,a,cell,array
%
%     >> str = strjoin([1,2,2],'_')
%     str =
%     1_2_2
%
%     >> str = strjoin({1,2,2,'string'},'\t')
%     str =
%     1 2 2 string
%

  if nargin < 2, separator = ','; end
  assert(ischar(separator), 'Invalid separator input: %s', class(separator));
  separator = strrep(separator, '%', '%%');

  output = '';
  if ~isempty(input)
    if ischar(input)
      input = cellstr(input);
    end
    if isnumeric(input) || islogical(input)
      output = [repmat(sprintf(['%.15g', separator], input(1:end-1)), ...
                       1, ~isscalar(input)), ...
                sprintf('%.15g', input(end))];
    elseif iscellstr(input)
      output = [repmat(sprintf(['%s', separator], input{1:end-1}), ...
                       1, ~isscalar(input)), ...
                sprintf('%s', input{end})];
    elseif iscell(input)
      output = strjoin(cellfun(@(x)strjoin(x, separator), input, ...
                               'UniformOutput', false), ...
                       separator);
    else
      error('strjoin:invalidInput', 'Unsupported input: %s', class(input));
    end
  end
end

author: Sable

repo: mcbench-benchmarks


text:
function [LRatio IsolDist ID] = plx_getClusterQuality(varargin)
% USAGE
%
% [LRatio IsolDist ID] = getClusterQuality(plxfile)
%
% INPUTS
%
%  plx - filename of a *.plx file in the current working directory (should
%        already be manually cluster cut)
%
% OUTPUTS
%
%  LRatio - L_Ratio.m
%  IsolDist - IsolationDistance.m
%  ID -
%
%
% pulls the spike waveforms out of a .plx file and returns cluster quality
% metrics

if nargin < 1
    plx = dir('*plx');
    for i=1:length(plx)
        if ~isempty(findstr(lower(plx(i).name),'cut'))
            plx = plx(i).name;
        end
    end
else
    plx = varargin{1};
end

if isempty(plx)
    LRatio=nan;
    IsolDist=nan;
    ID=nan;
    return
end
try
    [OpenedFileName, Version, Freq, Comment, Trodalness, NPW, PreTresh, ...
        SpikePeakV, SpikeADResBits, SlowPeakV, SlowADResBits, Duration, ...
        DateTime] = plx_information(plx);

    for i=1:Trodalness:48
        for j=1:50  % ..if you're recording > 50 units per *trode, you may want to do some merging...
        [n, npw, ts, wave{i,j}] = plx_waves_v(plx,i,j-1); % -1 includes 'unsorted spikes'
        [n, npw, ts, wave{i+1,j}] = plx_waves_v(plx,i+1,j-1);
        end
    end

    c=1;
    for j=1:Trodalness:48
        w=[];id=[];
        for i=1:50
            if wave{j,i}~=-1
                if size(wave{j,i},1) == size(wave{j+1,i},1)
                w =[w; wave{j,i}, wave{j+1,i}];
                id = [id;i*ones(size(wave{j,i},1),1)];
                end
            end
        end
        for i=1:50
            if ~isempty(find(id==i))
           IsolDist(c) = IsolationDistance(w, find(id==i));
           [L, LRatio(c), df]  = L_Ratio(w, find(id==i));
           if i == 1
               ID(c) = 0;
           else
               ID(c) = 1;
           end
           c=1+c;
            end
        end
        clear w id
    end
end


text:
function fh = decomp_reconst_full(im,Nsc,Nor,block,noise,parent,covariance,optim,sig);

% Decompose image into subbands, denoise, and recompose again.
%		fh = decomp_reconst(im,Nsc,Nor,block,noise,parent,covariance,optim,sig);
%       covariance:	 are we considering covariance or just variance?
%       optim:		 for choosing between BLS-GSM (optim = 1) and MAP-GSM (optim = 0)
%       sig:        standard deviation (scalar for uniform noise or matrix for spatially varying noise)
% Version using the Full steerable pyramid (2) (High pass residual
% splitted into orientations).

% JPM, Univ. de Granada, 5/02
% Last Revision: 11/04

if (block(1)/2==floor(block(1)/2))|(block(2)/2==floor(block(2)/2)),
   error('Spatial dimensions of neighborhood must be odd!');
end

if ~exist('parent'),
        parent = 1;
end

if ~exist('covariance'),
        covariance = 1;
end

if ~exist('optim'),
        optim = 1;
end

if ~exist('sig'),
        sig = sqrt(mean(noise.^2));
end

[pyr,pind] = buildFullSFpyr2(im,Nsc,Nor-1);
[pyrN,pind] = buildFullSFpyr2(noise,Nsc,Nor-1);
pyrh = real(pyr);
Nband = size(pind,1)-1;
for nband = 2:Nband, % everything except the low-pass residual
  fprintf('%d % ',round(100*(nband-1)/(Nband-1)))
  aux = pyrBand(pyr, pind, nband);
  auxn = pyrBand(pyrN, pind, nband);
  [Nsy,Nsx] = size(aux);
  prnt = parent & (nband < Nband-Nor);   % has the subband a parent?
  BL = zeros(size(aux,1),size(aux,2),1 + prnt);
  BLn = zeros(size(aux,1),size(aux,2),1 + prnt);
  BL(:,:,1) = aux;
  BLn(:,:,1) = auxn*sqrt(((Nsy-2)*(Nsx-2))/(Nsy*Nsx));     % because we are discarding 2 coefficients on every dimension
  if prnt,
  	aux = pyrBand(pyr, pind, nband+Nor);
    auxn = pyrBand(pyrN, pind, nband+Nor);
    if nband>Nor+1,     % resample 2x2 the parent if not in the high-pass oriented subbands.
	   aux = real(expand(aux,2));
       auxn = real(expand(auxn,2));
    end
  	BL(:,:,2) = aux;
    BLn(:,:,2) = auxn*sqrt(((Nsy-2)*(Nsx-2))/(Nsy*Nsx)); % because we are discarding 2 coefficients on every dimension
  end

  sy2 = mean2(BL(:,:,1).^2);
  sn2 = mean2(BLn(:,:,1).^2);
  if sy2>sn2,
     SNRin = 10*log10((sy2-sn2)/sn2);
  else
     disp('Signal is not detectable in noisy subband');
  end

  % main
  BL = denoi_BLS_GSM_band(BL,block,BLn,prnt,covariance,optim,sig);
  pyrh(pyrBandIndices(pind,nband)) = BL(:)';
end
fh = reconFullSFpyr2(pyrh,pind);

author: ricedsp


text:
% Serhat Selcuk Bucak, bucakser@msu.edu

% Lets assume that all samples are stored in data matrix V (each column is a sample)

%Execute NMF for the first n samples
%% These number are selected just for demonstration
V=rand(40,500);
n=100;
rdim=10;
maxiter=150;
[W, H, objhistory] = nmf(V(:,1:n), rdim, 0, maxiter);
% Now we can execute inmf on each new samples
maxiter=50;
A=V(:,1:n)*H';
B=H*H';
h=H(:,end); % Warm start for h
for i=n+1:size(V,2)
    i
    V_new=V(:,i);
    [W_new, h, A, B] = inmf( V_new, W, h, A, B, rdim, 0.9, 0.1, maxiter);
    H_store(:,i-n)=h; %Just for demonstration
end

author: andrewssobral

repo: lrslibrary


text:
%% Copyright (C) 2016, 2018-2019, 2022 Colin B. Macdonald
%%
%% 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
%% @defmethod  @@sym dilog (@var{z})
%% Symbolic dilogarithm function.
%%
%% Example:
%% @example
%% @group
%% syms z
%% dilog (z)
%%   @result{} ans = (sym) polylog(2, 1 - z)
%% @end group
%% @end example
%%
%% @seealso{@@sym/polylog}
%% @end defmethod

function L = dilog(z)
  if (nargin ~= 1)
    print_usage ();
  end

  L = polylog(2, 1 - z);
end

%!assert (isequal (dilog (sym(1)), sym(0)))
%!assert (isequal (dilog (sym(0)), sym(pi)^2/6))

%!assert (isequal (dilog (sym(2)), -sym(pi)^2/12))

%!assert (double(dilog(sym(-1))), pi^2/4 - pi*1i*log(2), eps)

%!test
%! % round-trip
%! syms x
%! f = dilog (x);
%! h = function_handle (f);
%! A = h (1.1);
%! B = dilog (1.1);
%! assert (A, B, -eps)

repo: octsympy


text:
function d = setoptparams(d, varargin)

%tstoolbox/@description/setoptparams
%   Syntax:
%     * d = setoptparams(d, nr, param)
%
%   set optional parameter number nr
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(2,3);

if nargin < 3
	nr = 1;
	param = varargin{1};
else
	nr = varargin{1};
	param = varargin{2};
end

d.optparam{nr} = param;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ out ] = proj_box( x,l,u )
%PROJ_BOX computes the orthogonal projection onto the box {x:l<=x<=u}
%
%  Usage:
%  out = PROJ_BOX(x,l,u)
%  ===========================================
%  Input:
%  x - point to be projected (vector/matrix)
%  l - lower bound (vector/matrix/scalar)
%  u - upper bound (vector/matrix/scalar)
%  ===========================================
%  Assumptions:
%  l<=u
%  ===========================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%reading the user x and
if (nargin < 3)
    error ('usage: proj_box( x,l,u )') ;
end

if any(any((l > u)))
    error('Set is infeasible') ;
end

out= min(max(l,x),u) ;
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
% eeg_miclust() - Calculate mutual information matrix between independent component
%                 activations and cluster the components using it. Opens several figures:
%                -  A figure with sorted silhouette values of each cluster
%                -  A figure showing clusters in 3-D, with brightness of
%                     each point proportional to its silhouette value
%                -  Multiple 2-D plots, each for one cluster. Background
%                     color shows interpolated silhouette value. A bright
%                     background means the component fits well in the cluster;
%                     a dark background means it belongs near equally
%                     to another cluster or to no cluster.
%                -  A dendrogram of componets if 'linkage' method is used
% Usage:
%              >> EEG = eeg_miclust(EEG);
%              >> EEG = eeg_miclust(EEG,N,components, clusterMethod);
% Inputs:
%
%   EEG        - EEG data structure
%   N          - number of clusters to produce {default: 5}
%   components - a vector containing component indices (rows of sim matrix)
%                to cluster. For example [1:10] uses only first 10 components
%                {Default: [1:120] or all}
%
% Optional Inputs:
%
%   clusterMethod - which clustering method to use, options are 'linkage'
%                   and 'kmeans' {defgault = 'linkage'}
% Output:
%
%   EEG        - input EEG structure containing mutual information
%                between specified components in field EEG.etc.miclust.mutual_info and indices of
%                these components in field EEG.etc.miclust.allcomponents Cluster information
%                is placed in EEG.etc.miclust field.
% Example:
%
%   % Cluster components 1 to 30 into 4 clusters using mutual information.
%   >> EEG = eeg_miclust(EEG,4,1:30);
%
% See also: %   getmiclusts(), showmiclusts(), mi_pairs()
%
% Author: Nima Bigdely Shamlo, SCCN/INC/UCSD, 2006

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

% Edit history:
%
% 5/26/06 added default for N, test for <120 comps, extended help msg, formatted -sm
% 5/30/06 use only first 120 comps -nb

function EEGout = eeg_miclust(EEG,n,comps, clusterMethod)

if nargin<2
   n = 5; % default to 5 clusters
end
if nargin<3
    if size(EEG.icaact,1) >= 120
      comps = 120;
    else
      comps = 1:size(EEG.icaact,1);
    end;
end;

if nargin<4
   clusterMethod = 'linkage' ;
else
    if ~strcmp(clusterMethod, 'kmeans') && ~strcmp(clusterMethod, 'linkage')
        fprintf('Clustering method not found.');
        return;
    end;
end

if ~(isfield(EEG,'etc') && isfield(EEG.etc,'miclust') && isfield(EEG.etc.miclust,'mutual_info')  ) || isempty(EEG.etc.miclust.allcomponents) || ~isfield(EEG.etc.miclust,'allcomponents') || length(EEG.etc.miclust.allcomponents) ~= length(comps) || sum(EEG.etc.miclust.allcomponents ~= comps) > 0
    fprintf('Calculating mutual information matrix (time consuming)...');
    mutual_info = mi_pairs(EEG.icaact(comps,:),'full',false);
else
    mutual_info = EEG.etc.miclust.mutual_info;
end

EEG.etc.miclust = getmiclusts(mutual_info,n, comps, clusterMethod);

EEG.etc.miclust.mutual_info = mutual_info;

plotTopo = @(x) topoplot(EEG.icawinv(:,x), EEG.chanlocs,'electrodes','off');

showmiclusts(EEG.etc.miclust,plotTopo);

if nargout>0
    EEGout = EEG;
end;

author: PatternRecognition

repo: OpenBMI


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

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

if SO3F.c0 ~= 0
  disp('  uniform component');
  disp(['  weight: ',xnum2str(SO3F.c0)]);

  disp(' ');
end

if ~isempty(SO3F.center)

  if length(SO3F.center)==1
    disp('  unimodal component');
  else
    disp('  multimodal components');
  end
  disp(['  kernel: ',char(SO3F.psi)]);
  if isa(SO3F.center,'SO3Grid')
    disp(['  center: ',char(SO3F.center)]);
    disp(['  weight: ',xnum2str(sum(SO3F.weights(:)))]);
    disp(' ');
  else
    disp(['  center: ',num2str(length(SO3F.center)), ' orientations']);
    s.weight = SO3F.weights(:);

    if length(SO3F.center) < 20 && ~isempty(SO3F.center)
       Euler(SO3F.center,s)
    elseif ~getMTEXpref('generatingHelpMode')
      disp(' ')
      s = setappdata(0,'data2beDisplayed',SO3F.center);
      setappdata(0,'data2beDisplayedWeights',s);
      disp(['  show centers of the components and corresponding weights'])
      disp(' ')
    end

  end
end

end

author: mtex-toolbox

repo: mtex

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


text:
function test_bug3461

% WALLTIME 00:20:00
% MEM 3gb
% DEPENDENCY ft_sourcestatistics ft_sourcegrandaverage

%%

% this contains allSourceDiff050
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3461.mat'));

%% Grand average

cfg = [];
cfg.parameter = 'coh';
cfg.keepindividual = 'yes';
avgSourceDiff050 = ft_sourcegrandaverage(cfg, allSourceDiff050{:});

%% Make null condition

cfg = [];
cfg.parameter = 'coh';
cfg.operation = 'multiply';
cfg.scalar = 0;
avgSourceZero = ft_math(cfg, avgSourceDiff050);

%% Run statistics

cfg = [];
cfg.parameter   = 'coh';
cfg.method      = 'montecarlo';
cfg.statistic   = 'depsamplesT';
cfg.alpha       = 0.05;
cfg.tail        = 0;
cfg.correcttail = 'alpha';
cfg.correctm    = 'cluster';
cfg.numrandomization = 1000;

Nsub = 37;
cfg.design(1,1:2*Nsub)  = [ones(1,Nsub) 2*ones(1,Nsub)];
cfg.design(2,1:2*Nsub)  = [1:Nsub 1:Nsub];
cfg.ivar                = 1; % the 1st row in cfg.design contains the independent variable
cfg.uvar                = 2; % the 2nd row in cfg.design contains the subject number

tval_diff_050 = ft_sourcestatistics(cfg, avgSourceDiff050, avgSourceZero);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [testsResults] = runTestSuite(suiteTag)
%runtests: Run tagged tests from all subdirectories
%   suiteTag: String matching tags from test classes in subdirectories
%
    import matlab.unittest.TestSuite;
    fullSuite = TestSuite.fromFolder(pwd, 'IncludingSubfolders', true);
    persistenceSuite = fullSuite.selectIf('Tag',suiteTag);
    testsResults = run(persistenceSuite)
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%% SIMPLE_ODE113 solves the three body problem using MATLAB's ODE113.
%
%  Discussion:
%
%    Three bodies, regarded as point masses, are constrained to lie in a plane.
%    The masses of each body are given, as are the positions and velocities
%    at a starting time T = 0.  The bodies move in accordance with the gravitational
%    force between them.
%
%    The force exerted on the 0-th body by the 1st body can be written:
%
%      F = - m0 m1 ( p0 - p1 ) / |p0 - p1|^3
%
%    assuming that units have been normalized to that the gravitational
%    coefficient is 1.  Newton's laws of motion can be written:
%
%      m0 p0'' = - m0 m1 ( p0 - p1 ) / |p0 - p1|^3
%                - m0 m2 ( p0 - p2 ) / |p0 - p2|^3
%
%      m1 p1'' = - m1 m0 ( p1 - p0 ) / |p1 - p0|^3
%                - m1 m2 ( p1 - p2 ) / |p1 - p2|^3
%
%      m2 p2'' = - m2 m0 ( p2 - p0 ) / |p2 - p0|^3
%                - m2 m1 ( p2 - p1 ) / |p2 - p1|^3
%
%    Letting
%
%      y1 = p0(x)
%      y2 = p0(y)
%      y3 = p0'(x)
%      y4 = p0'(y)
%
%    and using similar definitions for p1 and p2, the 3 second order vector
%    equations can be rewritten as 12 first order equations.  In particular,
%    the first four are:
%
%      y1' = y3
%      y2' = y4
%      y3' = - m1 ( y1 - y5  ) / |(y1,y2) - (y5,y6) |^3
%            - m2 ( y1 - y9  ) / |(y1,y2) - (y9,y10)|^3
%      y4' = - m1 ( y2 - y6  ) / |(y1,y2) - (y5,y6) |^3
%            - m2 ( y2 - y10 ) / |(y1,y2) - (y9,y10)|^3
%
%    and so on.
%
%    This first order system can be integrated by a standard ODE solver.
%
%    Note that when any two bodies come close together, the solution changes
%    very rapidly, and very small steps must be taken by the ODE solver.
%    For this system, the first near collision occurs around T=15.8299, and
%    the results produced by MATLAB's ode113 will not be very accurate after
%    that point.
%
%  Modified:
%
%    03 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dominik Gruntz, Joerg Waldvogel,
%    "Orbits in the Planar Three-Body Problem",
%    Walter Gander, Jiri Hrebicek,
%    Solving Problems in Scientific Computing using Maple and Matlab,
%    Springer, 1997,
%    ISBN: 3-540-61793-0,
%    LC: Q183.9.G36.
%
  addpath ( '../rkf45' );

  global m0 m1 m2

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIMPLE_RKF45:\n' );
  fprintf ( 1, '  A simple formulation of the planar three-body problem.\n' );
  fprintf ( 1, '  This program uses RKF45 for the ODE solver.\n' );
%
%  Set the masses.
%
  m0 = 5.0;
  m1 = 3.0;
  m2 = 4.0;
%
%  Set the time range.
%
  t_start = 0.0;
  t_stop = 63.0;
  step_num = 631;
%
%  Set input to RKF45.
%
  neqn = 12;
  y = [ 1.0; -1.0;  0.0;  0.0;
        1.0;  3.0;  0.0;  0.0;
       -2.0; -1.0;  0.0;  0.0 ];
  yp = simple_f ( t_start, y );
  relerr = 1.0E-10;
  abserr = 1.0E-10;
  flag = 1;
%
%  Set up the array that stores the data.
%
  T1 = zeros ( step_num + 1 );
  Y1 = zeros ( step_num + 1, neqn );
  step = 0;
  T1(step+1) = t_start;
  Y1(step+1,1:neqn) = y';

  for step = 1 : step_num

    t = ( ( step_num - step + 1 ) * t_start ...
        + (            step - 1 ) * t_stop ) ...
        / ( step_num              );

    t_out = ( ( step_num - step ) * t_start ...
            + (            step ) * t_stop ) ...
            / ( step_num          );

    [ y, yp, t, flag ] = r8_rkf45 ( @simple_f, neqn, y,  yp, t, t_out, ...
      relerr, abserr, flag );

    T1(step+1) = t;
    Y1(step+1,1:neqn) = y';

  end
%
%  Display the results.
%
  figure ( 1 )

  [ ~, i10 ] = min ( abs ( T1 - 10.0 ) );

  R1 = 1 : i10;

  plot ( Y1(R1,1), Y1(R1,2), 'b.', ...
         Y1(R1,5), Y1(R1,6), 'r.', ...
         Y1(R1,9), Y1(R1,10), 'g.' )
  title ( '0 <= T <= 10' )

  figure ( 2 )

  [ ~, i20 ] = min ( abs ( T1 - 20.0 ) );

  R2 = i10 : i20;

  plot ( Y1(R2,1), Y1(R2,2), 'b.', ...
         Y1(R2,5), Y1(R2,6), 'r.', ...
         Y1(R2,9), Y1(R2,10), 'g.' )
  title ( '10 <= T <= 20' )

  figure ( 3 )

  [ ~, i50 ] = min ( abs ( T1 - 50.0 ) );
  i63 = length ( T1 );

  R3 = i50 : i63;

  plot ( Y1(R3,1), Y1(R3,2), 'b.', ...
         Y1(R3,5), Y1(R3,6), 'r.', ...
         Y1(R3,9), Y1(R3,10), 'g.' )
  title ( '50 <= T <= 63' )

  figure ( 4 )

  R4 = 1 : i63;

  plot ( Y1(R4,1), Y1(R4,2), 'b.', ...
         Y1(R4,5), Y1(R4,6), 'r.', ...
         Y1(R4,9), Y1(R4,10), 'g.' )
  title ( '0 <= T <= 63' )

  rmpath ( '../rkf45' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIMPLE_RKF45:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

author: johannesgerer

repo: jburkardt-m


text:
function x=randcorr(n,R)
% RANDCORR Generates corremlated random variables
% Generates n vector valued variates with uniform marginals and correlation
% matrix R.
% Returns an nxk matrix, where k is the order of R.
  k=size(R,1);
  R=2*sin((pi/6)*R);
  x=normcdf(randn(n,k)*chol(R));

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function out=subsref(x,varargin)

out=builtin('subsref',x,varargin{:});

author: opencobra

repo: cobratoolbox


text:
function d=center(x,op);

%  CENTER  Delay calculation for Wavelet transform alignment.
%
%          CENTER (X, OP) calculates the delay for filter in X
%          according to the alignment method indicated in OP.
%          This delay is used by Wavelet transform functions.
%          The value of OP can be:
%              0 : First Absolute Maxima Location
%              1 : Zero delay in analysis (Full for synthesis).
%              2 : Mass center (sum(m*d)/sum(m))
%              3 : Energy center (sum(m^2 *d)/sum(m^2))
%
%          If no output argument is given, then the vector X will
%          be plotted in the current figure, and a color line will be
%          marking the result.(red: OP=0; green: OP=1; cyan: OP=2;
%          blue: OP=4)
%
%          See also: WTCENTER, WTMETHOD

%--------------------------------------------------------
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
%
%
% Uvi_Wave is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2, or (at your option) any
% later version.
%
% Uvi_Wave 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 Uvi_Wave; see the file COPYING.  If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%       Author: Sergio J. Garcia Galan
%       e-mail: Uvi_Wave@tsc.uvigo.es
%--------------------------------------------------------

lx=length(x);
l=1:lx;

if op==1
	d=0;
else
	if op==2
		xx=abs(x(:)');
		L=l;
	end
	if op==3
		xx=x(:)'.^2;
		L=l;
	end
	if op==0
		[mx,d]=max(abs(x));
	else

		d=sum(xx.*L)/sum(xx);
	end
end

if nargout==0,
  cad='rgcbk';
  plot(x)
  l=line([d,d],[min(x),max(x)]);
  set(l,'Color',cad(op+1));
end

author: alexandrebarachant


text:
function SwatheOpt = PlotSwatheOption
% =======================================================================
% Optional inputs for PlotSwathe
% =======================================================================
% SwatheOpt = PlotSwatheOption
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2015. Updated November 2020
% -----------------------------------------------------------------------

%% INITIALIZE VARIABLES
%color = bone(8);
SwatheOpt.swathecol  = [138, 178, 212]./255; % color(6,:);
SwatheOpt.linecol    = [13, 54, 84]./255;    % color(1,:);
SwatheOpt.do_dates   = 0;
SwatheOpt.frequency  = 'q';
SwatheOpt.swatheonly = 0;
SwatheOpt.marker     = 'none';
SwatheOpt.transp     = 0;   % Allow for transparency when plotting multiple swathes
SwatheOpt.alpha      = 1;   % Sets alpha for transparent swathes. Default is 1

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
%% Test suite for Minimum bounding objects: RECT, CIRCLE, SPHERE, ELLIPSE

%% Minimal area rectangle around random data

% Random points in the unit square
x = rand(1000,1);
y = rand(1000,1);

% should have a nearly square bounding rectangle
tic,[xrect,yrect] = minboundrect(x,y),toc

%
plot(x,y,'ro',xrect,yrect,'b-')

%% Minimal circle around (normally distributed) random data
x = randn(100000,1);
y = randn(100000,1);

% It should have a center of roughly (0,0).
tic,[center,radius] = minboundcircle(x,y),toc

theta = linspace(0,2*pi,100);
xc = center(1) + radius*cos(theta);
yc = center(2) + radius*sin(theta);

% The circle will generally pass exactly through
% three of the points.
plot(x,y,'ro',xc,yc,'b-')

%% Minimal circle around nearly linear data.
x = rand(5,1);
y = x + (rand(size(x))-.5)/10;

[center,radius] = minboundcircle(x,y)

theta = linspace(0,2*pi,100);
xc = center(1) + radius*cos(theta);
yc = center(2) + radius*sin(theta);

% Sometimes the minimal radius circle may pass
% through only two points.
plot(x,y,'ro',xc,yc,'b-')
axis equal

%% Degenerate cases - a minimal circle around exactly two points
x = [0 1];
y = x;

[center,radius] = minboundcircle(x,y)

%% Degenerate cases - a minimal circle around exactly one point
x = 1;
y = 2;

% The radius must be zero
[center,radius] = minboundcircle(x,y)

%% Minimum radius sphere around random data
xyz = rand(10000,3);

% The sphere should have approximate center [.5 .5 .5]
% and approximate radius sqrt(3)/2 = .86603...
[center,radius] = minboundsphere(xyz)

%% Minimum area ellipse around exactly circular data
theta = linspace(0,2*pi,100);
x = cos(theta);
y = sin(theta);

% will have center == [0,0], radius = 1, both to within
% double precision accuracy in Matlab.
[center,H,xe,ye] = minboundellipse(x,y)

plot(x,y,'ro',xe,ye,'b-')
axis equal

%% Minimum area ellipse around correlated data
xy = randn(100,2);
xy = xy*rand(2)^2;

% will have center [0,0], radius = 1, both to within
% double precision accuracy in Matlab.
[center,H,xe,ye] = minboundellipse(xy(:,1),xy(:,2))

plot(xy(:,1),xy(:,2),'ro',xe,ye,'b-')
axis equal

%%


text:
function b = hoffman_rhs ( )

%*****************************************************************************80
%
%% HOFFMAN_RHS returns the Hoffman right hand side.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real B(3,1), the right hand side vector.
%
  b = [ 1.0; ...
        0.0; ...
        0.0 ];

  return
end


text:
classdef CEC2017_F23 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{12}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 30
                obj.D   = 10;
                obj.Mat = Data{12}.M_10;
            elseif obj.D < 50
                obj.D   = 30;
                obj.Mat = Data{12}.M_30;
            elseif obj.D < 100
                obj.D   = 50;
                obj.Mat = Data{12}.M_50;
            else
                obj.D   = 100;
                obj.Mat = Data{12}.M_100;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Z = Z*obj.Mat';
            PopObj = -20*exp(-0.2*sqrt(mean(Z.^2,2))) + 20 - exp(mean(cos(2*pi*Z),2)) + exp(1);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Z = Z*obj.Mat';
            PopCon(:,1) = sum(Z(:,2:end).^2,2) + 1 - abs(Z(:,1));
            PopCon(:,2) = abs(sum(Z.^2,2)-4) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [AddedRxnsExtended] = postProcessGapFillSolutions(AddedRxns, model, BlockedRxns, IdentifyPW)
% This function uses the output from `fastGapFill` and signs generalized
% function to the identified solutions from the universal database (e.g.,
% 'Metabolic reaction', 'Transport reaction', 'Exchange reaction'). It also
% provides basic statistics for the solution. One can use the option
% 'IdentifyPW' to compute for each solved blocked reaction a flux vector,
% such that the flux through this reaction is maximized while the sum of
% all fluxes is minimized. This analysis can be helpful to put the solution
% reactions into the network context.
%
% USAGE:
%
%    [AddedRxnsExtended] = postProcessGapFillSolutions(AddedRxns, model, BlockedRxns, IdentifyPW)
%
% INPUTS:
%    AddedRxns:            Output from `fastGapFill.m`
%    model:                Model structure. Original model given as input into `prepareFastGapFill`
%    IdentifyPW:           Optional: identifies reactions that are need from
%                          consistent model and added reactions to have flux through previously
%                          blocked reaction, default: false
%
% OUTPUT:
%    AddedRxnsExtended:    Structure containing the information present in
%                          `AddedRxns` along with the statistics and if desired pathways containing
%                          the flux vectors.
%
% .. Author: - Ines Thiele, June 2013, http://thielelab.eu

if ~exist('IdentifyPW','var')
    IdentifyPW = 0;
end

% postprocess solutions
if length(AddedRxns.rxns) > 0
    [a,b] = strtok(AddedRxns.rxns,'\[');
    b = regexprep(b,'\]r','');
    b = regexprep(b,'\[','');
    Metabolic = 0;
    Transport = 0;
    Exchange = 0;
    for i = 1 : length(b)
        if isempty(char(b(i)))
            AddedRxns.subSystem{i}={'Metabolic reaction'};
            Metabolic = Metabolic + 1;
        elseif strcmp('e]',b(i))
            AddedRxns.subSystem{i}={'Exchange reaction'};
            Exchange = Exchange + 1;
        elseif strcmp('n',b(i))
            AddedRxns.subSystem{i}={'Nucleus <=> cytosol transport'};
            Transport = Transport + 1;
        elseif strcmp('g',b(i))
            AddedRxns.subSystem{i}={'Golgi <=> cytosol transport'};
            Transport = Transport + 1;
        elseif strcmp('m',b(i))
            AddedRxns.subSystem{i}={'Mitochondrion <=> cytosol transport'};
            Transport = Transport + 1;
        elseif strcmp('l',b(i))
            AddedRxns.subSystem{i}={'Lysosome <=> cytosol transport'};
            Transport = Transport + 1;
        elseif strcmp('r',b(i))
            AddedRxns.subSystem{i}={'ER <=> cytosol transport'};
            Transport = Transport + 1;
        elseif strcmp('e',b(i))
            AddedRxns.subSystem{i}={'Extracellular space <=> cytosol transport'};
            Transport = Transport + 1;
        elseif strcmp('x',b(i))
            AddedRxns.subSystem{i}={'Peroxisome <=> cytosol transport'};
            Transport = Transport + 1;
        elseif strcmp('p',b(i))
            AddedRxns.subSystem{i}={'Periplasm <=> cytosol transport'};
            Transport = Transport + 1;
        else
            AddedRxns.subSystem{i}={strcat(b{i},' <=> cytosol transport')};
            Transport = Transport + 1;
        end
    end
    AddedRxns.subSystem = AddedRxns.subSystem';
    AddedRxns.Stats.metabolicSol = Metabolic;
    AddedRxns.Stats.transportSol = Transport;
    AddedRxns.Stats.exchangeSol = Exchange;

    if IdentifyPW
        % model
        modelFilled = model;
        % Index = length(model.rxns)+1;
        % add found solutions to model
        for i = 1 : length(AddedRxns.rxns)
            modelFilled = addReaction(modelFilled,char(AddedRxns.rxns(i)),char(AddedRxns.rxnFormula(i)));
        end

        % find blocked reactions in modelFilled
        Rxns = find(ismember(modelFilled.rxns,BlockedRxns.solvableRxns));

        for i = 1 : length(BlockedRxns.solvableRxns)
            modelFilled.c = zeros(length(modelFilled.c),1);
            % set blocked reaction as objective
            modelFilled.c(Rxns(i)) = 1;
            % perform LP to reconstitute pathway
            % I am using taxi cab (Manhattan norm) as the shortest possible
            % solution
            FBA = optimizeCbModel(modelFilled,'max','one');
            if FBA.f == 0
                % minimize flux through reaction
                % perform LP to reconstitute pathway
                FBA = optimizeCbModel(modelFilled,'min','one');
            end
            % find reactions in FBA solution
            name= modelFilled.rxns{Rxns(i)};
            name = regexprep(name,'\(e\)','');
            name = regexprep(name,'^\d','x');
            if ~isempty(FBA.x)
                FBA.x(FBA.x>=-1e-8&FBA.x<=1e-8)=0;
                if ~isempty(find(FBA.x)~=0)
                    AddedRxns.solvableRxnsSol.(name).solution(:,1) = modelFilled.rxns(FBA.x~=0);
                    AddedRxns.solvableRxnsSol.(name).solution(:,2) = printRxnFormula(modelFilled,AddedRxns.solvableRxnsSol.(name).solution(:,1),false);
                    tmp = FBA.x(FBA.x~=0);
                    for k = 1 : length(tmp)
                        AddedRxns.solvableRxnsSol.(name).solution{k,3}=num2str(tmp(k));
                    end
                    clear tmp
                    % highlight which reaction have not been part of flux consistent
                    % network
                    tmp2 = find(ismember(AddedRxns.solvableRxnsSol.(name).solution(:,1),AddedRxns.rxns));
                    AddedRxns.solvableRxnsSol.(name).solution(tmp2,4)={'added'};
                    clear tmp2
                    tmp2 = find(ismember(AddedRxns.solvableRxnsSol.(name).solution(:,1),BlockedRxns.solvableRxns));
                    AddedRxns.solvableRxnsSol.(name).solution(tmp2,4)={'previously blocked'};

                    clear tmp2
                else
                    AddedRxns.Problem.(name) = BlockedRxns.solvableFormula(i);
                end
            else
                AddedRxns.Problem.(name) = BlockedRxns.solvableFormula(i);
            end
            clear FBA name
        end
    end
    AddedRxnsExtended = AddedRxns;
else
    AddedRxnsExtended = [];
end


text:
% orientation filtering function(Eq.(9)in the paper)
function K = angfilter(I,Theta,angle,sigma)

b = angle;
c = sigma;
[m,n] = size(Theta);

X = 0:180;
f = (1/(sqrt(2*pi)*c))*exp(-(((X-b).^2)/(2*(c^2))));

for i=1:m
    for j=1:n
        T(i,j)=f(Theta(i,j));
    end
end

K = I.*T;


text:
function [K,H,Pmax,Pmin] = surfature(X,Y,Z),
% SURFATURE -  COMPUTE GAUSSIAN AND MEAN CURVATURES OF A SURFACE
%   [K,H] = SURFATURE(X,Y,Z), WHERE X,Y,Z ARE 2D ARRAYS OF POINTS ON THE
%   SURFACE.  K AND H ARE THE GAUSSIAN AND MEAN CURVATURES, RESPECTIVELY.
%   SURFATURE RETURNS 2 ADDITIONAL ARGUEMENTS,
%   [K,H,Pmax,Pmin] = SURFATURE(...), WHERE Pmax AND Pmin ARE THE MINIMUM
%   AND MAXIMUM CURVATURES AT EACH POINT, RESPECTIVELY.

% First Derivatives
[Xu,Xv] = gradient(X);
[Yu,Yv] = gradient(Y);
[Zu,Zv] = gradient(Z);

% Second Derivatives
[Xuu,Xuv] = gradient(Xu);
[Yuu,Yuv] = gradient(Yu);
[Zuu,Zuv] = gradient(Zu);

[Xuv,Xvv] = gradient(Xv);
[Yuv,Yvv] = gradient(Yv);
[Zuv,Zvv] = gradient(Zv);

% Reshape 2D Arrays into Vectors
Xu = Xu(:);   Yu = Yu(:);   Zu = Zu(:);
Xv = Xv(:);   Yv = Yv(:);   Zv = Zv(:);
Xuu = Xuu(:); Yuu = Yuu(:); Zuu = Zuu(:);
Xuv = Xuv(:); Yuv = Yuv(:); Zuv = Zuv(:);
Xvv = Xvv(:); Yvv = Yvv(:); Zvv = Zvv(:);

Xu          =   [Xu Yu Zu];
Xv          =   [Xv Yv Zv];
Xuu         =   [Xuu Yuu Zuu];
Xuv         =   [Xuv Yuv Zuv];
Xvv         =   [Xvv Yvv Zvv];

% First fundamental Coeffecients of the surface (E,F,G)
E           =   dot(Xu,Xu,2);
F           =   dot(Xu,Xv,2);
G           =   dot(Xv,Xv,2);

m           =   cross(Xu,Xv,2);
p           =   sqrt(dot(m,m,2));
n           =   m./[p p p];

% Second fundamental Coeffecients of the surface (L,M,N)
L           =   dot(Xuu,n,2);
M           =   dot(Xuv,n,2);
N           =   dot(Xvv,n,2);

[s,t] = size(Z);

% Gaussian Curvature
K = (L.*N - M.^2)./(E.*G - F.^2);
K = reshape(K,s,t);

% Mean Curvature
H = (E.*N + G.*L - 2.*F.*M)./(2*(E.*G - F.^2));
H = reshape(H,s,t);

% Principal Curvatures
Pmax = H + sqrt(H.^2 - K);
Pmin = H - sqrt(H.^2 - K);


text:
function p=v_psychofunc(m,q,x,r)
%V_PSYCHOFUNC Calculate psychometric functions: trial success probability versus SNR
%
% Usage: p=v_psychofunc('',q,x)       % calculate probabilities
%        b=v_psychofunc('r',q,x)       % generate boolean variables with success prob p
%        p=v_psychofunc(m,q,x,r)     % Calculate likelihoods for observations r
%        x=v_psychofunc([m 'i'],q,p) % Calculate inverse
%
% Inputs:
%          m        mode string [may be omitted if not required]
%                      'n'   do not normalize likelihoods
%                      'f'   do not squeeze output arrays to remove singleton dimensions
%                      'i'   calculate inverse function
%                      'r'   calculate binary random variables with probability p
%                      ['s'   calculate sweet points for threshold and slope]
%                      ['d'   calculate partial derivatives with respect to q(1:5)]
%                      'g'   plot graph
%                      'G'   plot image
%                      'c'   include colourbar
%          q        model parameters. Either a column vector with a single model,
%                   a matrix with one model per column or a cell array with multiple values for
%                   some or all of the parameters
%                      1  probability at threshold [0.5]
%                      2  threshhold [0 dB]
%                      3  slope at threshold [0.1 prob/dB ]
%                      4  miss or lapse probability [0]
%                      5  guess probability   [0]
%                      6  psychometric function type [1]
%                          1 = logistic
%                          2 = cumulative Gaussian
%                          3 = Weibull
%                          [4 = reversed Weibull]
%                          [5 = Gumbell]
%                          [6 = reversed Gumbell]
%          x        vector of SNR values
%          r        test results (0 or 1) corresponding to x
%          p        vector of probabilities
%
% Outputs:
%          p        array of probabilities or random variates ('r' option).
%                   p is a squeezed 7-dimensional array
%                   whose dimensions correspond to x followed by the 6 model parameter entries.
%                   if q is a cell array, singleton dimensions are removed unless the 'f' option is given.
%          x        Inverse function gives SNR, x, as a function of p
%          b        array of boolean variables

%      Copyright (C) Mike Brookes 2009-2010
%      Version: $Id: v_psychofunc.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first sort out input arguments
minp=0.01;          % minimum probability to use for inverse function by default
qq=[0.5 0 0.1 0 0 1]';  % default values for q
if nargin<4
    r=[];
    if nargin<3
        x=[];
        if nargin<2
            q=[];
            if ~nargin
                m='';
            end
        end
    end
end
if ~ischar(m);      % mode argument is optional
    r=x;
    x=q;
    q=m;
    m='';
end
sq=size(q);
ckmod=0;
if iscell(q)
    nq=ones(1,6);
    qax=num2cell([0; qq]);  % used for plotting
    for i=1:min(numel(q),6)
        nq(i)=numel(q{i});
        if nq(i)>=1
            nr=size(qq,2);
            qax{i+1}=q{i};
            if i<=5             % do not replicate for multiple models
                qq=repmat(qq,1,nq(i));
                qq(i,:)=reshape(repmat(q{i}(:)',nr,1),1,nr*nq(i));
            else
                qq(i,:)=repmat(q{i}(1),1,nr);
            end
        end
    end
    nq=max(nq,1);
    nmod=nq(6);
    if nmod>1      % list of models to use
        modlist=q{6};
    else
        modlist=qq(6,1);  % default model
    end
else
    nq=sq(2);
    if nq
        ql=repmat(qq,1,nq);
        ql(1:sq(1),:)=q;
    else
        ql=qq;
        nq=1;
    end
    modlist=unique(ql(6,:));
    nmod=length(modlist);
    ckmod=nmod>1;   % need to check model list
    qq=ql;
end
% now perform the calculation
nx=numel(x);
npt=50;       % number of points
if any(m=='i') % doing inverse
    if ~nx
        nx=npt;
        xlim=[max(qq(5,:)),1-max(qq(4,:))]*[1-minp minp; minp 1-minp];
        x=linspace(xlim(1),xlim(2),nx)';
    end
    p=zeros([nx nq]);  % space for SNRs
    ia=0;
    for i=1:nmod % loop for each model type
        mod=modlist(i);
        if ckmod
            qq=ql(:,ql(6,:)==mod);
        end
        pscale=1-qq(4,:)-qq(5,:);
        pstd=(qq(1,:)-qq(5,:))./pscale; % prob target compensating for miss and lapse probs
        sstd=qq(3,:)./pscale; % slope compensating for miss and lapse probs
        px=x(:)*pscale.^(-1)-repmat(qq(5,:)./pscale,nx,1);  % adjust for miss and lapse probs
        switch mod
            case 1
                beta=sstd./(pstd.*(1-pstd));
%                 alpha=qq(2,:)+log((1-pstd)./pstd)./beta;
                px=repmat(qq(2,:)+log((1-pstd)./pstd)./beta,nx,1)-log(px.^(-1)-1).*repmat(beta.^(-1),nx,1);
            case 2   % cumulative Gaussian function
                xtstd=norminv(pstd); % x position of target in std measure
                sig=normpdf(xtstd)./sstd;
                px= repmat(qq(2,:)-sig.*xtstd,nx,1) + repmat(sig,nx,1).*norminv(px);
            case 3
                wlog=log(1-pstd);
                kbeta=sstd./((pstd-1).*wlog);
                alpha=qq(2,:)-log(-wlog)./kbeta;
                px=repmat(alpha,nx,1)+log(-log(1-px)).*repmat(kbeta.^(-1),nx,1);
            otherwise
                error('Invalid psychometric model index');
        end
        if ckmod
            p(:,ql(6,:)==i)=px;
        else
            ib=ia+numel(p)/nmod;
            p(ia+1:ib)=px(:);
            ia=ib;
        end
    end
else % doing forward mapping
    if ~nx
        ef=2;         % expansion factor
        nx=npt;
        x=linspace(min(qq(2,:)-ef*(qq(1,:)-qq(5,:))./qq(3,:)), ...
            max(qq(2,:)+ef*(1-qq(1,:)-qq(4,:))./qq(3,:)),nx)';
    end
    p=zeros([nx nq]);  % space for probabilities
    ia=0;
    for i=1:nmod % loop for each model type
        mod=modlist(i);
        if ckmod
            qq=ql(:,ql(6,:)==mod);
        end
        pscale=1-qq(4,:)-qq(5,:);  % prob range excluding miss and lapse probs
        pstd=(qq(1,:)-qq(5,:))./pscale; % prob target compensating for miss and lapse probs
        sstd=qq(3,:)./pscale; % slope compensating for miss and lapse probs
        switch mod
            case 1   % logistic function
                beta=sstd./(pstd.*(1-pstd));
%                 alpha=qq(2,:)+log((1-pstd)./pstd)./beta;
                px=(1+exp(repmat(beta.*qq(2,:)+log((1-pstd)./pstd),nx,1)-x(:)*beta)).^(-1);
            case 2   % cumulative Gaussian function
                xtstd=norminv(pstd); % x position of target in std measure
                sigi=sstd./normpdf(xtstd);
                px=normcdf(x(:)*sigi-repmat(qq(2,:).*sigi-xtstd,nx,1));
            case 3
                wlog=log(1-pstd);
                kbeta=sstd./((pstd-1).*wlog);
                alpha=qq(2,:)-log(-wlog)./kbeta;
                px=1-exp(-exp(x(:)*kbeta-repmat(alpha.*kbeta,nx,1)));
            otherwise
                error('Invalid psychometric model index');
        end
        px=repmat(qq(5,:),nx,1)+repmat(pscale,nx,1).*px;  % adjust for miss and lapse probs
        if ckmod
            p(:,ql(6,:)==i)=px;
        else
            ib=ia+numel(p)/nmod;
            p(ia+1:ib)=px(:);
            ia=ib;
        end
    end
    if numel(r)                 % we are calculating likelihoods
        mk=r(:)==0;
        p(mk,:)=1-p(mk,:);      % invert probability for results that are zero
        if nx>1
            if any(m=='n')
                p=prod(p,1);
            else
                p=sum(log(p),1);
                p=exp(p-max(p(:)));
                p=p/sum(p(:));     % normalize to equal 1
            end
            nx=1;
        end
    end

end
pg=p;       % save unsqueezed p for plotting
if ~any(m=='f') && iscell(q) % remove all singleton dimensions
    szp=size(p);
    szq=szp(szp>1);
    szq=[szq ones(1,max(0,2-numel(szq)))];
    p=reshape(p,szq);
end
if any(m=='r') && ~any(m=='i');
    p=rand(size(p))<p;
end

if ~nargout || any(lower(m)=='g')
    clf;
    szp=[nx nq];
    czp=sum(szp>1);
    if czp>0  % check if there is anything to plot
        if iscell(q)
            axlab={'Input SNR','Threshold prob','Threshold SNR','Threshold slope','Lapse prob','Guess prob','Sigmoid type'};
            [szs,izs]=sort(szp,'descend');
            pg=permute(pg,izs);
            qax{1}=x;
            if any(m=='G') || czp>2 % image
                ngr=prod(szs(3:end));
                ncol=ceil(sqrt(ngr));
                nrow=ceil(ngr/ncol);
                npix=szs(1)*szs(2);
                ia=0;
                for i=1:ngr
                    subplot(nrow,ncol,i);
                    ib=ia+npix;
                    imagesc(qax{izs(1)},qax{izs(2)},reshape(pg(ia+1:ib),szs(1:2))');
                    axis 'xy'
                    if any(m=='c')
                        colorbar;
                    end
                    if nrow*ncol-i<ncol
                        xlabel(axlab(izs(1)));
                    end
                    if rem(i-1,ncol)==0
                        ylabel(axlab(izs(2)));
                    end
                    ia=ib;
                end
            else                    % graph
                plot(qax{izs(1)},reshape(permute(pg,izs),szs(1:2)),'-');
                xlabel(axlab{izs(1)});
            end
        else
            if any(m=='G')  % image
                imagesc(pg');
                axis 'xy'
                if any(m=='c')
                    colorbar;
                end
                xlabel('Input SNR (dB)');
                ylabel('Model Index');
            else            % graph
                if nx>=nq
                    plot(x,pg,'-');
                    xlabel('Input SNR (dB)');
                else
                    plot(1:nq,pg','-');
                    xlabel('Model Index');
                end
            end
        end
    end
end

author: ImperialCollegeLondon

repo: sap-voicebox