text:
function varargout = implies(varargin)
%IMPLIES Logical implication
%
% IMPLIES(X,Y) creates a mixed integer representation of
% the constraint X --> Y, i.e. Y is true if X is true.
%
% Syntax
%   F = implies(X,Y,tol)
%
% Input
%   X : binary SDPVAR variable or a constraint
%   Y : binary SDPVAR variable or a constraint
%  tol: Optional threshhold for defining zero (see NOTE)
%
% Output
%   F : SET object
%
% Examples
%
%  binvar X Y; F = implies(X,Y);
%  binvar X;sdpvar Y; F = [implies(X,Y>=5), -10 <= Y <= 100];
%  binvar X;Y=sdpvar(3,1); F = [implies(X,[sum(Y);Y(2)]>=[5;0]), -1<= Y <= 10];
%
% Note
%
%  All variables in the expressions have to be explicitly bounded somewhere
%  in the model (implicit constraints are not sufficients such as [Y <= Z,...,Z<= 10]
%
%  Using implies with X non-binary is highly sensitive numerically.
%  The problem comes from the definition of 0 in a floating-point
%  environment, and precision in the solver. To account for this,
%  the user can supply a third argument to define a dead-zone around
%  zero, i.e Implies(X<=0,Y) will be replaced with IMPLIES(X<=-tol,Y)
%  Note, you typically need to tweak this number for your
%  application/solver. By default, YALMIP uses tol = 0, which means you
%  easily can get garbage... A positive number means YALMIP is cautious in
%  terms of activating the condition, while a negative number means YALMIP
%  will be aggressive  in activating the condition.
%
%   See also @SDPVAR/AND, @SDPVAR/OR, IFF

% There are some cases to take care of...
%
% X --> Y     binary/binary :                     Implemented
% X --> Y     binary/(LP,equality,sdp)            Implemented
% X --> Y     (LP,equality,sdp)/binary            Not implemented
% X --> Y     (LP,equality,sdp)/(LP,equality,sdp) Not implemented

X = varargin{1};
Y = varargin{2};

if isempty(X)
    varargout{1} = [];
end

switch class(X)

    case {'sdpvar','constraint','lmi'}
        varargout{1} = setupMeta(lmi([]), mfilename,varargin{:});

    case 'char'
        varargout{1} = implies_internal(varargin{3:end});

    case 'logical'
        if length(X)==1
            if X
                varargout{1} = Y;
            else
                varargout{1} = [];
            end
        else
            if length(X) == length(Y)
                i = find(X);
                if isempty(i)
                    varargout{1} = [];
                else
                    varargout{1} = Y(i);
                end
            else
                error('Size mismatch in input arguments');
            end
        end
end

author: shamilmamedov


text:
function [m, pot] = marginal_family(engine, i)
% MARGINAL_FAMILY Compute the marginal on i's family (global_inf_engine)
% [m, pot] = marginal_family(engine, i)
%

bnet = bnet_from_engine(engine);
[m, pot] = marginal_nodes(engine, family(bnet.dag, i));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function r8lib_test122 ( )

%*****************************************************************************80
%
%% R8LIB_TEST122 tests R8VEC_HISTOGRAM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  histo_num = 20;
  n = 1000;

  seed = 123456789;
  test_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8LIB_TEST122\n' );
  fprintf ( 1, '  R8VEC_HISTOGRAM histograms an integer vector.\n' );

  for test = 1 : test_num

    if ( test == 1 )

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Uniform data:\n' );

      a_lo =  0.0;
      a_hi = +1.0;
      [ a, seed ] = r8vec_uniform_01 ( n, seed );

    elseif ( test == 2 )

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Normal data:\n' );
      a_lo = -3.0;
      a_hi = +3.0;
      [ a, seed ] = r8vec_normal_01 ( n, seed );

    end

    histo_gram = r8vec_histogram ( n, a, a_lo, a_hi, histo_num );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Histogram of data:\n' );
    fprintf ( 1, '\n' );

    for i = 0 : histo_num+1

      if ( i == 0 )

        fprintf ( 1, '              %10f  %6d\n', a_lo, histo_gram(i+1) );

      elseif ( i <= histo_num )

        bin_lo = ( ( histo_num - i + 1 ) * a_lo   ...
                 + (             i - 1 ) * a_hi ) ...
                 / ( histo_num         );

        bin_hi = ( ( histo_num - i ) * a_lo   ...
                 + (             i ) * a_hi ) ...
                 / ( histo_num     );

        fprintf ( 1, '  %10f  %10f  %6d\n', bin_lo, bin_hi, histo_gram(i+1) );

      elseif ( i == histo_num+1 )

        fprintf ( 1, '  %10f              %6d\n', a_hi, histo_gram(i+1) );

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = lut_fun(x, lut_table, lut_signs, varargin)
    i0 = 1;
    if numel(varargin) > 0
        i0 = 1 + varargin{1};
    end
    y = lut_table(round(x * (10^lut_signs)) + i0);
    y = reshape(y, size(x));
end

repo: pTVreg


text:
function regModel = readBooleanRegModel(metModel, fileName)
% Reads Boolean regulatory network model
%
% USAGE:
%
%    regModel = readBooleanRegModel(metModel, fileName)
%
% INPUT:
%    metModel:    model
%
% OPTIONAL INPUT:
%    fileName:    file name
%
% OUTPUTS
%    regModel:    model containing the following fields:
%
%                   * .mets - Metabolite rules:
%
%                     * name - Metabolite/pool name (internal to the reg network model)
%                     * rule - Metabolite 'activation' rule
%                     * type - Metabolite type (extra/intracellular/pool)
%                     * excInd - Exchange flux indices corresponding to extracellular
%                       metabolites
%                     * icmRules - Intracellular metabolite 'activation' rules (based on
%                       a flux vector - `fluxVector`)
%                     * pool - Pool components
%                   * .regs - Regulator rules
%
%                     * name - Regulator name
%                     * rule - Regulator rule
%                     * comp - Regulator rule components (i.e. metabolites or other
%                       regulators that affect the state of this regulator)
%                     * ruleParsed - Rule in parsed format (based on metabolite state -
%                       `metState`, regulator state - `regState`)
%                   * .tars - Target rules
%
%                     * name - Target name
%                     * rule - Target rule
%                     * comp - Target rule components (i.e. metabolites or other
%                       regulators that affect the state of this regulator)
%                     * ruleParsed - Rule in parsed format (based on metabolite state -
%                       `metState`, regulator state - `regState`)
%
% .. Author: - Markus Herrgard 12/5/07

if (nargin < 2)
    fileName = 'iMHruletest.xls';
end

%% Read xls file
[tmp,textData] = xlsread(fileName);
nameList = textData(2:end,1);
ruleList = textData(2:end,2);
ruleType = textData(2:end,3);

%% Assign metabolite, regulator, and target rules
regModel = [];

metCnt = 0;
regCnt = 0;
tarCnt = 0;
nameList = strrep(nameList,'[','_');
nameList = strrep(nameList,']','');
nameList = deblank(nameList);
ruleList = strrep(ruleList,'[','_');
ruleList = strrep(ruleList,']','');
ruleList = deblank(ruleList);
for i = 1:length(ruleList)
    ruleList{i} = strrep(ruleList{i},'if ','');
    switch ruleType{i}
        case {'xcm','icm','xcp','icp'}
            metCnt = metCnt + 1;
            regModel.mets.name{metCnt} = nameList{i};
            regModel.mets.rule{metCnt} = ruleList{i};
            if (strcmp(ruleType{i},'xcp') | strcmp(ruleType{i},'icp'))
                regModel.mets.rule{metCnt} = splitString(ruleList{i},'+');
            end
            regModel.mets.type{metCnt} = ruleType{i};
        case {'reg','regc'}
            regCnt = regCnt + 1;
            regModel.regs.name{regCnt} = nameList{i};
            regModel.regs.rule{regCnt} = ruleList{i};
            [elements,newRule] = parseBoolean(ruleList{i});
            regModel.regs.comp{regCnt} = elements;
            regModel.regs.ruleParsed{regCnt} = newRule;
        case 'tar'
            tarCnt = tarCnt + 1;
            regModel.tars.name{tarCnt} = nameList{i};
            regModel.tars.rule{tarCnt} = ruleList{i};
            [elements,newRule] = parseBoolean(ruleList{i});
            regModel.tars.comp{tarCnt} = elements;
            regModel.tars.ruleParsed{tarCnt} = newRule;
    end
end

%% Parse metabolite rules

regModel.mets.name = columnVector(regModel.mets.name);
regModel.mets.rule = columnVector(regModel.mets.rule);
regModel.mets.type = columnVector(regModel.mets.type);

% Extracellular
selEcMets = find(strcmp(regModel.mets.type,'xcm'));
xcMets = strrep(regModel.mets.rule(selEcMets),'-','_');
regModel.mets.excInd = findRxnIDs(metModel,xcMets);

notInModel = xcMets(regModel.mets.excInd == 0);
for i = 1:length(notInModel)
    if (~strcmp(notInModel{i},'NA'))
        fprintf([regModel.mets.name{selEcMets(i)} ': exchange rxn for ' notInModel{i} ' not in metabolic model\n']);
    end
end

% Intracellular
selIcMets = find(strcmp(regModel.mets.type,'icm'));

for i = 1:length(selIcMets)
    [elements,newRule] = parseBoolean(regModel.mets.rule{selIcMets(i)});
    rxnInd = findRxnIDs(metModel,elements);
    for j = 1:length(rxnInd)
        if (rxnInd(j) > 0)
            newRule = strrep(newRule,['x(' num2str(j) ')'],['fluxVector(' num2str(rxnInd(j)) '_TMP_)']);
        else
            if (~strcmp(elements{j},'=0'))
                fprintf([regModel.mets.name{selIcMets(i)} ': rxn ' elements{j} ' not found in the metabolic model']);
            end
        end
    end
    regModel.mets.icmRules{i} = strrep(newRule,'_TMP_','');
end
regModel.mets.icmRules = columnVector(regModel.mets.icmRules);

% Pools

selPoolMets = find(strcmp(regModel.mets.type,'xcp') | strcmp(regModel.mets.type,'icp'));

for i = 1:length(selPoolMets)
    poolMets = regModel.mets.rule{selPoolMets(i)};
    [isInModel,metInd] = ismember(poolMets,regModel.mets.name);
    notInModel = poolMets(~isInModel);
    for j = 1:length(notInModel)
        fprintf([regModel.mets.name{selPoolMets(i)} ': metabolite ' notInModel{j} ' not in regulatory network model\n']);
    end
    regModel.mets.pool{i} = metInd(isInModel == 1);
end
regModel.mets.pool = columnVector(regModel.mets.pool);

%% Parse regulator rules

regModel.regs.name = columnVector(regModel.regs.name);
regModel.regs.rule = columnVector(regModel.regs.rule);
regModel.regs.comp = columnVector(regModel.regs.comp);
regModel.regs.ruleParsed = columnVector(regModel.regs.ruleParsed);

for i = 1:length(regModel.regs.name)
    elements = regModel.regs.comp{i};
    [isMet,metInd] = ismember(elements,regModel.mets.name);
    [isReg,regInd] = ismember(elements,regModel.regs.name);
    newRule = regModel.regs.ruleParsed{i};
    for j = 1:length(elements)
        if (metInd(j) > 0)
            newRule = strrep(newRule,['x(' num2str(j) ')'],['metState(' num2str(metInd(j)) '_TMP_)']);
        elseif (regInd(j) > 0)
            newRule = strrep(newRule,['x(' num2str(j) ')'],['regState(' num2str(regInd(j)) '_TMP_)']);
        else
            if (~strcmp(elements{j},'NA'))
                fprintf([regModel.regs.name{i} ': metabolite or regulator ' elements{j} ' not found in the model\n']);
            end
            newRule = '';
        end
    end
    regModel.regs.ruleParsed{i} = strrep(newRule,'_TMP_','');
end

%% Parse target rules

regModel.tars.name = columnVector(regModel.tars.name);
regModel.tars.rule = columnVector(regModel.tars.rule);
regModel.tars.comp = columnVector(regModel.tars.comp);
regModel.tars.ruleParsed = columnVector(regModel.tars.ruleParsed);

for i = 1:length(regModel.tars.name)
    elements = regModel.tars.comp{i};
    [isMet,metInd] = ismember(elements,regModel.mets.name);
    [isReg,regInd] = ismember(elements,regModel.regs.name);
    newRule = regModel.tars.ruleParsed{i};
    for j = 1:length(elements)
        if (metInd(j) > 0)
            newRule = strrep(newRule,['x(' num2str(j) ')'],['metState(' num2str(metInd(j)) '_TMP_)']);
        elseif (regInd(j) > 0)
            newRule = strrep(newRule,['x(' num2str(j) ')'],['regState(' num2str(regInd(j)) '_TMP_)']);
        else
            if (~strcmp(elements{j},'NA'))
                fprintf([regModel.tars.name{i} ': metabolite or regulator ' elements{j} ' not found in the model\n']);
            end
            newRule = '';
        end
    end
    regModel.tars.ruleParsed{i} = strrep(newRule,'_TMP_','');
end

author: opencobra

repo: cobratoolbox


text:
function M = spm_mesh_refine(M)
% Refine a triangle mesh
% FORMAT M = spm_mesh_refine(M)
% M        - a patch structure or gifti object
%__________________________________________________________________________
%
% See also:
%
% R.E. Bank, A.H. Sherman and A. Weiser. Refinement Algorithms and Data
% Structures for Regular Local Mesh Refinement. Scientific Computing
% (Applications of Mathematics and Computing to the Physical Sciences)
% (R. S. Stepleman, ed.), North-Holland (1983), 3-17.
% https://ccom.ucsd.edu/~reb/reports/a23.pdf.gz
%__________________________________________________________________________
% Copyright (C) 2018-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_refine.m 7514 2019-01-16 16:18:52Z guillaume $

V  = M.vertices;
F  = M.faces;

Nv = size(V,1);
Nf = size(F,1);

Vo = V;
Fo = zeros(4*Nf,3);
A  = spm_mesh_adjacency(M);
if isfield(M,'cdata')
    cdata = M.cdata;
else
    cdata = [];
end

for f=1:Nf
    T0 = F(f,:);
    T1 = T0([2 3 1]);
    v  = (V(T0,:) + V(T1,:)) / 2;
    if ~isempty(cdata), C = (cdata(T0,:) + cdata(T1,:)) / 2; end

    s = 1:3;
    b = [false false false];
    for j=1:3
        if A(T0(j),T1(j)) == 1
            s(j) = Nv + 1;
            A(T0(j),T1(j)) = s(j);
            A(T1(j),T0(j)) = s(j);
            Nv = s(j);
            b(j) = true;
        else
            s(j) = A(T0(j),T1(j));
        end
    end
    Vo(end+1:end+nnz(b),:) = v(b,:);
    if ~isempty(cdata), cdata(end+1:end+nnz(b),:) = C(b,:); end
    T0(4:6) = s;

    Fo(4*f+(-3:0),:) = T0([1 4 6;4 2 5;6 5 3;4 5 6]);
end

M.faces = Fo;
M.vertices = Vo;
if ~isempty(cdata), M.cdata = cdata; end

author: spm


text:
% Test file for chebtech/feval.m

function pass = test_feval(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebtech.techPref();
end

% Generate a few random points to use as test values.
seedRNG(7681);
x = 2 * rand(1000, 1) - 1;

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %%
    % Spot-check values for a couple of functions.  We can only expect
    % accuracy on % the order of the truncation level, so we use this as our
    % criterion.

    f = testclass.make(@(x) exp(x) - 1, [], pref);
    f_exact = @(x) exp(x) - 1;
    pass(n, 1) = (norm(feval(f, x) - f_exact(x), inf) < ...
        10*vscale(f)*eps);

    f = testclass.make(@(x) 1./(1 + x.^2), [], pref);
    f_exact = @(x) 1./(1 + x.^2);
    pass(n, 2) = (norm(feval(f, x) - f_exact(x), inf) < ...
        10*vscale(f)*eps);

    f = testclass.make(@(x) cos(1e4*x), [], pref);
    f_exact = @(x) cos(1e4*x);
    pass(n, 3) = (norm(feval(f, x) - f_exact(x), inf) < ...
        2e4*vscale(f)*eps);

    z = exp(2*pi*1i/6);
    f = testclass.make(@(t) sinh(t*z), [], pref);
    f_exact = @(t) sinh(t*z);
    pass(n, 4) = (norm(feval(f, x) - f_exact(x), inf) < ...
        10*vscale(f)*eps);

    %%
    % Check row vector and matrix input.

    err = feval(f, x.') - f_exact(x.');
    pass(n, 5) = (all(size(err) == [1 1000])) && (norm(err(:), inf) < ...
        10*vscale(f)*eps);

    x_mtx = reshape(x, [100 10]);
    err = feval(f, x_mtx) - f_exact(x_mtx);
    pass(n, 6) = (all(size(err) == [100 10])) && (norm(err(:), inf) < ...
        10*vscale(f)*eps);

    x_3mtx = reshape(x, [10 10 10]);
    err = feval(f, x_3mtx) - f_exact(x_3mtx);
    pass(n, 7) = (all(size(err) == [10 10 10])) && (norm(err(:), inf) < ...
        10*vscale(f)*eps);

    %%
    % Check operation for array-valued chebtech objects.

    f = testclass.make(@(x) [sin(x) x.^2 exp(1i*x)], [], pref);
    f_exact = @(x) [sin(x) x.^2 exp(1i*x)];
    err = feval(f, x) - f_exact(x);
    pass(n, 8) = all(max(abs(err)) < 10*max(vscale(f)*eps));

    %%
    % Test for evaluating array-valued chebtech objects at matrix arguments if
    % the operation makes sense.
    f = testclass.make(@(x) [sin(pi*x) cos(pi*x) exp(pi*x)], [], pref);
    x2 = [-1 0 1 ; .25 .5 .75];
    fx = feval(f, x2);
    f_exact = [0 0 0 -1 1 -1 exp(-pi) 1 exp(pi)
              [1 sqrt(2) 1 1 0 -1]/sqrt(2) exp(pi.*[.25 .5 .75])];
    pass(n, 9) = all(all(abs(fx - f_exact) < 10*max(vscale(f)*eps)));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

[c,p] = max(sigmoid(X*all_theta'),[],2);

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

end

author: Ayatans

repo: Machine-Learning-homework


text:
function [y,w] = idctn(y,DIM,w)

%IDCTN N-D inverse discrete cosine transform.
%   X = IDCTN(Y) inverts the N-D DCT transform, returning the original
%   array if Y was obtained using Y = DCTN(X).
%
%   IDCTN(X,DIM) applies the IDCTN operation across the dimension DIM.
%
%   Class Support
%   -------------
%   Input array can be numeric or logical. The returned array is of class
%   double.
%
%   Reference
%   ---------
%   Narasimha M. et al, On the computation of the discrete cosine
%   transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936.
%
%   Example
%   -------
%       RGB = imread('autumn.tif');
%       I = rgb2gray(RGB);
%       J = dctn(I);
%       imshow(log(abs(J)),[]), colormap(jet), colorbar
%
%   The commands below set values less than magnitude 10 in the DCT matrix
%   to zero, then reconstruct the image using the inverse DCT.
%
%       J(abs(J)<10) = 0;
%       K = idctn(J);
%       figure, imshow(I)
%       figure, imshow(K,[0 255])
%
%   See also DCTN, IDSTN, IDCT, IDCT2.
%
%   -- Damien Garcia -- 2009/04, revised 2009/11
%   website: <a
%   href="matlab:web('http://www.biomecardio.com')">www.BiomeCardio.com</a>

% ----------
%   [Y,W] = IDCTN(X,DIM,W) uses and returns the weights which are used by
%   the program. If IDCTN is required for several large arrays of same
%   size, the weights can be reused to make the algorithm faster. A typical
%   syntax is the following:
%      w = [];
%      for k = 1:10
%          [y{k},w] = idctn(x{k},[],w);
%      end
%   The weights (w) are calculated during the first call of IDCTN then
%   reused in the next calls.
% ----------

error(nargchk(1,3,nargin))

y = double(y);
sizy = size(y);

% Test DIM argument
if ~exist('DIM','var'), DIM = []; end
assert(~isempty(DIM) || ~isscalar(DIM),...
    'DIM must be a scalar or an empty array')
assert(isempty(DIM) || DIM==round(DIM) && DIM>0,...
    'Dimension argument must be a positive integer scalar within indexing range.')

% If DIM is empty, a DCT is performed across each dimension

if isempty(DIM), y = squeeze(y); end % Working across singleton dimensions is useless
dimy = ndims(y);

% Some modifications are required if Y is a vector
if isvector(y)
    dimy = 1;
    if size(y,1)==1
        if DIM==1, w = []; return
        elseif DIM==2, DIM=1;
        end
        y = y.';
    elseif DIM==2, w = []; return
    end
end

% Weighing vectors
if ~exist('w','var') || isempty(w)
    w = cell(1,dimy);
    for dim = 1:dimy
        if ~isempty(DIM) && dim~=DIM, continue, end
        n = (dimy==1)*numel(y) + (dimy>1)*sizy(dim);
        w{dim} = exp(1i*(0:n-1)'*pi/2/n);
    end
end

% --- IDCT algorithm ---
if ~isreal(y)
    y = complex(idctn(real(y),DIM,w),idctn(imag(y),DIM,w));
else
    for dim = 1:dimy
        if ~isempty(DIM) && dim~=DIM
            y = shiftdim(y,1);
            continue
        end
        siz = size(y);
        n = siz(1);
        y = reshape(y,n,[]);
        y = bsxfun(@times,y,w{dim});
        y(1,:) = y(1,:)/sqrt(2);
        y = ifft(y,[],1);
        y = real(y*sqrt(2*n));
        I = (1:n)*0.5+0.5;
        I(2:2:end) = n-I(1:2:end-1)+1;
        y = y(I,:);
        y = reshape(y,siz);
        y = shiftdim(y,1);
    end
end

y = reshape(y,sizy);

author: Sable

repo: mcbench-benchmarks


text:
function [xe,Pe,lm_seq] = update_std(xe,Pe,lm_seq,z,R)

lenz = length(find(z(3,:)>0));
lenx = size(xe,1);

global gDISTBEAR

nf = 0;

% % update
for i= 1:lenz
    % data association (based on landmark id)
    % TODO: use nearest neighbor to do the job
    is_exist = ~(lm_seq - z(3,i));
    idx = find(is_exist);

    ii = 2*i+(-1:0);

    % update: already in the state vecor
    if ~isempty(idx)
        nf = nf+1;
        jj = 2*nf+(-1:0);
        if gDISTBEAR
            [zhat,Hii] = measurement_model_std(xe,idx);   %Hekf=Hii
            H(jj,:) = Hii;
            r(jj,1) = [ z(1,i)-zhat(1,1);  pi_to_pi(z(2,i)-zhat(2,1)) ];
            Rf(jj,jj) = R(ii,ii);
        else
            [zhat,Hii] = measurement_model_std(xe,idx);   %Hekf=Hii
            H(jj,:) = Hii;
            r(jj,1) = [ z(1,i)-zhat(1,1);  (z(2,i)-zhat(2,1)) ];
            Rf(jj,jj) = R(ii,ii);
        end
    end
end

if nf~=0
    S = H*Pe*H'+ Rf;
    S = (S+S')*0.5;

    if isspd(S)

        K = Pe*H'/S;
        xe = xe + K*r;
        Pe = (eye(length(Pe)) - K*H) * Pe *(eye(length(Pe)) - K*H)' + K*Rf*K';

    end
end

% % augment
for i= 1:lenz
    % data association (known)
    is_exist = ~(lm_seq - z(3,i));
    idx = find(is_exist);

    lenx= size(xe,1);
    ii = 2*i + (-1:0);

    % add the new landmark into the state vector
    if isempty(idx)
        lm_seq = [lm_seq; z(3,i)];

        if gDISTBEAR
            % augment state
            d = z(1,i);
            th = z(2,i);
            k_xL = [d*cos(th); d*sin(th)];

            x_L = xe(1:2,1) + [ d*cos(th+xe(3)); d*sin(th+xe(3)) ];
            xe = [xe; x_L];

            % jacobians
            J = [0 -1; 1 0];
            C = [cos(xe(3))  -sin(xe(3));   sin(xe(3))   cos(xe(3)) ];

            H_Lk = [ 1/norm(k_xL)*k_xL'; 1/norm(k_xL)^2*k_xL'*J' ];

            HR = - H_Lk* C'*[eye(2)  J*(x_L-xe(1:2,1))];
            HL = H_Lk* C';

        else
            % augment state
            k_xL = z(1:2,i);

            C = [cos(xe(3)) -sin(xe(3));  sin(xe(3)) cos(xe(3)) ];
            x_L = xe(1:2,1) + C*k_xL;
            xe = [xe; x_L];

            % jacobians
            J = [0 -1; 1 0];
            C = [cos(xe(3))  -sin(xe(3));   sin(xe(3))   cos(xe(3)) ];

            HR = - C'*[eye(2)  J*(x_L-xe(1:2,1))];
            HL = C';
        end

        V_std = null([HR,HL]);

        % augment covariance
        rng= lenx+1:lenx+2;
        Pe(rng,rng)= inv(HL)*HR*Pe(1:3,1:3)*HR'*inv(HL)' + inv(HL)*R(ii,ii)*inv(HL)'; % landmark cov
        Pe(rng,1:3)= -inv(HL)*HR*Pe(1:3,1:3); % landmark-robot xcorr
        Pe(1:3,rng)= Pe(rng,1:3)';
        if lenx>3
            rnm= 4:lenx;
            Pe(rng,rnm)= -inv(HL)*HR*Pe(1:3,rnm);
            Pe(rnm,rng)= Pe(rng,rnm)';
        end
        Pe = 0.5*(Pe+Pe');
    end
end

author: rpng

repo: ocekf-slam

save_path: github-repos/MATLAB/rpng-ocekf-slam


text:
function sandia_sgmgg_test01 ( )

%*****************************************************************************80
%
%% SANDIA_SGMGG_TEST01 demonstrates the naive coefficient calculations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 August 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SANDIA_SGMGG_TEST01:\n' );
  fprintf ( 1, '  Demonstrate naive coefficient calculations.\n' );
%
%  Isotropic grid in 2D.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  1) Isotropic grid in 2D\n' );

  dim_num = 2;
  point_num = 7;

  sparse_index = [ ...
    0, 2; ...
    0, 3; ...
    1, 1; ...
    1, 2; ...
    2, 0; ...
    2, 1; ...
    3, 0 ]';

  sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index );
%
%  Isotropic grid in 3D.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  2) Isotropic grid in 3D\n' );

  dim_num = 3;
  point_num = 19;

  sparse_index = [ ...
    0, 1, 0; ...
    0, 2, 0; ...
    0, 3, 0; ...
    1, 0, 0; ...
    1, 1, 0; ...
    1, 2, 0; ...
    2, 0, 0; ...
    2, 1, 0; ...
    3, 0, 0; ...
    0, 0, 1; ...
    0, 1, 1; ...
    0, 2, 1; ...
    1, 0, 1; ...
    1, 1, 1; ...
    2, 0, 1; ...
    0, 0, 2; ...
    0, 1, 2; ...
    1, 0, 2; ...
    0, 0, 3 ]';

  sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index );
%
%  Anisotropic grid in 2D.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  3) Anisotropic grid in 2D\n' );

  dim_num = 2;
  point_num = 8;

  sparse_index = [ ...
    0, 2; ...
    1, 1; ...
    1, 2; ...
    2, 1; ...
    3, 0; ...
    3, 1; ...
    4, 0; ...
    5, 0 ]';

  sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index );
%
%  Generalized grid in 2D.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  4) Generalized grid in 2D\n' );

  dim_num = 2;
  point_num = 8;

  sparse_index = [ ...
    0, 0; ...
    0, 1; ...
    0, 2; ...
    0, 3; ...
    1, 0; ...
    1, 1; ...
    2, 0; ...
    3, 0 ]';

  sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index );

  return
