text:
function [X, y, XTest, yTest] = mappingLoadData(dataset, seedVal)

% MAPPINGLOADDATA Load a regression or classification dataset.

if nargin < 2
  seedVal = 1e5;
end
randn('seed', seedVal)
rand('seed', seedVal)
XTest = [];
yTest = [];
switch dataset
 %/~
 case 'pumadynSeeger'
  data = load('Dataset.data');
  ind = randperm(size(data, 1));
  indTr = ind(1:7168);
  indTe = ind(7169:end);
  X = data(indTr, 1:end-1);
  y = data(indTr, end);
  XTest = data(indTe, 1:end-1);
  yTest = data(indTe, end);
  Xscale = sqrt(var(X));
  for j = 1:length(Xscale);
    X(:, j) = X(:, j)/Xscale(j);
    XTest(:, j) = XTest(:, j)/Xscale(j);
  end
  augX = [X ones(size(X, 1), 1)];
  w = inv(augX'*augX)*augX'*y;
  y = y -augX*w;
  augTestX = [XTest ones(size(XTest, 1), 1)];
  yTest = yTest - augTestX*w;
  yscale = sqrt(var(y));
  y = y/yscale;
  yTest = yTest/yscale;
 case 'pumadyn'

  % Data is variance 1, no need to normalise.
  data = load('Dataset.data');
  ind = randperm(size(data, 1));
  indTr = ind(1:7168);
  indTe = ind(7169:end);
  X = data(indTr, 1:end-1);
  y = data(indTr, end);
  XTest = data(indTe, 1:end-1);
  yTest = data(indTe, end);
  %~/
 case 'usps'
  load usps_train
  X = ALL_DATA;
  range =  min(ALL_T):max(ALL_T);
  for i = 1:length(range)
    y(:, i) = (ALL_T == range(i))*2 - 1;
  end
  if nargout > 2
    load usps_test
    XTest = ALL_DATA;
    range =  min(ALL_T):max(ALL_T);
    for i = 1:length(range)
      yTest(:, i) = (ALL_T == range(i))*2 - 1;
    end
  end

 case {'usps0', 'usps1', 'usps2', 'usps3', 'usps4', 'usps5', 'usps6', 'usps7', 'usps8', 'usps9'}
  digitNo = str2num(dataset(end));
  load usps_train
  X = ALL_DATA;
  range =  min(ALL_T):max(ALL_T);
  for i = 1:length(range)
    y(:, i) = (ALL_T == range(i))*2 - 1;
  end
  if nargout > 2
    load usps_test
    XTest = ALL_DATA;
    range =  min(ALL_T):max(ALL_T);
    for i = 1:length(range)
      yTest(:, i) = (ALL_T == range(i))*2 - 1;
    end
  end
  y = y(:, digitNo+1);
  yTest = yTest(:, digitNo+1);

 case 'regressionOne'
  try
    load regressionOneData.mat
  catch

    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile')
      numIn = 2;
      N = 500;
      X = zeros(N, numIn);
      X(1:floor(N/2), :) = ...
          randn(floor(N/2), numIn)*.5+1;
      X(floor(N/2)+1:end, :) = ...
          randn(ceil(N/2), numIn)*.5-1;
      kern = kernCreate(X, 'rbfard');
      kern.variance = 1;
      kern.inverseWidth = 20;
      kern.inputScales = [0 0.999];

      K = kernCompute(kern, X);
      y = real(gaussSamp(K, 1)') + randn(N, 1)*0.01;

      save('regressionOneData.mat', 'numIn', 'N', 'X', 'y')
    else
      error(lasterr);
    end

  end

 case 'regressionTwo'
  try
    load regressionTwoData.mat
  catch

    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile')
      numIn = 2;
      N = 500;
      X = zeros(N, numIn);
      X(1:floor(N/2), :) = ...
          randn(floor(N/2), numIn)*.5+1;
      X(floor(N/2)+1:end, :) = ...
          randn(ceil(N/2), numIn)*.5-1;
      kern = kernCreate(X, 'rbfard');
      kern.variance = 1;
      kern.inverseWidth = 20;
      kern.inputScales = [0.999 .2];

      K = kernCompute(kern, X);
      y = real(gaussSamp(K, 1)') + randn(N, 1)*0.01;

      save('regressionTwoData.mat', 'numIn', 'N', 'X', 'y')
    else
      error(lasterr)
    end

  end

 case 'regressionThree'
  try
    load regressionThreeData.mat
  catch

    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile')
      numIn = 2;
      N = 500;
      X = zeros(N, numIn);
      X(1:floor(N/2), :) = ...
          randn(floor(N/2), numIn)*.5+1;
      X(floor(N/2)+1:end, :) = ...
          randn(ceil(N/2), numIn)*.5-1;
      kern = kernCreate(X, 'lin');
      kern.variance = 1;

      K = kernCompute(kern, X);
      y = real(gaussSamp(K, 1)') + randn(N, 1)*0.01;
      save('regressionThreeData.mat', 'numIn', 'N', 'X', 'y')
    else
      error(lasterr);
    end
  end

 case 'regressionFour'
  try
    load regressionFourData.mat
  catch

    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile')
      numIn = 2;
      N = 500;
      X = zeros(N, numIn);
      X(1:floor(N/2), :) = ...
          randn(floor(N/2), numIn)*.5+1;
      X(floor(N/2)+1:end, :) = ...
          randn(ceil(N/2), numIn)*.5-1;
      kern = kernCreate(X, 'mlp');
      kern.variance = 1;
      kern.weightVariance = 1;
      kern.biasVariance = 1;
      K = kernCompute(kern, X);
      y = real(gaussSamp(K, 1)') + randn(N, 1)*0.01;
      save('regressionFourData.mat', 'numIn', 'N', 'X', 'y')
    else
      error(lasterr);
    end
  end
 case 'classificationOne'
  try
    load classificationOneData.mat
  catch
    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile')
      X = [randn(100,2)-[zeros(100, 1) 6*ones(100, 1)]; randn(100,2)+[zeros(100, 1) 6*ones(100, 1)]; randn(100, 2)];
      y = [ones(200, 1); -ones(100, 1)];
      save('classificationOneData.mat', 'X', 'y')
    else
      error(lasterr);
    end
  end
 case 'classificationTwo'

  try
    load classificationTwoData.mat
  catch
    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile')
      numIn = 2;
      N = 500;

      X = zeros(N, numIn);
      X = rand(N, numIn);

      kern = kernCreate(X, 'rbf');
      kern.variance = 10;
      kern.inverseWidth = 10;

      K = kernCompute(kern, X);
      u = real(gaussSamp(K, 1)');

      p = cumGaussian(u);
      y = 2*(rand(size(u))>p)-1;
      save('classificationTwoData.mat', 'numIn', 'N', 'X', 'u', 'y')
    else
      error(lasterr);
    end

  end

 case 'classificationThree'

  try
    load classificationThreeData.mat
  catch
    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile')
      numIn = 2;
      N = 500;

      X = zeros(N, numIn);
      X = rand(N, numIn);

      kern = kernCreate(X, 'rbf');
      kern.variance = 10;
      kern.inverseWidth = 10;

      K = kernCompute(kern, X);
      u = real(gaussSamp(K, 1)');
      a = 3;
      pMinus = cumGaussian(u-a/2);
      pPlus = cumGaussian(-u-a/2);
      p =rand(size(u));
      indMinus = find(p<pMinus);
      indPlus = find(p>pMinus & p<pMinus+pPlus);
      indNone = find(p>pMinus+pPlus);
      y = zeros(N, 1);
      y(indPlus) = 1;
      y(indMinus) = -1;
      y(indNone, :) = [];
      X(indNone, :) = [];
      save('classificationThreeData.mat', 'numIn', 'N', 'X', 'u', 'y')
    else
      error(lasterr);
    end

  end

 case 'orderedOne'
  dataPerCat = 30;
  spacing = 3;

  % Generate a toy data-set of (linear) ordered categories.
  X = [randn(dataPerCat,2)-[zeros(dataPerCat, 1) repmat(3*spacing, dataPerCat, 1)]; ...
       randn(dataPerCat,2)-[zeros(dataPerCat, 1) repmat(2*spacing, dataPerCat, 1)]; ...
       randn(dataPerCat,2)-[zeros(dataPerCat, 1) repmat(spacing, dataPerCat, 1)]; ...
       randn(dataPerCat, 2); ...
       randn(dataPerCat,2)+[zeros(dataPerCat, 1) repmat(spacing, dataPerCat, 1)]; ...
       randn(dataPerCat,2)+[zeros(dataPerCat, 1) repmat(2*spacing, dataPerCat, 1)]; ...
       randn(dataPerCat,2)+[zeros(dataPerCat, 1) repmat(3*spacing, dataPerCat, 1)]];
  y = [zeros(dataPerCat, 1); ...
       repmat(1, dataPerCat, 1); repmat(2, dataPerCat, 1); ...
       repmat(3, dataPerCat, 1); repmat(4, dataPerCat, 1); ...
       repmat(5, dataPerCat, 1); repmat(6, dataPerCat, 1)];

 case 'orderedTwo'
  dataPerCat = 30;
  spacing = 3;

  % Generate a toy data-set of (linear) ordered categories.
  thetaR = [randn(dataPerCat,2)-[zeros(dataPerCat, 1) repmat(3*spacing, dataPerCat, 1)]; ...
            randn(dataPerCat,2)-[zeros(dataPerCat, 1) repmat(2*spacing, dataPerCat, 1)]; ...
            randn(dataPerCat,2)-[zeros(dataPerCat, 1) repmat(spacing, dataPerCat, 1)]; ...
            randn(dataPerCat, 2); ...
            randn(dataPerCat,2)+[zeros(dataPerCat, 1) repmat(spacing, dataPerCat, 1)]; ...
            randn(dataPerCat,2)+[zeros(dataPerCat, 1) repmat(2*spacing, dataPerCat, 1)]; ...
            randn(dataPerCat,2)+[zeros(dataPerCat, 1) repmat(3*spacing, dataPerCat, 1)]];
  thetaR(:, 1) = thetaR(:, 1);
  thetaR(:, 2) = thetaR(:, 2) + 15;
  X = [sin(thetaR(:, 1)).*thetaR(:, 2) cos(thetaR(:, 1)).*thetaR(:, 2)];
  y = [zeros(dataPerCat, 1); ...
       repmat(1, dataPerCat, 1); repmat(2, dataPerCat, 1); ...
       repmat(3, dataPerCat, 1); repmat(4, dataPerCat, 1); ...
       repmat(5, dataPerCat, 1); repmat(6, dataPerCat, 1)];

end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function functionStruct = agilentMSO6014A(visaRsrcString)

iObj = visa('agilent',visaRsrcString);
set(iObj,'InputBufferSize',2048);
scopeObj = icdevice('AGNMSO6014A',iObj);
connect(scopeObj);
fprintf(scopeObj.Interface,'*RST');
functionStruct.setTimeBase = @setTimeBase;
functionStruct. getTimeBase = @getTimeBase;
functionStruct.setVerticalScale = @setVerticalScale;
functionStruct.getVerticalScale = @getVerticalScale;
functionStruct.getWaveform = @getWaveform;
functionStruct.autoscale = @autoscale;
functionStruct.close = @close;

    function setTimeBase(timebase)
        set(scopeObj,'timebase',timebase);
    end

    function timebase = getTimeBase
        timebase = get(scopeObj,'timebase');
    end

    function setVerticalScale(scale)
        set(scopeObj,'verticalscale',scale);
    end

    function scale = getVerticalScale
        scale = get(scopeObj,'verticalscale');
    end

    function [yData xData] = getWaveform
        [yData xData] = invoke(scopeObj,'getwaveform');
    end

    function autoscale
        invoke(scopeObj,'autoscale');
    end

    function close
        disconnect(scopeObj);
        delete(scopeObj);
        delete(iObj);
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function agm_values_test ( )

%*****************************************************************************80
%
%% AGM_VALUE_TEST demonstrates the use of AGM_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'AGM_VALUE_TEST:\n' );
  fprintf ( 1, '  AGM_VALUES stores values of \n' );
  fprintf ( 1, '  the arithmetic geometric mean function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A           B         AGM(A,B)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, fx ] = agm_values ( n_data );

    if ( n_data == 0 )
      break
    end

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

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function EEG = hlp_laplaceRef(EEG,varargin)
% Perform laplacian re-referecing using the spherical spline method of Perrin et al 1989
% This computes scalp current density or current source density
%
% This function requires the CSD toolbox by Jorgen Kayser
% (http://psychophysiology.cpmc.columbia.edu/Software/CSDtoolbox)
%
% Author: Tim Mullen, SCCN/INC/UCSD Dec, 2013

% check if CSD toolbox exists
if ~exist('CSD','file') || ~exist('GetGH.m','file')
    error('CSD toolbox not installed.  Please download from [http://psychophysiology.cpmc.columbia.edu/Software/CSDtoolbox] and place in Matlab path');
end

% parse inputs
g = finputcheck(varargin, ...
            {'sphSplineFlex' 'real' {} 4    ...
             'headRadius',   'real' {} 10   ...
             'smoothLambda', 'real' {} 1.0e-5 ...
            },'ignore','quiet');

% set up paths
locsdir = [tempdir() 'tmplocs.locs'];
csddir  = [tempdir() 'tmplocs.csd'];

% write the chanlocs to a temp location
writelocs(EEG.chanlocs,locsdir,'filetype','loc');

% convert chanlocs to CSD toolbox format
evalstr = sprintf(['ConvertLocations(' hlp_tostring(locsdir) ',' hlp_tostring(csddir) ');']);
evalc(evalstr);

% read locations
[~,M] = evalc(['ExtractMontage(' hlp_tostring(csddir) ',' hlp_tostring({EEG.chanlocs.labels}') ');']);

% clean up
delete(locsdir);
delete(csddir);

% make G and H matrices
[G H] = GetGH(M,g.sphSplineFlex);

% compute current source density
EEG.data = double(CSD(single(EEG.data),G,H));

repo: REST


text:
function str = regexpsplit(varargin)
% mimics regexp(a,b,'split') as introduced in MATLAB version

[a,b] = regexp(varargin{:});

a = [a-1,length(varargin{1})];
b = [1,b+1];

str = arrayfun(@(i) varargin{1}(b(i):a(i)),1:length(a),'uniformoutput',0);
ind = cellfun('isempty',str);
str(ind) = {''};
% if any(ind),  str{ind} = ''; end

author: mtex-toolbox

repo: mtex

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


text:
function x = r83t_cg ( n, a, b, x )

%*****************************************************************************80
%
%% R83T_CG uses the conjugate gradient method on an R83T system.
%
%  Discussion:
%
%    The R83T storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1:N-1,3), the diagonal in
%    entries (1:N,2), and the subdiagonal in (2:N,1).  Thus, the
%    original matrix is "collapsed" horizontally into the array.
%
%    The matrix A must be a positive definite symmetric band matrix.
%
%    The method is designed to reach the solution after N computational
%    steps.  However, roundoff may introduce unacceptably large errors for
%    some problems.  In such a case, calling the routine again, using
%    the computed solution as the new starting estimate, should improve
%    the results.
%
%  Example:
%
%    Here is how an R83T matrix of order 5 would be stored:
%
%       *  A11 A12
%      A21 A22 A23
%      A32 A33 A34
%      A43 A44 A45
%      A54 A55  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Frank Beckman,
%    The Solution of Linear Equations by the Conjugate Gradient Method,
%    in Mathematical Methods for Digital Computers,
%    edited by John Ralston, Herbert Wilf,
%    Wiley, 1967,
%    ISBN: 0471706892,
%    LC: QA76.5.R3.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A(N,3), the matrix.
%
%    Input, real B(N), the right hand side vector.
%
%    Input/output, real X(N).
%    On input, an estimate for the solution, which may be 0.
%    On output, the approximate solution vector.
%
  b = b(:);
  x = x(:);
%
%  Initialize
%    AP = A * x,
%    R  = b - A * x,
%    P  = b - A * x.
%
  ap = r83t_mv ( n, n, a, x );

  r(1:n,1) = b(1:n,1) - ap(1:n,1);
  p(1:n,1) = b(1:n,1) - ap(1:n,1);
%
%  Do the N steps of the conjugate gradient method.
%
  for it = 1 : n
%
%  Compute the matrix*vector product AP=A*P.
%
    ap = r83t_mv ( n, n, a, p );
%
%  Compute the dot products
%    PAP = P*AP,
%    PR  = P*R
%  Set
%    ALPHA = PR / PAP.
%
    pap = p' * ap;
    pr = p' * r;

    if ( pap == 0.0 )
      return
    end

    alpha = pr / pap;
%
%  Set
%    X = X + ALPHA * P
%    R = R - ALPHA * AP.
%
    x(1:n,1) = x(1:n,1) + alpha * p(1:n,1);
    r(1:n,1) = r(1:n,1) - alpha * ap(1:n,1);
%
%  Compute the vector dot product
%    RAP = R*AP
%  Set
%    BETA = - RAP / PAP.
%
    rap = r' * ap;

    beta = - rap / pap;
%
%  Update the perturbation vector
%    P = R + BETA * P.
%
    p(1:n,1) = r(1:n,1) + beta * p(1:n,1);

  end

  return
end


text:
function [logp, yhat, res] = tapas_unitsq_sgm_mu3(r, infStates, ptrans)
% Calculates the log-probability of response y=1 under the unit-square sigmoid model
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Initialize returned log-probabilities as NaNs so that NaN is
% returned for all irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Weed irregular trials out from inferred states and responses
mu1hat = infStates(:,1,1);
mu1hat(r.irr) = [];
mu3hat = infStates(:,3,1);
mu3hat(r.irr) = [];
y = r.y(:,1);
y(r.irr) = [];

% Decision temperature is exponential of log-volatility
% (i.e., inverse decision temperature is exponential of negative log-volatility)
ze = exp(-mu3hat);

% Avoid any numerical problems when taking logarithms close to 1
x = mu1hat;
logx = log(x);
log1pxm1 = log1p(x-1);
logx(1-x<1e-4) = log1pxm1(1-x<1e-4);
log1mx = log(1-x);
log1pmx = log1p(-x);
log1mx(x<1e-4) = log1pmx(x<1e-4);

% Calculate log-probabilities for non-irregular trials
reg = ~ismember(1:n,r.irr);
logp(reg) = y.*ze.*(logx -log1mx) +ze.*log1mx -log((1-x).^ze +x.^ze);
yhat(reg) = x;
res(reg) = (y-x)./sqrt(x.*(1-x));

return;

author: translationalneuromodeling

repo: tapas


text:
function [gKern1, gKern2, gVarmeans, gVarcovars, gInd] = kernkernVardistPsi2Gradient(kern1, kern2, vardist, Z, covGrad, learnInducing)

% KERNKERNVARDISTPSI2GRADIENT Description

% VARGPLVM

if nargin < 6
    learnInducing = true;
end

if strcmp(kern2.type, kern1.type)
    error('Kern1 and Kern2 cannot be the same')
    return;
elseif ~strcmp(kern1.type, 'white') &  ~strcmp(kern2.type, 'white')
    %
    swap = 0;
    if strcmp(kern2.type, 'rbfard2')
       swap = 1;
    elseif strcmp(kern2.type, 'linard2') & ~strcmp(kern1.type, 'rbfard2')
       swap = 1;
    end
    %
    if swap == 0
       fhandle = str2func([kern1.type kern2.type 'VardistPsi2Gradient']);
       [gKern1, gKern2, gVarmeans, gVarcovars, gInd] = fhandle(kern1, kern2, vardist, Z, covGrad, learnInducing);
    else
       fhandle = str2func([kern2.type kern1.type 'VardistPsi2Gradient']);
       [gKern2, gKern1, gVarmeans, gVarcovars, gInd] = fhandle(kern2, kern1, vardist, Z, covGrad, learnInducing);
    end
    % !!! Do not take transformation here: they applied in the kernkernVardistPsi2Gradient
    %gKern1 = paramTransformPsi2(kern1, gKern1);
    %gKern2 = paramTransformPsi2(kern2, gKern2);
    % variational variances are positive
    %gVarcovars = (gVarcovars(:).*vardist.covars(:))';
   %
else
   %
    % the white kernel gives zeros
    gKern1 = zeros(1,kern1.nParams);
    gKern2 = zeros(1,kern2.nParams);
    gVarmeans = zeros(1,prod(size(vardist.means)));
    gInd = zeros(1,prod(size(Z)));
    gVarcovars = zeros(1,prod(size(vardist.covars)));
   %
end

%%
%-----------------------------------------------------
% This applies transformations
% This must be done similar to kernGradient at some point
function gKern = paramTransformPsi2(kern, gKern)
%
%

if strcmp(kern.type,'rbfard2')
   gKern(1) = gKern(1)*kern.variance;
   gKern(2:end) = gKern(2:end).*kern.inputScales;
elseif strcmp(kern.type,'linard2')
   gKern(1:end) = gKern(1:end).*kern.inputScales;
elseif strcmp(kern.type,'bias')
   gKern = gKern*kern.variance;
end

repo: vargplvm


text:
function y = lnCumGaussSum(u1, u2, w1, w2)

% LNCUMGAUSSSUM The log of the weighted sum of two cumulative Gaussians.
% FORMAT
% DESC returns the logarithm of the weighted sum of two cumulative
% Gaussians.
% ARG u1 : argument of the first cumulative Gaussian.
% ARG u2 : argument of the second cumulative Gaussian.
% ARG w1 : weight of the first cumulative Gaussian.
% ARG w2 : weight of the second cumulative Gaussian.
%
% SEEALSO : cumGaussian, lnCumGaussian, lnDiffCumGaussian
%
% COPYRIGHT : Neil D. Lawrence, 2004

% NDLUTIL

y = zeros(size(u1));
safeCond = u1 > 0 & u2 > 0;
index = find(safeCond);
if ~isempty(index)
  y(index) = log(w1.*cumGaussian(u1(index)) ...
                 + w2*cumGaussian(u2(index)));
end
index = find(~safeCond & u1>u2);
if ~isempty(index)
  y(index) = log(w1) + lnCumGaussian(u1(index))...
      + log(1 + w2/w1*exp(lnCumGaussian(u2(index))...
                          -lnCumGaussian(u1(index))));
end
index = find(~safeCond & u2>=u1);
if ~isempty(index)
  y(index) = log(w2) + lnCumGaussian(u2(index))...
      + log(1 + w1/w2*exp(lnCumGaussian(u1(index))...
                          -lnCumGaussian(u2(index))));
end


text:
classdef FLEA < ALGORITHM
% <multi/many> <real> <large/none>
% Fast sampling based evolutionary algorithm

%------------------------------- Reference --------------------------------
% L. Li, C. He, R. Cheng, H. Li, L. Pan, and Y. Jin, A fast sampling based
% evolutionary algorithm for million-dimensional multiobjective
% optimization, Swarm and Evolutionary Computation, 2022, 75: 101181.
%------------------------------- 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".
%--------------------------------------------------------------------------

% This function is written by Cheng He & Lianghao Li

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            Population = Problem.Initialization();
            RN         = ceil(sqrt(Problem.N));	% Number of reference solutions

            %% Optimization
            while Algorithm.NotTerminated(Population)
                FrontNo          = NDSort(Population.objs,Population.cons,inf);
                [Ref,Population] = RefSelection(Problem,Population,RN,(Problem.FE/Problem.maxFE)^2);	% APD_Based_Selection
                [id_obj,id_dec]  = Neighborhood_Association(Population,Ref);
                cm               = 0.1*floor(10*(Problem.FE/Problem.maxFE));	% Ratio of convergence population
                dirV_obj         = Direction_Calculation(Problem,Population,Ref,id_obj,FrontNo);
                dirV_dec         = Direction_Calculation(Problem,Population,Ref,id_dec,FrontNo);
                Offspring        = Reproduction(Problem,Ref,cm,dirV_obj,dirV_dec);
                [Population,~,~] = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

function [Ref,Population] = RefSelection(Problem,Population,RN,theta)
% The environmental selection of RVEA
    [V,~] = UniformPoint(RN,Problem.M);
    V(1:RN,:) = V.*repmat(max(Population.objs,[],1)-min(Population.objs,[],1)+eps,size(V,1),1);
    PopObj = Population.objs;
    [N,M]  = size(PopObj);
    NV     = size(V,1);
    PopObj = PopObj - repmat(min(PopObj,[],1),N,1);
    cosine = 1 - pdist2(V,V,'cosine');
    cosine(logical(eye(length(cosine)))) = 0;
    gamma  = min(acos(cosine),[],2);

    %% Associate each solution to a reference vector
    Angle  = acos(1-pdist2(PopObj,V,'cosine'));
    Next   = zeros(1,NV);
    for i  = 1 : NV
        APD = (1+M*theta*Angle(:,i)/gamma(i)).*sqrt(sum(PopObj.^2,2));
        [~,Next(i)] = min(APD);
        Angle(Next(i),:) = inf;
    end
    % Population for next generation
    Ref    = Population(Next);
    [~,order1] = sort(Ref.objs);
    Ref    = Ref(order1(:,1));
end

function [OCid,DCid] = Neighborhood_Association(Population,Ref)
    % Neighborhood in objective space
    RefObj   = Ref.objs; PopObj = Population.objs;
    Distance_Obj   = pdist2(RefObj,PopObj,'chebychev');
    [~,OCid] = min(Distance_Obj);
    % Neighborhood in decision space
    RefDec = Ref.decs; PopDec = Population.decs;
    Distance_Dec = pdist2(RefDec,PopDec,'chebychev');
    [~,DCid] = min(Distance_Dec);
end

function dirV = Direction_Calculation(Problem,Population,Ref,id,FrontNo)
    %Direction of convergence
    RefDec   = Ref.decs;
    RN = length(Ref);
    dirV = zeros(size(RefDec));
    for i = 1:RN
        Neighbor = Population(id==i);
        F_Neighbor = FrontNo(id==i);
        if ~isempty(Neighbor) && min(F_Neighbor)>1
            dirV(i,:) = mean(Neighbor(F_Neighbor==min(F_Neighbor)).decs,1)-RefDec(i,:);
        elseif rand>0.5
            dirV(i,:)  = (RefDec(i,:)-Problem.lower)./Problem.D;
        else
            dirV(i,:)  = (Problem.upper-RefDec(i,:))./Problem.D;
        end
    end
end

function Offspring = Reproduction(Problem,Ref,cm,dirV_obj,dirV_dec)
    RefDec   = Ref.decs;
    RN = length(Ref);
    CN = ceil(cm*0.5*Problem.N/RN); %Size of convergence population
    DN = Problem.N-CN*RN*2; %Size of diversity population
    OffDec = zeros(CN*RN*2+DN,Problem.D);
    for i = 1:RN
        mu = repmat(RefDec(i,:),CN,1);
        %% Convergence offsprings
        OffDec(1+(i-1)*CN:i*CN,:) = mu+repmat(randn(1,CN)',1,Problem.D).*repmat(dirV_obj(i,:),CN,1);
        OffDec(CN*RN+1+(i-1)*CN:CN*RN+i*CN,:) = mu+repmat(randn(1,CN)',1,Problem.D).*repmat(dirV_dec(i,:),CN,1);
    end
    %% Diversity offsprings
    x_1 = randi(RN,1,DN);x_2 = randi(RN,1,DN);
    dirV  = Ref(x_1).decs-Ref(x_2).decs;
    OffDec(Problem.N-DN+1:end,:) = Ref(x_1).decs + randn(1,DN)'.*dirV;
    %% Evaluation
    Offspring = Problem.Evaluation(OffDec);
end

function [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Population,N)
% The environmental selection of NSGA-II
    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function pixelMeasurements = genFeatureMeasurements(T_wCam_GT, landmarks_w, K, simSetup)
%Will return -1 if the feature cannot be seen

camRes = simSetup.cameraResolution;

pixelMeasurements = zeros(2, size(T_wCam_GT, 3));

% Extract all viewable measurements
for step_i = 1:size(T_wCam_GT,3)

        %Transform and project landmarks into the camera frame
        landmarks_cam = homo2cart(inv(T_wCam_GT(:,:,step_i))*cart2homo(landmarks_w));
        pixels = homo2cart(K*landmarks_cam);

        %Determine which landmarks are viewable
        isInFieldOfView = pixels(1,:) > 0 & pixels(1,:) < camRes(2) & pixels(2,:) > 0 & pixels(2,:) < camRes(1);
        isInFieldOfView = isInFieldOfView & landmarks_cam(3,:) > 0;

        %viewableLandmarkIds = find(viewableLandmarksIdx);
        %viewableLandmarkIds = viewableLandmarkIds(:);
        %viewableLandmarks = landmarks_cam(:, viewableLandmarksIdx);
        if ~isInFieldOfView
            pixels = [-1; -1];
            n_p = [0; 0];
        else
            n_p = simSetup.pixelNoiseStd*randn(2,size(pixels,2)); %image noise
        end

        pixelMeasurements(:, step_i) =  pixels + n_p;

end

end

repo: MSCKF


text:
function variance = circular_normal_01_variance ( )

%*****************************************************************************80
%
%% CIRCULAR_NORMAL_01_VARIANCE returns the variance of the Circular Normal 01 PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VARIANCE(2), the variance of the PDF.
%
  variance(1:2) = 1.0;

  return
end


text:
% HOMOTRANS - 2D homogeneous transformation of points/lines.
%
% Function to perform a transformation on 2D homogeneous points/lines
% The resulting points/lines are normalised to lie in the z = 1 plane
%
% Usage:
%           t = homoTrans(P,v);
%
% Arguments:
%           P  - 3 x 3 transformation matrix
%           v  - 3 x n matrix of points/lines

%  Peter Kovesi
%  School of Computer Science & Software Engineering
%  The University of Western Australia
%  pk @ csse uwa edu au
%  http://www.csse.uwa.edu.au/~pk
%
%  April 2000

function t = homoTrans(P,v);

t = P*v;
t(1,:) = t(1,:)./t(3,:);   %  Now normalise
t(2,:) = t(2,:)./t(3,:);
t(3,:) = ones(1,size(v,2));

author: christianwengert


text:
function crfw = crfPairwiseTrain(vlab, hlab, pvSP, phSP, adjlist, pE, params, edgelen)

% get f1(y1, y2) = log(P(y1|x)) + log(P(y2|x))
% and f2(y1, y2) = I(y1=y2)*log(P(y1=y2|x)) + I(y1~=y2)*log(P(y1~=y2|x))

npair = 0;
for f = 1:numel(adjlist)
    npair = npair + size(adjlist{f}, 1);
end

f1 = zeros(npair, 49);
f2 = zeros(npair, 49);
f3 = zeros(npair, 49);
lab = zeros(npair, 1);
%nbpair = zeros(npair, 2);

c = 0;
for f = 1:numel(adjlist)
    npairf = size(adjlist{f}, 1);
    pg = [pvSP{f}(:, 1)  repmat(pvSP{f}(:, 2), 1, 5).*phSP{f}  pvSP{f}(:, 3)];
    labf = (vlab{f}==1)*1 + (vlab{f}==2).*(hlab{f}>0).*(hlab{f}+1) + (vlab{f}==3)*7;

    s1 = adjlist{f}(:, 1);
    s2 = adjlist{f}(:, 2);

    edgeperc = zeros(numel(s1), 2);
    nb = zeros(size(pg, 1), 1);
    for k = 1:numel(s1)
        nb(s1(k)) = nb(s1(k))+1;
        nb(s2(k)) = nb(s2(k))+1;
        edgeperc(k, 1) = edgelen{f}(s1(k), s2(k)) / (sum(edgelen{f}(s1(k), :))+sum(edgelen{f}(:, s1(k))));
        edgeperc(k, 2) = edgelen{f}(s1(k), s2(k)) / (sum(edgelen{f}(s2(k), :))+sum(edgelen{f}(:, s2(k))));
    end
    %nbpair(c+1:c+npairf, :) = [nb(s1) nb(s2)];
%    min(1./edgeperc(:, 1) + 1./edgeperc(:, 2))
%    max(1./edgeperc(:, 1) + 1./edgeperc(:, 2))
    % unaries normalize for number of edges, and pairwise normalize for
    % edge length
    for k1 = 1:7
        for k2 = 1:7
            f1(c+1:c+npairf, (k1-1)*7+k2) = log(pg(s1, k1))./nb(s1) + log(pg(s2, k2))./nb(s2);
            kp = params(k1 + (k2-1)*7, :);
            if k1==k2
                f2(c+1:c+npairf, (k1-1)*7+k2) = ...
                    log(kp(3)./(1+exp(-kp(1)-kp(2)*(log(pE{f})-log(1-pE{f}))))).*...
                    (edgeperc(:, 1) + edgeperc(:, 2));
            else
                f3(c+1:c+npairf, (k1-1)*7+k2) = ...
                    log(kp(3)./(1+exp(-kp(1)-kp(2)*(log(pE{f})-log(1-pE{f}))))).*...
                    (edgeperc(:, 1) + edgeperc(:, 2));
            end
        end
    end

    lab(c+1:c+npairf) = (labf(s1)>0 & labf(s2)>0).*((labf(s1)-1)*7 + labf(s2));

    c = c + npairf;
end

ind = find(lab==0);
lab(ind) = [];
f1(ind, :) = [];
f2(ind, :) = [];
f3(ind, :) = [];
%nbpair(ind, :) = [];

f2 = f2 - repmat(log(sum(exp(f2), 2)), [1 size(f2, 2)]);
f3 = f3 - repmat(log(sum(exp(f3), 2)), [1 size(f3, 2)]);

fnc{1} = f1;
fnc{2} = f2;
fnc{3} = f3;
%disp(numel(lab))

[p, err] = avep([1 0 0], fnc, lab);
disp(num2str([p err]))

crfw = fminsearch(@(x) objective(x, fnc, lab), [1 1 1]);

[p, err] = avep(crfw, fnc, lab);
disp(num2str([p err]))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function lpw = objective(w, f, lab)

correctind = [1:numel(lab)]' + (lab-1)*numel(lab);

lpw = 0;
tmpz = 0;
for k = 1:numel(f)
    %if k==1
        lpw = lpw + w(k)*f{k}(correctind);
        tmpz = tmpz + w(k)*f{k};
%     else
%         lpw = lpw + w(k)*f{k}(correctind)./nb(:, 1).^w(end) + w(k)*f{k}(correctind)./nb(:, 2).^w(end);
%         tmpz = tmpz + w(k)*f{k}./repmat(nb(:, 1).^w(end), [1 size(f{k},2)]) + w(k)*f{k}./repmat(nb(:, 2).^w(end), [1 size(f{k},2)]);
    %end
end
z = log(sum(exp(tmpz), 2));

lpw = -(sum(lpw - z));

%disp(num2str([lpw w]))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [p err] = avep(w, f, lab)

correctind = [1:numel(lab)]' + (lab-1)*numel(lab);

lpw = 0;
tmpz = 0;
for k = 1:numel(f)
    %if k==1
        lpw = lpw + w(k)*f{k}(correctind);
        tmpz = tmpz + w(k)*f{k};
%     else
%         lpw = lpw + w(k)*f{k}(correctind)./nb(:, 1).^w(end) + w(k)*f{k}(correctind)./nb(:, 2).^w(end);
%         tmpz = tmpz + w(k)*f{k}./repmat(nb(:, 1).^w(end), [1 size(f{k},2)]) + w(k)*f{k}./repmat(nb(:, 2).^w(end), [1 size(f{k},2)]);
    %end
end
z = log(sum(exp(tmpz), 2));

p = mean(exp(lpw-z));

[tmp, guess] = max(tmpz, [], 2);
err = mean(guess~=lab);

author: Cloud-CV

repo: object-proposals


text:
function out = iimg_reslice(matchto, reslicethis, varargin)
% :Usage:
% ::
%
%     out = iimg_reslice(matchto, reslicethis, varargin)
%
%     out = iimg_reslice(matchto, reslicethis, 'write', 'outname', 'myimg.img')
%
% ..
%    tor wager
%    nov. 06
% ..

    % default flags
    % interp = 2 is trilinear
    flags = struct('interp', 2, 'vox', NaN, 'bb', NaN, 'wrap', [0 0 0], 'preserve', 0);
    % spm options for interpolation and wrapping
    d  = [flags.interp*[1 1 1]' flags.wrap(:)];

    write_vol = 0;
    if any(strcmp(varargin, 'write'))
        write_vol = 1;
    end

    % get output image name
    outname = 'resliced_image.img';
    argno = find(strcmp(varargin, 'outname'));
    if argno
        outname = varargin{argno + 1};
    end

    % get info of image to match to
    volInfo = iimg_read_img(matchto);

    % grid in space of target image
    [x1, x2] = ndgrid(1:volInfo.dim(1), 1:volInfo.dim(2));

    out = zeros(volInfo.dim(1:3));

    for i = 1	% for each image
        % get info for image to reslice
        volInfo2 = spm_vol(reslicethis);
        %[volInfo2, reslicethis] = iimg_read_img(reslicethis);

        if write_vol
            VO = volInfo2;
            VO.fname   = outname;

            switch(spm('Ver'))
                case 'SPM2'
                    VO.dim     = [volInfo.dim(1:3) volInfo2.dim(4)];
                case {'SPM5', 'SPM8'}
                    VO.dt = volInfo.dt;
                    VO.private.dat.fname = outname;

                otherwise
                    error('Unknown SPM version "%s": neuroscientists of the future, fix me!', spm('Ver'));
            end

            VO.mat     = volInfo.mat;
            VO.descrip = 'iimg_reslice - resliced - trilinear';
        end

        % final affine matrix
        % mat is affine mtx of image
        % Target\Object maps target to object
        affinemat = inv(volInfo.mat\volInfo2.mat);

        C = spm_bsplinc(volInfo2, d);

        % for each slice

        for x3 = 1:volInfo.dim(3)
            [tmp, y1, y2, y3] = getmask(affinemat, x1, x2, x3, volInfo2.dim(1:3), flags.wrap);
            out(:, :, x3)              = spm_bsplins(C, y1, y2, y3, d);
        end

        if write_vol
            spm_write_vol(VO, out);
        end
    end
end

function [Mask, y1, y2, y3] = getmask(M, x1, x2, x3, dim, wrp)
    tiny = 5e-2; % From spm_vol_utils.c
    y1   = M(1, 1)*x1+M(1, 2)*x2+(M(1, 3)*x3+M(1, 4));
    y2   = M(2, 1)*x1+M(2, 2)*x2+(M(2, 3)*x3+M(2, 4));
    y3   = M(3, 1)*x1+M(3, 2)*x2+(M(3, 3)*x3+M(3, 4));
    Mask = logical(ones(size(y1)));
    if ~wrp(1), Mask = Mask & (y1 >= (1-tiny) & y1 <= (dim(1)+tiny)); end;
    if ~wrp(2), Mask = Mask & (y2 >= (1-tiny) & y2 <= (dim(2)+tiny)); end;
    if ~wrp(3), Mask = Mask & (y3 >= (1-tiny) & y3 <= (dim(3)+tiny)); end;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [Y, optinf] = bpdngrp(D, S, lambda, mu, g, opt)

% bpdngrp -- Basis Pursuit DeNoising with l2,1 group sparsity
%
%         argmin_x (1/2)||D*x - s||_2^2 + lambda*||x||_1 +
%                  mu * \sum_l ||G_l(x)||_2
%
%         The solution is computed using the ADMM approach (see
%         boyd-2010-distributed for details).
%
% Usage:
%       [Y, optinf] = bpdngrp(D, S, lambda, mu, g, opt)
%
% Input:
%       D           Dictionary matrix
%       S           Signal vector (or matrix)
%       lambda      Regularization parameter
%       mu          l2,1 regularization parameter
%       g           Vector containing index values indicating the
%                   group number for each dictionary element. The
%                   first group index is 1 (0 indicates no group).
%                   Numbers must be contiguous. Overlapping groups
%                   are not supported.
%       opt         Options/algorithm parameters structure (see below)
%
% Output:
%       Y           Dictionary coefficient vector (or matrix)
%       optinf      Details of optimisation
%
%
% Options structure fields:
%   Verbose           Flag determining whether iteration status is displayed.
%                     Fields are iteration number, functional value,
%                     data fidelity term, l1 regularisation term, l2,1
%                     regularisation term, and primal and dual residuals
%                     (see Sec. 3.3 of boyd-2010-distributed). The value of
%                     rho is also displayed if options request that it is
%                     automatically adjusted.
%   MaxMainIter       Maximum main iterations
%   AbsStopTol        Absolute convergence tolerance (see Sec. 3.3.1 of
%                     boyd-2010-distributed)
%   RelStopTol        Relative convergence tolerance (see Sec. 3.3.1 of
%                     boyd-2010-distributed)
%   Y0                Initial value for Y
%   U0                Initial value for U
%   rho               ADMM penalty parameter
%   AutoRho           Flag determining whether rho is automatically updated
%                     (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoRhoPeriod     Iteration period on which rho is updated
%   RhoRsdlRatio      Primal/dual residual ratio in rho update test
%   RhoScaling        Multiplier applied to rho when updated
%   AutoRhoScaling    Flag determining whether RhoScaling value is
%                     adaptively determined (see wohlberg-2015-adaptive). If
%                     enabled, RhoScaling specifies a maximum allowed
%                     multiplier instead of a fixed multiplier.
%   RhoRsdlTarget     Residual ratio targeted by auto rho update policy.
%   StdResiduals      Flag determining whether standard residual definitions
%                     (see Sec 3.3 of boyd-2010-distributed) are used instead
%                     of normalised residuals (see wohlberg-2015-adaptive)
%   RelaxParam        Relaxation parameter (see Sec. 3.4.3 of
%                     boyd-2010-distributed)
%   AuxVarObj         Flag determining whether objective function is computed
%                     using the auxiliary (split) variable
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-07-10
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 6,
  opt = [];
end
checkopt(opt, defaultopts([]));
opt = defaultopts(opt);

% Set up status display for verbose operation
hstr = 'Itn   Fnc       DFid      l1        l2,1      r         s      ';
sfms = '%4d %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e';
nsep = 64;
if opt.AutoRho,
  hstr = [hstr '   rho   '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.Verbose && opt.MaxMainIter > 0,
  disp(hstr);
  disp(char('-' * ones(1,nsep)));
end

% Start timer
tstart = tic;

% Set up algorithm parameters and initialise variables
rho = opt.rho;
if isempty(rho), rho = 50*lambda+1; end;
[Nr, Nc] = size(D);
Nm = size(S,2);
Nx = Nc*Nm;
Ng = max(g);
DTS = D'*S;
[luL, luU] = factorise(D, rho);
optinf = struct('itstat', [], 'opt', opt);
r = Inf;
s = Inf;
epri = 0;
edua = 0;

% Initialise main working variables
X = [];
if isempty(opt.Y0),
  Y = zeros(Nc,Nm);
else
  Y = opt.Y0;
end
Yprv = Y;
if isempty(opt.U0),
  if isempty(opt.Y0),
    U = zeros(Nc,Nm);
  else
    U = (lambda/rho)*sign(Y);
  end
else
  U = opt.U0;
end

% Main loop
k = 1;
while k <= opt.MaxMainIter && (r > epri | s > edua),

  % Solve X subproblem
  X = linsolve(D, rho, luL, luU, DTS + rho*(Y - U));

  % See pg. 21 of boyd-2010-distributed
  if opt.RelaxParam == 1,
    Xr = X;
  else
    Xr = opt.RelaxParam*X + (1-opt.RelaxParam)*Y;
  end

  % Solve Y subproblem
  Y = shrink_groups(Xr + U, g, Ng, lambda/rho, mu/rho);

  % Update dual variable
  U = U + Xr - Y;

  % Objective function and convergence measures
  if opt.AuxVarObj,
    Jdf = sum(vec(abs(D*Y - S).^2))/2;
    Jl1 = sum(abs(vec(Y)));
    Jl21 = norm21(Y, g, Ng);
  else
    Jdf = sum(vec(abs(D*X - S).^2))/2;
    Jl1 = sum(abs(vec(X)));
    Jl21 = norm21(X, g, Ng);
  end
  Jfn = Jdf + lambda*Jl1 + mu*Jl21;

  nX = norm(X(:)); nY = norm(Y(:)); nU = norm(U(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    r = norm(vec(X - Y));
    s = norm(vec(rho*(Yprv - Y)));
    epri = sqrt(Nx)*opt.AbsStopTol+max(nX,nY)*opt.RelStopTol;
    edua = sqrt(Nx)*opt.AbsStopTol+rho*nU*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    r = norm(vec(X - Y))/max(nX,nY);
    s = norm(vec(Yprv - Y))/nU;
    epri = sqrt(Nx)*opt.AbsStopTol/max(nX,nY)+opt.RelStopTol;
    edua = sqrt(Nx)*opt.AbsStopTol/(rho*nU)+opt.RelStopTol;
  end

  % Record and display iteration details
  tk = toc(tstart);
  optinf.itstat = [optinf.itstat;[k Jfn Jdf Jl1 Jl21 r s epri edua rho tk]];
  if opt.Verbose,
    if opt.AutoRho,
      disp(sprintf(sfms, k, Jfn, Jdf, Jl1, Jl21, r, s, rho));
    else
      disp(sprintf(sfms, k, Jfn, Jdf, Jl1, Jl21, r, s));
    end
  end

  % See wohlberg-2015-adaptive and pp. 20-21 of boyd-2010-distributed
  if opt.AutoRho,
    if k ~= 1 && mod(k, opt.AutoRhoPeriod) == 0,
      if opt.AutoRhoScaling,
        rhomlt = sqrt(r/(s*opt.RhoRsdlTarget));
        if rhomlt < 1, rhomlt = 1/rhomlt; end
        if rhomlt > opt.RhoScaling, rhomlt = opt.RhoScaling; end
      else
        rhomlt = opt.RhoScaling;
      end
      rsf = 1;
      if r > opt.RhoRsdlTarget*opt.RhoRsdlRatio*s, rsf = rhomlt; end
      if s > (opt.RhoRsdlRatio/opt.RhoRsdlTarget)*r, rsf = 1/rhomlt; end
      rho = rsf*rho;
      U = U/rsf;
      if rsf ~= 1,
        [luL, luU] = factorise(D, rho);
      end
    end
  end

  Yprv = Y;
  k = k + 1;

end

% Record run time and working variables
optinf.runtime = toc(tstart);
optinf.X = X;
optinf.Y = Y;
optinf.U = U;
optinf.lambda = lambda;
optinf.rho = rho;

% End status display for verbose operation
if opt.Verbose && opt.MaxMainIter > 0,
  disp(char('-' * ones(1,nsep)));
end

return

function u = vec(v)

  u = v(:);

return

function u = shrink1(v, a)

  u = sign(v).*max(0, abs(v) - a);

return

function U = shrink2_col_vec(V, a)

  % Additional complexity here allows simultaenous shrinkage of a
  % set of column vectors
  n2v = sqrt(sum(V.^2,1));
  n2v(n2v == 0) = 1;
  U = bsxfun(@times, V, max(0, n2v - a)./n2v);

return

function U = shrink21(V, a, b)

  % See wohlberg-2012-local and chartrand-2013-nonconvex
  U = shrink2_col_vec(shrink1(V,  a), b);

return

function U = shrink_groups(V, g, Ng, a, b)

  U = zeros(size(V));
  U(g==0,:) = shrink1(V(g==0,:), a);
  for l = 1:Ng,
    U(g==l,:) = shrink21(V(g==l,:), a, b);
  end

return

function x = norm21(u, g, Ng)

  x = 0;
  for l = 1:Ng,
    x = x + sqrt(sum(u(g==l,:).^2, 1));
  end
  x = sum(x); % In case u is a matrix (i.e. not a column vector)

return

function [L,U] = factorise(A, c)

  [N,M] = size(A);
  % If N < M it is cheaper to factorise A*A' + cI and then use the
  % matrix inversion lemma to compute the inverse of A'*A + cI
  if N >= M,
    [L,U] = lu(A'*A + c*eye(M,M));
  else
    [L,U] = lu(A*A' + c*eye(N,N));
  end

return

function x = linsolve(A, c, L, U, b)

  [N,M] = size(A);
  if N >= M,
    x = U \ (L \ b);
  else
    x = (b - A'*(U \ (L \ (A*b))))/c;
  end

return

function opt = defaultopts(opt)

  if ~isfield(opt,'Verbose'),
    opt.Verbose = 0;
  end
  if ~isfield(opt,'MaxMainIter'),
    opt.MaxMainIter = 1000;
  end
  if ~isfield(opt,'AbsStopTol'),
    opt.AbsStopTol = 0;
  end
  if ~isfield(opt,'RelStopTol'),
    opt.RelStopTol = 1e-4;
  end
  if ~isfield(opt,'Y0'),
    opt.Y0 = [];
  end
  if ~isfield(opt,'U0'),
    opt.U0 = [];
  end
  if ~isfield(opt,'rho'),
    opt.rho = [];
  end
  if ~isfield(opt,'AutoRho'),
    opt.AutoRho = 1;
  end
  if ~isfield(opt,'AutoRhoPeriod'),
    opt.AutoRhoPeriod = 10;
  end
  if ~isfield(opt,'RhoRsdlRatio'),
    opt.RhoRsdlRatio = 1.2;
  end
  if ~isfield(opt,'RhoScaling'),
    opt.RhoScaling = 100;
  end
  if ~isfield(opt,'AutoRhoScaling'),
    opt.AutoRhoScaling = 1;
  end
  if ~isfield(opt,'RhoRsdlTarget'),
    opt.RhoRsdlTarget = 1;
  end
  if ~isfield(opt,'StdResiduals'),
    opt.StdResiduals = 0;
  end
  if ~isfield(opt,'RelaxParam'),
    opt.RelaxParam = 1.8;
  end
  if ~isfield(opt,'AuxVarObj'),
    opt.AuxVarObj = 1;
  end

return

author: thfylsty


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   demo script for surface smoothing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% preparation
% user must add the path of iso2mesh to matlab path list
% addpath('../');

% user need to add the full path to .../iso2mesh/bin directory
% to windows/Linux/Unix PATH environment variable

%% load the sample data
load rat_head.mat

% volimage is a volumetric image such as an X-ray or MRI image
% A,b are registration matrix and vector, respectively
%% perform mesh generation

[node,face]=v2s(volimage,0.5,2,'cgalmesh');

face=face(:,1:3);

p0=min(node);
p1=max(node);

rownum=3;
colnum=4;
figure;
subplot(rownum,colnum,1);

plotmesh(node,face(:,1:3));
if(~isoctavemesh)
	title({'Laplacian+HC Smoothing Test','no smoothing'});
else
	title('Laplacian+HC - no smoothing');
end
axis equal;
set(gca,'xlim',[p0(1),p1(1)],'ylim',[p0(2),p1(2)],'zlim',[p0(3),p1(3)])

%=========================================================
% apply Laplacian+HC smoothing
%=========================================================

n1=node;
for i=1:rownum*colnum-1
  n1=sms(n1,face(:,1:3),1,0.5); % apply Laplacian+HC mesh smoothing
  subplot(rownum,colnum,i+1);
  plotmesh(n1,face(:,1:3));
  title(['iter=' num2str(i)]);
  axis equal;
  set(gca,'xlim',[p0(1),p1(1)],'ylim',[p0(2),p1(2)],'zlim',[p0(3),p1(3)])
end

%=========================================================
% apply Laplacian smoothing
%=========================================================

figure;
subplot(rownum,colnum,1);

plotmesh(node,face(:,1:3));
if(~isoctavemesh)
        title({'Laplacian Smoothing Test','no smoothing'});
else
        title('Laplacian - no smoothing');
end
axis equal;
set(gca,'xlim',[p0(1),p1(1)],'ylim',[p0(2),p1(2)],'zlim',[p0(3),p1(3)])

conn=meshconn(face(:,1:3),size(node,1));

n1=node;
for i=1:rownum*colnum-1
  n1=smoothsurf(n1,[],conn,1,0.5,'laplacian');
  subplot(rownum,colnum,i+1);
  plotmesh(n1,face(:,1:3));
  title(['iter=' num2str(i)]);
  axis equal;
  set(gca,'xlim',[p0(1),p1(1)],'ylim',[p0(2),p1(2)],'zlim',[p0(3),p1(3)])
end

%=========================================================
% apply Low-pass smoothing
%=========================================================

figure;
subplot(rownum,colnum,1);

plotmesh(node,face(:,1:3));
if(~isoctavemesh)
        title({'Low-pass Smoothing Test','no smoothing'});
else
        title('Low-pass - no smoothing');
end
axis equal;
set(gca,'xlim',[p0(1),p1(1)],'ylim',[p0(2),p1(2)],'zlim',[p0(3),p1(3)])

conn=meshconn(face(:,1:3),size(node,1));

n1=node;
for i=1:rownum*colnum-1
  n1=smoothsurf(n1,[],conn,1,0.5,'lowpass');
  subplot(rownum,colnum,i+1);
  plotmesh(n1,face(:,1:3));
  title(['iter=' num2str(i)]);
  axis equal;
  set(gca,'xlim',[p0(1),p1(1)],'ylim',[p0(2),p1(2)],'zlim',[p0(3),p1(3)])
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function r = hessianinit(v)
%HESSIANINIT  Initialization of hessian variable
%
%  x = hessianinit(v)
%
%The independent variable(s) x are identified and initialized with value(s) v.
%  v may be scalar for one independent variable, or vector (or matrix)
%  for several independent variables.
%
%The call
%
%  hessianinit
%
%without input parameters sets the number of independent variables to zero.
%
%Hessians are frequently sparsely populated. Therefore, one can choose to store
%first and second derivative dense or sparse, see "help sparsehessian".

%The hessian of arrays is stored in the 'next dimension'. So y.hx is
%  3-dimensional for a column vector y, and y.dx is 3-dimensional for
%  a gradient row vector y.
%Since Matlab does not support multi-dimensional sparse arrays, y.dx
%  and y.hx is not accessible in those cases for sparse y.
%
%As a simple example of hessians
%
%  u = hessianinit([ -3 ; 3+4i ])
%
%initializes the hessian package to have two independent variables. The
%variable u, a column vector, with values u.x(1) = -3 and u.x(2) = 3+4i
%has first derivatives u(1).dx = [1 0] and u(2).dx = [0 1] and second
%derivatives u(1).hx = u(2).hx = zeros(2).
%If after that the statement
%
%  v = hessian( intval('3.14159_') )
%
%is executed, the value v.x is the interval with left bound 3.14158 and
%right bound 3.14160 (correctly rounded), the first derivative v.dx = [0 0]
%and second derivative v.hx = zeros(2) (v is treated like a constant). Similarly,
%
%  u(2) = -4711
%
%produces u.x = [ -3 ; -4711 ], u.dx = [ 1 0 ; 0 0 ] and does not change u.hx.
%
%Hessians frequently depend only on few variables. Therefore, if not specified
%  otherwise, the first and second derivative are stored sparse for eight and
%  more unknowns (see above).
%To change only the display you may use "full(u)" or "sparse(u)".
%
%For some examples, see demohessian.
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/03/12     S.M. Rump  SlopeSparseArrayDeriv and see removed
%

  if nargin==0
    setappdata(0,'INTLAB_HESSIAN_NUMVAR',0);
    return
  end

  if ~( isa(v,'double') | isa(v,'intval') )
    error('invalid initialization of hessian')
  end
  setappdata(0,'INTLAB_HESSIAN_NUMVAR',prod(size(v)));

  dummy.init = v;
  r = hessian( dummy , 'hessianinit' );

repo: OpenMAS


text:
% MHGIBBSTRANS
%
%  MCMC Metropolis-Hastings transition function that
%  utilizes the Gibbs sampling distribution for proposals.
%  A - The current joint assignment.  This should be
%      updated to be the next assignment
%  G - The network
%  F - List of all factors
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function A = MHGibbsTrans(A, G, F)

% Draw proposed new state from Gibbs Transition distribution
A_prop = GibbsTrans(A, G, F);

% Compute acceptance probability
p_acceptance = 1.0;
lnum = LogProbOfJointAssignment(F,A_prop);
lden = LogProbOfJointAssignment(F,A);

p_acceptance = min(1,exp(lnum-lden));
% Accept or reject proposal
if rand() < p_acceptance
    A = A_prop;
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function value = alnorm ( x, upper )

%*****************************************************************************80
%
%% ALNORM computes the cumulative density of the standard normal distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by David Hill
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    David Hill,
%    Algorithm AS 66:
%    The Normal Integral,
%    Applied Statistics,
%    Volume 22, Number 3, 1973, pages 424-427.
%
%  Parameters:
%
%    Input, real X, is one endpoint of the semi-infinite interval
%    over which the integration takes place.
%
%    Input, logical UPPER, determines whether the upper or lower
%    interval is to be integrated:
%    1 => integrate from X to + Infinity;
%    0 => integrate from - Infinity to X.
%
%    Output, real VALUE, the integral of the standard normal
%    distribution over the desired interval.
%
  a1 = 5.75885480458;
  a2 = 2.62433121679;
  a3 = 5.92885724438;
  b1 = -29.8213557807;
  b2 = 48.6959930692;
  c1 = -0.000000038052;
  c2 = 0.000398064794;
  c3 = -0.151679116635;
  c4 = 4.8385912808;
  c5 = 0.742380924027;
  c6 = 3.99019417011;
  con = 1.28;
  d1 = 1.00000615302;
  d2 = 1.98615381364;
  d3 = 5.29330324926;
  d4 = -15.1508972451;
  d5 = 30.789933034;
  ltone = 7.0;
  p = 0.39894228044;
  q = 0.39990348504;
  r = 0.398942280385;
  utzero = 18.66;

  up = upper;
  z = x;

  if ( z < 0.0 )
    if ( up )
      up = 0;
    else
      up = 1;
    end
    z = - z;
  end

  if ( ltone < z & ( ( ~up ) | utzero < z ) )

    if ( up )
      value = 0.0;
    else
      value = 1.0;
    end

    return

  end

  y = 0.5  * z * z;

  if ( z <= con )

    value = 0.5  - z * ( p - q * y ...
      / ( y + a1 + b1 ...
      / ( y + a2 + b2 ...
      / ( y + a3 ))));

  else

    value = r * exp ( - y ) ...
      / ( z + c1 + d1 ...
      / ( z + c2 + d2 ...
      / ( z + c3 + d3 ...
      / ( z + c4 + d4 ...
      / ( z + c5 + d5 ...
      / ( z + c6 ))))));

  end

  if ( ~up )
    value = 1.0  - value;
  end

  return
end


function [total_boxes points] = detect_face(img,minsize,PNet,RNet,ONet,threshold,fastresize,factor)
	%im: input image
	%minsize: minimum of faces' size
	%pnet, rnet, onet: caffemodel
	%threshold: threshold=[th1 th2 th3], th1-3 are three steps's threshold
	%fastresize: resize img from last scale (using in high-resolution images) if fastresize==true
	factor_count=0;
	total_boxes=[];
	points=[];
	h=size(img,1);
	w=size(img,2);
	minl=min([w h]);
    img=single(img);
	if fastresize
		im_data=(single(img)-127.5)*0.0078125;
    end
    m=12/minsize;
	minl=minl*m;
	%creat scale pyramid
    scales=[];
	while (minl>=12)
		scales=[scales m*factor^(factor_count)];
		minl=minl*factor;
		factor_count=factor_count+1;
	end
	%first stage
	for j = 1:size(scales,2)
		scale=scales(j);
		hs=ceil(h*scale);
		ws=ceil(w*scale);
		if fastresize
			im_data=imResample(im_data,[hs ws],'bilinear');
		else
			im_data=(imResample(img,[hs ws],'bilinear')-127.5)*0.0078125;
		end
		PNet.blobs('data').reshape([hs ws 3 1]);
		out=PNet.forward({im_data});
		boxes=generateBoundingBox(out{2}(:,:,2),out{1},scale,threshold(1));
		%inter-scale nms
		pick=nms(boxes,0.5,'Union');
		boxes=boxes(pick,:);
		if ~isempty(boxes)
			total_boxes=[total_boxes;boxes];
		end
	end
	numbox=size(total_boxes,1);
	if ~isempty(total_boxes)
		pick=nms(total_boxes,0.7,'Union');
		total_boxes=total_boxes(pick,:);
		regw=total_boxes(:,3)-total_boxes(:,1);
		regh=total_boxes(:,4)-total_boxes(:,2);
		total_boxes=[total_boxes(:,1)+total_boxes(:,6).*regw total_boxes(:,2)+total_boxes(:,7).*regh total_boxes(:,3)+total_boxes(:,8).*regw total_boxes(:,4)+total_boxes(:,9).*regh total_boxes(:,5)];
		total_boxes=rerec(total_boxes);
		total_boxes(:,1:4)=fix(total_boxes(:,1:4));
		[dy edy dx edx y ey x ex tmpw tmph]=pad(total_boxes,w,h);
	end
	numbox=size(total_boxes,1);
	if numbox>0
		%second stage
 		tempimg=zeros(24,24,3,numbox);
		for k=1:numbox
			tmp=zeros(tmph(k),tmpw(k),3);
			tmp(dy(k):edy(k),dx(k):edx(k),:)=img(y(k):ey(k),x(k):ex(k),:);
			tempimg(:,:,:,k)=imResample(tmp,[24 24],'bilinear');
		end
        tempimg=(tempimg-127.5)*0.0078125;
		RNet.blobs('data').reshape([24 24 3 numbox]);
		out=RNet.forward({tempimg});
		score=squeeze(out{2}(2,:));
		pass=find(score>threshold(2));
		total_boxes=[total_boxes(pass,1:4) score(pass)'];
		mv=out{1}(:,pass);
		if size(total_boxes,1)>0
			pick=nms(total_boxes,0.7,'Union');
			total_boxes=total_boxes(pick,:);
            total_boxes=bbreg(total_boxes,mv(:,pick)');
            total_boxes=rerec(total_boxes);
		end
		numbox=size(total_boxes,1);
		if numbox>0
			%third stage
			total_boxes=fix(total_boxes);
			[dy edy dx edx y ey x ex tmpw tmph]=pad(total_boxes,w,h);
            tempimg=zeros(48,48,3,numbox);
			for k=1:numbox
				tmp=zeros(tmph(k),tmpw(k),3);
				tmp(dy(k):edy(k),dx(k):edx(k),:)=img(y(k):ey(k),x(k):ex(k),:);
				tempimg(:,:,:,k)=imResample(tmp,[48 48],'bilinear');
			end
			tempimg=(tempimg-127.5)*0.0078125;
			ONet.blobs('data').reshape([48 48 3 numbox]);
			out=ONet.forward({tempimg});
			score=squeeze(out{3}(2,:));
			points=out{2};
			pass=find(score>threshold(3));
			points=points(:,pass);
			total_boxes=[total_boxes(pass,1:4) score(pass)'];
			mv=out{1}(:,pass);
			w=total_boxes(:,3)-total_boxes(:,1)+1;
            h=total_boxes(:,4)-total_boxes(:,2)+1;
            points(1:5,:)=repmat(w',[5 1]).*points(1:5,:)+repmat(total_boxes(:,1)',[5 1])-1;
            points(6:10,:)=repmat(h',[5 1]).*points(6:10,:)+repmat(total_boxes(:,2)',[5 1])-1;
			if size(total_boxes,1)>0
				total_boxes=bbreg(total_boxes,mv(:,:)');
                pick=nms(total_boxes,0.7,'Min');
				total_boxes=total_boxes(pick,:);
                points=points(:,pick);
			end
		end
    end
end


text:
function [tightInset,figTightInset] = calcTightInset(mtexFig)
% determine tight inset for each axis

tightInset = zeros(1,4);
figTightInset = mtexFig.outerPlotSpacing * [1,1,1,1];

%return

if isempty(mtexFig.children), return; end
ax = mtexFig.children(1);

if strcmpi(get(ax,'visible'),'off') || strcmpi(get(ax,'XColor'),'none')

  xtl = get(ax,'xTickLabel');
  ytl = get(ax,'yTickLabel');
  %xl = get(ax,'xLabel');
  %yl = get(ax,'yLabel');
  set(ax,'xTickLabel',[],'yTickLabel',[],'units','pixel');
  tightInset = get(ax,'tightInset');
  %set(ax,'xTickLabel',xtl,'yTickLabel',ytl,'xlabel',xl,'ylabel',yl);
  set(ax,'xTickLabel',xtl,'yTickLabel',ytl);

  % consider text labels
  txt = findall(ax,'type','text','unit','data');
  s = get(txt,'string');
  if ~iscell(s), s = {s}; end
  ind = cellfun(@isempty,s);
  txt = txt(~ind);
  if ~isempty(txt)
    pos = ensurecell(get(txt,'position'));
    set(txt,'unit','pixel')
    ext = cell2mat(ensurecell(get(txt,'extent')));
    %set(txt,'units','data');
    for i=1:length(txt), set(txt(i),'units','data','position',pos{i}); end
    pos = get(ax,'position');
    tightInset(1:2) = max([tightInset(1:2);-ext(:,1:2)]);
    tightInset(3:4) = max([tightInset(3:4);ext(:,1:2)+ext(:,3:4)-repmat(pos(3:4),size(ext,1),1)]);
  end

elseif all(get(ax,'ticklength') == 0)

  xt = get(ax,'xtick');
  yt = get(ax,'ytick');
  set(ax,'xtick',[],'ytick',[]);
  %axis(ax,'normal');
  tightInset = get(ax,'tightInset');
  %axis(ax,'equal','tight');
  set(ax,'xtick',xt,'ytick',yt);

else %if strcmpi(get(ax,'PlotBoxAspectRatioMode'),'auto')

  tightInset = get(ax,'tightInset');

end

 % consider colorbar
 if ~isempty(mtexFig.cBarAxis)

   pos = get(mtexFig.cBarAxis(1),'position');
   pos = pos(3:4);
   pos(pos==max(pos)) = 0;
   fs = get(mtexFig.cBarAxis(1),'FontSize');

   try
     tiPos = get(mtexFig.cBarAxis(1),'tightInset');
     tiPos = tiPos(1:2) + tiPos(3:4);

   catch
     l = get(mtexFig.cBarAxis(1),'Label');
     s = get(l,'String');
     if isempty(s)
       tiPos = [3.5,1.5] * fs;
     elseif get(l,'Rotation') == 0
       tiPos = [3.5,3.5] * fs;
     else
       tiPos = [5.5,1.5] * fs;
     end
   end
   mtexFig.cBarShift = tiPos(pos>0) + fs / 2;

   pos(pos>0) = pos(pos>0) + mtexFig.cBarShift + fs / 2;

   if numel(mtexFig.cBarAxis) == numel(mtexFig.children)
     tightInset = tightInset + [0,pos(2),pos(1),0];
   else
     figTightInset = figTightInset + [0,pos(2),pos(1),0];
   end
 end

 % consider sgtitle
 h = findobj(gcf,'Type','subplottext');
 if ~isempty(h)
   figTightInset(4) = figTightInset(4) + 2.5*h(1).FontSize;
 end

end


text:
function n = norm(T)
%NORM Frobenius norm of a tensor.
%
%   NORM(X) returns the Frobenius norm of a tensor.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

v = reshape(T.data, numel(T.data), 1);
n = norm(v);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function value = year_length_roman ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_ROMAN returns the number of days in a Roman year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of days
%    in the year.
%
  if ( year_is_leap_roman ( y ) )
    value = 366;
  else
    value = 365;
  end

  return
end


text:
function df = ipoptgp_callback_df(x,model)

global latest_x_f
global latest_df

x = x(:);
if isequal(x,latest_x_f)
    df = latest_df;
else
    [f,df] = fmincon_fungp(x,model);
    latest_df = df;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function pass = test_sample( )
% Test spherefun sample() command

tol = 100*chebfunpref().cheb2Prefs.chebfun2eps;

% Function to test
f = spherefun(@(x,y,z) sin(pi*x.*y));

% Ensure the matrix of sampled values is correct.
[m,n] = length(f);
[nn,mm] = size(sample(f));
pass(1) = (m == mm) && (n == nn);

% Sample on fixed grids of various sizes to make sure the right size output
% is given.
m = 120;
n = 121;
[nn,mm] = size(sample(f, m, n));
pass(2) = (m == mm) && (n == nn);

m = 121;
n = 120;
[nn,mm] = size(sample(f, m, n));
pass(3) = (m == mm) && (n == nn);

% Check samples are correct.
% m and n even
m = 30;
n = 20;
[lam,th] = meshgrid(trigpts(m, [-pi, pi]), linspace(0, pi, n));
F = f(lam, th);
G = sample(f, m, n);
pass(4) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, n);
G = U * D * V.';
pass(5) = norm(F(:) - G(:), inf) < tol;

% m even and n odd
m = 30;
n = 21;
[lam, th] = meshgrid(trigpts(m, [-pi, pi]), linspace(0, pi, n));
F = f(lam, th);
G = sample(f, m, n);
pass(6) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, n);
G = U * D * V.';
pass(7) = norm(F(:) - G(:), inf) < tol;

% m odd and n even
m = 31;
n = 20;
[lam, th] = meshgrid(trigpts(m, [-pi, pi]), linspace(0, pi, n));
F = f(lam, th);
G = sample(f, m, n);
pass(8) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, n);
G = U * D * V.';
pass(9) = norm(F(:) - G(:), inf) < tol;

% m odd and n odd
m = 31;
n= 21;
[lam, th] = meshgrid(trigpts(m, [-pi, pi]), linspace(0, pi, n));
F = f(lam, th);
G = sample(f, m, n);
pass(10) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, n);
G = U * D * V.';
pass(11) = norm(F(:) - G(:), inf) < tol;

% Sample should return all ones for the function 1.
f = spherefun(@(x,y,z) 1 + 0*x);
F = sample(f, 128, 128);
pass(12) = norm(F(:) - 1, inf) < tol;

% Check that errors are caught
try
    F = sample(f, 0, 20);
    pass(13) = false;
catch ME
    pass(13) = strcmp(ME.identifier, 'CHEBFUN:SPHEREFUN:sample:inputs');
end

try
    F = sample(f, 20, 0);
    pass(14) = false;
catch ME
    pass(14) = strcmp(ME.identifier, 'CHEBFUN:SPHEREFUN:sample:inputs');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [output_stack, IRR, ORR, ARR] = evalfisx(input_stack, fis)
%EVALFISX Evaluation of a fuzzy inference system.
%
%   See evalfis1 for syntax and explanation.
%
%   It is completely based on MATLAB's evalfis1
%   with some modifications, so it's compatible
%   with extendent fuzzy rule structure.
%
%   Compare also code of this function
%   with code of original the evalfis1.
%
%   Example:
%   load carsmall;
%   fis = genfis4([Weight, Displacement], Acceleration, 'mamdani');
%   out = evalfisx([2000 100; 2000 200; 2000 300], fis);

%   Per Konstantin A. Sidelnikov, 2009.

persistent CURRENT_FIS;
persistent FIS_NAME FIS_TYPE IN_N OUT_N IN_MF_N OUT_MF_N RULE_N;
persistent AND_METHOD OR_METHOD IMP_METHOD AGG_METHOD DEFUZZ_METHOD;
persistent BOUND IN_MF_TYPE OUT_MF_TYPE;
persistent IN_RULE_LIST OUT_RULE_LIST AND_OR RULE_WEIGHT;
persistent IN_PARAM OUT_PARAM;
persistent OUT_TEMPLATE_MF OUT_MF QUALIFIED_OUT_MF OVERALL_OUT_MF;

point_n = 101;
mf_para_n = 4;

initialization = 1;
% Check if initialization necessary.
if isequal(fis, CURRENT_FIS)
    initialization = 0;
end

% Unpack data and initialize global variables.
if initialization
    CURRENT_FIS = fis;
    FIS_NAME = fis.name;
    FIS_TYPE = fis.type;
    IN_N = length(fis.input);
    OUT_N = length(fis.output);

    for i = 1 : IN_N
        IN_MF_N(i) = length(fis.input(i).mf);
    end
    for i = 1 : OUT_N
        OUT_MF_N(i) = length(fis.output(i).mf);
    end

    RULE_N = length(fis.rule);

    AND_METHOD = fis.andMethod;
    OR_METHOD =  fis.orMethod;
    IMP_METHOD = fis.impMethod;
    AGG_METHOD = fis.aggMethod;
    DEFUZZ_METHOD = fis.defuzzMethod;

    IN_MF_TYPE = [];
    for i = 1 : IN_N
        IN_MF_TYPE = [IN_MF_TYPE; {fis.input(i).mf.type}'];
    end
    OUT_MF_TYPE = [];
    for i = 1 : OUT_N
        OUT_MF_TYPE = [OUT_MF_TYPE; {fis.output(i).mf.type}'];
    end

    in_bound = reshape([fis.input.range], IN_N, 2);
    out_bound = reshape([fis.output.range], OUT_N, 2);
    BOUND = [in_bound; out_bound];

    RULE_WEIGHT = [fis.rule.weight]';
    AND_OR = [fis.rule.connection]';
    IN_RULE_LIST = {fis.rule.antecedent}';
    OUT_RULE_LIST = reshape([fis.rule.consequent], RULE_N, OUT_N);

    k = 1;
    totalInputMFs = sum(IN_MF_N);
    totalOutputMFs = sum(OUT_MF_N);
    IN_PARAM = zeros(totalInputMFs, 4);
    for i = 1 : IN_N
        for j = 1 : length(fis.input(i).mf)
            tmp = fis.input(i).mf(j).params;
            IN_PARAM(k, 1 : length(tmp)) = tmp;
            k = k + 1;
        end
    end
    k = 1;
    OUT_PARAM = zeros(totalOutputMFs, 4);
    for i = 1 : OUT_N
        for j = 1 : length(fis.output(i).mf)
            tmp = fis.output(i).mf(j).params;
            OUT_PARAM(k, 1 : length(tmp)) = tmp;
            k = k + 1;
        end
    end

    if strcmp(FIS_TYPE, 'sugeno')
        OUT_PARAM = OUT_PARAM(:, 1 : IN_N + 1);
    elseif strcmp(FIS_TYPE, 'mamdani')
        OUT_PARAM = OUT_PARAM(:, 1 : mf_para_n);
    else
        error('Unknown FIS type!');
    end

    if strcmp(FIS_TYPE, 'mamdani')
        % Compute OUT_TEMPLATE_MF
        OUT_TEMPLATE_MF = zeros(sum(OUT_MF_N), point_n);
        cum_mf = cumsum(OUT_MF_N);
        for i = 1 : sum(OUT_MF_N)
            % index for output
            output_index = find((cum_mf - i) >= 0, 1, 'first');
            tmp = linspace(BOUND(IN_N + output_index, 1), ...
                BOUND(IN_N + output_index, 2), point_n);
            OUT_TEMPLATE_MF(i, :) = ...
                evalmf(tmp, OUT_PARAM(i, :), OUT_MF_TYPE{i});
        end

        % Reorder to fill OUT_MF, an (RULE_N X point_n * OUT_N) matrix.
        OUT_MF = zeros(RULE_N, point_n * OUT_N);
        for i = 1 : RULE_N
            for j = 1 : OUT_N
                mf_index = OUT_RULE_LIST(i, j);
                index = sum(OUT_MF_N(1 : j - 1)) + abs(mf_index);
                tmp = (j - 1) * point_n + 1 : j * point_n;
                if mf_index > 0    % regular MF
                    OUT_MF(i, tmp) = OUT_TEMPLATE_MF(index, :);
                elseif mf_index < 0 % Linguistic hedge "NOT"
                    OUT_MF(i, tmp) = 1 - OUT_TEMPLATE_MF(index, :);
                else                % Don't care (MF index == 0)
                    OUT_MF(i, tmp) = ones(1, point_n);
                end
            end
        end

        % Allocate other matrices
        QUALIFIED_OUT_MF = zeros(RULE_N, point_n * OUT_N);
        OVERALL_OUT_MF = zeros(1, point_n * OUT_N);
    end
    %   fprintf('Global variables for %s FIS are initialized\n', FIS_NAME);
end
% End of initialization

% Error checking for input stack
m = size(input_stack, 1);
n = size(input_stack, 2);
if ~((n == IN_N) || ((n == 1) && (m == IN_N)))
    fprintf('The input stack is of size %dx%d,', m, n);
    fprintf('while expected input vector size is %d.\n', IN_N);
    error('Exiting ...');
end
if (n == 1) && (m == IN_N)
    data_n = 1;
    input_stack = input_stack';
else
    data_n = m;
end

% Allocate output stack
output_stack = zeros(data_n, OUT_N);

% Iteration through each row of input stack
for kkk = 1 : data_n
    input = input_stack(kkk, :);

    % Find in_template_mf_value
    in_template_mf_value = zeros(sum(IN_MF_N), 1);
    cum_mf = cumsum(IN_MF_N);
    for i = 1 : sum(IN_MF_N)
        input_index = find((cum_mf - i) >= 0, 1, 'first');
        in_template_mf_value(i) = ...
            evalmf(input(input_index), IN_PARAM(i, :), IN_MF_TYPE{i});
    end

    % Reordering to fill in_mf_value, which is an (RULE_N X 1) cell matrix.
    tmp = cumsum([0, IN_MF_N(1 : IN_N - 1)]);
    in_mf_value = cell(RULE_N, 1);
    for ind = 1 : RULE_N
        index = tmp(IN_RULE_LIST{ind}(1, :)) + abs(IN_RULE_LIST{ind}(2, :));
        in_mf_value{ind} = in_template_mf_value(index)';
        % Take care of linguistic hedge NOT (MF index is negative)
        neg_index = find(IN_RULE_LIST{ind}(2, :) < 0);
        in_mf_value{ind}(neg_index) = 1 - in_mf_value{ind}(neg_index);
    end

    % Find the firing strengths
    % AND_METHOD = 'min' or 'prod'; which is used as function name too
    % OR_METHOD = 'max' or 'probor'; which is used as function name too
    firing_strength = zeros(RULE_N, 1);
    and_index = find(AND_OR == 1);
    or_index = find(AND_OR == 2);
    if IN_N ~= 1
        firing_strength(and_index) = ...
            cellfun(@(x) feval(AND_METHOD, x'), in_mf_value(and_index, :))';
        firing_strength(or_index) = ...
            cellfun(@(x) feval(OR_METHOD, x'), in_mf_value(or_index, :))';
    else
        firing_strength = [in_mf_value{:}]';
    end

    % Recalculate firing strengths scaled by rule weights
    firing_strength = firing_strength .* RULE_WEIGHT;

    % Find output
    if strcmp(FIS_TYPE, 'sugeno')
        template_output = OUT_PARAM * [input(:); 1]; % Output for template

        % Reordering according to the output part of RULE_LIST;
        % Negative MF index will becomes positive
        tmp = cumsum([0, OUT_MF_N(1 : OUT_N - 1)]);
        index = repmat(tmp, RULE_N, 1) + abs(OUT_RULE_LIST);
        index(index == 0) = 1;   % temp. setting for easy indexing
        rule_output = template_output(index);
        rule_output(OUT_RULE_LIST == 0) = 0;  % take care of zero index
        sum_firing_strength = sum(firing_strength);

        if sum_firing_strength == 0
            fprintf('input = [');
            for i = 1 : IN_N
                fprintf('%f ', input(i));
            end
            fprintf(']\n');
            error('Total firing strength is zero!');
        end

        switch DEFUZZ_METHOD
            case 'wtaver'
                output_stack(kkk, :) = firing_strength' * rule_output / ...
                    sum_firing_strength;
            case 'wtsum'
                output_stack(kkk, :) = firing_strength' * rule_output;
            otherwise
                error('Unknown defuzzification method!');
        end
    elseif strcmp(FIS_TYPE, 'mamdani')
        % Transform OUT_MF to QUALIFIED_OUT_MF
        % Duplicate firing_strength.
        tmp = firing_strength(:, ones(1, point_n * OUT_N));

        if strcmp(IMP_METHOD, 'prod')      % IMP_METHOD == 'prod'
            QUALIFIED_OUT_MF = tmp .* OUT_MF;
        elseif strcmp(IMP_METHOD, 'min')   % IMP_METHOD == 'min'
            QUALIFIED_OUT_MF = feval(IMP_METHOD, tmp, OUT_MF);
        else    % IMP_METHOD is user-defined
            tmp1 = feval(IMP_METHOD, [tmp(:)'; OUT_MF(:)']);
            QUALIFIED_OUT_MF = reshape(tmp1, RULE_N, point_n * OUT_N);
        end

        % AGG_METHOD = 'sum' or 'max' or 'probor' or user-defined
        OVERALL_OUT_MF = feval(AGG_METHOD, QUALIFIED_OUT_MF);

        for i = 1 : OUT_N
            tmp = linspace(BOUND(IN_N + i, 1), BOUND(IN_N + i, 2), point_n);
            tmp1 = (i - 1) * point_n + 1 : i * point_n;
            output_stack(kkk, i) = ...
                defuzz(tmp, OVERALL_OUT_MF(1, tmp1), DEFUZZ_METHOD);
        end
    else
        fprintf('fis_type = %d\n', FIS_TYPE);
        error('Unknown FIS type!');
    end
end

if nargout >= 2
    IRR = in_mf_value;
end

if strcmp(FIS_TYPE, 'sugeno')
    if nargout >= 3
        ORR = rule_output;
    end
    if nargout >= 4
        ARR = firing_strength(:, ones(1, OUT_N));
    end
else
    if nargout >= 3
        ORR = [];
        for iii = 1 : OUT_N
            tmp = (iii - 1) * point_n + 1 : iii * point_n;
            ORR = [ORR; QUALIFIED_OUT_MF(:, tmp)];
        end
        ORR = ORR';
    end
    if nargout >= 4
        ARR = reshape(OVERALL_OUT_MF, point_n, OUT_N);
    end
end


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
function ok=initJava
ok=false;
try
    edu.stanford.facs.swing.StochasticGradientDescent.EPOCH_REPORTS;
catch
    jar=fullfile(fileparts(mfilename('fullpath')), 'umap.jar');
    javaaddpath(jar);
end
try
    edu.stanford.facs.swing.StochasticGradientDescent.EPOCH_REPORTS;
    ok=true;
catch
end


text:
function tapas_rw_binary_dual_plotTraj(r)
% Plots the estimated or generated trajectories for the binary HGF perceptual model for multi-armed
% bandit situations.
%
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_rw_binary_dual_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Optional plotting of responses (true or false)
ploty = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.7*scrsz(4),0.8*scrsz(3),0.3*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Set up colors
colors = [1 0 0; 0.67 0 1; 0 0.67 1; 0.67 1 0];

% Number of bandits
b = 2;

% Number of trials
n = size(r.u,1);

% Time axis
t = ones(1,n);
ts = cumsum(t);
ts = [0, ts];

% Plot
for j=1:b
    plot(ts, [r.p_prc.v_0(j); r.traj.v(:,j)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.v_0(j), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0 0]); % inputs
if (ploty == true) && ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    y = r.y(:,1);
    if ~isempty(find(strcmp(fieldnames(r),'irr')))
        y(r.irr) = NaN; % weed out irregular responses
        plot(ts(r.irr),  1.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses
    end
    for j=1:b
        plot(find(y==j), 1.08*ones([1 length(find(y==j))]), '.', 'Color', colors(j,:)); % responses
    end
    title(['Response y, input u (black), and posterior expectation of reward v ', ...
           'for \alpha=', num2str(r.p_prc.al)], ...
      'FontWeight', 'bold');
    ylabel('y, u, s(\mu_2)');
    axis([0 ts(end) -0.15 1.15]);
else
    title(['Input u (black) and posterior expectation of input s(\mu_2) ', ...
           'for \alpha=', num2str(r.p_prc.al)], ...
      'FontWeight', 'bold');
    ylabel('u, s(\mu_2)');
    axis([0 ts(end) -0.1 1.1]);
end
plot(ts(2:end), 0.5, 'k');
xlabel('Trial number');
hold off;


text:
%% Copyright (C) 2014, 2016, 2018 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 conj (@var{z})
%% Symbolic conjugate function.
%%
%% Examples:
%% @example
%% @group
%% syms z
%% conj(z)
%%   @result{} ans = (sym)
%%       _
%%       z
%%
%% @end group
%% @group
%% syms x real
%% conj(x)
%%   @result{} ans = (sym) x
%%
%% conj(sym(pi) + 6i)
%%   @result{} ans = (sym) π - 6⋅ⅈ
%% @end group
%% @end example
%%
%% Unlike @ref{@@sym/ctranspose}, this command does not transpose
%% a matrix:
%% @example
%% @group
%% A = [1 z x; sym(4) 5 6+7i]
%%   @result{} A = (sym 2×3 matrix)
%%       ⎡1  z     x   ⎤
%%       ⎢             ⎥
%%       ⎣4  5  6 + 7⋅ⅈ⎦
%% conj(A)
%%   @result{} ans = (sym 2×3 matrix)
%%       ⎡   _         ⎤
%%       ⎢1  z     x   ⎥
%%       ⎢             ⎥
%%       ⎣4  5  6 - 7⋅ⅈ⎦
%% @end group
%% @end example
%%
%% @seealso{@@sym/ctranspose, @@sym/real, @@sym/imag}
%% @end defmethod

function z = conj(x)

  if (nargin ~= 1)
    print_usage ();
  end

  % can just use .conjugate on matrix but avoids S.true err
  sf = { 'def _op(x):'
         '    if x in (S.true, S.false):'
         '        return x'
         '    return x.conjugate()' };

  z = elementwise_op (sf, x);

end

%!test
%! a = sym(6);
%! b = sym(5i);
%! assert (isequal (conj(a), a))
%! assert (isequal (conj(b), -b))
%! assert (isequal (conj(a+b), a-b))

%!test
%! syms x
%! assert (isequal (conj(conj(x)), x))

%!test
%! syms x real
%! assert (isequal (conj(x), x))

%!test
%! % array
%! syms x
%! A = [x 6+1i; sym(1) x+2i];
%! B = [conj(x) 6-1i; sym(1) conj(x)-2i];
%! assert (isequal (conj(A), B))

%!test
%! % true/false
%! t = sym(true);
%! f = sym(false);
%! assert (isequal ( conj(t), t))
%! assert (isequal ( conj(f), f))

%!test
%! % round trip
%! syms x
%! d = 3 - 5i;
%! f = conj (x);
%! A = conj (d);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B)

repo: octsympy


text:
function [indexes] = FindCornerSolutions(front)

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

% This function is written by Annibale Panichella

[m,n] = size(front);

%% let's normalize the objectives
if m<=n
  indexes = 1:m;
  return
end

%% let's define the axes of the n-dimensional spaces
W = eye(n);
[r,~]= size(W);
indexes = zeros(1,n);
for i=1:r
   [~, index] = min(Point2LineDistance(front, zeros(1,n), W(i,:)));
   indexes(i) = index;
end

end


text:
% Fig. 8.13   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

w=logspace(-1,2,100);
l=length(w);

% s-plane
num=5;
den=[1 5];
[magC,phaseC]=bode(num,den,w);
subplot(2,1,1),
loglog(w,magC,'-')
axis([.1  100  .01 1]);
ylabel('Magnitude')
%xlabel('Frequency (rad/sec)')
title('Fig. 8.13  (a) 100 rad/sec  T=1/15 sec')
nicegrid
hold on

%              T = 1/15 sec
T = 1/15;

%  MPZ

num=.143*[1 1];
den=[1 -.715];
[mag1,phase]=dbode(num,den,T,w);

%  MMPZ

num=.285;
den=[1 -.715];
[mag2,phase]=dbode(num,den,T,w);

% Tustins

num=.143*[1 1];
den=[1 -.713];
[mag3,phase]=dbode(num,den,T,w);

loglog(w,mag1,'--',w,mag2,'--',w,mag3,'-.')
hold off

%              T = 1/3 sec

subplot(2,1,2),loglog(w,magC,'-')
ylabel('Magnitude')
xlabel('Frequency (rad/sec)')
title('(b) 20 rad/sec  T=1/3 sec')
nicegrid
hold on

T = 1/3;

%  MPZ

num=.405*[1 1];
den=[1 -.189];
wp=w;
wp=wp(wp < 24);
[mag1,phase]=dbode(num,den,T,wp);

%  MMPZ

num=.811;
den=[1 -.189];
[mag2,phase]=dbode(num,den,T,wp);

% Tustins

num=.454*[1 1];
den=[1 -.0914];
[mag3,phase]=dbode(num,den,T,wp);

loglog(wp,mag1,'--',wp,mag2,'--',wp,mag3,'-.')

hold off

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_KUKA_KR5_SIXX_R650(robot, T)
%   Solves the inverse kinematic problem for the KUKA KR5_SIXX_R650 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_KUKA_KR5_SIXX_R650 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('kuka', 'KR5_sixx_R650');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_kuka_kr5_sixx_r650(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=abs(d(6));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist(robot, q(:,i), T, 1,'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'geometric'); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));
A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi - eta;
q3(2) = pi - phi + eta;

repo: ARTE


text:
function [trx,fieldnames,data,units] = ReadStructSVMTrxFile(trxfile)

MAXSTRINGLENGTH = 10000;
fid = fopen(trxfile,'rb');
version = fread(fid,1,'double'); %#ok<NASGU>
moviename_length = fread(fid,1,'int');
if moviename_length > MAXSTRINGLENGTH,
  error('Movie name length = %d > %d\n',moviename_length,MAXSTRINGLENGTH);
end
moviename = fread(fid,moviename_length,'char=>char')';
matname_length = fread(fid,1,'int');
if matname_length > MAXSTRINGLENGTH,
  error('Mat name length = %d > %d\n',matname_length,MAXSTRINGLENGTH);
end
matname = fread(fid,matname_length,'char=>char')'; %#ok<NASGU>
nflies = fread(fid,1,'int'); %#ok<NASGU>
fly_ids = fread(fid,1,'int');
firstframe = fread(fid,1,'int');
firstframe = firstframe + 1;
lastframe = fread(fid,1,'int');
lastframe = lastframe + 1;
sex = fread(fid,1,'char=>char');
fps = fread(fid,1,'double');
nfields = fread(fid,1,'int');

fieldnames = cell(1,nfields);
data = cell(1,nfields);
units = cell(1,nfields);
nframes = lastframe-firstframe+1;

for i = 1:nfields,
  fieldname_length = fread(fid,1,'int');
  if fieldname_length > MAXSTRINGLENGTH,
    error('Field name length = %d > %d\n',fieldname_length,MAXSTRINGLENGTH);
  end
  fieldname = fread(fid,fieldname_length,'char=>char')';
  fieldnames{i} = fieldname;
  units_numerator_length = fread(fid,1,'int');
  if units_numerator_length > MAXSTRINGLENGTH,
    error('Units numerator length = %d > %d\n',units_numerator_length,MAXSTRINGLENGTH);
  end
  units_numerator = fread(fid,units_numerator_length,'char=>char')';
  units_numerator = regexp(strtrim(units_numerator),'\s','split');
  idx = strcmp(units_numerator,'1');
  units_numerator(idx) = [];

  units_denominator_length = fread(fid,1,'int');
  if units_denominator_length > MAXSTRINGLENGTH,
    error('Units denominator length = %d > %d\n',units_denominator_length,MAXSTRINGLENGTH);
  end
  units_denominator = fread(fid,units_denominator_length,'char=>char')';
  units_denominator = regexp(strtrim(units_denominator),'\s','split');
  idx = strcmp(units_denominator,'1');
  units_denominator(idx) = [];

  units{i} = struct('numerator',{units_numerator},'denominator',{units_denominator});
  datacurr = fread(fid,nframes,'double');
  data{i} = datacurr(:)';
end

fclose(fid);

trxfns = {'x','y','a','b','theta','x_mm','y_mm','a_mm','b_mm','theta_mm','timestamps'};

trx = struct;
for i = 1:numel(trxfns),
  fn = trxfns{i};
  j = find(strcmp(fn,fieldnames),1);
  if isempty(j), continue; end
  trx.(fn) = data{j};
end

if isfield(trx,'theta') && ~isfield(trx,'theta_mm') && isfield(trx,'x_mm'),
  trx.theta_mm = trx.theta;
end

j = find(strcmp('timestamp',fieldnames),1);
if ~isfield(trx,'timestamps') && ~isempty(j),
  trx.timestamps = data{j};
end
if ~isfield(trx,'timestamps'),
  trx.dt = diff(trx.timestamps);
end
if ~isfield(trx,'dt'),
  trx.dt = repmat(1/fps,[1,nframes-1]);
end
if ~isfield(trx,'timestamps'),
  trx.timestamps = cumsum([(firstframe-1)/fps,trx.dt]);
end

trx.firstframe = firstframe;
trx.off = 1-trx.firstframe;
trx.endframe = lastframe;
trx.nframes = nframes;
trx.id = fly_ids;
trx.moviename = moviename;
trx.sex = sex;
trx.fps = fps;

author: kristinbranson

repo: JAABA


text:
%L2depths Compute depths of PRMM from basis L.
%
%  No known depths are exploited. These can be different because of noise
%  etc.
%
%  Parameters:
%    opt.verbose(1) .. whether display info

function [Mdepths, lambda] = L2depths(L, M, Idepths, opt)

if nargin < 4, opt = []; end
if ~isfield(opt,'verbose')
  opt.verbose = 1; end

if opt.verbose, fprintf('Computing depths...'); tic; end

Mdepths = M;

[m n] = size(M); m = m/3;
lambda(m, n) = 0;  % memory allocation

for j = 1:n
  full     = find(~isnan(M(1:3:end,j)));
  mis_rows = intersect(find(Idepths(:,j)==0),full);
  if length(mis_rows) > 0
    submatrix = spread_depths_col(M(k2i(full),j),Idepths(full,j));

    % We want submatrix to be in the space L ->
    % we search for combination of columns of the base of L i.e.
    % L(k2i(full),:)*res(1:4)-submatrix*[1 res(5:length(res))] = 0
    right = submatrix(:,1);
    A     = [ L(k2i(full),:) -(submatrix(:,2:size(submatrix,2))) ];
    if rank(A) < size(A, 2)  % depths cannot be computed => kill the data
      kill = full(~Idepths(full,j));
      Mdepths(k2i(kill),j) = NaN; lambda(kill,j) = NaN;
    else
      res   = A \ right;

      %test: er should be near to zero
      %er=L(k2i(full),:)*res(1:4)-submatrix*[1 res(5:length(res))']'

      % depth corresponding to right is/are set to 1
      i = full(find(right(1:3:end))); lambda(i,j) = 1;
      Mdepths(k2i(i),j) = M(k2i(i),j);

      for ii = 1:size(submatrix,2)-1
        i = full(find(submatrix(1:3:end,1+ii))); lambda(i,j) = res(4+ii);
        Mdepths(k2i(i),j) = M(k2i(i),j)*lambda(i,j);
      end
    end
  end
end

if opt.verbose, disp(['(' num2str(toc) ' sec)']); end

author: strawlab

repo: MultiCamSelfCal


text:
%ESTIMATEPOSEBOARD  Pose estimation for a board of markers
%
%     [rvec, tvec, num] = cv.estimatePoseBoard(corners, ids, board, cameraMatrix, distCoeffs)
%     [rvec, tvec, num] = cv.estimatePoseBoard(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __corners__ cell array of already detected markers corners. For each
%   marker, its four corners are provided, (e.g `{{[x,y],..}, ..}`). The order
%   of the corners should be clockwise.
% * __ids__ list of identifiers for each marker in `corners` (0-based).
% * __board__ layout of markers in the board. The layout is composed by the
%   marker identifiers and the positions of each marker corner in the board
%   reference system. You can specify the board as a cell-array that starts
%   with the type name followed by option arguments `{Type, ...}`. There are
%   three types of boards available:
%   * __Board__ `{'Board', objPoints, dictionary, ids}`.
%     Creates a board of markers.
%   * __GridBoard__ `{'GridBoard', markersX, markersY, markerLength, markerSeparation, dictionary, 'FirstMarker',firstMarker}`.
%     Creates a a GridBoard object given the number of markers in each
%     direction and the marker size and marker separation.
%   * __CharucoBoard__ `{'GridBoard', squaresX, squaresY, squareLength, markerLength, dictionary}`.
%     Creates a CharucoBoard object given the number of squares in each
%     direction and the size of the markers and chessboard squares.
% * __cameraMatrix__ input 3x3 floating-point camera matrix
%   `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
% * __distCoeffs__ vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4]` of 4, 5, 8 or 12 elements.
%
% ## Output
% * __rvec__ Output vector `[x,y,z]` corresponding to the rotation vector of
%   the board.
% * __tvec__ Output vector `[x,y,z]` corresponding to the translation vector
%   of the board.
% * __num__ The number of markers from the input employed for the board pose
%   estimation. Note that returning a 0 means the pose has not been estimated.
%
% ## Options
% * __Rvec__, __Tvec__ Initial `rvec` and `tvec`. Used as initial guess if not
%   empty. The function uses the provided values as initial approximations of
%   the rotation and translation vectors, respectively, and further optimizes
%   them. Not set by default.
% * __UseExtrinsicGuess__ defines whether initial guess for `rvec` and `tvec`
%   will be used or not. default false.
%
% ## Inputs for Board
% * __objPoints__ array of object points of all the marker corners in the
%   board, i.e. their coordinates with respect to the board system. Each
%   marker include its 4 corners in CCW order
%   `{{[x1,y1,z1],[x2,y2,z2],[x3,y3,z3],[x4,y4,z4]}, ..}`. Usually object
%   points are planar with Z=0.
% * __dictionary__ the dictionary of markers employed for this board. This is
%   specified in the same format described in cv.detectMarkers.
% * __ids__ vector of the identifiers of the markers in the board (same size
%   as `objPoints`). The identifiers refers to the board dictionary (0-based).
%
% ## Inputs for GridBoard
% * __markersX__ number of markers in X direction.
% * __markersY__ number of markers in Y direction.
% * __markerLength__ marker side length (normally in meters).
% * __markerSeparation__ separation between two markers in the grid (same unit
%   as `markerLength`).
% * __dictionary__ dictionary of markers indicating the type of markers. This
%   is specified in the same format described in cv.detectMarkers.
% * __FirstMarker__ optional 0-based id of first marker in dictionary to use
%   on board. default 0
%
% ## Inputs for CharucoBoard
% * __squaresX__ number of chessboard squares in X direction.
% * __squaresY__ number of chessboard squares in Y direction.
% * __squareLength__ chessboard square side length (normally in meters).
% * __markerLength__ marker side length (same unit as `squareLength`)
% * __dictionary__ dictionary of markers indicating the type of markers. The
%   first markers in the dictionary are used to fill the white chessboard
%   squares. This is specified in the same format described in
%   cv.detectMarkers.
%
% The cv.estimatePoseBoard function receives the detected markers and returns
% the pose of a marker board composed by those markers. A board of marker has
% a single world coordinate system which is defined by the board layout. The
% returned transformation is the one that transforms points from the board
% coordinate system to the camera coordinate system. Input markers that are
% not included in the board layout are ignored.
%
% # Board of markers
%
% A board is a set of markers in the 3D space with a common cordinate system.
% The common form of a board of marker is a planar (2D) board, however any 3D
% layout can be used.
%
% # Grid Board
%
% A GridBoard is a special case of planar boards with grid arrangement of
% markers. It is the most common type of board. All markers are placed in the
% same plane in a grid arrangment.
%
% # ChArUco board
%
% Specific class for ChArUco boards. A ChArUco board is a planar board where
% the markers are placed inside the white squares of a chessboard. The
% benefits of ChArUco boards is that they provide both, ArUco markers
% versatility and chessboard corner precision, which is important for
% calibration and pose estimation.
%
% See also: cv.detectMarkers, cv.estimatePoseSingleMarkers, cv.Rodrigues,
%  cv.solvePnP
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv