text:
function test_bug1249

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_componentanalysis ft_rejectcomponent

load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf275.mat'));

cfg = [];
cfg.method = 'fastica';
cfg.numcomponent = 20;
cfg.channel = 'MEG';
comp = ft_componentanalysis(cfg, data);

cfg = [];
cfg.component = 1:5;
datanew = ft_rejectcomponent(cfg, comp);

if size(datanew.grad.chanpos,1)~=numel(datanew.grad.label)
  error('labels in the grad structure are inconsistent with the chanpos');
end


text:
function [colors, allRoiVertInds, dataMask, vw, newColors] = meshDrawROIs(vw, colors, v2g, perimThickness, dataMask, prefs)
%
% Overlay ROI colors on top of the data color overlay
%
% NOTE: an ROI named 'mask' (case-insensitive) will NOT be drawn- it will
% be used to mask away data outside the mask ROI.
%
% NOTE: the following code is slow and should be optimized.
%
% HISTORY:
%  broken off by ras, sometime in 2006.
%  ras 03/07: updated ROI outline parsing to sync with the new way of
%             specifying what ROIs to show (see roiSetOptions).
%  2007.08.23 RFD: fixed broken data masking code.
%  ras, 11/07: streamlined code (broke off two subfunctions), fixed a bug
%				in which perimeters of multiple ROIs were interfering with
%				one another.
%  2009.12.30 JW: allow filling of single ROIs (with name containing
%               string 'fill'), even if other ROIs are drawn as perimeter.
%               Useful for masking out artifacts.
%

newColors = zeros(size(colors));
% If we want the ROIs to be translucent patches, then make the color of the
% ROIs a weighted sum of the ROI color and the underlay color. If we are
% not drawing the ROIs as translucent patches, then make the weight 1 such
% that we are painting the ROI voxels (presumably just the boundaries) only
% the color of the ROI, with no underlay.
if strcmpi('patches', viewGet(vw, 'roidrawmethod'))
    w = .3;
else
    w = 1;
end

% get list of ROIs to show
roiList = viewGet(vw, 'ROIsToDisplay');

% roiList can (should) not be 0
if(roiList == 0)
    roiList = [];
end

% Put any mask ROIs at the end of the list.
if ~isempty(roiList)
    ROInames = viewGet(vw, 'ROI names');
    maskRoiIndex = strcmpi('mask',ROInames(roiList));
    if any(maskRoiIndex)
        maskRoiIndex = ismember(roiList, maskRoiIndex);
        [~,inds] = sort(maskRoiIndex);
        roiList = roiList(inds);
    end
end

if notDefined('perimThickness')
    % Option to display only the ROI perimeter.
    % '1' is a good number to use here.
    % ras 01/06: inherits from the view's settings
    if checkfields(vw, 'ui', 'roiDrawMethod')
        switch lower(viewGet(vw, 'roi draw method'))
            case 'boxes', perimThickness = 0;  % filled
            case 'perimeter', perimThickness = 1; % thin outlines
            case 'filled perimeter', perimThickness = 2; % thick outlines
            case 'patches', perimThickness = 0; % filled, translucent patches
        end
    end
end

if ~exist('prefs', 'var'),
    % this takes a LONG time!
    prefs = mrmPreferences;
end

% we'll need some standard colors for mapping labels
% (e.g. 'y') to vals (e.g. [1 1 0]):
[stdColVals, stdColLabs]  = meshStandardColors;

%%%%%%Main part:
% Substitute the ROI color for appropriate nearest neighbor vertices
allRoiVertInds = [];
if ~isempty(viewGet(vw, 'ROIs'))

    %% loop across ROIs
    for iR = roiList
        roi = viewGet(vw, 'ROIs', iR);
        %% find indices of those mesh vertices belonging to this ROI:
        % Note this is a new method of getting roiVertInds. We now use a
        % viewGet. We can also use a viewSet in order to store the indices
        % so that they do not need to be calculated each time. See, e.g,
        % roiSetVertInds.m
        roiVertInds = viewGet(vw, 'roiVertInds', iR, prefs);

        %% modify mesh colors using roiVertInds
		% If this is a mask ROI, don't show the ROI, just mask the data.
		if strfind(lower(roi.name), 'mask')
			tmp = false(size(dataMask));
			tmp(roiVertInds) = true;
			dataMask = dataMask & tmp;
			clear tmp;
            % The lines below
            % roiMaskIndices = roiVertInds;
            % roiMaskIndices = adjustPerimeter(roiMaskIndices, [], vw, prefs);
            % dataMask(roiMaskIndices) = 1;
			continue
		end

		% find ROI perimeter / dilate ROI, depending on the settings
        %   if the roi name contains the string 'fill', we display the ROI
        %   as a solid color instead of a perimeter
        if strfind(lower(roi.name), 'fill'), p = []; else p = perimThickness; end
        roiVertInds = adjustPerimeter(roiVertInds, p, vw, prefs);

		% expand dataMask so that the entire ROI will be drawn
        if exist('roiMaskIndices', 'var'),
    		dataMask(intersect(roiVertInds,roiMaskIndices)) = 1;
        else
            dataMask(roiVertInds) = 1;
        end

		%% modify the colors for this ROI's vertices
        if (ischar(roi.color))
            colVal = stdColVals(:, strfind(stdColLabs, roi.color));
        else
            % Check to see if we have a 3x1 vector.
            colVal = roi.color(:);
            if (length(colVal)~=3)
                error('Bad ROI color');
            end
        end

		newColors(1:3,roiVertInds) = repmat(colVal*255, 1, length(roiVertInds));

        newColors(newColors>255) = 255;
        newColors(newColors<0) = 0;

        colors(1:3,roiVertInds) = w*newColors(1:3,roiVertInds) + (1-w)*colors(1:3,roiVertInds);
		%% store this ROI's vertices in a list of all ROI vertex locations
		allRoiVertInds = [allRoiVertInds(:); roiVertInds(:)];
    end
end

colors(colors>255) = 255;
colors(colors<0) = 0;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% cross_validation on the MOT benchmark
function MOT_cross_validation

% set is_train to 0 if testing trained trackers only
is_train = 1;
opt = globals();

mot2d_train_seqs = {'TUD-Stadtmitte', 'TUD-Campus', 'PETS09-S2L1', ...
   'ETH-Bahnhof', 'ETH-Sunnyday', 'ETH-Pedcross2', 'ADL-Rundle-6', ...
   'ADL-Rundle-8', 'KITTI-13', 'KITTI-17', 'Venice-2'};

% training and testing pairs
seq_idx_train = {{1}, {4},    {7},     {9}};
seq_idx_test  = {{2}, {5, 6}, {8, 11}, {10}};

seq_set_test = 'train';
N = numel(seq_idx_train);

% for each training-testing pair
for i = 1:N
    % training
    idx_train = seq_idx_train{i};

    if is_train
        % number of training sequences
        num = numel(idx_train);
        tracker = [];

        % online training
        for j = 1:num
            fprintf('Online training on sequence: %s\n', mot2d_train_seqs{idx_train{j}});
            tracker = MDP_train(idx_train{j}, tracker);
        end
        fprintf('%d training examples after online training\n', size(tracker.f_occluded, 1));

    else
        % load tracker from file
        filename = sprintf('%s/%s_tracker.mat', opt.results, mot2d_train_seqs{idx_train{end}});
        object = load(filename);
        tracker = object.tracker;
        fprintf('load tracker from file %s\n', filename);
    end

    % testing
    idx_test = seq_idx_test{i};
    % number of testing sequences
    num = numel(idx_test);
    for j = 1:num
        fprintf('Testing on sequence: %s\n', mot2d_train_seqs{idx_test{j}});
        MDP_test(idx_test{j}, seq_set_test, tracker);
    end
end

% evaluation for all test sequences
benchmark_dir = fullfile(opt.mot, opt.mot2d, seq_set_test, filesep);
seqs = {'TUD-Campus', 'ETH-Sunnyday', 'ETH-Pedcross2', ...
   'ADL-Rundle-8', 'Venice-2', 'KITTI-17'};
evaluateTracking(seqs, opt.results, benchmark_dir);

author: yuxng


text:
function H= plot_scalpPatternsPlusChannel(erp, mnt, clab, ival, varargin)
%PLOT_SCALPPATTERNSPLUSCHANNEL - Display Classwise topographies and one channel
%
%Synposis:
% H= plot_scalpPatternsPlusChannel(ERP, MNT, CLAB, IVAL, <OPTS>)
%
%Input:
% ERP:  struct of epoched EEG data
% MNT:  struct defining an electrode montage
% CLAB: label of the channel(s) which are to be displayed in the
%        ERP plot.
% IVAL: time interval for which scalp topographies are to be plotted.
%       May be either one interval for all Classes, or specific
%       intervals for each Class. In the latter case the k-th row of IVAL
%       defines the interval for the k-th Class.
% OPTS: struct or property/value list of optional fields/properties:
%  .LegendPos - specifies the position of the legend in the ERP plot,
%               default 'best' (see man page of legend for choices).
%  .MarkIval  - When true, the time interval is marked in the channel plot.
%
%The opts struct is passed to plot_scalpPattern
%
%Output:
% H: Handle to several graphical objects
%
%See also plot_scalpPatterns, plot_scalpEvolutionPlusChannel, plot_scalp.

% Author(s): Benjamin Blankertz, Jan 2005

props = {'LineWidth',           3,                  'DOUBLE';
         'IvalColor',           0.85*[1 1 1],       'DOUBLE[3]';
         'ColorOrder',          [0 0 0],            'DOUBLE[3]';
         'MarkIval',            0,                  'BOOL';
         'XGrid',               'on',               'CHAR';
         'YGrid',               'on',               'CHAR';
         'XUnit',               '[ms]',             'CHAR';
         'PlotChannel',         1,                  'BOOL';
         'ScalePos',            'vert',             'CHAR';
         'LegendPos',           'best',             'CHAR';
         'NewColormap',         0,                  'BOOL';
         'MarkPatterns',        [],                 'DOUBLE|CHAR';
         'MarkStyle',           {'LineWidth',3},    'CELL';
         'Subplot',             [],                 'DOUBLE'};

props_scalpPattern= plot_scalpPattern;
props_channel= plotutil_channel1D;

if nargin==0,
  H= opt_catProps(props, props_scalpPattern, props_channel); return
end

opt= opt_proplistToStruct(varargin{:});
[opt,isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props, props_scalpPattern, props_channel);
[opt,isdefault] = opt_overrideIfDefault(opt,isdefault, 'MarkIval',size(ival,1)==1,'PlotChannel',~isempty(clab));

opt_scalpPattern= opt_substruct(opt, props_scalpPattern(:,1));
opt_channel= opt_substruct(opt, props_channel(:,1));

misc_checkType(erp,'!STRUCT(x y className)');
misc_checkType(mnt,'!STRUCT(x)');

if nargin<4 || isempty(ival),
  ival= [NaN NaN]; %erp.t([1 end]);
elseif size(ival,2)==1,
  ival= [ival ival];
elseif size(ival,2)>2,
  error('IVAL must be sized [N 2]');
end

if isfield(erp, 'XUnit'),
  [opt,isdefault]= opt_overrideIfDefault(opt, isdefault, 'XUnit', erp.xUnit);
end

if max(sum(erp.y,2))>1,
  erp= proc_average(erp);
end
nClasses= length(erp.className);

% Diese Abfrage ist nicht korrekt, da alle Felder schon gesetzt sein
% muessen
if ~isfield(opt, 'CLim') && all(erp.x>=0),
  opt.CLim= 'range';
  warning('automatically set opt.CLim= ''range''.');
end