end


text:
% VL_DEMO_KMEANS_ANN_SPEED   Compares Lloyd's, Elkan, and ANN k-means

numCenters = 200 ;
numTrials = 3 ;
maxNumIterations = 10 ;
initialization = 'randsel' ;
%initialization = 'plusplus';
distance = 'l2' ;

%% Create an example dataset

dimension = 32 ;
numData = 50000 ;

X = randn(dimension,numData);

%% Run various k-means algorithms on the data
algorithms = {'Lloyd','Elkan','ANN 1/4','ANN 1/10','ANN 1/50' } ;
options = {{'Algorithm', 'Lloyd'}, ...
           {'Algorithm', 'Elkan'}, ...
           {'Algorithm', 'ANN', 'MaxNumComparisons', ceil(numCenters / 4)}, ...
           {'Algorithm', 'ANN', 'MaxNumComparisons', ceil(numCenters / 10)}, ...
           {'Algorithm', 'ANN', 'MaxNumComparisons', ceil(numCenters / 50)}} ;
numCpus = [1 0] ;

clear time energy ;
for n = 1:2
  for a = 1:numel(algorithms)
    for t = 1:numTrials
      vl_threads(numCpus(n)) ;
      start = tic ;
      [C, A, E] = vl_kmeans(X, ...
                            numCenters, 'Verbose', ...
                            'Distance', distance, ...
                            'MaxNumIterations', maxNumIterations, ...
                            options{a}{:}) ;
      if vl_isoctave()
        time(t,a,n) = (tic() - start) / 1e6 ;
      else
        time(t,a,n) = toc(start) ;
      end
      energy(t,a,n) = E ;
    end
  end
end

% average over tirals
time = squeeze(mean(time,1)) ;
energy = squeeze(mean(energy,1)) ;

figure(1) ; clf ;
for n=1:2
  if n == 1
    str = 'Serial' ;
  else
    str = 'Parallel' ;
  end

  subplot(3,2,(n-1)+1) ;
  bar(time(:,n)) ;
  set(gca,'XTickLabel',algorithms);
  set(gca,'FontSize',8),
  xlabel('Algorithm');
  ylabel('Time [s]');
  title(str) ;

  subplot(3,2,(n-1)+3) ;
  bar(energy(:,n));
  set(gca,'XTickLabel',algorithms);
  set(gca,'FontSize',8),
  xlabel('Algorithm');
  ylabel('Energy');
  title(str) ;

  subplot(3,2,(n-1)+5) ;
  bar(time(1,1)./time(:,n)) ;
  set(gca,'XTickLabel',algorithms);
  set(gca,'FontSize',8),
  xlabel('Algorithm');
  ylabel('Speedup');
  title(str) ;
end

vl_demo_print('kmeans_speed',1);

author: jianxiongxiao

repo: ProfXkit


text:
function [startingvals,nu,lambda,LLs,outputParameters]=igarch_starting_values(startingvals,epsilon,fepsilon,p,q,T,errorType,igarchType,constant)
% Perform a grid search to find decent starting values for IGARCH(P,Q)
% esimtation.  If starting values are user supplied (and thus nonempty), reformats
% starting values depending on ERRORTYPE.
%
% USAGE:
%   [STARTINGVALS,NU,LAMBDA,LLS,OUTPUTPARAMETERS] = ...
%        igarch_starting_values(STARTINGVALS,EPSILON,EPSILON2,P,Q,T,ERRORTYPE,IGARCHTYPE,CONSTANT)
%
% INPUTS:
%   STARTINGVALS     - A vector of starting values or empty to perform a grid search
%   EPSILON          - A column of mean zero data
%   FEPSILON         - Either abs(EPSILON) or EPSILON.^2, depending on IGARCHTYPE
%   P                - Positive, scalar integer representing the number of
%                      symmetric innovations
%   Q                - Non-negative, scalar integer representing the number
%                      of lags of conditional variance (0 for ARCH)
%   T                - Length of EPSILON
%   ERRORTYPE        - [OPTIONAL] The error distribution used, valid types are:
%                       'NORMAL'    - Gaussian Innovations [DEFAULT]
%                       'STUDENTST' - T distributed errors
%                       'GED'       - Generalized Error Distribution
%                       'SKEWT'     - Skewed T distribution
%   IGARCHTYPE       - [OPTIONAL] The type of variance process, either
%                        1 - Model evolves in absolute values
%                        2 - Model evolves in squares [DEFAULT]
%   CONSTANT         - [OPTIONAL] Logical value indicating whether model
%                       should include a constant.  Default is true (include).
%
% OUTPUTS:
%   STARTINGVALS     - A vector of starting values (CONSTANT+p+q) by 1
%   NU               - Distribution kurtosis parameter, empty if not applicable
%   LAMBDA           - Distribution asymmetry parameter, empty if not applicable
%   LLS              - A vector of log likelihoods corresponding to OUTPUTPARAMETERS
%   OUTPUTPARAMETERS - A matrix of alternative starting values, sorted by log likelihood
%
% COMMENTS:
%   See also IGARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 7/12/2009

%Initialize variables
LLs=[];
outputParameters=[];

%No starting values provided
if isempty(startingvals)
    nu=[];
    lambda=[];

    %Procedure is to find best starting values, using a grid search find values for normal, then

    %Possible starting values based on commonly estimated values
    a=[.05 .1 .2];
    la=length(a);
    ab=1;
    lb=length(ab);

    %Many outputParameters and LLs
    outputParameters=zeros(la*lb,1+p+q-1);
    LLs=zeros(la,1);

    %Adjustment is needed to the intercept.  Assumes normality
    if igarchType==1
        adjFactor=sqrt(2/pi);
        backCast=mean(abs(epsilon));
    else
        adjFactor=1;
        backCast=cov(epsilon);
    end

    %Use an index to count
    index=1;

    for i=1:la
        %Loop over a
        alpha=a(i);
        tempAlpha=alpha*ones(p,1)/p;
        %Beta must satisfy the unit root
        beta=1-sum(tempAlpha);
        %Pick omega to match the unconditional
        if constant
            omega=backCast*.01*adjFactor;
        else
            omega = [];
        end
        %Build the parameter vector

        if q==0
            parameters=[omega; tempAlpha];
        else
            parameters=[omega; tempAlpha; beta*ones(q-1,1)/q];
        end
        %Set the output parameters
        outputParameters(index,:)=parameters';
        %Set the log likelihoods
        LLs(index)=igarch_likelihood(parameters, epsilon, fepsilon, p, q, 1, igarchType, constant, backCast, T, false);
        %Increment
        index=index+1;
    end
    %Sort the LLs so the best (lowest, since we minimize the -1*LL)
    [LLs,index]=sort(LLs);
    %Order the starting values
    startingvals=outputParameters(index(1),:)';
    %Order the ouputParameters
    outputParameters=outputParameters(index,:);
    %Use generic values for nu and lambda if needed
    if errorType==2
        nu=8;
    elseif errorType==3
        nu=1.9;
    elseif errorType==4
        nu=8;
        lambda=-.1;
    end
else
    %Values provided, only parse them
    nu=[];
    lambda=[];
    if errorType==2 || errorType==3 || errorType==4
        nu=startingvals(constant+p+q);
    end
    if errorType==4
        lambda=startingvals(constant+p+q+1);
    end
    startingvals=startingvals(1:constant+p+q-1);
end

author: bashtage

repo: mfe-toolbox


text:
% H = XLABELS(TEXT,...)
% H = XLABELS(HAX,TEXT,...)

% Last modified 2010-06-11
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function h = xlabels(text,varargin)

options = struct('Separation', []);

if isnumeric(text)
  hax = text;
  text = varargin{1};
  [options, errmsg, remopts] = argparse(options, varargin{2:end});
  error(errmsg);
else
  hax = gca;
  [options, errmsg, remopts] = argparse(options, varargin{:});
  error(errmsg);
end

h = zeros(size(hax));
for n=1:numel(hax)

  % Get axes position
  pos_ax = get(hax(n), 'Position');

  % Add a colorbar
  if iscell(text)
    h(n) = xlabel(hax(n), text{n}, remopts{:});
  else
    h(n) = xlabel(hax(n), text, remopts{:});
  end

  % Use the same units
  units = get(hax(n), 'Units');
  set(h(n), 'Units', units);

  pos = get(h(n), 'Position');
  pos(2) = -options.Separation;
  set(h(n), 'Position', pos);

end

if nargout < 1
  clear hcb;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% GeneralizedFibonacci.m by David Terr, Raytheon, 5-11-04

% Given integers n, a, and b, compute the nth generalized Fibonacci number
% U_n, defined by the recurrence U_n = a U_{n-1} + b U_{n-2} and the
% initial conditions U_0 = 0 and U_1 = 1.

function gfib = GeneralizedFibonacci(n,a,b)

% Make sure arguments are integers.
if ( n ~= floor(n) )
    error('First argument must be a nonnegative integer.');
    return;
end

if ~isreal(n) || n < 0
    error('First argument must be a nonnegative integer.');
    return;
end

if size(n,1) ~= 1 || size(n,2) ~= 1
    error('First argument must be a nonnegative integer.');
    return;
end

if ( a ~= floor(a) )
    error('Second argument must be an integer.');
    return;
end

if ~isreal(a)
    error('Second argument must be an integer.');
    return;
end

if size(a,1) ~= 1 || size(a,2) ~= 1
    error('Second argument must be an integer.');
    return;
end

if ( b ~= floor(b) )
    error('Third argument must be an integer.');
    return;
end

if ~isreal(b)
    error('Third argument must be an integer.');
    return;
end

if size(b,1) ~= 1 || size(b,2) ~= 1
    error('Third argument must be an integer.');
    return;
end

if ( n == 0 )
    gfib = 0;
    return;
end

D = a^2 + 4*b;  % discriminant

if ( D <= 0 )
    error('Discriminant must be positive.');
    return;
end

alpha = (a + sqrt(D))/2;
beta = (a - sqrt(D))/2;
gfib = round( ( alpha^n - beta^n ) / sqrt(D) );