if ischar(opt.MarkPatterns),
  if ~strcmpi(opt.MarkPatterns, 'all'),
    warning('unknown value for opt.MarkPatterns ignored');
  end
  opt.MarkPatterns= 1:nPat;
end

[AxesStyle, lineStyle]= opt_extractPlotStyles(opt);
nIvals= size(ival,1);
if nIvals<nClasses,
  ival= repmat(ival, [nClasses, 1]);
end

if opt.NewColormap,
  acm= visutil_addColormap(opt.Colormap);
end
if isempty(opt.Subplot),
  clf;
end
for cc= 1:nClasses,
  if isempty(opt.Subplot),
    subplotxl(1, nClasses+opt.PlotChannel, cc+opt.PlotChannel, ...
              0.07, [0.07 0.02 0.1]);
  else
    axis_getQuietly(opt.Subplot(cc));
  end
  opt_scalpPattern= setfield(opt_scalpPattern, 'ScalePos','none');
  if opt.NewColormap,
    opt_scalpPattern= rmfield(opt_scalpPattern, 'Colormap');
    opt_scalpPattern= setfield(opt_scalpPattern, 'NewColormap',0);
  end
  opt_scalpPattern.Class= cc;
  clscol= opt.ColorOrder(min(cc,size(opt.ColorOrder,1)),:);
  opt_scalpPattern= opt_setDefaults(opt_scalpPattern, 'Linespec', {'linewidth',2, 'Color',clscol});
  h.H_topo(cc)= plot_scalpPattern(erp, mnt, ival(cc,:), opt_scalpPattern);
  if ismember(cc, opt.MarkPatterns,'legacy'),
    set([h.H_topo(cc).head h.H_topo(cc).nose], opt.MarkStyle{:});
  end
  yLim= get(gca, 'yLim');
  h.text(cc)= text(mean(xlim), yLim(2)+0.06*diff(yLim), erp.className{cc});
  set(h.text(cc), 'Color',clscol);
  if ~any(isnan(ival(cc,:))),
    ival_str= sprintf('%g - %g %s', ival(cc,:), opt.XUnit);
  else
    ival_str= '';
  end
  h.text_ival(cc)= text(mean(xlim), yLim(1)-0.04*diff(yLim), ival_str);
%  axis_aspectRatioToPosition;   %% makes Colorbar appear in correct size
end
set(h.text, 'horizontalAli','center', ...
            'Visible','on', ...
            'FontSize',12, 'fontWeight','bold');
set(h.text_ival, 'verticalAli','top', 'horizontalAli','center', ...
                 'Visible','on');
if ismember(opt.ScalePos, {'horiz','vert'},'legacy'),
  h.cb= plotutil_colorbarAside(opt.ScalePos);
  % hack to fix a matlab bug
  ud= get(h.cb, 'UserData');
  ud.orientation= opt.ScalePos;
  set(h.cb, 'UserData',ud);
  % put YUnit on top of Colorbar
  if isfield(erp, 'yUnit'),
    axes(h.cb);
    yl= get(h.cb, 'YLim');
    h.YUnit= text(mean(xlim), yl(2), erp.yUnit);
    set(h.YUnit, 'horizontalAli','center', 'verticalAli','bottom');
  end
  if nClasses>1,
%    visutil_unifyCLim([h.H_topo.ax], [zeros(1,nClasses-1) h.cb]);
    visutil_unifyCLim([h.H_topo.ax]);
  end
elseif nClasses>1,
  visutil_unifyCLim([h.H_topo.ax], [zeros(1,nClasses)]);
end
if opt.NewColormap,
  visutil_acmAdaptCLim(acm, [h.H_topo.ax]);
end

if opt.PlotChannel,
  if isempty(opt.Subplot),
    h.ax_erp= subplotxl(1, nClasses+1, 1, 0.1, [0.06 0.02 0.1]);
  else
    h.ax_erp= opt.Subplot(nClasses+1);
    axes(h.ax_erp);
  end
  topopos= get(h.H_topo(end).ax, 'position');
  pos= get(h.ax_erp, 'position');
  pos([2 4])= topopos([2 4]);
  set(h.ax_erp, 'position',pos);
  if ~isempty(AxesStyle),
    set(h.ax_erp, AxesStyle{:});
  end
  hold on;   %% otherwise axis properties like ColorOrder are lost
  H.channel= plot_channel(erp, clab, opt_channel, 'legend',0);
  if opt.MarkIval,
    for cc= 1:nIvals,
      grid_markInterval(ival(cc,:), clab, ...
                    opt.IvalColor(min(cc,size(opt.IvalColor,1)),:));
    end
    axis_redrawFrame(h.ax_erp);
  end
  set(get(h.ax_erp, 'title'), 'FontSize',12, 'fontWeight','bold');
  if ~isequal(opt.LegendPos, 'none'),
    h.leg= legend(erp.className, 'Location',opt.LegendPos);
  end
end

if nargout<1,
  clear h
end

author: bbci


text:
function [grad, grad_W, grad_b, validFrameMask]= B_affine_transform(input_layers, curr_layer, future_layers, skip_grad)
if length(input_layers)==1
    input = input_layers{1}.a;
    W = curr_layer.W;
    b = curr_layer.b;
elseif length(input_layers)==2
    W = input_layers{1}.a;
    b = zeros(size(W,1),1);
    input = input_layers{2}.a;
end
validFrameMask = [];

if curr_layer.update==0
    grad_W = [];
    grad_b = [];
end

future_grad = GetFutureGrad(future_layers, curr_layer);

[n1,n2,n3] = size(future_grad);
if n3>1     % reshape the matrix to 2D
    [validFrameMask, variableLength] = getValidFrameMask(input_layers{length(input_layers)});
    if variableLength
        future_grad = PadShortTrajectory(future_grad, validFrameMask, 0);
%         input = PadShortTrajectory(input, validFrameMask, 0);
%         future_grad = future_grad(:,validFrameMask==0);
%         input = input(:,validFrameMask==0);
    end
    future_grad = reshape(future_grad, n1,n2*n3);
    input = reshape(input, size(input,1), n2*n3);
end