text:
function skeletonLabeler( objTypes, vidNm, annNm )
% Video bound box (vbb) Labeler.
%
% USAGE
%  vbbLabeler( [objTypes], [vidNm], [annNm] )
%
% INPUTS
%  objTypes - [{'object'}] list of object types to annotate
%  imgDir   - [] initial video to load
%  resDir   - [] initial annotation to load
%

    % defaults
    if(nargin<1 || isempty(objTypes)), objTypes={'object'}; end
    if(nargin<2 || isempty(vidNm)), vidNm=''; end
    if(nargin<3 || isempty(annNm)), annNm=''; end

    % handles to gui objects / other globals
    [hFig, pSet, pCp, pMp] = deal([]);
    [curInd,curSeg, setApi,dispApi] = deal([]);

    makeLayout();
    setApi = setMakeApi(); % set the parameters
    dispApi  = dispMakeApi(); % display & update
    setApi.closeVid();
    set(hFig,'Visible','on');
    drawnow;

    % optionally load default data
    if(~isempty(vidNm)), setApi.openVid(vidNm); end
    if(~isempty(annNm)), setApi.openAnn(annNm); end

    function makeLayout()
        % hFig: main figure
        set(0,'Units','pixels'); ss=get(0,'ScreenSize');
        if( ss(3)<800 || ss(4)<600 ), error('screen too small'); end
        pos = [(ss(3)-580-400)/2 (ss(4)-650-200)/2 580+400 680+200];
        hFig = figure( 'Name','Labeler', 'NumberTitle','off', ...
          'Toolbar','auto', 'MenuBar','none', 'Resize','on', ...
          'Color','k', 'Visible','off', 'DeleteFcn',@exitLb, 'Position',pos, ...
          'keyPressFcn',@keyPress );

        % pSet: top panel containing setting
        pSet.h=uipanel('Units','pixels', 'BackgroundColor',[.1 .1 .1],'BorderType','none','Position', [8 645+200 580 30], 'Parent',hFig);
        pSet.hOpen   = uicontrol(pSet.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [5  3 90 20],'String','Open');
        pSet.nImgLbl = uicontrol(pSet.h,'Style','text','FontSize',8,'BackgroundColor',[.1 .1 .1],'ForegroundColor','w','Position',[185 0 50 20],'HorizontalAlignment','Left', 'String', '/');
        pSet.IndLbl  = uicontrol(pSet.h,'Style','edit','FontSize',8,'BackgroundColor',[.1 .1 .1],'ForegroundColor','w','Position',[135 3 50 20],'HorizontalAlignment','Right');
        pSet.hImageName = uicontrol(pSet.h,'Style','text','FontSize',8,'BackgroundColor',[.1 .1 .1],'ForegroundColor','w','Position',[240 0 90 20],'HorizontalAlignment','Left');

        % pCp: control pannel
        pCp.h=uipanel('Units','pixels', 'BackgroundColor',[.1 .1 .1],'BorderType','none','Position', [8 613+200 580 30],'Parent',hFig);
        pCp.hBtPrvImg    = uicontrol(pCp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [5  3 90 25],'String','PreviousImage');
        pCp.hBtShowImg   = uicontrol(pCp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [100  3 90 25],'String','ShowImage');
        pCp.hBtShowMask  = uicontrol(pCp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [195 3 80 25],'String','ShowMask');
        pCp.hBtShowSkel  = uicontrol(pCp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [280 3 90 25],'String','ShowSkeleton');
        pCp.hBtNextImg   = uicontrol(pCp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [375 3 80 25],'String','NextImage');
        pCp.hBtAnnSave   = uicontrol(pCp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [460 3 95 25],'String','AnnotationSave');

        % pMp main panel
        pMp.h=uipanel('Units','pixels', 'BackgroundColor',[.1 .1 .1],'BorderType','none','Position', [8 8 580+400 604+200],'Parent',hFig);
        pMp.hAx=axes('Units','pixels', 'Position', [5 40 560+400 560+200], 'Parent',pMp.h, 'Color', [1 1 1]);
        pMp.hBtPrvObj  = uicontrol(pMp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position', [5  3 90 25],'String','PreviousObject');
        pMp.hBtOriMask = uicontrol(pMp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position',[100  3 90 25],'String','OriginalMask');
        pMp.hBtExdDisc = uicontrol(pMp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position',[195 3 80 25],'String','ExtendDisc');
        pMp.hBtExdPoly = uicontrol(pMp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position',[280 3 80 25],'String','ExtendPoly');
        pMp.hBtAnn     = uicontrol(pMp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position',[365 3 80 25],'String','Annotation');
        pMp.hBtNxtObj  = uicontrol(pMp.h,'BackgroundColor',[.7 .7 .7],'FontSize', 10, 'Position',[450 3 80 25],'String','NextObject');

        function exitLb( h, e ) %#ok<INUSD>
        	setApi.closeVid();
        end

        function keyPress( h, evnt ) %#ok<INUSL>
            char=int8(evnt.Character); if(isempty(char)), char=0; end;
            if( char==28 ), dispApi.setPlay(-inf); end
            if( char==29 ), dispApi.setPlay(+inf); end
            if( char==31 ), dispApi.setPlay(0); end
        end
    end

    function api = setMakeApi()
        % variables
        [imgPath, segPath, imgLst] = deal([]);
        [imgROI, seg, segIdx, skeGt,skeInst, extIdx] = deal([]);
        [segMask, extMask, skelFlag, annFlag] = deal([]);
        [h, w] = deal([]);

        imgPath     = [VOC2011Path '\VOCdevkit\VOC2011\JPEGImages';
        segPath     = [VOC2011Path '\VOCdevkit\VOC2011\SegmentationObject';
        imgLstPath  = [VOC2011Path '\VOCdevkit\VOC2011\ImageSets\Segmentation\train.txt';
        f=fopen(imgLstPath); imgLst=textscan(f,'%s %*s'); imgLst=imgLst{1}; fclose(f);
        set(pSet.nImgLbl, 'String', ['\' num2str(length(imgLst))]);

        % callbacks
        set(pSet.hOpen, 'callback', @(h,e) imgOpen());
        set(pCp.hBtNextImg, 'callback', @(h,e) nextImg());
        set(pCp.hBtPrvImg, 'callback', @(h,e) previousImg());
        set(pCp.hBtShowImg, 'callback', @(h,e) showImg());
        set(pCp.hBtShowMask, 'callback', @(h,e) showMask());
        set(pCp.hBtShowSkel, 'callback', @(h,e) showSkel());
        set(pCp.hBtAnnSave, 'callback', @(h,e) annSave());

        set(pMp.hBtPrvObj, 'callback', @(h,e) previousObj());
        set(pMp.hBtNxtObj, 'callback', @(h,e) nextObj());
        set(pMp.hBtOriMask, 'callback', @(h,e) oriMask());
        set(pMp.hBtExdDisc, 'callback', @(h,e) exdMaskDisc());
        set(pMp.hBtExdPoly, 'callback', @(h,e) exdMaskPoly());
        set(pMp.hBtAnn, 'callback', @(h,e) annotation());

        % create api
        api=struct('imgSet',@imgSet, ...
          'imgOpen',@imgOpen, 'closeVid', @closeVid);

        function imgShow()
            curSeg = 0;
            set(pSet.IndLbl, 'String', num2str(curInd));
            set(pSet.hImageName, 'String', [imgLst{curInd} '.jpg']);
            set(pCp.hBtShowSkel, 'Enable','off');
            set(pCp.hBtAnnSave, 'Enable','off');

            cla(pMp.hAx); set(pMp.hAx,'XTick',[],'YTick',[]);
            set( hFig, 'CurrentAxes', pMp.hAx );
            imgROI = imread([imgPath '\' imgLst{curInd} '.jpg']);
            [h,w,~] =size(imgROI);
            set( pMp.hAx, 'Position', [5 40 980 804]);
            img = uint8(50*ones(804,980,3));
            img(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,1) = imgROI(:,:,1);
            img(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,2) = imgROI(:,:,2);
            img(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,3) = imgROI(:,:,3);
            image(img); axis off;
            seg = imread([segPath '\' imgLst{curInd} '.png']);
            segIdx = unique(seg);
            skeGt = false(h,w);
            extMask = cell(50,1);
            skeInst = cell(length(segIdx)-2,1);
            skelFlag= false(length(segIdx)-2,1);
            segMask = false(804,980);
        end

        function segShow(segmask)
            tmp = uint8(50*ones(804,980,3));
            tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,1) = imgROI(:,:,1);
            tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,2) = imgROI(:,:,2);
            tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,3) = imgROI(:,:,3);
            tmp(:,:,1) = tmp(:,:,1).*uint8(~segmask);
            tmp(:,:,2) = tmp(:,:,2).*uint8(~segmask);
            tmp(:,:,3) = tmp(:,:,3).*uint8(~segmask);
            image(tmp);
            axis off;
        end

        function imgOpen()
            num = get(pSet.IndLbl, 'String');
            if(isempty(num))
                curInd = 1;
            else
                curInd = uint16(str2num(num));
            end
            imgShow();
        end

        function nextImg()
            if(curInd < 1112)
                curInd = curInd+1;
            end
            imgShow();
        end

        function previousImg()
            if(curInd>1)
                curInd = curInd-1;
            end
            imgShow();
        end

        function showImg()
            imgROI = imread([imgPath '\' imgLst{curInd} '.jpg']);
            [h,w,~] =size(imgROI);
            set( pMp.hAx, 'Position', [5 40 980 804]);
            img = uint8(50*ones(804,980,3));
            img(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,1) = imgROI(:,:,1);
            img(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,2) = imgROI(:,:,2);
            img(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,3) = imgROI(:,:,3);
            image(img); axis off;
        end

        function showMask()
            segMaskROI = seg == segIdx(length(segIdx));
            segMask(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w) = segMaskROI;
            segShow(segMask);
        end

        function showSkel()
            for i= 1:length(segIdx)-2
                skeGt = skeGt|skeInst{i};
            end

            tmp = uint8(50*ones(804,980,3));
            tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,1) = imgROI(:,:,1);
            tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,2) = imgROI(:,:,2);
            tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,3) = imgROI(:,:,3);
            tmpMask = false(804,980);
            tmpMask(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w)= seg>0;
            tmpGt = false(804,980);
            tmpGt(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w) = skeGt;
            tmpMask = uint8(~tmpMask);
            tmpMask(tmpGt) = 255;

            tmp(:,:,1) = tmp(:,:,1).*tmpMask;
            tmp(:,:,2) = tmp(:,:,2).*tmpMask;
            tmp(:,:,3) = tmp(:,:,3).*tmpMask;
            image(tmp);axis off;

            set(pCp.hBtAnnSave, 'Enable','on');
        end

        function annSave()
            sym = skeGt;
            save(['skeGt/' imgLst{curInd} '.mat'], 'sym');
            msgbox('save sucessed.');
        end

        function previousObj()
            if(curSeg>1)
                curSeg = curSeg-1;
            else
               curSeg = 1;
            end
            annFlag = true;
            segMaskROI = seg == segIdx(curSeg+1);
            segMask(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w) = segMaskROI;
            segShow(segMask);
            extIdx = 1;
        end

        function nextObj()
            if(curSeg<length(segIdx)-2)
                curSeg = curSeg+1;
            end
            annFlag = true;
            segMaskROI = seg == segIdx(curSeg+1);
            segMask(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w) = segMaskROI;
            segShow(segMask);
            extIdx = 1;
        end

        function oriMask()
            extIdx = 1;
            annFlag = true;
            segShow(segMask);
        end

        function exdMaskDisc()
            if curSeg >0
                [x,y] = ginput(2);
                if(length(x)<2)
                    return;
                end

                tmpMask = false(size(segMask));
                % extend by disc
                dot1 = [x(1) y(1)];
                dot2 = [x(2) y(2)];
                tmpMask = extendDisc(tmpMask, dot1, dot2);

                tmpMask = imfill(tmpMask, 'holes');

                if extIdx==1
                    extMask{extIdx} = tmpMask;
                else
                    extMask{extIdx} = extMask{extIdx-1}|tmpMask;
                end
                tmpMask = segMask|extMask{extIdx};
                segShow(tmpMask);
                if (extIdx<51)
                    extIdx = extIdx+1;
                end
            end
        end

        function exdMaskPoly()
            if curSeg >0
                [x,y] = ginput(4);
                if(length(x)<2)
                    return;
                end
                tmpMask = false(size(segMask));
                % extend by polygon
                x = [x(1) x(2) x(3) x(4)];
                y = [y(1) y(2) y(3) y(4)];
                tmpMask = extendPolygon(tmpMask, x, y);
                tmpMask = imfill(tmpMask, 'holes');
                if extIdx==1
                    extMask{extIdx} = tmpMask;
                else
                    extMask{extIdx} = extMask{extIdx-1}|tmpMask;
                end
                tmpMask = segMask|extMask{extIdx};
                segShow(tmpMask);
                if (extIdx<51)
                    extIdx = extIdx+1;
                end
            end
        end

        function annotation()
            if curSeg >0
                if(extIdx>1)
                    tmpMask = extMask{extIdx-1};
                    tmpMask = segMask|tmpMask;
                else
                    tmpMask = segMask;
                end

                try

                    epsilon = 0.8;
                    skel = skel_pruning_bpr(tmpMask, epsilon, 9);
                    skel = skel>0;

                    tmp = uint8(50*ones(804,980,3));
                    tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,1) = imgROI(:,:,1);
                    tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,2) = imgROI(:,:,2);
                    tmp(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w,3) = imgROI(:,:,3);

                    if annFlag
                        tmpMask = uint8(~tmpMask);
                    else
                        skel = skel & segMask;
                        tmpMask = uint8(~segMask);
                    end

                    tmpMask(skel) = 255;

                    tmp(:,:,1) = tmp(:,:,1).*tmpMask;
                    tmp(:,:,2) = tmp(:,:,2).*tmpMask;
                    tmp(:,:,3) = tmp(:,:,3).*tmpMask;
                    image(tmp);axis off;
                    annFlag = ~annFlag;

                    skelFlag(curSeg) = true;
                    skeInst{curSeg} = skel(420-floor(h/2):419-floor(h/2)+h,485-floor(w/2):484-floor(w/2)+w);

                    if all(skelFlag==true)
                        set(pCp.hBtShowSkel, 'Enable','on');
                    end
                end
            end
        end

        function mask = extendDisc(extMask, dot1, dot2)
            dot1 = double(dot1);
            dot2 = double(dot2);
            mask = extMask;
            x0 = (dot1(1)+dot2(1))/2;
            y0 = (dot1(2)+dot2(2))/2;
            R0 = sqrt((dot1(1)-dot2(1))^2+(dot1(2)-dot2(2))^2)/2;
            dy = dot2(2)-dot1(2);
            dx = dot2(1)-dot1(1);
            theta0 = atan((dot2(2)-dot1(2))/(dot2(1)-dot1(1)));
            if dot2(2)-dot1(2)>=0 && dot2(1)-dot1(1)<0
                theta0 = atan((dot2(2)-dot1(2))/(dot2(1)-dot1(1)))-pi;
            elseif dot2(2)-dot1(2)>=0 && dot2(1)-dot1(1)>0
                theta0 = atan((dot2(2)-dot1(2))/(dot2(1)-dot1(1)));
            elseif dot2(2)-dot1(2)<0 && dot2(1)-dot1(1)<0
                theta0 = atan((dot2(2)-dot1(2))/(dot2(1)-dot1(1)))-pi;
            elseif dot2(2)-dot1(2)<0 && dot2(1)-dot1(1)>0
                theta0 = atan((dot2(2)-dot1(2))/(dot2(1)-dot1(1)));
            elseif dot2(2)-dot1(2)<0 && dot2(1)-dot1(1)==0
                theta0 = -pi/2;
            elseif dot2(2)-dot1(2)>0 && dot2(1)-dot1(1)==0
                theta0 = pi/2;
            end

            Theta=linspace(theta0,theta0+pi,1000);
            x = x0+ cos(Theta)*R0;
            y = y0 + sin(Theta)*R0;
            for ii=1:1000
                mask(round(y(ii)),round(x(ii)))=true;
            end
            x = linspace(dot1(1),dot2(1),500);
            y = linspace(dot1(2),dot2(2),500);
            for ii=1:500
                mask(round(y(ii)),round(x(ii)))=true;
            end

        end

        function mask = extendPolygon(extMask, x, y)
            mask = extMask;
            a = [linspace(x(1),x(2),500) linspace(x(2),x(3),500) linspace(x(3),x(4),500) linspace(x(4),x(1),500)];
            b = [linspace(y(1),y(2),500) linspace(y(2),y(3),500) linspace(y(3),y(4),500) linspace(y(4),y(1),500)];
            for ii=1:2000
                mask(round(b(ii)),round(a(ii)))=true;
            end
        end

        function closeVid()
%             dispApi.init([]);
        end
    end

    function api = dispMakeApi()
        api=struct( 'showImage',@showImage, 'imgSet',@imgSet, ...
          'imgOpen',@imgOpen);

    end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% removeNotIndexed
% Below is a demonstration of the features of the |removeNotIndexed| function

%%
clear; close all; clc;

%% Syntax
% |[F,V]=removeNotIndexed(F,V);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function a = i4mat_flip_rows ( m, n, a )

%*****************************************************************************80
%
%% I4MAT_FLIP_ROWS swaps the rows of an I4MAT.
%
%  Discussion:
%
%    An I4MAT is an integer matrix.
%
%    To "flip" the rows of an IMAT is to start with something like
%
%      11 12 13 14 15
%      21 22 23 24 25
%      31 32 33 34 35
%      41 42 43 44 45
%      51 52 53 54 55
%
%    and return
%
%      51 52 53 54 55
%      41 42 43 44 45
%      31 32 33 34 35
%      21 22 23 24 25
%      11 12 13 14 15
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A[M,N], the matrix whose rows are to be flipped.
%
%    Output, integer A[M,N], the flipped matrix.
%
  for i = 1 : floor ( m / 2 )
    b(1,    1:n) = a(    i,1:n);
    a(    i,1:n) = a(m+1-i,1:n);
    a(m+1-i,1:n) = b(1,    1:n);
  end

  return
end


function blur_img = synthesize_nonuniform_blur(img, info_data, gyro_data)

    img= im2double(img);

    bd = 20; %% tmp
    img_pad = padarray(img, [bd, bd], 'replicate', 'both');

    num_row = size(img_pad, 1);
    num_col = size(img_pad, 2);

    time_interval = 0.01;

    % set up start and end time
    st_time = info_data(1, 3) + (-8) * time_interval - info_data(1, 2);
    ed_time = info_data(1, 3) + (-8) * time_interval;

    valid_row = 1:size(gyro_data, 1);%find(gyro_data(:, end) < ed_time + time_interval);

    % interpolate
    portion = 3;%ceil(1/infoData(1,2));
    intervalAfterTime = time_interval/portion;

    % linear interpolation
    tempInd = 1:portion:1+(length(valid_row)-1)*portion;
    x = tempInd;
    v = gyro_data(valid_row,1:end);
    xq = 1:1+(length(valid_row)-1)*portion;
    vq = interp1(x,v,xq);
    gyroInterData = vq;

    % compute angular velocity
    avelInter = gyroInterData(2:end,4:6).*repmat((gyroInterData(2:end,end)- ...
        gyroInterData(1:end-1,end)),[1,3]);

    avelInter = cumsum(avelInter)/1;

    % switch the tx ty for headleft
    tmp = gyroInterData(:,1);
    gyroInterData(:,1) = -gyroInterData(:,2);
    gyroInterData(:,2) = tmp;

    % compute velocity
    comDrift = [0,0.00,0];
    velInter = (gyroInterData(2:end,1:3)-repmat(comDrift,[size(gyroInterData,1)-1,1]))...
        .*repmat((gyroInterData(2:end,end)- gyroInterData(1:end-1,end)),[1,3]);

    velInter = cumsum(velInter*9.8);

    % initial speed
    initialSpeed = [0,0,0];
    velInter = velInter + repmat(initialSpeed,[size(velInter,1),1]);

    % compute the position
    activeRow = intersect(find(gyroInterData(:,end)<ed_time+time_interval),...
        find(gyroInterData(:,end)>=st_time));
    gyroInterData = gyroInterData(activeRow,:);
    activeRow(end)=[];

    rotInter = avelInter(activeRow,:).*repmat((gyroInterData(2:end,end)...
        - gyroInterData(1:end-1,end)),[1,3]);
    rotInter = cumsum(rotInter);

    tranInter = velInter(activeRow,:).*repmat((gyroInterData(2:end,end)...
        - gyroInterData(1:end-1,end)),[1,3]);
    tranInter = cumsum(tranInter)/1;

    pixelPitch = 1e-5;
    poseInter = [tranInter,avelInter(activeRow,:)];
    poseInter = poseInter - repmat(poseInter(1,:),[size(poseInter,1),1]);

    if(1)
        factor = [1,1,1,-1,-1,-1];
    else
        factor = [1,1,1,1,1,1];
    end
    poseInter = poseInter.*repmat(factor,[size(poseInter,1),1]);

    ox = 0;
    oy = 0;
    rotCen = [0,0,0,ox,oy,0];
    poseInterOff = poseInter + repmat(rotCen,[size(poseInter,1),1]);

    fprintf('Generate local PSF...\n');

    data = poseInterOff;
    N_pose = size(data,1);
    weight = 1/N_pose;

    m = num_row;
    n = num_col;
    K = sparse([],[],[],m*n,m*n);

    for j = 1:N_pose
        pose = [data(j,4),data(j,5),data(j,6),data(j,1),data(j,2),data(j,3)];

        K_mat = compute_Kpose_6D(m,n,pose,1);

        K = K + weight * K_mat;
    end

    blur_img = img_pad;
    for c = 1:size(img_pad, 3)

        img_vector = img_pad(:, :, c);
        img_vector = K * img_vector(:);
        blur_img(:, :, c) = reshape(img_vector, num_row, num_col);

    end

    blur_img = blur_img(1+bd:end-bd, 1+bd:end-bd, :);

end


text:
%% t_glmComputeContrastMap
%
% Illustrates how to compute a contrast map for later projection onto a 3D
% mesh.
%
% Before proceeding with tutorial:
%	* Run t_glmRun.
%
% See also T_GLMRUN.
%
% Tested 01/04/2011 - MATLAB r2008a, Fedora 12, Current Repos
%
% Stanford VISTA
%

%% Initialize the key variables and data path:
% Data directory (where the mrSession file is located)
dataDir = fullfile(mrvDataRootPath,'functional','vwfaLoc');
dataType = 'MotionComp';

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

%% Retrieve data structure and set data type:
vw_ip = initHiddenInplane();
vw_ip = viewSet(vw_ip, 'currentDataType', dataType);

%% Get information re: the experiment/trials:
stimuli = er_concatParfiles(vw_ip);
nConds = length(stimuli.condNums);

%% Print condition numbers and names for input into contrast fxn:
fprintf('[##] - Condition Name\n');
fprintf('---------------------\n');
for i = 1:nConds
    fprintf('[%02d] - %s\n', stimuli.condNums(i), stimuli.condNames{i});
end

%% Choose active and control conditions:
activeConds     = [0]; % Fixation
% versus
controlConds    = [1 2]; % Word & WordScramble

% Choose a name for the contrast - left empty to assign default
contrastName    = [];

vw_ip = viewSet(vw_ip, 'currentDataType', 'GLMs');
%% Compute the contrast map for view on mesh:
computeContrastMap2(vw_ip, activeConds, controlConds, contrastName);

%% Return to the original directory.
cd(curDir);

%% END

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = T2quaternion2(T)
%   Returns the quaternion corresponding to an homogeneous transformation
%   matrix T. Only the 3x3 rotation matrix in T is used.
%
%   See also QPROD, QUATERNION2T.
% The method implemented here was extracted from:
% Accurate Computation of Quaternions from Rotation Matrices.
% Soheil Sarabandi and Federico Thomas
%  http://www.iri.upc.edu/files/scidoc/2068-Accurate-Computation-of-Quaternions-from-Rotation-Matrices.pdf
%   Author: Arturo Gil. Universidad Miguel Hern�ndez de Elche. email:
%   arturo.gil@umh.es date:   21/04/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 = T2quaternion2(T)
% use only the orientation from T
R = T(1:3, 1:3);

if R(1,1) + R(2,2) + R(3,3) > 0
    Q(1) = 0.5*sqrt(1+R(1,1) + R(2,2) + R(3,3));
else
    num = (R(3,2)-R(2,3))^2 + (R(1,3)-R(3,1))^2 + (R(2,1)-R(1,2))^2;
    den = 3 - R(1,1) - R(2,2) - R(3,3);
    Q(1) = 0.5*sqrt(num/den);
end

if R(1,1) - R(2,2) - R(3,3) > 0
    Q(2) = 0.5*sqrt(1 + R(1,1) - R(2,2) - R(3,3));
else
    num = (R(3,2)-R(2,3))^2 + (R(1,3)+R(3,1))^2 + (R(2,1)+R(1,2))^2;
    den = 3 - R(1,1) + R(2,2) + R(3,3);
    Q(2) = 0.5*sqrt(num/den);
end

if -R(1,1) + R(2,2) - R(3,3) > 0
    Q(3) = 0.5*sqrt(1 - R(1,1) + R(2,2) - R(3,3));
else
    num = (R(3,2)+R(2,3))^2 + (R(1,3)-R(3,1))^2 + (R(2,1)+R(1,2))^2;
    den = 3 + R(1,1) - R(2,2) + R(3,3);
    Q(3) = 0.5*sqrt(num/den);
end

if -R(1,1) - R(2,2) + R(3,3) > 0
    Q(4) = 0.5*sqrt(1 - R(1,1) - R(2,2) + R(3,3));
else
    num = (R(3,2)+R(2,3))^2 + (R(1,3)+R(3,1))^2 + (R(2,1)-R(1,2))^2;
    den = 3 + R(1,1) + R(2,2) + R(3,3);
    Q(4) = 0.5*sqrt(num/den);
end

repo: ARTE


text:
function value = r8_gamma_01_sample ( a )

%*****************************************************************************80
%
%% R8_GAMMA_01_SAMPLE samples the standard Gamma distribution.
%
%  Discussion:
%
%    This procedure corresponds to algorithm GD in the reference.
%
%    pdf ( a; x ) = 1/gamma(a) * x^(a-1) * exp ( - x )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    Original FORTRAN77 version by Barry Brown, James Lovato.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joachim Ahrens, Ulrich Dieter,
%    Generating Gamma Variates by a Modified Rejection Technique,
%    Communications of the ACM,
%    Volume 25, Number 1, January 1982, pages 47-54.
%
%  Parameters:
%
%    Input, real A, the parameter of the standard gamma
%    distribution.  0.0 < A < 1.0.
%
%    Output, real VALUE, a random deviate from the distribution.
%
  a1 =  0.3333333;
  a2 = -0.2500030;
  a3 =  0.2000062;
  a4 = -0.1662921;
  a5 =  0.1423657;
  a6 = -0.1367177;
  a7 =  0.1233795;

  e1 = 1.0;
  e2 = 0.4999897;
  e3 = 0.1668290;
  e4 = 0.0407753;
  e5 = 0.0102930;

  q1 =  0.04166669;
  q2 =  0.02083148;
  q3 =  0.00801191;
  q4 =  0.00144121;
  q5 = -0.00007388;
  q6 =  0.00024511;
  q7 =  0.00024240;

  sqrt32 = 5.656854;

  if ( 1.0 <= a )

    s2 = a - 0.5;
    s = sqrt ( s2 );
    d = sqrt32 - 12.0 * s;
%
%  Immediate acceptance.
%
    t = r8_normal_01_sample ( );
    x = s + 0.5 * t;
    value = x * x;

    if ( 0.0 <= t )
      return
    end
%
%  Squeeze acceptance.
%
    u = r8_uniform_01_sample ( );
    if ( d * u <= t * t * t )
      return
    end

    r = 1.0 / a;
    q0 = (((((( q7 ...
      * r + q6 ) ...
      * r + q5 ) ...
      * r + q4 ) ...
      * r + q3 ) ...
      * r + q2 ) ...
      * r + q1 ) ...
      * r;
%
%  Approximation depending on size of parameter A.
%
    if ( 13.022 < a )
      b = 1.77;
      si = 0.75;
      c = 0.1515 / s;
    elseif ( 3.686 < a )
      b = 1.654 + 0.0076 * s2;
      si = 1.68 / s + 0.275;
      c = 0.062 / s + 0.024;
    else
      b = 0.463 + s + 0.178 * s2;
      si = 1.235;
      c = 0.195 / s - 0.079 + 0.16 * s;
    end
%
%  Quotient test.
%
    if ( 0.0 < x )

      v = 0.5 * t / s;

      if ( 0.25 < abs ( v ) )
        q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v );
      else
        q = q0 + 0.5 * t * t * (((((( a7 ...
          * v + a6 ) ...
          * v + a5 ) ...
          * v + a4 ) ...
          * v + a3 ) ...
          * v + a2 ) ...
          * v + a1 ) ...
          * v;
      end

      if ( log ( 1.0 - u ) <= q )
        return
      end

    end

    while ( 1 )

      e = r8_exponential_01_sample ( );
      u = 2.0 * r8_uniform_01_sample ( ) - 1.0;

      if ( 0.0 <= u )
        t = b + abs ( si * e );
      else
        t = b - abs ( si * e );
      end
%
%  Possible rejection.
%
      if ( t < -0.7187449 )
        continue
      end
%
%  Calculate V and quotient Q.
%
      v = 0.5 * t / s;

      if ( 0.25 < abs ( v ) )
        q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v );
      else
        q = q0 + 0.5 * t * t * (((((( a7 ...
          * v + a6 ) ...
          * v + a5 ) ...
          * v + a4 ) ...
          * v + a3 ) ...
          * v + a2 ) ...
          * v + a1 ) ...
          *  v;
      end
%
%  Hat acceptance.
%
      if ( q <= 0.0 )
        continue
      end

      if ( 0.5 < q )
        w = exp ( q ) - 1.0;
      else
        w = (((( e5 * q + e4 ) * q + e3 ) * q + e2 ) * q + e1 ) * q;
      end
%
%  May have to sample again.
%
      if ( c * abs ( u ) <= w * exp ( e - 0.5 * t * t ) )
        break
      end

    end

    x = s + 0.5 * t;
    value = x * x;

    return
%
%  Method for A < 1.
%
  else

    b = 1.0 + 0.3678794 * a;

    while ( 1 )

      p = b * r8_uniform_01_sample ( );

      if ( p < 1.0 )

        value = exp ( log ( p ) / a );

        if ( value <= r8_exponential_01_sample ( ) )
          return
        end

        continue

      end

      value = - log ( ( b - p ) / a );

      if ( ( 1.0 - a ) * log ( value ) <= r8_exponential_01_sample ( ) )
        break
      end

    end

  end

  return
end


text:
function tparameters=scalar_vt_vech_transform(parameters,p,o,q,kappa)
% SCALAR_VT_VECH(P,Q) parameter transformation.  Used to map parameters
% from a scalar MVGARCH process to the real line. Used in the estimation of SCALAR_VT_VECH.
%
% USAGE:
%   [TPARAMETERS]=scalar_vt_vech_transform(PARAMETERS,P,O,Q,KAPPA)
%
% INPUTS:
%   PARAMETERS       - Column parameter vector
%   P                - Positive, scalar integer representing the number of symmetric innovations
%   Q                - Non-negative, scalar integer representing the number of lags of conditional variance
%
% OUTPUTS:
%   TPARAMETERS      - A 1+p+q column vector of transformed parameters corresponding to
%                      [alpha(1),...,alpha(p), beta1 ... beta(q)]'
%
% COMMENTS:
%   Input parameters must satisfy:
%    (1) alpha(i) >= 0 for i = 1,2,...,p
%    (2) beta(i)  >= 0 for i = 1,2,...,q
%    (3) sum(alpha) + sum(beta) < 1
%
% See also SCALAR_VT_VECH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 9/1/2005

if size(parameters,2)>size(parameters,1)
   parameters = parameters';
end
%Upper bound to keep it a bit away from 1
UB=.999998;

alpha=parameters(1:p);
gamma = parameters(p+1:p+o);
beta=parameters(p+o+1:p+o+q);
%Check that the parameters satisfy the necessary constraints
if  any(alpha<0) || any(beta<0) || any(gamma<0) || (sum(alpha)+sum(gamma)/kappa+sum(beta))>=UB
    error('These do not conform to the necessary set of restrictions to be transformed.')
end

%Alpha, beta cannot be exactly zero or there will be problems with log()
alpha(alpha<1e-8)=1e-8;
gamma(gamma<1e-8)=1e-8;
beta(beta<1e-8)=1e-8;
%Up the upper bound a small amount to make sure it is satisfied
UB=UB+1e-8*(p+o+q);

%Set the scale
scale=UB;
%Initialize the transformed parameters
parameters=[alpha;gamma;beta];
tparameters=[alpha;gamma;beta];
for i=1:(p)
    %Scale the parameters
    tparameters(i)=tparameters(i)/scale;
    %Use an inverse logistic
    tparameters(i)=log(tparameters(i)/(1-tparameters(i)));
    %Update the scale
    scale=scale-parameters(i);
end
for i=p+1:p+o
    %Scale the parameters
    tparameters(i)=tparameters(i)/scale/kappa;
    %Use an inverse logistic
    tparameters(i)=log(tparameters(i)/(1-tparameters(i)));
    %Update the scale
    scale=scale-parameters(i)/kappa;
end

for i=p+o+1:p+o+q
    %Scale the parameters
    tparameters(i)=tparameters(i)/scale;
    %Use an inverse logistic
    tparameters(i)=log(tparameters(i)/(1-tparameters(i)));
    %Update the scale
    scale=scale-parameters(i);
end


text:
%% GSA_GetSy: calculate the Sobol' sensitivity indices
%
% Usage:
%   [S eS pro] = GSA_GetSy(pro, iset, verbose)
%
% Inputs:
%    pro                project structure
%    iset               cell array or array of inputs of the considered set, they can be selected
%                       by index (1,2,3 ...) or by name ('in1','x',..) or
%                       mixed
%    verbose            if not empty, it shows the time (in hours) for
%                       finishing
%
% Output:
%    S                  sensitivity coefficient
%    eS                 error of sensitivity coefficient
%    pro                project structure
%
% ------------------------------------------------------------------------
% See also
%
% Author : Flavio Cannavo'
% e-mail: flavio(dot)cannavo(at)gmail(dot)com
% Release: 1.0
% Date   : 15-02-2011
%
% History:
% 1.0  15-04-2011  Added verbose parameter
% 1.0  15-01-2011  First release.
%%

function [S eS pro] = GSA_GetSy(pro, iset, verbose)

if ~exist('verbose','var')
    verbose = 0;
else
    verbose = ~isempty(verbose) && verbose;
end

index = fnc_SelectInput(pro, iset);

if isempty(index)
    S = 0;
    eS = 0;
else
    S = 0;
    eS = 0;
    n = length(index);
    L = 2^n;

    if verbose
        tic
    end
    for i=1:(L-1)
        ii = fnc_GetInputs(i);
        si = fnc_GetIndex(index(ii));
        if isnan(pro.GSA.GSI(si))

            %-------
            if isnan(pro.GSA.Di(si))

                ixi = fnc_GetInputs(si);
                s = length(ixi);
                l = 2^s - 1;

                %======
                if isnan(pro.GSA.Dmi(si))
                    n = length(pro.Inputs.pdfs);
                    N = size(pro.SampleSets.E,1);
                    H = pro.SampleSets.E(:,:);
                    cii = fnc_GetComplementayInputs(si, n);
                    H(:,cii) = pro.SampleSets.T(:,cii);
                    ff = nan(N,1);

                    for j=1:N
                        ff(j) = pro.GSA.fE(j)*(pro.Model.handle(H(j,:))-pro.GSA.mfE);
                    end

                    pro.GSA.Dmi(si)  = nanmean(ff);
                    pro.GSA.eDmi(si) = 0.6745*sqrt((nanmean(ff.^2) - pro.GSA.Dmi(si)^2)/sum(~isnan(ff)));
                end
                %=======

                Di = pro.GSA.Dmi(si);
                eDi = pro.GSA.eDmi(si)^2;

                for j=1:(l-1)
                    sii = fnc_GetInputs(j);
                    k = fnc_GetIndex(ixi(sii));
                    s_r = s - length(sii);
                    Di = Di + pro.GSA.Dmi(k)*((-1)^s_r);
                    eDi = eDi + pro.GSA.eDmi(k)^2;
                end

                pro.GSA.Di(si) = Di + (pro.GSA.f0^2)*((-1)^s);
                pro.GSA.eDi(si) = sqrt(eDi + 2*(pro.GSA.ef0^2));

            end
            %------
            pro.GSA.GSI(si) = pro.GSA.Di(si)/pro.GSA.D;
            pro.GSA.eGSI(si) = pro.GSA.GSI(si)*pro.GSA.eDi(si)/pro.GSA.D;
        end
        S = S + pro.GSA.GSI(si);
        eS = eS + pro.GSA.eGSI(si);

        if verbose
            timelapse = toc;
            disp(timelapse*(L-1-i)/i/60/60);
        end
    end
end


text:
function varargout = chebkind(kind)
%CHEBKIND   Set the default Chebyshev grid type.
%   CHEBKIND 1 or CHEBKIND 2 globally sets the default grid used by the CHERBFUN
%   constructor to be Chebyshev points of the first or second kind,
%   respectively, and is equivalent to setting the default 'tech' preference in
%   CHEBFUNPREF to be either @chebtech1 or @chebtech2.
%
%   CHEBKIND by itself returns the current default Chebyshev grid type. If the
%   default technology is not Chebyshev, then zero is returned.
%
% See also CHEBTECH1, CHEBTECH2, CHEBFUNPREF, CHEBFUNPREF.SETDEFAULTS.

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

if ( nargin == 0 )
    tech = feval(chebfunpref().tech);
    if ( isa(tech, 'chebtech') )
        tech = class(tech);
        kind = str2double(tech(end));
    else
        kind = 0;
    end
elseif ( (kind == 1) || any(strcmp(kind, {'1', '1st'})) )
    chebfunpref.setDefaults('tech', @chebtech1);
    kind = 1;
elseif ( (kind == 2) || any(strcmp(kind, {'2', '2nd'})) )
    chebfunpref.setDefaults('tech', @chebtech2);
    kind = 2;
else
    error('CHEBFUN:chebkind:unknown', ...
        'Unknown input to CHEBKIND(). Valid options are 1 or 2.');
end

if ( nargout > 0 )
    varargout{1} = kind;
end

end


text:
function [gp, indA] = passgp(gp, x, y, varargin)
%PASSGP   Optimize active set and hyperparameters of PASS-GP
%
%  Description
%    [GP, INDA] = PASSGP(GP, X, Y, OPTIONS)
%    Returns GP structure with hyperparameters optimized according to
%    PASS-GP routine (Henao & Winther, 2012) and active set indices
%    INDA for X and Y
%
%   OPTIONS is optional parameter-value pair
%      npass - Number of passes the algorithm takes over the whole training
%              data set
%      ninit - Initial active set size
%      nsub  - Number of subsets we process at each pass (in how many
%              parts we divide the whole data)
%      pinc  - Predictive density threshold for inclusion in active set.
%      pdel  - LOO-predictive density threshold for deletion from active
%              set
%      pexc  - Exchange proportion for fixed PASS-GP
%      opt   - Options structure for optimizer
%      fixed - Whether we use fixed size of active set or not. Default
%               'off'
%      display - Whether to display additional info or not. Default 'off'
%      optimn - Whether to optimize always or only after every nth
%               deletion/addition to active set. Default 1 (every time). If
%               given e.g. value 3, optimizes after every 3rd
%               addition/deletion.
%
%  See also
%    GP_SET, LIK_*
%
%  Reference:
%    Ricardo Henao & Ole Winther (2012). Preditive active set selection
%    methods for Gaussian processes. Neurocomputing 80 (2012), 10-18.

% Copyright (c) 2013 Ville Tolvanen

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

  ip=inputParser;
  ip.FunctionName = 'PASSGP';
  ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('npass', 3, @(x) isscalar(x) && x > 0)
  ip.addParamValue('ninit', 100, @(x) isscalar(x) && x > 0)
  ip.addParamValue('nsub', 10, @(x) isscalar(x) && x > 0)
  ip.addParamValue('pinc', 0.5, @(x) isscalar(x) && x > 0)
  ip.addParamValue('pdel', 0.99, @(x) isscalar(x) && x > 0)
  ip.addParamValue('pexc', 0.1, @(x) isscalar(x) && x > 0)
  ip.addParamValue('opt', [], @(x) isstruct(x))
  ip.addParamValue('fixed', 'off', @(x) ismember(x,{'on','off'}))
  ip.addParamValue('display', 'off', @(x) ismember(x,{'on','off'}))
  ip.addParamValue('optimn', 1', @(x) isscalar(x) && x > 0 && rem(10*x,2)==0)
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.parse(gp, x, y, varargin{:});
  opt=ip.Results.opt;
  fixed=ip.Results.fixed;
  options.z = ip.Results.z;
  npass=ip.Results.npass;
  ninit=ip.Results.ninit;
  nsub=ip.Results.nsub;
  pinc=ip.Results.pinc;
  pdel=ip.Results.pdel;
  pexc=ip.Results.pexc;
  display=ip.Results.display;
  optimn=ip.Results.optimn;

  if isequal(display,'on')
    display=1;
  else
    display=0;
  end
  if isequal(fixed, 'on')
    fixed=1;
  else
    fixed=0;
  end

  if ninit > size(x,1)
    error('Initial active set must be subset of original data');
  end
  if nsub > (size(x,1) - ninit)
    error('nsub must be lower than size(x,1) - ninit');
  end

  [n,nin]=size(x);

  % Initial active set

  indA=sort(randperm(n, ninit),'ascend');
  % Inclusions/deletions per iteration for fixed pass-gp
  nexc=floor(ninit*pexc);
  iter=optimn-1;
  for i=1:npass
    if display
      fprintf('Pass %d / %d.\n', i, npass)
    end
    [tmp,indSub]=cvit(n, nsub, floor(10*rand(1)));
    for j=1:nsub
      iter=iter+1;
      inds=indSub{j};
      % Remove indices that are already in active set
      inds(ismember(inds,indA))=[];

      if iter==optimn
        % Optimize hyperparameters
        gp = gp_optim(gp, x(indA,:), y(indA), 'opt', opt, options);
        iter=0;
      end

      % Calculate weights for active set inputs (loo predictive densities)
      [tmp,tmp,lpyt]=gp_loopred(gp,x(indA,:),y(indA), options);

      % Remove active set indices according to removal rule
      if ~fixed
        indA(find(exp(lpyt)>pdel))=[];
      else
        [tmp,ii]=sort(lpyt, 'descend');
        indA(ii(1:nexc))=[];
      end
      % Calculate weights for inputs not in active set (predictive density)
      [tmp,tmp,lpyt]=gp_pred(gp, x(indA,:), y(indA), x(inds,:), 'yt', y(inds), options);

      % Add indices to active set according to addition rule
      if ~fixed
        ind=find(exp(lpyt)<pinc)';
      else
        [tmp,ii]=sort(lpyt, 'ascend');
        ind=ii(1:nexc);
      end
      indA=[indA inds(ind)];

    end

  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function CPD = hhmmF_CPD(bnet, self, Qself, Fbelow, varargin)
% HHMMF_CPD Make the CPD for an F node in a hierarchical HMM
% CPD = hhmmF_CPD(bnet, self, Qself,  Fbelow, ...)
%
%        Qps
%          \
%           \
%           Fself
%         /   |
%        /    |
%       Qself Fbelow
%
% We assume nodes are ordered (numbered) as follows: Qps, Q, Fbelow, F
% All nodes numbers should be from slice 1.
%
% If Fbelow if missing, this becomes a regular tabular_CPD.
% Qps may be omitted.
%
% optional args [defaults]
%
% Qps - node numbers.
% termprob - termprob(k,i,2) = prob finishing given Q(d)=i and Q(1:d-1)=k [ finish in last state wp 0.9]
%
% hhmmF_CPD is a subclass of tabular_CPD so we inherit inference methods like CPD_to_pot, etc.
%
% We create an isolated tabular_CPD with no F parent to learn termprob
% so we can avail of e.g., entropic or Dirichlet priors.
%
% For details, see "Linear-time inference in hierarchical HMMs", Murphy and Paskin, NIPS'01.

Qps = [];
% get parents
for i=1:2:length(varargin)
  switch varargin{i},
   case 'Qps', Qps = varargin{i+1};
  end
end

ns = bnet.node_sizes(:);
Qsz = ns(Qself);
Qpsz = prod(ns(Qps));
CPD.Qsz = Qsz;
CPD.Qpsz = Qpsz;

ps = parents(bnet.dag, self);
CPD.Fbelow_ndx = find_equiv_posns(Fbelow, ps);
CPD.Qps_ndx = find_equiv_posns(Qps, ps);
CPD.Qself_ndx = find_equiv_posns(Qself, ps);

% set default arguments
p = 0.9;
%termprob(k,i,t) Might terminate if i=Qsz; will not terminate if i<Qsz
termprob = zeros(Qpsz, Qsz, 2);
termprob(:, Qsz, 2) = p;
termprob(:, Qsz, 1) = 1-p;
termprob(:, 1:(Qsz-1), 1) = 1;

for i=1:2:length(varargin)
  switch varargin{i},
   case 'termprob', termprob = varargin{i+1};
  end
end

CPD.sub_CPD_term = mk_isolated_tabular_CPD([Qpsz Qsz 2], {'CPT', termprob});
S = struct(CPD.sub_CPD_term);
CPD.termprob = S.CPT;

CPD = class(CPD, 'hhmmF_CPD', tabular_CPD(bnet, self));

CPD = update_CPT(CPD);


text:
classdef FieldPrinter < handle

    properties (Access = protected)
        fileID
        fieldValues
        fieldName
        iter
        fieldPosition
        simulationStr
        fieldRepresenter
    end

    properties (Access = protected, Abstract)
       fieldType
    end

    methods (Access = protected)

        function printValuesLine(obj)
            iD = obj.fileID;
            fprintf(iD,'Values\n');
        end

        function printEndValuesLine(obj)
            iD = obj.fileID;
            fprintf(iD,'End Values\n');
        end

        function printFieldLines(obj)
            obj.fieldRepresenter.printFieldLines()
        end

        function printResultsLineHeader(obj)
           obj.fieldRepresenter.printResultsLineHeader()
        end

    end

 end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function PlottinMaxWithPNorm

p = 64;

%nIteration = 498;
%fCase = 'ExperimentingPlot';
%folder = ['/home/alex/git-repos/Swan/Output/',fCase];

%nIteration = 545;
%fCase = 'LatticeExperimentInputCantileverSymmetricMeshSuperEllipsePDE';
%folder = '/media/alex/My Passport/LatticeResults/CantileverSymmetricMeshSuperEllipsePDEEpsilonEhGradient2000';

nIteration = 235;
fCase = 'LatticeExperimentInputCantileverSymmetricMeshSuperEllipsePDE';
folder = '/media/alex/My Passport/LatticeResults/CantileverSymmetricMeshSuperEllipsePDEGradientEpsilonH';

%nIteration = 2362;
%fCase = 'LatticeExperimentInputCantileverSymmetricMeshSuperEllipsePDE';
%folder = '/media/alex/My Passport/LatticeResults/CantileverSymmetricMeshSuperEllipsePDE10000iteration';
for iter = 1:nIteration

    s.fileName = [fCase,num2str(iter)];
    s.folderPath = fullfile(folder);

    wM = WrapperMshResFiles(s);
    wM.compute();

    mesh = wM.mesh;

    quad = Quadrature.set(mesh.type);
    quad.computeQuadrature('CONSTANT');

    dvolum = mesh.computeDvolume(quad);
    sl2Norm = wM.dataRes.StressNormGauss;
    sl2NormP = (sl2Norm).^(p);
    int = sl2NormP.*dvolum';
    intOpt = int(:);
    stressNorml2LpNorm = sum(intOpt);
    stresLpNorm(iter)  = stressNorml2LpNorm.^(1/p);
    stressMax(iter)    = max(abs(sl2Norm));
    iter

    if iter == 1 || iter/10 == floor(iter/10)
        print(folder,fCase,stresLpNorm,stressMax,iter)
    end

end
print(folder,fCase,stresLpNorm,stressMax,iter)
end

function print(folder,fCase,stresLpNorm,stressMax,nIteration)

fid = fopen(fullfile(folder,[fCase,'.txt']));
tLine = textscan(fid,'%s','delimiter','\n', 'headerlines',0);

a = tLine{1};
costA = a(end-9,1);
costAs = split(costA);
costAt = str2double(costAs(2));

a = tLine{1};
costR = a(end,1);
costRs = split(costR);
costRt = str2double(costRs(2));

const = costRt/costAt;

fNameMon = fullfile(folder,'Monitoring.fig');
h = openfig(fNameMon);
handles = findobj(h,'Type','line');
iterC = get(handles(5),'Xdata');
cost = get(handles(5),'Ydata');
close all
f = figure();
pN{1} = plot(iterC,cost*const);
hold on
pN{2} = plot(1:nIteration,stresLpNorm);
hold on
pN{3} = plot(1:nIteration,stressMax);

leg = legend({'Cost function','L^{64}(\Omega) norm','max norm'});
set(leg);

savefig(f,fullfile(folder,[fCase,'CostVsLinf.fig']))

p = plotPrinter(f,pN);
p.print(fullfile(folder,[fCase,'Cost']));

end


text:
function datToNPY(inFilename, outFilename, dataType, shape, varargin)
% function datToNPY(inFilename, outFilename, shape, dataType, [fortranOrder, littleEndian])
%
% make a NPY file from a flat binary file, given that you know the shape,
% dataType, ordering, and endianness of the flat binary file.
%
% The point here is you don't want to read in all the data from the
% existing binary file - instead you can just create the appropriate header
% and then concatenate it with the data.
%
% ** completely untested

if ~isempty(varargin)
    fortranOrder = varargin{1}; % must be true/false
    littleEndian = varargin{2}; % must be true/false
else
    fortranOrder = true;
    littleEndian = true;
end

header = constructNPYheader(dataType, shape, fortranOrder, littleEndian);

% ** TODO: need to put the header into a temp file instead, in case the
% outFilename is the same as the inFilename (and then delete the temp file
% later)
fid = fopen(tempFilename, 'w');
fwrite(fid, header, 'uint8');
fclose(fid)

str = computer;
switch str
    case {'PCWIN', 'PCWIN64'}
        [~,~] = system(sprintf('copy /b %s+%s %s', tempFilename, inFilename, outFilename));
    case {'GLNXA64', 'MACI64'}
        [~,~] = system(sprintf('cat %s %s > %s', tempFilename, inFilename, outFilename));

    otherwise
        fprintf(1, 'I don''t know how to concatenate files for your OS, but you can finish making the NPY youself by concatenating %s with %s.\n', tempFilename, inFilename);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function det = sgbdi ( abd, lda, n, ml, mu, ipvt )

%*****************************************************************************80
%
%% SGBDI computes the determinant of a band matrix factored by SGBCO or SGBFA.
%
%  Discussion:
%
%    If the inverse is needed, use SGBSL N times.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real ABD(LDA,N), the output from SGBCO or SGBFA.
%
%    Input, integer LDA, the leading dimension of the array ABD.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, MU, the number of diagonals below and above the
%    main diagonal.  0 <= ML < N, 0 <= MU < N.
%
%    Input, integer IPVT(N), the pivot vector from SGBCO or SGBFA.
%
%    Output, real DET(2), the determinant of the original matrix.
%      determinant = DET(1) * 10.0**DET(2)
%    with  1.0 <= abs ( DET(1) ) < 10.0 or DET(1) = 0.0.
%
  ten = 10.0;

  m = ml + mu + 1;
  det(1) = 1.0;
  det(2) = 0.0;

  for i = 1 : n

    if ( ipvt(i) ~= i )
      det(1) = -det(1);
    end

    det(1) = abd(m,i) * det(1);

    if ( det(1) == 0.0 )
      return
    end

    while ( abs ( det(1) ) < 1.0 )
      det(1) = ten * det(1);
      det(2) = det(2) - 1.0;
    end

    while ( ten <= abs ( det(1) ) )
      det(1) = det(1) / ten;
      det(2) = det(2) + 1.0;
    end

  end

  return
end

repo: jburkardt-m


text:
function time=str2time(str)

global gls
time=gls.gtime;

ep(1:6)=0;
for i=1:6
    [tmp,str]=strtok(str); %#ok
    if isempty(tmp)
        break;
    end
    ep(i)=str2double(tmp);
end

idx=(ep==NaN); %#ok
if any(idx),return;end

if ep(1)<100
    if ep(1)<80
        ep(1)=ep(1)+2000;
    else
        ep(1)=ep(1)+1900;
    end
end

time=epoch2time(ep);

return

repo: GINav


text:
% Function to get all files only under specific directory
% (including subdirectories)
%
% ref: stackoverflow
%      2652630/how-to-get-all-files-under-a-specific-directory-in-matlab

function fileList = getAllFiles(dirName,fileExtension)

  dirData = dir([dirName '/' fileExtension]);      %# Get the data for the current directory
  dirIndex = [dirData.isdir];  %# Find the index for directories
  fileList = {dirData(~dirIndex).name}';  %'# Get a list of the files
  if ~isempty(fileList)
    fileList = cellfun(@(x) fullfile(dirName,x),...  %# Prepend path to files
                       fileList,'UniformOutput',false);
  end
  subDirs = {dirData(dirIndex).name};  %# Get a list of the subdirectories
  validIndex = ~ismember(subDirs,{'.','..'});  %# Find index of subdirectories
                                               %#   that are not '.' or '..'
  for iDir = find(validIndex)                  %# Loop over valid subdirectories
    nextDir = fullfile(dirName,subDirs{iDir});    %# Get the subdirectory path
    fileList = [fileList; getAllFiles(nextDir,fileExtension)];  %# Recursively call getAllFiles
  end

end


text:
% This demo loads two color images and efficiently computes the emd_hat
% between them. The ground distance is a thresholded linear combination of
% the spatial and color distance between pixels. This is similar to the
% distance I used for color images in my ICCV paper. Image sizes do not
% need to be the same. As a color distance I used CIEDE2000.
% emd_hat is described in the paper:
%  A Linear Time Histogram Metric for Improved SIFT Matching
%  Ofir Pele, Michael Werman
%  ECCV 2008
% The efficient algorithm is described in the paper:
%  Fast and Robust Earth Mover's Distances
%  Ofir Pele, Michael Werman
%  ICCV 2009
% CIEDE2000 is described in the papers:
% The development of the CIE 2000 colour-difference formula: CIEDE2000
%  M. R. Luo, G. Cui, B. Rigg
%  CRA 2001
% The CIEDE2000 color-difference formula: Implementation notes, supplementary test data, and mathematical observations
%  Gaurav Sharma, Wencheng Wu, Edul N. Dalal
%  CRA 2004

clc; close all; clear all;

addpath ../

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Distance parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In the paper I used threshold=20, but I didn't use alpha_color and
% 1-alpha_color for the spatial. I just added them. So using threshold=10
% is equivalent to what I used in the paper.
threshold= 10;
alpha_color= 1/2;
% Center. In the ICCV paper center and identity were equivalent as the image
% sizes were the same.
coordinates_transformation= 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
im1= imresize( imread('1.jpg') , 1/30);
im2= imresize( imread('2.jpg') , 1/30);
% 3.jpg is more similar to 1.jpg and thus the running time will be longer.
%im2= imresize( imread('3.jpg') , 1/30);
im1_Y= size(im1,1);
im1_X= size(im1,2);
im1_N= im1_Y*im1_X;
im2_Y= size(im2,1);
im2_X= size(im2,2);
im2_N= im2_Y*im2_X;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EMD input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P= [ ones(im1_N,1)  ;  zeros(im2_N,1) ];
Q= [ zeros(im1_N,1) ;  ones(im2_N,1)  ];

im1= double(im1)./255;
im2= double(im2)./255;
cform = makecform('srgb2lab');
im1_lab= applycform(im1, cform);
im2_lab= applycform(im2, cform);

% Creating the ground distance matrix.
% Loops in Matlab are very time consuming, so I use mex.
% Matlab corresponding code is commented out below.
tic
[ground_distance_matrix]= color_spatial_EMD_ground_distance(im1_lab,im2_lab,...
                                                  alpha_color,threshold, ...
                                                  coordinates_transformation);
fprintf(1,'computing the ground distance matrix, time in seconds: %f\n', ...
        toc);

% In ICCV paper extra_mass_penalty did not matter as image sizes were the same.
extra_mass_penalty= -1;
flowType= 3;

tic
[emd_hat_mex_val_with_flow F]= emd_hat_mex(P,Q,ground_distance_matrix,extra_mass_penalty,flowType);
fprintf(1,'Note that the ground distance here is not a metric.\n');
fprintf(1,'emd_hat_mex, time in seconds: %f\n',toc);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (c) 2009-2012, Ofir Pele
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%    * Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
%    * Neither the name of the The Hebrew University of Jerusalem nor the
%    names of its contributors may be used to endorse or promote products
%    derived from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
%--------------------------------------------------------------------------
%   [yn,ref_dn,W,en] = lms2(xn,dn,M,mu)
%--------------------------------------------------------------------------
%   功能：
%   LMS自适应滤波器 带通道时延
%--------------------------------------------------------------------------
%   输入:
%           xn  输入信号的序列 列向量
%           dn  所期望的响应序列 列向量
%           M   滤波器阶数
%           mu  收敛因子(步长)
%   输出:
%           W   滤波器的权值矩阵     大小为M * iter
%           en  误差序列（itr * 1） 列向量
%           yn  实际输出序列        列向量
%--------------------------------------------------------------------------
%   例子：
%   暂无
%--------------------------------------------------------------------------
function [yn,ref_dn,W,en] = lms2(xn,dn,M,mu)
if rem(M,2)==0
    disp('请将滤波器阶数为奇数')
    return
end

w_ref = [zeros((M-1)/2,1);1;zeros((M-1)/2,1)];                              %正常延迟信号 延迟(M-1)/2个点
xn = xn(:);
dn =conv(w_ref,dn(:));                                                      %正常信号延迟
ref_dn = dn;
itr = length(xn);
en = zeros(itr,1);                                                          %误差序列，en(k)表示第k次迭代时预期输出与实际输出的误差
W = zeros(M,itr);                                                           %滤波器的权值矩阵（矩阵）
%--------------------------------------------------------------------------
%   迭代计算
%--------------------------------------------------------------------------
for idx = M:itr                                                             %第k次迭代
    x = xn(idx:-1:idx-M+1);                                                 %滤波器M个抽头的输入
    y = W(:,idx-1).'*x;                                                     %滤波器的输出
    en(idx) = dn(idx) - y;                                                  %第k次迭代的误差
    %----------------------------------------------------------------------
    %   滤波器权值计算的迭代式
    %----------------------------------------------------------------------
    W(:,idx) = W(:,idx-1) + 2*mu*en(idx)*x;
end

%--------------------------------------------------------------------------
%   求最优时滤波器的输出序列，r如果没有yn的返回参数可以不要下面的
%--------------------------------------------------------------------------
yn = conv(W(:,end),xn);

repo: SP


text:
function [angle, check] = WardSimpleRotAux(img1, img2, rect)
%
%
%       rot = WardSimpleRotAux(img1, img2, rect)
%
%       This function computes the Ward's MTB.
%
%       Input:
%           -img1: the target image
%           -img2: the image that needs to be aligned to img1
%
%       Output:
%           -rot: rotation angle (degree) for aligning img2 into img1.
%
%     Copyright (C) 2013-15 Francesco Banterle. A big thank to Greg J. Ward
%     for help during the implementation.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

[r, c, ~] = size(img1);

maxDivergence = 0.005;

%First block
r_img1  = img1(rect(1):rect(2),rect(3):rect(4),:);
r_img2  = img2(rect(1):rect(2),rect(3):rect(4),:);
r1_shift = WardGetExpShift(r_img1, r_img2);

%Mirror block
rect_mirror(1) = r - rect(2) + 1;
rect_mirror(2) = r - rect(1) + 1;
rect_mirror(3) = c - rect(4) + 1;
rect_mirror(4) = c - rect(3) + 1;

r_img1  = img1(rect_mirror(1):rect_mirror(2),rect_mirror(3):rect_mirror(4),:);
r_img2  = img2(rect_mirror(1):rect_mirror(2),rect_mirror(3):rect_mirror(4),:);
r2_shift = WardGetExpShift(r_img1, r_img2);

dx = rect_mirror(3) - rect(3);
dy = rect_mirror(1) - rect(1);

dxr = dx + 0.5 * (r2_shift(1) - r1_shift(1));
dyr = dy + 0.5 * (r2_shift(2) - r1_shift(2));

value = abs(sqrt((dxr * dxr + dyr * dyr) / (dx * dx + dy * dy)) - 1.0);

if(value <= maxDivergence)
    angle = atan2(dyr, dxr) - atan2(dy, dx);
    angle = (angle * 180.0) / pi;
    check = 1;
else
    angle = 0.0;
    check = 0;
end

end

author: banterle


text:
function [y,deriv] = affineTrans(w,affineMap,linMap,transMap)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% Applies affine transform y = affineMap(w). It needs also needs
% linMap, the linear part of the mapping, as well as transMap, the
% transpose of linMap. All of affineMap, linMap and transMap are function
% handles.
%
% Note, linMap(x) =  J*x where J is the Jacobian of affineMap; and
% transMap(y) = J'y.

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w)affineTrans(w,affineMap,linMap,transMap);
    return;
end

if isa(w,'function_handle')
    outer = affineTrans([],affineMap,linMap,transMap);
    y = compose_mv(outer,w,[]);
    return;
end

y = affineMap(w);
y = y(:);

deriv = @(g2) deriv_this(g2,linMap,transMap);

function [g,hess,linear] = deriv_this(g2,linMap,transMap)
g = transMap(g2);
g = g(:);
%linear = false;  % use this to test linearity of affineMap, if in doubt
linear = true;
hess = @(d) hess_this(linMap,d);

function [h,Jd] = hess_this(linMap,d)
%h=zeros(size(d)); % use this to test linearity of affineMap, if in doubt
h = [];
if nargout>1
    Jd = linMap(d);
    Jd = Jd(:);
end

function test_this()
A = randn(4,5);
k = randn(4,1);
affineMap = @(w) A*w+k;
linMap = @(w) A*w;
transMap = @(y) (y.'*A).'; % faster than A'*y, if A is big
f = affineTrans([],affineMap,linMap,transMap);
test_MV2DF(f,randn(5,1));

author: nesl