if curr_layer.update
    if issparse(input)
        visible_nonzero_idx = find(sum(abs(input),2)>0);
        visible_nonzero = full(input(visible_nonzero_idx,:));
        grad_W = sparse(size(W,1),size(W,2));
        grad_W(:,visible_nonzero_idx) = future_grad * visible_nonzero';
    else
        grad_W = conj(future_grad * input');
    end
    if isfield(curr_layer, 'mask')
        grad_W = grad_W .* curr_layer.mask;
    end
    if isfield(curr_layer, 'updateBias') && ~curr_layer.updateBias  % sometimes, we don't use bias
        grad_b = future_grad(:,1)*0;
    else
        grad_b = sum(future_grad,2);
    end
end
if isfield(curr_layer, 'skip_grad')
    skip_grad = skip_grad + curr_layer.skip_grad;
end
if skip_grad==0
    grad = W' * future_grad;
    if n3>1
        grad = reshape(grad, size(grad,1), n2, n3);
    end
else
    grad = [];
end

end

author: singaxiong

repo: SignalGraph


text:
function [x, infos] = rank2nmf(V, in_options)
% Rank-two NMF computes a rank-two NMF of X
%
%       If X >= 0 and rank(X) = 1, take W as any non-zero column of X and compute
%       H accordingly
%
%       If X >= 0 and rank(X) = 2, this is Algorithm 4.1 in Chapter 4,
%       with alpha1 = alpha2 = 0. This is an exact algorithm.
%
%       Otherwise, the input matrix has is replaced by max(0,X2) where X2 is its
%       best rank-two approximation; cf. Chapter 6.
%
% Inputs:
%       matrix      V
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       Nicolas Gillis,
%       "Nonnegative Matrix Factorization,"
%       SIAM, 2020.
%
%
% This file is part of NMFLibrary.
%
% This file has been ported from
%       Rank2NMF.m at https://gitlab.com/ngillis/nmfbook/-/tree/master/algorithms
%       by Nicolas Gillis (nicolas.gillis@umons.ac.be)
%
% Change log:
%
%       June 13, 2022 (Mitsuhiko Horie and Hiroyuki Kasai): Ported initial version
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.delta = 1e-6;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    if options.verbose > 0
        fprintf('# Rank2-NMF: started ...\n');
    end

    if min(V(:)) < 0
        if options.verbose > 0
            fprintf('  * Rank2-NMF: The input matrix is not nonnegative: V <-- max(V,0).\n');
        end
        V = max(0,V);
    end
    [u, s, v] = svds(V, 3);

    % Case 0. rank(V) = 0
    if s(1, 1) <= eps
        if options.verbose > 0
            fprintf('  * Rank2-NMF: The input matrix is (close to) zero.');
        end
        [m,n] = size(V);
        W = zeros(m,1);
        H = zeros(1,n);
        return;
    end

    % Case 1. rank(V) = 1
    if s(2, 2) < 1e-9 * s(1, 1)
        if options.verbose > 0
            fprintf('  * Rank2-NMF: The input matrix has rank one, an optimal solution is returned.\n');
        end
        sumV = sum(V);
        [~,j] = max(sumV);
        W = V(:,j);
        H = sumV/sum(W);
        return;
    end

    % Case 2. rank(V) >= 2
    if s(3, 3) > 1e-9 * s(2, 2)
        if options.verbose > 0
	        fprintf('  * Rank2-NMF: The input matrix does not have rank two.\n');
            fprintf('  * Rank2-NMF: V <-- rank(V2,0) where V2 is the rank-2 truncated SVD of V.\n');
            fprintf('  * Rank2-NMF: The solution computed might not be optimal.\n');
        end
        [u, s, v] = svds(V, 2);
        Vj = max(u*s*v', 0);
    else
        if options.verbose > 0
            fprintf('  * Rank2-NMF: The input matrix has rank two, an optimal solution is returned.\n');
        end
        Vj = V;
    end

    % remove zero columns
    sumVj = sum(Vj);
    J = find(sumVj > 0);
    Vj = Vj(:, J);

    % normalize columns
    [m, n] = size(Vj);
    for i = 1 : n
        Vj(:, i) = Vj(:, i) / sum(Vj(:, i));
    end

    % compute extreme points of conv(Vj) hence a feasible solution W
    % given that rank(V) = 2.
    W = zeros(m, 2);
    [~, j1] = max(sum(Vj.^2));
    W(:, 1) = Vj(:, j1);
    [~, j2] = max(sum( (Vj-repmat(Vj(:, j1), 1, n)).^2));
    W(:, 2) = Vj(:, j2);

    W = full(W);
    WtW = W' * W;
    WtV = W' * V;

    % If no initial matrices are provided, H is initialized as follows:
    if ~isfield(options,'init') || isempty(options.init)
        H = nnls_init_nmflibrary(V, W, WtW, WtV);
    else
        H = options.init;
    end

    %% perform main routine
    options.alg_name = 'Rank2-NMF';
    options.V = V;
    options.W = W;
    options.main_routine_mode = true;
    options.eps0 = 0;
    options.eps = 1;

    [H, WtW, WtV, infos] = nnls_solver(V, W, options);

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function zoomfixticks_rev(fig,method,varargin)
%ZOOMFIXTICKS allows zoom to work on a figure whose AXES tick labels are manually set
%   ZOOMFIXTICKS(FIG) operates on the figure with handle FIG.
%   ZOOMFIXTICKS(AX) operates on the figure with child AXES handle AX.
%   ZOOMFIXTICKS with no arguments operates on the current figure.
%
%   Note that when using this function, you will need to modify the FIXTICKS
%   subfunction to meet your application needs.  Alternatively, you can comment
%   out the function and supply it as an M-file function on the MATLAB path.
%   If you supply your own, make sure that the M-file is named "fixticks.m", and
%   it takes a single input argument, the handle to an AXES object.
%
%   Example:
%       hf = figure;
%       xdata = [0:.01:1];
%       ydata = 100000*sin(2*pi*5*xdata);
%       plot(xdata,ydata);
%       zoomfixticks(hf)

%   Greg Aloe
%   8-8-2002
%   Copyright 2002-2004 The MathWorks, Inc.

%   Edited by Alex Taylor, 09-01-2004, to handle change to zoom button tag
%   between version 6 and 7

% Set up the figure if there are less than 2 arguments
if nargin < 2
    % Make sure the input is a handle to an AXES or FIGURE
    if nargin==0
        % If no arguments, assume the current figure
        fig = gcf;
    elseif ishandle(fig) & ~strcmp(lower(get(fig,'type')),'figure')
        % If FIG is a handle, but not to a figure, assume it's an axes and get its parent
        fig = get(fig,'parent');
    end

    % Error out if the input is not a handle, or not a handle to a figure or axes
    if ~ishandle(fig) | (ishandle(fig) & ~strcmp(lower(get(fig,'type')),'figure'))
        error('Input argument must be a handle to an AXES or FIGURE object.')
    end

    % Initialize the setup of the figure
    init(fig)
else
    % if 2 or more arguments, perform the desired operation
    feval(method,varargin{:})
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function init(fig)
% This function sets up the figure so that zooms will work even when the AXES
% ticklabels are set manually

%detect version
v=version;

%Get proper handles to zoom in and zoom out uitoggletool buttons.  Account
%for difference in tag names between version 6 and version 7

if (str2num(v(1))<7)
    hin = findall(fig,'tag','figToolZoomIn');
    hout = findall(fig,'tag','figToolZoomOut');
else
    hin = findall(fig,'tag','Exploration.ZoomIn');
    hout = findall(fig,'tag','Exploration.ZoomOut');
end

% Change the callback for the "Zoom In" toolbar button
set(hin,'ClickedCallback','zoomfixticks(gcbf,''zoomInCallback'',gcbo,gcbf)')

% Change the callback for the "Zoom Out" toolbar button
set(hout,'ClickedCallback','zoomfixticks(gcbf,''zoomOutCallback'',gcbo,gcbf)')

% Fix the ticks on initialization
fixticks(gca(fig))

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function zoomInCallback(hcbo,fig)
% Callback for the "Zoom In" button

% First call what is normally called during a zoom in
putdowntext('zoomin',hcbo)

% If the button is depressed, force WindowButtonDownFcn to call the DOWNFCN subfunction
if strcmp(get(hcbo,'state'),'on')
    set(fig,'WindowButtonDownFcn','zoomfixticks(gcbf,''downFcn'',gcbf)')
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function zoomOutCallback(hcbo,fig)
% Callback for the "Zoom Out" button

% First call what is normally called during a zoom out
putdowntext('zoomout',hcbo)

% If the button is depressed, force WindowButtonDownFcn to call the DOWNFCN subfunction
if strcmp(get(hcbo,'state'),'on')
    set(fig,'WindowButtonDownFcn','zoomfixticks(gcbf,''downFcn'',gcbf)')
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function downFcn(fig)
% Since the button was pressed, tell the figure that we are starting a zoom
zoom(fig,'down');

% After the zoom, fix the ticks
fixticks(gca(fig))

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function fixticks(ax)
% This function will perform the desired ticklabel fixing
% This is the function you will want to modify to meet your application needs
% Alternatively, you can comment out this function to call an M-file function
% in its place.

% EXAMPLE:
% First update the axes so the ticks are as expected
drawnow
% Get the ticks
tick=get(ax,'ytick');
% Convert ticks to a string of the desired format
tickstr=num2str(tick',7);
% Reset the labels to the new format
set(ax,'yticklabel',tickstr);

author: Sable

repo: mcbench-benchmarks


text:
function a_inverse = r8ge_inverse ( n, a_lu, pivot )

%*****************************************************************************80
%
%% R8GE_INVERSE computes the inverse of a matrix factored by R8GE_FA.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%    R8GE_INVERSE is a simplified standalone version of the LINPACK routine
%    R8GEDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix A.
%
%    Input, real A_LU(N,N), the factor information computed by R8GE_FA.
%
%    Input, integer PIVOT(N), the pivot vector from R8GE_FA.
%
%    Output, real A_INVERSE(N,N), the inverse matrix.
%
  a_inverse(1:n,1:n) = a_lu(1:n,1:n);
%
%  Compute Inverse(U).
%
  for k = 1 : n

    a_inverse(k,k) = 1.0E+00 / a_inverse(k,k);
    a_inverse(1:k-1,k) = -a_inverse(1:k-1,k) * a_inverse(k,k);

    for j = k + 1 : n

      temp = a_inverse(k,j);
      a_inverse(k,j) = 0.0E+00;
      a_inverse(1:k,j) = a_inverse(1:k,j) + a_inverse(1:k,k) * temp;

    end

  end
%
%  Form Inverse(U) * Inverse(L).
%
  for k = n - 1 : -1 : 1

    work(k+1:n) = a_inverse(k+1:n,k);
    a_inverse(k+1:n,k) = 0.0E+00;

    for j = k + 1 : n
      a_inverse(1:n,k) = a_inverse(1:n,k) + a_inverse(1:n,j) * work(j);
    end

    if ( pivot(k) ~= k )

      for i = 1 : n
        t                     = a_inverse(i,k);
        a_inverse(i,k)        = a_inverse(i,pivot(k));
        a_inverse(i,pivot(k)) = t;
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gap = evalgap(fnc, fspec, dnorm, ww, uu, A, B, lambda)

[ff,gg]=evalloss(fnc,ww,uu,A,B);

fval = ff+lambda*sum(fspec(ww));
dval = evaldual(fnc,dnorm,-gg,A,B,lambda);

gap = (fval+dval)/fval;

function [fval,gg]=evalloss(fnc, ww, uu, A, B)

if ~isempty(uu)
  zz=A*ww+B*uu;
else
  zz=A*ww;
end

[fval, gg] =fnc.p(zz, fnc.args{:});

function dval = evaldual(fnc, dnorm, aa, A, B, lambda)

mm=length(aa);

if ~isempty(B)
  aa=aa-B*((B'*B)\(B'*aa));
end

[dnm,ishard] =dnorm(A'*aa);

if ishard && dnm>0
  aa  = min(1, lambda/dnm)*aa;
  dnm = 0;
end

dval = fnc.d(aa, fnc.args{:})+dnm;

repo: REST


text:
function [X,proc]=genMarkedPointProc(lambda,dt,ts,tf,X0,mark,eventDist,timeDist,staycell)
%%GENMARKEDPOINTPROC A function to simulate a marked point process with
%                    rate parameter lambda, initial value X0, and given
%                    time parameters. Using the defaults for the last 3
%                    parameters will generate a compound Poisson process
%                    provided the mark distribution is properly defined
%                    (see note below). The null mark is assumed to be zero.
%                    See [1] for more details on general marked point
%                    processes.
%
%INPUTS: lambda The unit rate parameter for the default Poisson
%               distribution. This may be arbitrarily defined if eventDist
%               is given.
%            dt Time step for event draws.
%            ts Start time for process.
%            tf End time for process.
%            X0 The initial value for the process. Defaults to 0 if
%               omitted.
%          mark An optional function handle f(x,t,j,proc) to be used for
%               sampling mark values. This must accept three scalar inputs
%               and a structure input:
%               x The process value at the beginning of a dt time step.
%               t The elapsed time at the beginning of the dt time step,
%                 starting from ts.
%               j The index corresponding to the ordinal number of the
%                 jump since the beginning of the current time step.
%               proc A structure containing information about the process as
%                 currently constructed. See output information below for
%                 more details.
%               If mark is not given or is empty, it defaults to the unity
%               function for all potential values (jumps will have mark 1
%               regardless of the variables).
%     eventDist An optional function handle for the distribution of events.
%               This must output a scalar value in the interval [0,inf). If
%               not given, it defaults to a Poisson distribution with rate
%               parameter lambda*dt. This function handle must take as
%               input:
%               x The process value at the beginning of a dt time step.
%               t The elapsed time at the beginning of the dt time step,
%                 starting from ts.
%               i The index corresponding to the ordinal number of the
%                 time step since the beginning of the process at time
%                 ts.
%               proc A structure containing information about the process as
%                 currently constructed. See output information below for
%                 more details.
%     timeDist An optional function handle for the distribution of
%              times during a step. This must output a vector of values in
%              the interval (0,1). If not given, it defaults to a uniform
%              random distribution on (0,1). This function handle must take
%              as input:
%              x The process value at the beginning of a dt time step.
%              t The elapsed time at the beginning of the dt time step,
%                 starting from ts.
%              i The index corresponding to the ordinal number of the
%                 time step since the beginning of the process at time
%                 ts.
%              proc A structure containing information about the process as
%                 currently constructed. See output information below for
%                 more details.
%     staycell A logical value. If this is equivalent to true, then X is
%              returned as a cell array. Otherwise, X is returned as a
%              vector. Defaults to false.
%
%OUTPUTS: X Defaults to a 1Xnsteps+1 vector of the process values after
%           each step. X(1) = X0. If the staycell variable evaluates to
%           true, then X is returned as a 1Xnsteps+1 cell array.
%      proc A structure containing information which can be used to
%           reconstruct the exact process X. Members are:
%           nEvents A 1Xnsteps+1 vector of the number of jumps which
%                   occured during each step prior to the index. So
%                   nEvents(1) = 0 and nEvents(i) will be the number of
%                   jumps which occured during the (i-1) time step
%                   corresponding to the change X{i}-X{i-1}.
%               mks A 1Xnsteps cell array where each cell contains a row
%                   vector with the mark values generated for a step.
%              times A 1X((tf-ts)/dt) cell array with each cell
%                   containing a row vector containing the jump times for
%                   the process at the corresponding time step.
%
%Note: If the mark distribution is defined such that all marks are
%independently and identically distributed and the events are Poisson
%distributed, then the generated process will be a compound Poisson
%process. For more on compound Poisson processes, see pages 10-12 of [2].
%
%WARNING: No check is performed to ensure that eventDist and timeDist
%output reasonable values. This function could technically work with marks
%which are not numbers, provided some addition operation is defined for
%them. However, read  the documentation above to ensure your inputs will
%generate a proper process as intended.
%
%EXAMPLE 1: Generating several realizations of a time dependent marked
%           Poisson process and comparing it with the expected path.
% %Set path parameters
% ts = 40; % Start time
% tf = 50; % End time
% dt = 1; % Time step
% numpaths = 100; % Number of simulated paths
% t = ts:dt:tf;
% lambda1 = 5; % Poisson rate parameter for unit time
% lambda2 = 6; % Exponential rate parameter
% X0 = 100; % Initial value for the process
%
% close all
% figure
% hold on
%
% pathends = zeros([1,numpaths]);
% %Generate paths
% for i = 1:numpaths
%     [cp,proc] = genMarkedPointProc(lambda1,dt,ts,tf,X0,@(x,t,j,proc)(tf-t)*ExponentialD.rand(1,lambda2));
%     plot(t,cp)
%     pathends(i) = cp(end);
% end
% avgend = mean(pathends);
% stdend = std(pathends);
%
% %Compute the expected path
% expected = X0+[0,cumsum(lambda1*dt*(tf-(t(1:end-1)+t(2:end))/2)/lambda2)];
% expectedstd = sqrt(cumsum(lambda1*dt*((tf-(t(1:end-1)+t(2:end))/2)).^2*2/(lambda2)^2));
% plot(t,expected,'--k','LineWidth',5)
% xlabel('time (t)')
%
% %Compare simulated path statistics to the expected value
% fprintf("The average path value at t = %0.0f is : %0.5f\n",tf,avgend)
% fprintf("The path standard deviation at t = %0.0f is : %0.5f\n",tf,stdend)
% fprintf("The expected path value at t = %0.0f is : %0.5f\n",tf,expected(end))
% fprintf("The expected standard deviation at t = %0.0f is : %0.5f\n",tf,expectedstd(end))
%
%REFERENCES:
%[1] Stover, Christopher. "Marked Point Process." From MathWorld--A
%    Wolfram Web Resource, created by Eric W. Weisstein.
%    http://mathworld.wolfram.com/MarkedPointProcess.html
%[2] Platen, Eckhard, and Nicola Bruti-Liberati. Numerical solution of
%    stochastic differential equations with jumps in finance. Vol. 64.
%    Springer Science & Business Media, 2010.
%
%July 2019 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(~exist('X0','var')||isempty(X0))
    X0 = 0;
end
if(~exist('mark','var')||isempty(mark))
    mark = @(x,t,j,proc) 1;
end
if(~exist('eventDist','var')||isempty(eventDist))
    eventDist = @(x,t,i,proc) PoissonD.rand(1,lambda*dt);
end
if(~exist('timeDist','var')||isempty(timeDist))
    timeDist = @(x,t,j,proc) rand([1,proc.nEvents(j)]);
end
if(~exist('staycell','var'))
    staycell = false;
elseif(~islogical(staycell))
    error('The value in staycell is not a logical value.')
end

nsteps = floor((tf-ts)/dt);
X = cell(1,nsteps+1);
X{1} = X0;

proc.nEvents = zeros([1,length(X)]);
proc.mks = cell([1,length(X)]);
proc.times = cell([1,length(X)]);

for i = 1:nsteps
    t = (i-1)*dt+ts;
    proc.nEvents(i) = eventDist(X{i},t,i,proc);
    stepmks = cell([1,proc.nEvents(i)]);
    steptimes = sort(timeDist(X{i},t,i,proc)*dt+t);
    proc.times{i} = steptimes;

    stepIntensity = 0;
    for j = 1:proc.nEvents(i)
        stepmks{j} = mark(X{i},steptimes(j),j,proc);
        proc.mks{i}{j} = stepmks{j};
        stepIntensity = stepIntensity+stepmks{j};
    end
    X{i+1} = X{i}+stepIntensity;
end

if(~staycell)
    X = cell2mat(X);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function val= froNormMatn(x,y)

val =sqrt(sum( (x(:)-y(:)).^2));

author: thomaskuestner


text:
function [peakOutputStat] = computePeakStatistics(inputSignals,varargin)
	% Get slope ratio, the average trace from detected peaks, and other peak-related statistics
	% Biafra Ahanonu
	% started: 2013.12.09
	% inputs
		% inputSignals = [n m] matrix where n={1....N}
	% outputs
		% slopeRatio
		% traceErr
		% fwhmSignal
		% avgPeakAmplitude
		% spikeCenterTrace
		% pwelchPxx
		% pwelchf

	[peakOutputStat] = ciapkg.signal_processing.computePeakStatistics(inputSignals,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function varargout = type(varargin)
% TYPE M-file for type.fig
%      TYPE, by itself, creates a new TYPE or raises the existing
%      singleton*.
%
%      H = TYPE returns the handle to a new TYPE or the handle to
%      the existing singleton*.
%
%      TYPE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in TYPE.M with the given input arguments.
%
%      TYPE('Property','Value',...) creates a new TYPE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before type_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to type_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help type

% Last Modified by GUIDE v2.5 14-Sep-2005 17:55:01

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @type_OpeningFcn, ...
                   'gui_OutputFcn',  @type_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before type is made visible.
function type_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to type (see VARARGIN)

% Choose default command line output for type
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes type wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = type_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close main;
run('..\temp\main')

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close main
run('..\pre\main')

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close main
run('..\level\main')

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close main
run('..\level\main')


text:
function value = ball_volume_3d ( r )

%*****************************************************************************80
%
%% BALL_VOLUME_3D computes the volume of a ball in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that
%
%      X**2 + Y**2 + Z**2 <= R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the ball.
%
%    Output, real BALL_VOLUME_3D, the volume of the ball.
%
  value = ( 4.0E+00 / 3.0E+00 ) * pi * r * r * r;

  return
end


text:
function restricedRowBool = getCorrespondingRows(S, rowBool, colBool, mode)
% Returns a boolean vector that is true for a subset of the true rows in
% `rowBool` according to whether the rows 'exclusive', 'inclusive', or
% 'partial' -ly correspond to true entries in `colBool`
%
% USAGE:
%
%    restricedRowBool = getCorrespondingRows(S, rowBool, colBool, mode)
%
% INPUTS:
%    S:                     `m x n` stoichiometric matrix
%    rowBool:               `m x 1` boolean vector
%    colBool:               `n x 1` boolean vector
%    mode:                  'exclusive' , 'inclusive' or 'partial'
%
% OUTPUT:
%    restrictedRowBool:     `m x 1` boolean vector
%
% EXAMPLE:
%
%    S =
%        -1     0     0     0     0
%         2    -3     0     0     0
%         0     4    -5     0     0
%         0     0     6    -7     0
%         0     0     0     0     0
%
%    rowBool = [1; 1; 1; 1; 1];
%    colBool = [1; 1; 1; 0; 0];
%
%    % Therefore, the subset of rows and columns considered for inclusion are
%    %    -1     0     0
%    %     2    -3     0
%    %     0     4    -5
%    %     0     0     6
%    %     0     0     0
%
%    % If mode = 'exclusive' then restrictedRowBool corresponds to this subset
%    %    -1     0     0
%    %     2    -3     0
%    %     0     4    -5
%    % i.e subset of rowBool metabolites exclusively involved in the colBool
%    % reactions
%
%    % If mode = 'inclusive' then restrictedRowBool corresponds to this subset
%    %    -1     0     0
%    %     2    -3     0
%    %     0     4    -5
%    %     0     0     6
%    % i.e subset of rowBool metabolites involved in colBool reactions
%
%    % If mode ='partial' then restrictedRowBool corresponds to the extra rows
%    % with inclusive that are not present with exclusive.
%
% .. Author: - Ronan Fleming, July 2016

if ~islogical(rowBool)
    error('rowBool must be a logical vector')
end
if ~islogical(colBool)
    error('colBool must be a logical vector')
end

[mlt,nlt]=size(S);

if length(rowBool)~=mlt
    error('length of rowBool must equal size(S,1)')
end

if length(colBool)~=nlt
    error('length of rowBool must equal size(S,2)')
end

restricedRowBool=false(mlt,1);
switch mode
    case 'exclusive'
        %metatbolites exclusively involved in certain reactions
        restricedRowBool(rowBool)=     any(S(rowBool, colBool),2)...
                                    & ~any(S(rowBool,~colBool),2);
    case 'inclusive'
        %corresponding reactions involving certain metabolites
        restricedRowBool(rowBool) = any(S(rowBool, colBool),2);
    case 'partial'
        %metatbolites exclusively involved in certain reactions
        restricedRowBool(rowBool)=     any(S(rowBool, colBool),2)...
                                    & ~any(S(rowBool,~colBool),2);

        %corresponding reactions involving certain metabolites
        restricedRowBool2=false(mlt,1);
        restricedRowBool2(rowBool) = any(S(rowBool, colBool),2);
        %difference
        restricedRowBool= restricedRowBool2  & ~restricedRowBool;
   otherwise
        error(['Did not recognise mode: ' mode])
end

author: opencobra

repo: cobratoolbox


text:
function [U,N,Y] = exterior_skinning_modes(V,D,H,M,k)
  % EXTERIOR_SKINNING_MODES Perform modal analysis on H,M within the subspace
  % where vertices in "linear blend skinning handle" selected by D of a mesh with
  % vertices V deform affinely.
  %
  % [U,N,Y,H,M] = exterior_skinning_modes(V,D,H,M,k)
  %
  % Inputs:
  %   V  #V by dim
  %   D  #V by #H bone-handle incidence matrix D(i,j) = 1 means vertex i belongs
  %     to bone j
  %   H  #V*dim by #V*dim stiffness matrix
  %   M  #V*dim by #V*dim mass matrix
  %   k  number of eigen modes to compute
  % Outputs:
  %   U  #V*3 by k list of eigen modes in the "maximal" space
  %   N  #V*3 by #N subspace reduction matrix
  %   Y  #N by k list of subspace modes so that U = N * Y
  %
  % Example:
  %   H = arap_hessian(V,F);
  %   M = repdiag(massmatrix(V,F),3);
  %   [~,N,Y] = exterior_skinning_modes(V,D,H,M,20);
  %

  assert(max( sum(D,2) == 1));
  assert(all( any(D,1) ));
  n = size(V,1);
  dim = size(V,2);
  assert(size(H,1) == n*dim);
  assert(size(M,1) == n*dim);
  assert(k <= n);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Build the null-space so that V = N * Y
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  int = find(~any(D,2));
  b = find(any(D,2));
  ni = numel(int);
  m = size(D,2);
  S = sparse(lbs_matrix(V(b,:),D(b,:)));
  NI = speye(dim*ni,dim*ni+size(S,2));
  vec = @(X) X(:);
  b3 = vec(b+(0:dim-1)*n);
  int3 = vec(int+(0:dim-1)*n);
  N = sparse(int3,1:numel(int)*dim,1,n*dim,dim*ni+size(S,2));
  N(b3,:) = [sparse(dim*(n-ni),dim*ni) S];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Test that applying identity transformations to handles produces V
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Astack = repmat([eye(3,3) zeros(3,1)],[1 1 m]);
  %% collect transformations into column
  %A = reshape(permute(Astack,[3 1 2]),m*3*(3+1),1);
  %Y = [reshape(V(int,:),[],1);A];
  %tsurf(F,reshape(N*Y,[],3),'CData',1*any(D,2),fphong)
  %axis equal;
  %view(90,0);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Handle constrained Eigen Decomposition
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % wow, this symmetrization is important...
  smash = @(A,P) 0.5*(P'*A*P + (P'*A*P)');
  [Y,YD] = eigs(-0.5*smash(H,N),0.5*smash(M,N),k,'sm');

  U = N*Y;
end

author: alecjacobson

repo: gptoolbox


text:
function v = lift_catreg(r,t,c,theta,cBound,tBound)
%% Lifting for category registration
%% Heng Yang, July 05, 2021
r           = r(:);
t           = t(:);
c           = c(:);
K           = length(c);
theta       = theta(:);
cBoundSq    = cBound^2;
tBoundSq    = tBound^2;

cr          = kron(c,r);
x           = [r;t];
v1          = [1;x;c;cr;theta;kron(theta,x);kron(theta,cr)];

if tBoundSq < t'*t
    v2      = [1;theta] * 0;
else
    v2      = [1;theta] * sqrt(tBoundSq - t'*t);
end

if cBoundSq < c'*c
    tmp     = [1;r] * 0;
else
    tmp     = [1;r] * sqrt(cBoundSq - c'*c);
end

v           = {v1;v2;tmp};

for k = 1:K
    if c(k) < 0
        tmp = [1;r] * 0;
    else
        tmp = [1;r] * sqrt(c(k));
    end
    v       = [v;{tmp}];
end

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function bf10 = corrbfrep(r,n)
%
% bf10 = corrbfrep(r,n)
%
% Calculates replication correlation Bayes Factor using a uniform prior
% between 0-1 (one sided test) as used by Boekel et al, 2014, Cortex.
%

% Function to be integrated
F = @(rho,r,n) ((1-rho.^2).^(1./2*(n-1))) ./ (1-rho*r).^(n-3./2);

% Bayes factor calculation
bf10 = integral(@(rho) F(rho,r,n), 0, 1);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function prevstate = mcolonops(state)
% overload MATLAB colon operator with multi-colon capability
%
% mcolonops(STATE) where STATE is
%   'colon' 'off' 'matlab' : Set to Matlab original COLON mode.
% 	'mcolon' 'on' 'compatible': mcolon "compatible" mode, similar results
%       to Matlab colon returned when invoking with scalar inputs.
% 	'fullmcolon': full mcolon mode, arrays will be generated by MCOLON
%        engine in all configurations.
%
% INQUIRING: prevstate = mcolonops(...) return the current state before
%                                       the new state is set.
%
% See also: colon, mcolon, mcolsettol
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Last update: 31/December/2010

persistent OVERLOADCOLON_STATE
persistent OVERLOAD_PATH

if isempty(OVERLOADCOLON_STATE)
    OVERLOADCOLON_STATE = 'colon';
end
if nargout>=1 || nargin<1
    prevstate = OVERLOADCOLON_STATE; % Return the previous state
end

if nargin>=1
    if ~isempty(OVERLOAD_PATH)
        try
            rmpath(OVERLOAD_PATH);
        catch %#ok
            warning('mcolon:rmpath', ...
                    'mcolonops can''t remove path %s', OVERLOAD_PATH);
        end
        OVERLOADCOLON_STATE = 'colon';
        OVERLOAD_PATH = '';
    end

    bsxopsroot = fileparts(mfilename('fullpath'));
    switch lower(state)
        case {'c' 'colon' 'off' 'matlab'}
        case {'m' 'mc' 'mcolon' 'on' 'compatible'}
            OVERLOADCOLON_STATE = 'mcolon';
            OVERLOAD_PATH = [bsxopsroot '/overload/compatible'];
            addpath(OVERLOAD_PATH);
        case {'f' 'fullmcolon'}
            OVERLOADCOLON_STATE = 'fullmcolon';
            OVERLOAD_PATH = [bsxopsroot '/overload/fullmcolon'];
            addpath(OVERLOAD_PATH);
        otherwise
            error('mcolonops: unknown state <%s>', state);
    end
end

% Prevent the persistent variable to be cleared
mlock;

end


text:
classdef CarlaEnvironment < matlab.System & matlab.system.mixin.Propagates
    % Untitled Add summary here
    %
    % This template includes the minimum set of functions required
    % to define a System object with discrete state.

    % Public, tunable properties
    properties

    end

    properties(DiscreteState)

    end

    % Pre-computed constants
    properties(Access = private)
        car;
    end

    methods(Access = protected)
        function setupImpl(obj)
            % Perform one-time calculations, such as computing constants
            port = int16(2000);
            client = py.carla.Client('localhost', port);
            client.set_timeout(10.0);
            world = client.get_world();

            % Spawn Vehicle
            blueprint_library = world.get_blueprint_library();
            car_list = py.list(blueprint_library.filter("model3"));
            car_bp = car_list{1};
            spawn_point = py.random.choice(world.get_map().get_spawn_points());
            obj.car = world.spawn_actor(car_bp, spawn_point);
            obj.car.set_autopilot(true);

        end

        function [x_position, x_velocity] = stepImpl(obj)
            %
            pause(0.001);
            x_position = obj.car.get_location().x;
            x_velocity = obj.car.get_velocity().x;
        end

        function [distance, velocity] = isOutputComplexImpl(~)
            distance = false;
            velocity = false;
        end

        function [distance, velocity] = getOutputSizeImpl(~)
            distance = [1 1];
            velocity = [1 1];
        end

        function [distance, velocity] = getOutputDataTypeImpl(~)
            distance = 'double';
            velocity = 'double';
        end

        function [distance, velocity] = isOutputFixedSizeImpl(~)
            distance = true;
            velocity = true;
        end

        function resetImpl(~)
            % Initialize / reset discrete-state properties
        end
    end

    methods(Access= public)
        function delete(obj)
            % Delete the car from the Carla world
            if ~isempty(obj.car)
                obj.car.destroy();
            end
        end
    end
end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
function y = gauss(mu, covar, x)
%GAUSS	Evaluate a Gaussian distribution.
%
%	Description
%
%	Y = GAUSS(MU, COVAR, X) evaluates a multi-variate Gaussian  density
%	in D-dimensions at a set of points given by the rows of the matrix X.
%	The Gaussian density has mean vector MU and covariance matrix COVAR.
%
%	See also
%	GSAMP, DEMGAUSS
%

%	Copyright (c) Ian T Nabney (1996-2001)

[n, d] = size(x);

[j, k] = size(covar);

% Check that the covariance matrix is the correct dimension
if ((j ~= d) | (k ~=d))
  error('Dimension of the covariance matrix and data should match');
end

invcov = inv(covar);
mu = reshape(mu, 1, d);    % Ensure that mu is a row vector

x = x - ones(n, 1)*mu;
fact = sum(((x*invcov).*x), 2);

y = exp(-0.5*fact);

y = y./sqrt((2*pi)^d*det(covar));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%TEST_SmoothBnd
%
% This function plots the output of SmoothBnd for a few sample values of
% the smoothing parameter alphaover the domain [-1,1]
%
% Written by Matthew Kelly
% October 2013
% Cornell University
%

t = linspace(-1,2,1000);
Bnd = [0,1];
alpha = [0.01,0.05,0.2];
N=length(alpha);

figure(102); clf; hold on;
for i=1:N
    x1 = SmoothBnd(t,alpha(i),[0,1]);

    subplot(N,1,i); hold on;
    plot(t,x1,'b-','LineWidth',2)
    c1 = t<Bnd(1); c2 = t>=Bnd(2); c3 = ~c1&~c2;
    plot(t(c1),Bnd(1),'k:')
    plot(t(c2),Bnd(2),'k:')
    plot(t(c3),t(c3),'k:')
    title(['Alpha = ' num2str(alpha(i))],'FontSize',14);
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [translatedModel] = translateDraftReconstruction(model)
% This function translates reaction and metabolite IDs in KBase draft
% reconstructions to VMH nomenclature so the test suite can be performed on
% them. It does not make changes to the reaction content otherwise.
%
% USAGE: [translatedModel] = translateDraftReconstruction(model)
%
% INPUT
% model                     COBRA model structure generated by KBase
%                           reconstruction pipeline
% OUTPUT
% translatedModel:          COBRA model structure with metabolites and
%                           reactions translated to VMH identifiers (but
%                           not ran through refinement pipeline)
%
% AUTHOR
%       - Almut Heinken 09/2020

tol=0.0000001;

% Load metabolite and reaction translation tables
database=loadVMHDatabase;

translateRxns = readtable('ReactionTranslationTable.txt', 'Delimiter', '\t');
translateRxns = [translateRxns.Properties.VariableDescriptions;table2cell(translateRxns)];
translateRxns(find(~strncmp(translateRxns(:,1),'EX_',3)),:)=[];
translateMets = readtable('MetaboliteTranslationTable.txt', 'Delimiter', '\t');
translateMets = [translateMets.Properties.VariableDescriptions;table2cell(translateMets)];

model.mets=strrep(model.mets,'[c0]','[c]');
model.mets=strrep(model.mets,'[e0]','[e]');
model.mets=strrep(model.mets,'_c0','[c]');
model.mets=strrep(model.mets,'_e0','[e]');

for k=1:length(translateMets)
    model.mets = strrep(model.mets,translateMets{k,1}, translateMets{k,2});
end

% special case: two biomass reactions in model (bug in KBase?)
if length(intersect(model.rxns,{'bio1','bio2'}))==2
    model=removeRxns(model,'bio2');
end

% add some reactions that need to be in every reconstruction
essentialRxns={'DM_atp_c_','sink_PGPm1[c]'};

for i=1:length(essentialRxns)
    model = addReaction(model, essentialRxns{i}, 'reactionFormula', database.reactions{find(ismember(database.reactions(:, 1), essentialRxns{i})), 3});
end

% translate while keeping track of original reaction
rxns = strrep(model.rxns, '_c0', '');
rxns = strrep(rxns, '_e0', '');
%     rxns = strrep(rxns, 'R_', '');

for k=1:length(translateRxns)
    rxns = strrep(rxns,translateRxns{k,1}, translateRxns{k,2});
end

% rebuild model to get reactions in order-will get "reactions not
% sorted" error otherwise
[rxnsSorted,oldInd]=sort(rxns);

newModel=struct;
newModel.rxns={};
newModel.rxnNames={};
newModel.genes={};
newModel.grRules={};
newModel.rules={};
newModel.subSystems={};
newModel.lb=zeros(1);
newModel.ub=zeros(1);
newModel.S=[];
newModel.mets={};
newModel.metFormulas={};
newModel.metCharges=[];

for k=1:length(rxnsSorted)
    form=printRxnFormula(model,model.rxns{oldInd(k),1});
    name=model.rxnNames{oldInd(k)};
    lb=model.lb(oldInd(k));
    ub=model.ub(oldInd(k));
    geneRule=model.grRules{oldInd(k)};
    newModel=addReaction(newModel,rxnsSorted{k},'reactionName',name,'reactionFormula',...
        form{1},'lowerBound',lb,'upperBound',ub,'geneRule',geneRule,'subSystem','');
end
model=newModel;

% add metabolite formulas and charges
for k=1:length(model.mets)
    met=strrep(model.mets{k},'[c]','');
    met=strrep(met,'[e]','');
    if ~isempty(find(strcmp(database.metabolites(:,1),met)))
        model.metFormulas{k,1}=database.metabolites{find(strcmp(database.metabolites(:,1),met)),4};
        model.metCharges(k,1)=str2double(database.metabolites{find(strcmp(database.metabolites(:,1),met)),5});
    else
        model.metFormulas{k,1}='X';
        model.metCharges(k,1)=0;
    end
end
model=convertOldStyleModel(model);
model=changeObjective(model,model.rxns{find(strncmp(model.rxns,'bio',3))});

% Delete sink for petidoglycan if not needed
modelTest = removeRxns(model, 'sink_PGPm1[c]');
FBA = optimizeCbModel(modelTest, 'max');
if FBA.f>tol
    model= modelTest;
end

translatedModel=model;

end


text:
function [pathstr,name,ext] = filepartsVector(filename)

pathstr = [];
name = [];
ext = [];
for i=1:length(filename)
    [pathstr{i},name{i},ext{i}] = fileparts(filename{i});
end

author: jianxiongxiao

repo: ProfXkit


text:
function popselected = select(pop, fitness, nselected, method)
% SELECT
% popselected = SELECT(pop, fitness, nselected, method) select the fittest
% individuals to survive to the next generation.
%

popSize = size(pop,1);

switch method

    case 'roulette'
        p=fitness/sum(fitness); % probabilities of select
        cump=cumsum(p);         % cumulative sum of probabilities
        I = interp1([0 cump],1:(popSize+1),rand(1,nselected),'linear');
        % random numbers from 1:nselected according probabilities
        I = floor(I);

    case 'competition'
        % randomly generated two sets of population
        i1 = ceil( popSize*rand(1,nselected) );
        i2 = ceil( popSize*rand(1,nselected) );
        % compare the fitness and select the fitter
        I = i1.*( fitness(i1)>=fitness(i2) ) + ...
            i2.*( fitness(i1)< fitness(i2) );

end

popselected=pop(I,:);

author: HuangCongQing


text:
function [mse, rmse] = RMSE2(signal1, signal2)

originalRowSize = size(signal1,1);
originalColSize = size(signal1,2);

signal1 = signal1(:);
signal2 = signal2(:);

mse = sum((signal1 - signal2).^2)./(originalRowSize*originalColSize);
rmse = sqrt(mse);

end


text:
function [binaryStateSpace,multiStateSpace] = prtRvUtilNChooseKStateSpace(n,k)
%[b,m] = nChooseKStateSpace(n,k)  N-Choose-K state space
%   [b,m] = nChooseKStateSpace(n,k) outputs a boolean matrix of
%   size(nchoosek(n,k),n) where each row is a unique set of k "true" values.
%   Each row of m contains the indices of the corresponding row b which are
%   true.
%
%   Inputs:
%       n - (integer, 1 <= n <= 36)
%       k - (integer, vector, 0 <= k(i) <= n) - if k is a vector,
%           nChooseKStateSpace returns a concatenation of
%           nChooseKStateSpace(n,k(i)) \forall i.  Note, when k is a
%           vector, multiStateSpace is a cell array of matrices of indices
%
%           Note: If k is zero, binaryStateSpace is a row vector of zeros,
%           and multiStateSpace is empty

if k == 0
    binaryStateSpace = zeros(1,n);
    multiStateSpace = [];
    return
end

if numel(k) > 1
    binaryStateSpace = [];
    for i = 1:numel(k)
        [bss,multiStateSpace{i}] = prtRvUtilNChooseKStateSpace(n,k(i));
        binaryStateSpace = cat(1,binaryStateSpace,bss);
    end
    return;
end

if k > n
    error('k (%d) must be <= n (%d)',n,k);
end
if n <= 0 || k < 0
    error('k (%d) must be > 0 and n (%d) must be >= 0',n,k);
end
if nchoosek(n,k) > 1e6
    error('Maximum matrix size exceeded. nchoosek(n,k) must be below 10^6.')
end

multiStateSpace = prtRvUtilMultinomialStateSpace(k,n);
multiStateSpace = sort(multiStateSpace,2);

multiStateSpace = unique(multiStateSpace,'rows');  %remove duplicates, only one of [1,2], [2,1] should be in there
multiStateSpace = multiStateSpace + 1;             % 0-->1, 1-->2, etc.

keepIndices = ~any(diff(multiStateSpace,1,2) == 0,2); %remove rows with same indices [3 3], [1 2 2] etc.
multiStateSpace = multiStateSpace(keepIndices,:);

binaryStateSpace = false(size(multiStateSpace,1),n);
for i = 1:size(binaryStateSpace,1)
    binaryStateSpace(i,multiStateSpace(i,:)) = true;  %make the binaryStateSpace matrix have ones where multiStateSpace says to put em
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function s=sprintcpx(z,f)
%SPRINTCPX  format a complex number for printing S=(Z,F)
%
% Usage: fprintf('%s',sprintcpx(z));
%
%  Inputs: z   a complex number to print
%          f   optional formatting string as in fprintf e.g. '0.2f' [default: 'g']
%              may also include 'i' or 'j' [default] to control sqrt(-1) symbol.
%
% Outputs: s   formatted output string

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

if nargin<2 || ~numel(f)
    f='g';
end
if any(f=='i')
    ij='i';
else
    ij='j';
end
f((f=='i')|(f=='j'))=[]; % remove i and j specifiers
if ~numel(f)
    f='g';
end
if any(f=='+')
    pl='';
else
    pl='+';
end
f=['%' f];
a=real(z);
b=imag(z);
jx=[1 3 2 4 3 4 1 3 2];
ix=jx(3*sign(a)+sign(b)+5);
switch(ix)
    case 1
        s=sprintf([f f ij],a,b);
    case 2
        s=sprintf([f pl f ij],a,b);
    case 3
        s=sprintf(f,a);
    case 4
        s=sprintf([f ij],b);
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%TPOLY Generate scalar polynomial trajectory
%
% [S,SD,SDD] = TPOLY(S0, SF, M) is a scalar trajectory (Mx1) that varies
% smoothly from S0 to SF in M steps using a quintic (5th order) polynomial.
% Velocity and acceleration can be optionally returned as SD (Mx1) and SDD
% (Mx1) respectively.
%
% TPOLY(S0, SF, M) as above but plots S, SD and SDD versus time in a single
% figure.
%
% [S,SD,SDD] = TPOLY(S0, SF, T) as above but the trajectory is computed at
% each point in the time vector T (Mx1).
%
% [S,SD,SDD] = TPOLY(S0, SF, T, QD0, QD1) as above but also specifies the
% initial and final velocity of the trajectory.
%
% Notes::
% - If M is given
%   - Velocity is in units of distance per trajectory step, not per second.
%   - Acceleration is in units of distance per trajectory step squared, not
%     per second squared.
% - If T is given then results are scaled to units of time.
% - The time vector T is assumed to be monotonically increasing, and time
%   scaling is based on the first and last element.
%
% Reference:
%  Robotics, Vision & Control
%  Chap 3
%  Springer 2011
%
% See also LSPB, JTRAJ.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% [S,SD,SDD] = TPOLY(S0, SF, N, SD0, SDF) as above but specifies initial
% and final joint velocity for the trajectory.
%
% [S,SD,SDD] = TPOLY(S0, SF, T, SD0, SDF) as above but specifies initial
% and final joint velocity for the trajectory and time vector T.
%
% Notes::
% - In all cases if no output arguments are specified S, SD, and SDD are plotted
%   against time.
%
% See also LSPB, JTRAJ.

function [s,sd,sdd] = tpoly(q0, qf, t, qd0, qdf)

    t0 = t;
    if isscalar(t)
		t = (0:t-1)';
    else
        t = t(:);
    end
    if nargin < 4
        qd0 = 0;
    end
    if nargin < 5
        qdf = 0;
    end

    plotsargs = {'.-', 'Markersize', 16};

    tf = max(t);
    % solve for the polynomial coefficients using least squares
    X = [
        0           0           0         0       0   1
        tf^5        tf^4        tf^3      tf^2    tf  1
        0           0           0         0       1   0
        5*tf^4      4*tf^3      3*tf^2    2*tf    1   0
        0           0           0         2       0   0
        20*tf^3     12*tf^2     6*tf      2       0   0
    ];
    coeffs = (X \ [q0 qf qd0 qdf 0 0]')';

    % coefficients of derivatives
    coeffs_d = coeffs(1:5) .* (5:-1:1);
    coeffs_dd = coeffs_d(1:4) .* (4:-1:1);

    % evaluate the polynomials
    p = polyval(coeffs, t);
    pd = polyval(coeffs_d, t);
    pdd = polyval(coeffs_dd, t);

    switch nargout
        case 0
            if isscalar(t0)
                % for scalar time steps, axis is labeled 1 .. M
                xt = t+1;
            else
                % for vector time steps, axis is labeled by vector M
                xt = t;
            end

            clf
            subplot(311)
            plot(xt, p, plotsargs{:}); grid; ylabel('$s$', 'FontSize', 16, 'Interpreter','latex');

            subplot(312)
            plot(xt, pd, plotsargs{:}); grid;
            if isscalar(t0)
                ylabel('$ds/dk$', 'FontSize', 16, 'Interpreter','latex');
            else
                ylabel('$ds/dt$', 'FontSize', 16, 'Interpreter','latex');
            end

            subplot(313)
            plot(xt, pdd, plotsargs{:}); grid;
            if isscalar(t0)
                ylabel('$ds^2/dk^2$', 'FontSize', 16, 'Interpreter','latex');
            else
                ylabel('$ds^2/dt^2$', 'FontSize', 16, 'Interpreter','latex');
            end

            if ~isscalar(t0)
                xlabel('t (seconds)')
            else
                xlabel('k (step)');
                for c=findobj(gcf, 'Type', 'axes')
                    set(c, 'XLim', [1 t0]);
                end
            end
            shg
        case 1
            s = p;
        case 2
            s = p;
            sd = pd;
        case 3
            s = p;
            sd = pd;
            sdd = pdd;
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ x, ierror ] = congruence ( a, b, c )

%*****************************************************************************80
%
%% CONGRUENCE solves a congruence of the form A * X = C ( mod B ).
%
%  Discussion:
%
%    A, B and C are given integers.  The equation is solvable if and only
%    if the greatest common divisor of A and B also divides C.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eric Weisstein, editor,
%    CRC Concise Encylopedia of Mathematics,
%    CRC Press, 1998, page 446.
%
%  Parameters:
%
%    Input, integer A, B, C, the coefficients of the Diophantine equation.
%
%    Output, integer X, the solution of the Diophantine equation.
%    X will be between 0 and B-1.
%
%    Output, integer IERROR, error flag.
%    0, no error, X was computed.
%    1, A = B = 0, C is nonzero.
%    2, A = 0, B and C nonzero, but C is not a multiple of B.
%    3, A nonzero, B zero, C nonzero, but C is not a multiple of A.
%    4, A, B, C nonzero, but GCD of A and B does not divide C.
%    5, algorithm ran out of internal space.
%
  nmax = 100;
%
%  Defaults for output parameters.
%
  ierror = 0;
  x = 0;
  y = 0;
%
%  Special cases.
%
  if ( a == 0 & b == 0 & c == 0 )
    x = 0;
    return
  elseif ( a == 0 & b == 0 & c ~= 0 )
    ierror = 1;
    x = 0;
    return
  elseif ( a == 0 & b ~= 0 & c == 0 )
    x = 0;
    return
  elseif ( a == 0 & b ~= 0 & c ~= 0 )
    x = 0;
    if ( mod ( c, b ) ~= 0 )
      ierror = 2;
    end
    return
  elseif ( a ~= 0 & b == 0 & c == 0 )
    x = 0;
    return
  elseif ( a ~= 0 & b == 0 & c ~= 0 )
    x = c / a;
    if ( mod ( c, a ) ~= 0 )
      ierror = 3;
    end
    return
  elseif ( a ~= 0 & b ~= 0 & c == 0 )
%   g = i4_gcd ( a, b );
%   x = b / g;
    x = 0;
    return
  end
%
%  Now handle the "general" case: A, B and C are nonzero.
%
%  Step 1: Compute the GCD of A and B, which must also divide C.
%
  g = i4_gcd ( a, b );

  if ( mod ( c, g ) ~= 0 )
    ierror = 4;
    return
  end

  a_copy = a / g;
  b_copy = b / g;
  c_copy = c / g;
%
%  Step 2: Split A and B into sign and magnitude.
%
  a_mag = abs ( a_copy );
  a_sign = i4_sign ( a_copy );
  b_mag = abs ( b_copy );
  b_sign = i4_sign ( b_copy );
%
%  Another special case, A_MAG = 1 or B_MAG = 1.
%
  if ( a_mag == 1 )
    x = a_sign * c_copy;
    return
  elseif ( b_mag == 1 )
    x = 0;
    return
  end
%
%  Step 3: Produce the Euclidean remainder sequence.
%
  if ( b_mag <= a_mag )

    swap = 0;
    q(1) = a_mag;
    q(2) = b_mag;

  else

    swap = 1;
    q(1) = b_mag;
    q(2) = a_mag;

  end

  n = 3;

  while ( 1 )

    q(n) = mod ( q(n-2), q(n-1) );

    if ( q(n) == 1 )
      break;
    end

    n = n + 1;

    if ( nmax < n )
      ierror = 5;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CONGRUENCE - Fatal error!\n' );
      fprintf ( 1, '  Exceeded number of iterations.\n' );
      error ( 'CONGRUENCE - Fatal error!' );
    end

  end
%
%  Step 4: Now go backwards to solve X * A_MAG + Y * B_MAG = 1.
%
  y = 0;
  for k = n : -1 : 2
    x = y;
    y = ( 1 - x * q(k-1) ) / q(k);
  end
%
%  Step 5: Undo the swapping.
%
  if ( swap )
    z = x;
    x = y;
    y = z;
  end
%
%  Step 6: Now apply signs to X and Y so that X * A + Y * B = 1.
%
  x = x * a_sign;
%
%  Step 7: Multiply by C, so that X * A + Y * B = C.
%
  x = x * c_copy;
%
%  Step 8: Now force 0 <= X < B.
%
  x = mod ( x, b );
%
%  Step 9: Force positivity.
%
  if ( x < 0 )
    x = x + b;
  end

  return
end


text:
function installDevTools(repoName)
% checks if the MATLAB.devTools are already present and installs the
% MATLAB.devTools next to the cloned folder cobratoolbox (if not present yet)
%
% USAGE:
%    installDevTools(repoName)
%
% OPTIONAL INPUT:
%    repoName:     Name of the repository (default: 'opencobra/cobratoolbox')
%
% NOTE:
%    In order to install the devTools, please make sure that the SSH key is
%    configured properly as explained `here <https://opencobra.github.io/MATLAB.devTools/stable/installation.html>`__.
%    The full documentation of the MATLAB.devTools is `here <https://opencobra.github.io/MATLAB.devTools>`__.
%
% .. Author: Laurent Heirendt, June 2017

    global ENV_VARS
    global CBTDIR

    % make sure that The COBRA Toolbox is properly installed;
    % git and curl are tested here
    if isempty(ENV_VARS)
        ENV_VARS.printLevel = false;
        initCobraToolbox(false); %Don't update the toolbox automatically
        ENV_VARS.printLevel = true;
    end

    currentDir = pwd;

    % change to the local directory
    localDir = [CBTDIR filesep '..'];
    cd(localDir);

    % set the default repoName
    if ~exist('repoName', 'var')
        repoName = 'opencobra/cobratoolbox';
    end

    % add the public key from github.com to the known hosts
    addKeyToKnownHosts();

    installFlag = false;

    % check if the devTools already exist
    if exist([CBTDIR filesep '..' filesep 'MATLAB.devTools'], 'dir') == 7

        reply = input([' > There is already a copy of the MATLAB.devTools installed in the default location.\n   Do you want to enter another location to install the MATLAB.devTools? Y/N [Y]: '], 's');
        localDir = '';

        % enter another location for the MATLAB.devTools
        if isempty(reply) || strcmpi(reply, 'y') || strcmpi(reply, 'yes')
            dirReply = '';
            while isempty(dirReply)
                dirReply = input(['\n -> Please define the installation path of the MATLAB.devTools\n    Enter the path: '], 's');

                if exist([dirReply filesep 'MATLAB.devTools'], 'dir') == 7
                    fprintf(['\n   The directory ', [dirReply filesep 'MATLAB.devTools'], ' already exists. Please enter a different directory.\n']);
                    dirReply = '';
                end

                % if the entered directory is not empty, try to add it to the path
                if ~isempty(dirReply)
                    addpath(dirReply);
                    if exist(dirReply, 'dir') == 7
                        localDir = dirReply;
                        break;
                    else
                        fprintf('\n   -> The directory does not exist. Please try again.\n');
                        dirReply = '';
                    end
                end
            end

            installFlag = true;
        else
            installFlag = false;
            fprintf(' > You cannot install the MATLAB.devTools. Please delete the folder %s\n', [CBTDIR filesep '..' filesep 'MATLAB.devTools']);
        end
    else
        installFlag = true;
    end

    % install a new copy of the MATLAB.devTools
    if installFlag

        % changing to the local directory
        cd(localDir);

        % install the devTools properly speaking
        fprintf(['\n > Installing the MATLAB.devTools into ' localDir ' (might take some time) ... ']);

        [status_gitClone, result_gitClone] = system(['git clone -c http.sslVerify=false --depth=1 git@github.com:opencobra/MATLAB.devTools.git']);

        if status_gitClone == 0

            if ~isempty(localDir)
                devToolsDir = localDir;
            end

            fprintf(['Done.\n > Location of the MATLAB.devTools: ', strrep(devToolsDir, '\', '\\'), '\n\n']);

            % add the path to the devTools
            addpath(genpath(devToolsDir));

            % removing .git folder
            rmpath(genpath([devToolsDir filesep '.git']));

            fprintf([' > In order to start contributing, you can now run:\n\n >> contribute(''' repoName '''); \n\n']);
        else
            fprintf(result_gitClone);
            error(['The MATLAB.devTools could not be cloned. Please make sure that your SSH key is configured correctly.']);
        end
    end

    % change back to the current directory
    cd(currentDir);
end


text:
disp('Test: Mutual information between two images')
load mri
A=D(:,:,8);
B=D(:,:,9);
mi(A,A),mi(A,B)

disp('Test: Mutual information between two signals')
load garchdata
nasdaq = price2ret(NASDAQ);
nyse = price2ret(NYSE);
mi(nasdaq,nasdaq), mi(nasdaq,nyse)

author: Sable


text:
%  Figure 7.16      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%
% script to generate Fig. 7.16
% fig7_16.m
%
clf;
f=[0 1;-1 0];
g=[0;1];
h=[1 0];
K=[3 4];
fc=f-g*K;
gc=[4*g [1;0]];
hc=[h;[0 1]; -K/4];
jc=[0 0;0 0; 1 0];
t=0:.1:7;
sys=ss(fc,gc,hc,jc);
[y,t]=step(sys,t);
plot(t,y(:,:,1));
xlabel('Time (sec)');
ylabel('Amplitude');
text(1.5,.9,'x_1');
text(1.6,.3,'x_2');
text(1.5,.05,'u/4');
text(6.5,.3,'u_{ss}');
text(6.5,.95,'x_{ss}');
title('Fig. 7.16 Step response of the oscillator to a reference input');
nicegrid;


text:
function [s, lam] = shcovft(x, shrinkvar)
% Shrinkage estimate of a covariance matrix, using optimal shrinkage coefficient.
% INPUT:
% x is n*p data matrix
% shrinkvar :
% 0: corshrink (default)
% 1: varshrink
%
% OUTPUT:
% s is the posdef p*p cov matrix
% lam is the shrinkage coefficient
%
% See  J. Schaefer and K. Strimmer.  2005.  A shrinkage approach to
%   large-scale covariance matrix estimation and implications
%   for functional genomics. Statist. Appl. Genet. Mol. Biol. 4:32.
% This code is based on their original code http://strimmerlab.org/software.html
% but has been vectorized and simplified by Kevin Murphy.
% Adapted by Marcel van Gerven

if nargin < 2, shrinkvar = 0; end

[n p] = size(x);
if p==1, s=var(x); return; end

switch num2str(shrinkvar)

  case '1' % Eqns 10 and 11 of Opgen-Rhein and Strimmer (2007)
    [v, lam] = varshrink(x);
    dsv = diag(sqrt(v));
    r = corshrink(x);
    s = dsv*r*dsv;

  otherwise % case 'D' of Schafer and Strimmer
    v = var(x);
    dsv = diag(sqrt(v));
    [r, lam] = corshrink(x);
    s = dsv*r*dsv;

end

%%%%%%%%

function [sv, lambda] = varshrink (x)
% Eqns 10 and 11 of Opgen-Rhein and Strimmer (2007)
[v, vv] = varcov(x);
v = diag(v); vv = diag(vv);
vtarget = median(v);
numerator = sum(vv);
denominator = sum((v-vtarget).^2);
lambda = numerator/denominator;
lambda = min(lambda, 1); lambda = max(lambda, 0);
sv = (1-lambda)*v + lambda*vtarget;

function [Rhat, lambda] = corshrink(x)
% Eqn on p4 of Schafer and Strimmer 2005
[n, p] = size(x);
sx = makeMeanZero(x); sx = makeStdOne(sx); % convert S to R
[r, vr] = varcov(sx);
offdiagsumrij2 = sum(sum(tril(r,-1).^2));
offdiagsumvrij = sum(sum(tril(vr,-1)));
lambda = offdiagsumvrij/offdiagsumrij2;
lambda = min(lambda, 1); lambda = max(lambda, 0);
Rhat = (1-lambda)*r;
Rhat(logical(eye(p))) = 1;

function [S, VS] = varcov(x)
% s(i,j) = cov X(i,j)
% vs(i,j) = est var s(i,j)
[n,p] = size(x);
xc = makeMeanZero(x);
S = cov(xc);
XC1 = repmat(reshape(xc', [p 1 n]), [1 p 1]); % size p*p*n !
XC2 = repmat(reshape(xc', [1 p n]),  [p 1 1]); % size p*p*n !
VS = var(XC1 .* XC2, 0,  3) * n/((n-1)^2);

function xc = makeMeanZero(x)
% make column means zero
[n,p] = size(x);
m = mean(x);
xc = x - ones(n, 1)*m;

function xc = makeStdOne(x)
% make column  variances one
[n,p] = size(x);
sd = ones(n, 1)*std(x);
xc = x ./ sd;

author: alexandrebarachant

repo: covariancetoolbox


text:
function result = glmnetCoef( object, s )

%--------------------------------------------------------------------------
% glmnetCoef.m: print coefficients from a "glmnet" object
%--------------------------------------------------------------------------
%
% USAGE:
%    glmnetCoef(fit);
%    glmnetCoef(fit, s);
%
% DETAILS:
%    glmnetCoef(fit, s) is equivalent to glmnetPredict(fit, "coefficients", [], s)
%    See glmnetPredict for more details.
%
% LICENSE: GPL-2
%
% DATE: 14 Jul 2009
%
% AUTHORS:
%    Algorithm was designed by Jerome Friedman, Trevor Hastie and Rob Tibshirani
%    Fortran code was written by Jerome Friedman
%    R wrapper (from which the MATLAB wrapper was adapted) was written by Trevor Hasite
%    MATLAB wrapper was written and maintained by Hui Jiang, jiangh@stanford.edu
%    Department of Statistics, Stanford University, Stanford, California, USA.
%
% REFERENCES:
%    Friedman, J., Hastie, T. and Tibshirani, R. (2009)
%    Regularization Paths for Generalized Linear Models via Coordinate Descent.
%    Journal of Statistical Software, 33(1), 2010
%
% SEE ALSO:
%    glmnet, glmnetSet, glmnetPrint, glmnetPredict and glmnetPlot methods.
%
% EXAMPLES:
%    x=randn(100,20);
%    y=randn(100,1);
%    fit1=glmnet(x,y);
%    glmnetCoef(fit1,0.01) % extract coefficients at a single value of lambda
%
% DEVELOPMENT: 14 Jul 2009: Original version of glmnet.m written.

%

if nargin < 2
    s = object.lambda;
end

result = glmnetPredict(object, 'coefficients', [], s);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = randomAngle3d(varargin)
%RANDOMANGLE3D Return a 3D angle uniformly distributed on unit sphere.
%
%   usage
%   [THETA PHI] = randomAngle3d
%   Generate an angle unformly distributed on the surface of the unit
%   sphere.
%
%   "Mathematical" convention is used: theta is the colatitude (angle with
%   vertical axis, 0 for north pole, +pi for south pole, pi/2 for points at
%   equator) with z=0.
%   phi is the same as matlab cart2sph: angle from Ox axis, counted
%   positively counter-clockwise.
%
%   [THETA PHI] = randomAngle3d(N)
%   generates N random angles (N is a scalar). The result is a N-by-2
%   array.
%
%   Example:
%     % Draw some points on the surface of a sphere
%     figure;
%     drawSphere; hold on;
%     drawPoint3d(pts, '.');
%     axis equal;
%
%   See also
%   angles3d, sph2cart2, cart2sph2

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2005-02-18
% Copyright 2005-2022 INRA - Cepia Software platform

N = 1;
if ~isempty(varargin)
    N = varargin{1};
end

phi = 2*pi*rand(N, 1);
theta = asin(2*rand(N, 1)-1) + pi/2;

if nargout<2
    var = [theta phi];
    varargout{1} = var;
else
    varargout{1} = theta;
    varargout{2} = phi;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
close all
clear all
path(path,'..\..\..\FUZZCLUST')
%the data
load motorcycle.txt
data.X = motorcycle(:,[1 2]);

[N,n]=size(data.X);

%data normalization
data = clust_normalize(data,'range');
plot(data.X(:,1),data.X(:,2),'.')
hold on
%parameters
param.c=4;
param.vis=1;
param.val=2;
%clustering
result=kmedoid(data,param);
%validation
result = validity(result,data,param);
result.validity

repo: mcbench-benchmarks


text:
function [handles,msg,lattice,msize,dim,normalization,comps]=vis_som_show_data(p,f)

% VIS_SOM_SHOW_DATA Checks and returns UserData and subplot handles stored by SOM_SHOW
%
% [handles,msg,lattice,msize,dim,normalization,comps] = vis_som_show_data(p, f)
%
%  Input and output arguments ([]'s are optional):
%   [p]           (vector) subplot numbers
%                 (string) 'all' to process all subplots, this is default
%                          'comp' to process only subplots which have
%                          component planes
%   [f]           (double) figure handle, default is current figure
%
%   handles       (vector) handles of requested subplots
%   msg           (string) error message or empty string (no error)
%   lattice       (string) map lattice: 'hexa' or 'rect'
%   msize         (vector) map grid size in figure
%   dim           (scalar) map data dimension in figure
%   normalization (struct) normalization struct used in the map in figure
%   comps         (vector) the component indexes in figure
%
% This function gets the handles of component planes and u-matrices in
% subplots p from figure f. SOM_SHOW writes the handles into the
% UserData field of the figure where their order won't be mixed
% up. This function reads the data according to the vector p. If the
% figure has been manipulated (original planes are missing) the function
% warns user or returns error string.
%
% The main purpose for this is to be a subfuncion for SOM_SHOW_ADD,
% SOM_SHOW_CLEAR and SOM_RECOLORBAR functions, but it may be used on
% command line in the followong manner:
%
%  % plots text on the fifth plane
%  axes(vis_som_show_data(5)); hold on; text(1,3,'I am here');
%
% See also SOM_SHOW, SOM_SHOW_ADD.

% Copyright (c) 1997-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta Johan 201099 juuso 160600

%% Check input args %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

error(nargchk(0, 2, nargin))  % check no. of input args

%% Init %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

handles=[];                                     % initialize output
normalize=[];
comps=[];
dim=[];
msize=[];
lattice=[];
msg=[];

cr=sprintf('\n');                               % carriage return

if nargin < 2 || isempty(f)
  f=gcf;                                        % default figure
end

if nargin < 1 || isempty(p)                      % default p
  p= 'all';
end

% Find component planes and u-matrices from the figure and get the
% UserData field where the handles for the components are
% in the original order.
% If the fields are corrupted, return an error message.

h_real = [findobj(f, 'Tag', 'Cplane'); ...
    findobj(f, 'Tag', 'Uplane'); ...
    findobj(f,'Tag','CplaneI'); ...
    findobj(f,'Tag','UplaneI')];
eval( 'h_stored=getfield(get(f,''UserData''),''subplotorder'');' , ...
    'msg=[ msg cr '' Missing SOM_SHOW.subplotorder''];');
eval( 'normalization=getfield(get(f,''UserData''),''comp_norm'');' , ...
    'msg=[msg cr '' Missing SOM_SHOW.comp_norm''];');
eval( 'comps=getfield(get(f,''UserData''),''comps'');' , ...
    'msg=[msg cr '' Missing SOM_SHOW.comps''];');
eval( 'msize=getfield(get(f,''UserData''),''msize'');' , ...
    'msg=[msg cr '' Missing SOM_SHOW.msize''];');
eval( 'dim=getfield(get(f,''UserData''),''dim'');' , ...
    'msg=[msg cr '' Missing SOM_SHOW.dim''];');
eval( 'lattice=getfield(get(f,''UserData''),''lattice'');' , ...
    'msg=[msg cr '' Missing SOM_SHOW.lattice''];');
if ~isempty(msg),
  msg=['The figure does not contain SOM_SHOW visualization or is corrupted.'...
	cr msg cr cr ...
	'This command may be applied only to a SOM_SHOW visualization.'];
  return;
end

%% Check arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

index=ismember(h_stored, h_real);  % the original order for plot axes

if ~prod(double(index))                    % missing planes?!
                                           % double added by kr 1.10.02
  l1= 'Some of the original planes seems to be missing.';
  l2= 'Subplot numbers now refer to the existing ones.';
  warning([l1 cr l2]);
end

if ~prod(double(ismember(h_real, h_stored))) % extra planes?!
                                             % double added by kr 5.9.02
  warning('There seems to be new planes. Subplot numbers refer to the old ones.');
end

h_stored=h_stored(index);          % existing original plots in original order

if ischar(p)                       % check if p is 'all'
  switch(p)
   case 'all'
    p=1:size(h_stored,1);          % all original subplots
   case 'comp'
    p=find(comps>0);
   otherwise
    msg= 'String value for subplot number vector has to be ''all''!';
    return;
  end
end

if ~vis_valuetype(p,{ '1xn','nx1'}) % check the size
  msg= 'Subplot numbers (argument p in help text) have to be in a vector!';
  return
end

if min(p) < 1                      % check for invalid values
  msg= 'Subplot numbers (argument p in help text) must be at least 1!';
  return
end

%% p is too large

if max(p) > size(h_stored,1)
  l1= 'There are not so many existing subplots created by SOM_SHOW in the';
  l2= 'figure as you are trying to refer with subplot numbers.';
  l3= 'This is probably caused by having a too large subplot number.';
  l4= 'However, the reason may be invalid manipulation of';
  l5= 'this figure object or a program failure, too.';
  msg=([l1 cr l2 cr l3 cr cr l4 cr l5]);
  return;
end

%% Action and building output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

handles=h_stored(p);
comps=comps(p);

author: ilarinieminen

repo: SOM-Toolbox


text:
function a = svm(hyper)
%=============================================================================
% SVM Support Vector Machine object
%=============================================================================
% a=svm(hyperParam)
%
% Generates a svm object with given hyperparameters.
%
%
%   Hyperparameters (with defaults)
%   child=kernel         -- the kernel is stored as a member called "child"
%   C=Inf                -- the soft margin C parameter
%   ridge=1e-13          -- a ridge on the kernel
%   balanced_ridge=0     -- for unbalanced data
%   nu = 0               -- Schoelkopf's nu svm parameter
%   optimizer='default'  -- other choices={andre,quadprog,svmlight,
%                                          libsvm,svmtorch(linux only)}
%                           For "libsvm" optimizer you can specify the used cache size
%                           by the global variable "libsvm_cachesize".
%   alpha_cutoff=-1;     -- keep alphas with abs(a_i)>alpha_cutoff
%                           default keeps all alphas, another
%                           reasonable choice is e.g alpha_cutoff=1e-5 to remove
%                           zero alphas (i.e non-SVs) to speed up computations.
%
%   Model
%    alpha               -- the weights
%    b0                  -- the threshold
%    Xsv                 -- the Support Vectors
%
% Methods:
%  train, test, get_w
%
% Example:
%
%  d=gen(spiral({'m=200','n=2','noise=0.35'}));
%  [r,a]=train(cv(svm({kernel('rbf',1),'optimizer="andre"'})),d)
%  plot(a{1})
%
%=============================================================================
% Reference : A Tutorial on Support Vector Machines for Pattern Recognition
% Author    : Christopher J. C. Burges
% Link      : http://citeseer.ist.psu.edu/burges98tutorial.html
%=============================================================================

  %<<------hyperparam initialisation------------->>
  a.child=kernel;
  a.C=Inf;
  a.ridge=1e-13;
  a.balanced_ridge=0;
  a.nu = 0;
  a.optimizer='default';
  a.alpha_cutoff=-1;

  % <<-------------model----------------->>
  a.alpha=[];
  a.b0=0;
  a.Xsv=[];
  a.nob=0;

  algoType=algorithm('svm');
  a= class(a,'svm',algoType);

  a.algorithm.alias={'kern','child'}; % kernel aliases

 if nargin==1,
    eval_hyper;
 end;


text:
function y = normpdfg(x,scale,shape,center,max)
% y = normpdfg(x,scale,shape,center,max)
%
% Evalute the generalized gaussian pdf at point(s) x
% the pdf is defined by the given scale (2*variance) and shape and mean
% (center)
%
% Note the variance of the distribution is scale/2
%
% As shape --> Inf we get a (step function) uniform distribution over [center-scale, center+scale]
% As shape --> 0 we get a delta function
% Beta = 1 is the Laplace
% Beta = 2 is gaussian
%
% (C) Tim Mullen, 2011. SCCN/INC UCSD

A = shape./(2.*scale.*gamma(1/shape));
C = max./A;

y = C.*A.*exp(-((abs(x-center)/scale).^shape));

repo: REST


text:
function wval = lagrange_factor ( npol, xpol, nval, xval )

%*****************************************************************************80
%
%% LAGRANGE_FACTOR evaluates the polynomial Lagrange factor at a point.
%
%  Formula:
%
%    W(X) = Product ( 1 <= I <= NPOL ) ( X - XPOL(I) )
%
%  Discussion:
%
%    Suppose F(X) is at least N times continuously differentiable in the
%    interval [A,B].  Pick NPOL distinct points XPOL(I) in [A,B] and compute
%    the interpolating polynomial P(X) of order NPOL ( and degree NPOL-1)
%    which passes through all the points ( XPOL(I), F(XPOL(I)) ).
%    Then in the interval [A,B], the maximum error
%
%      abs ( F(X) - P(X) )
%
%    is bounded by:
%
%      C * FNMAX * W(X)
%
%    where
%
%      C is a constant,
%      FNMAX is the maximum value of the NPOL-th derivative of F in [A,B],
%      W(X) is the Lagrange factor.
%
%    Thus, the value of W(X) is useful as part of an estimated bound
%    for the interpolation error.
%
%    Note that the Chebyshev abscissas have the property that they minimize
%    the value of W(X) over the interval [A,B].  Hence, if the abscissas may
%    be chosen arbitrarily, the Chebyshev abscissas have this advantage over
%    other choices.
%
%    For a set of points XPOL(I), 1 <= I <= NPOL, the IPOL-th Lagrange basis
%    polynomial L(IPOL)(X), has the property:
%
%      L(IPOL)( XPOL(J) ) = delta ( IPOL, J )
%
%    and may be expressed as:
%
%      L(IPOL)(X) = W(X) / ( ( X - XPOL(IPOL) ) * W'(XPOL(IPOL)) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NPOL, the number of abscissas.
%    NPOL must be at least 1.
%
%    Input, real XPOL(NPOL), the abscissas, which should
%    be distinct.
%
%    Input, integer NVAL, the number of evaluation points.
%
%    Input, real XVA(NVAL), the points at which the Lagrange
%    factor is to be evaluated.
%
%    Output, real WVAL(NVAL), the value of the Lagrange factor at XVAL.
%
  wval = zeros ( nval, 1 );

  for ival = 1 : nval

    wval(ival) = prod ( xval(ival) - xpol(1:npol) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m