text:
function test_suite = test_dim_find()
% tests for cosmo_dim_find
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_dim_find_basics()
    prefixes='sf';
    for dim=1:2
        if dim==1
            transposer=@transpose;
        else
            transposer=@(x)x;
        end

        prefix=prefixes(dim);
        attr_name=[prefix 'a'];
        dim_name=[prefix 'dim'];

        dim_labels={'i','j','k'};

        ds=cosmo_synthetic_dataset();
        if dim==1
            ds=cosmo_dim_transpose(ds,{'i','j','k'},1);
        end

        combis=cell2mat(cosmo_cartprod(repmat({[false,true]},3,1)));

        for j=1:size(combis,1);
            msk=combis(j,:);

            use_string_name=sum(msk)==1 && msk(2);

            if use_string_name
                name=dim_labels{msk};
            else
                name=dim_labels(msk);
            end

            [d,i,an,dn,vs]=cosmo_dim_find(ds,name,false);

            if ~any(msk)
                assertEqual({d,i,an,dn,vs},{[],[],[],[],[]});

            else
                assertEqual(d,dim);
                assertEqual(i,transposer(find(msk)'));
                assertEqual(an,attr_name);
                assertEqual(dn,dim_name);

                values=ds.a.(dim_name).values;

                if use_string_name
                    assertEqual(vs,values{msk});
                else
                    assertEqual(vs,values(transposer(msk)));
                end
            end
        end
    end

function test_dim_find_exceptions()
    aet=@(varargin)assertExceptionThrown(@()...
                cosmo_dim_find(varargin{:}),'');

    ds=cosmo_synthetic_dataset();

    ds=cosmo_dim_transpose(ds,'i',1);
    aet(ds,{'i','j'});
    assertEqual(cosmo_dim_find(ds,{'i','j'},false),[]);
    aet(ds,{'foo'});
    assertEqual(cosmo_dim_find(ds,{'foo'},false),[]);
    aet(ds,'foo');
    aet(ds,[]);
    assertEqual(cosmo_dim_find(ds,'foo',false),[]);
    ds=cosmo_dim_transpose(ds,'i',2);
    cosmo_dim_find(ds,'i');
    ds.a.fdim.labels{1}='k';
    aet(ds,'k',true);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function C = nancorr( X,Y )
% NANCORR calculates the sample correlation coefficient
%    for the series with NaNs expected.
%    X is the one series, Y is another.
X=X(:);
Y=Y(:);
L1=length(X);
L2=length(Y);

if L1 ~= L2
    error('The samples must be of the same length')
end

for i=1:L1,
    if isnan(X(i)),
        Y(i)=nan;
    end
    if isnan(Y(i)),
        X(i)=nan;
    end
end

Xm=nanmean(X);
Ym=nanmean(Y);
C=nansum((X-Xm).*(Y-Ym))/sqrt((nansum((X-Xm).^2))*(nansum((Y-Ym).^2)));
end

author: Sable

repo: mcbench-benchmarks


text:
function net = cnn_imagenet_init(varargin)
% CNN_IMAGENET_INIT  Initialize a standard CNN for ImageNet

opts.scale = 1 ;
opts.initBias = 0.1 ;
opts.weightDecay = 1 ;
%opts.weightInitMethod = 'xavierimproved' ;
opts.weightInitMethod = 'gaussian' ;
opts.model = 'alexnet' ;
opts.batchNormalization = false ;
opts.networkType = 'simplenn' ;
opts.cudnnWorkspaceLimit = 1024*1024*1204 ; % 1GB
opts.classNames = {} ;
opts.classDescriptions = {} ;
opts.averageImage = zeros(3,1) ;
opts.colorDeviation = zeros(3) ;
opts = vl_argparse(opts, varargin) ;

% Define layers
switch opts.model
  case 'alexnet'
    net.meta.normalization.imageSize = [227, 227, 3] ;
    net = alexnet(net, opts) ;
    bs = 256 ;
  case 'vgg-f'
    net.meta.normalization.imageSize = [224, 224, 3] ;
    net = vgg_f(net, opts) ;
    bs = 256 ;
  case 'vgg-m'
    net.meta.normalization.imageSize = [224, 224, 3] ;
    net = vgg_m(net, opts) ;
    bs = 196 ;
  case 'vgg-s'
    net.meta.normalization.imageSize = [224, 224, 3] ;
    net = vgg_s(net, opts) ;
    bs = 128 ;
  case 'vgg-vd-16'
    net.meta.normalization.imageSize = [224, 224, 3] ;
    net = vgg_vd(net, opts) ;
    bs = 32 ;
  case 'vgg-vd-19'
    net.meta.normalization.imageSize = [224, 224, 3] ;
    net = vgg_vd(net, opts) ;
    bs = 24 ;
  otherwise
    error('Unknown model ''%s''', opts.model) ;
end

% final touches
switch lower(opts.weightInitMethod)
  case {'xavier', 'xavierimproved'}
    net.layers{end}.weights{1} = net.layers{end}.weights{1} / 10 ;
end
net.layers{end+1} = struct('type', 'softmaxloss', 'name', 'loss') ;

% Meta parameters
net.meta.inputSize = [net.meta.normalization.imageSize, 32] ;
net.meta.normalization.cropSize = net.meta.normalization.imageSize(1) / 256 ;
net.meta.normalization.averageImage = opts.averageImage ;
net.meta.classes.name = opts.classNames ;
net.meta.classes.description = opts.classDescriptions;
net.meta.augmentation.jitterLocation = true ;
net.meta.augmentation.jitterFlip = true ;
net.meta.augmentation.jitterBrightness = double(0.1 * opts.colorDeviation) ;
net.meta.augmentation.jitterAspect = [3/4, 4/3] ;
%net.meta.augmentation.jitterContrast = 0.4 ;
%net.meta.augmentation.jitterSaturation = 0.4 ;
%net.meta.augmentation.jitterScale = [0.9, 1] ;

if ~opts.batchNormalization
  lr = logspace(-2, -4, 60) ;
else
  lr = logspace(-1, -4, 20) ;
end

net.meta.trainOpts.learningRate = lr ;
net.meta.trainOpts.numEpochs = numel(lr) ;
net.meta.trainOpts.batchSize = bs ;
net.meta.trainOpts.weightDecay = 0.0005 ;

% Fill in default values
net = vl_simplenn_tidy(net) ;

% Switch to DagNN if requested
switch lower(opts.networkType)
  case 'simplenn'
    % done
  case 'dagnn'
    net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;
    net.addLayer('top1err', dagnn.Loss('loss', 'classerror'), ...
                 {'prediction','label'}, 'top1err') ;
    net.addLayer('top5err', dagnn.Loss('loss', 'topkerror', ...
                                       'opts', {'topK',5}), ...
                 {'prediction','label'}, 'top5err') ;
  otherwise
    assert(false) ;
end

% --------------------------------------------------------------------
function net = add_block(net, opts, id, h, w, in, out, stride, pad)
% --------------------------------------------------------------------
info = vl_simplenn_display(net) ;
fc = (h == info.dataSize(1,end) && w == info.dataSize(2,end)) ;
if fc
  name = 'fc' ;
else
  name = 'conv' ;
end
convOpts = {'CudnnWorkspaceLimit', opts.cudnnWorkspaceLimit} ;
net.layers{end+1} = struct('type', 'conv', 'name', sprintf('%s%s', name, id), ...
                           'weights', {{init_weight(opts, h, w, in, out, 'single'), ...
                             ones(out, 1, 'single')*opts.initBias}}, ...
                           'stride', stride, ...
                           'pad', pad, ...
                           'learningRate', [1 2], ...
                           'weightDecay', [opts.weightDecay 0], ...
                           'opts', {convOpts}) ;
if opts.batchNormalization
  net.layers{end+1} = struct('type', 'bnorm', 'name', sprintf('bn%s',id), ...
                             'weights', {{ones(out, 1, 'single'), zeros(out, 1, 'single'), ...
                               zeros(out, 2, 'single')}}, ...
                             'learningRate', [2 1 0.3], ...
                             'weightDecay', [0 0]) ;
end
net.layers{end+1} = struct('type', 'relu', 'name', sprintf('relu%s',id)) ;

% -------------------------------------------------------------------------
function weights = init_weight(opts, h, w, in, out, type)
% -------------------------------------------------------------------------
% See K. He, X. Zhang, S. Ren, and J. Sun. Delving deep into
% rectifiers: Surpassing human-level performance on imagenet
% classification. CoRR, (arXiv:1502.01852v1), 2015.

switch lower(opts.weightInitMethod)
  case 'gaussian'
    sc = 0.01/opts.scale ;
    weights = randn(h, w, in, out, type)*sc;
  case 'xavier'
    sc = sqrt(3/(h*w*in)) ;
    weights = (rand(h, w, in, out, type)*2 - 1)*sc ;
  case 'xavierimproved'
    sc = sqrt(2/(h*w*out)) ;
    weights = randn(h, w, in, out, type)*sc ;
  otherwise
    error('Unknown weight initialization method''%s''', opts.weightInitMethod) ;
end

% --------------------------------------------------------------------
function net = add_norm(net, opts, id)
% --------------------------------------------------------------------
if ~opts.batchNormalization
  net.layers{end+1} = struct('type', 'normalize', ...
                             'name', sprintf('norm%s', id), ...
                             'param', [5 1 0.0001/5 0.75]) ;
end

% --------------------------------------------------------------------
function net = add_dropout(net, opts, id)
% --------------------------------------------------------------------
if ~opts.batchNormalization
  net.layers{end+1} = struct('type', 'dropout', ...
                             'name', sprintf('dropout%s', id), ...
                             'rate', 0.5) ;
end

% --------------------------------------------------------------------
function net = alexnet(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;

net = add_block(net, opts, '1', 11, 11, 3, 96, 4, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '2', 5, 5, 48, 256, 1, 2) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '3', 3, 3, 256, 384, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 192, 384, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 192, 256, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 6, 6, 256, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_s(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1', 7, 7, 3, 96, 2, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 3, ...
                           'pad', [0 2 0 2]) ;

net = add_block(net, opts, '2', 5, 5, 96, 256, 1, 0) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '3', 3, 3, 256, 512, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 512, 512, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 3, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '6', 6, 6, 512, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_m(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1', 7, 7, 3, 96, 2, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '2', 5, 5, 96, 256, 2, 1) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '3', 3, 3, 256, 512, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 512, 512, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 6, 6, 512, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_f(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1', 11, 11, 3, 64, 4, 0) ;
net = add_norm(net, opts, '1') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', [0 1 0 1]) ;

net = add_block(net, opts, '2', 5, 5, 64, 256, 1, 2) ;
net = add_norm(net, opts, '2') ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '3', 3, 3, 256, 256, 1, 1) ;
net = add_block(net, opts, '4', 3, 3, 256, 256, 1, 1) ;
net = add_block(net, opts, '5', 3, 3, 256, 256, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 6, 6, 256, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

% --------------------------------------------------------------------
function net = vgg_vd(net, opts)
% --------------------------------------------------------------------

net.layers = {} ;
net = add_block(net, opts, '1_1', 3, 3, 3, 64, 1, 1) ;
net = add_block(net, opts, '1_2', 3, 3, 64, 64, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool1', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '2_1', 3, 3, 64, 128, 1, 1) ;
net = add_block(net, opts, '2_2', 3, 3, 128, 128, 1, 1) ;
net.layers{end+1} = struct('type', 'pool', 'name', 'pool2', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '3_1', 3, 3, 128, 256, 1, 1) ;
net = add_block(net, opts, '3_2', 3, 3, 256, 256, 1, 1) ;
net = add_block(net, opts, '3_3', 3, 3, 256, 256, 1, 1) ;
if strcmp(opts.model, 'vgg-vd-19')
  net = add_block(net, opts, '3_4', 3, 3, 256, 256, 1, 1) ;
end
net.layers{end+1} = struct('type', 'pool', 'name', 'pool3', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '4_1', 3, 3, 256, 512, 1, 1) ;
net = add_block(net, opts, '4_2', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '4_3', 3, 3, 512, 512, 1, 1) ;
if strcmp(opts.model, 'vgg-vd-19')
  net = add_block(net, opts, '4_4', 3, 3, 512, 512, 1, 1) ;
end
net.layers{end+1} = struct('type', 'pool', 'name', 'pool4', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '5_1', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5_2', 3, 3, 512, 512, 1, 1) ;
net = add_block(net, opts, '5_3', 3, 3, 512, 512, 1, 1) ;
if strcmp(opts.model, 'vgg-vd-19')
  net = add_block(net, opts, '5_4', 3, 3, 512, 512, 1, 1) ;
end
net.layers{end+1} = struct('type', 'pool', 'name', 'pool5', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;

net = add_block(net, opts, '6', 7, 7, 512, 4096, 1, 0) ;
net = add_dropout(net, opts, '6') ;

net = add_block(net, opts, '7', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '7') ;

net = add_block(net, opts, '8', 1, 1, 4096, 1000, 1, 0) ;
net.layers(end) = [] ;
if opts.batchNormalization, net.layers(end) = [] ; end

repo: CREST-Release


text:
function [ p1, p2, p3 ] = plane_imp2exp_3d ( a, b, c, d )

%*****************************************************************************80
%
%% PLANE_IMP2EXP_3D converts an implicit plane to explicit form in 3D.
%
%  Discussion:
%
%    The implicit form of a plane in 3D is
%
%      A * X + B * Y + C * Z + D = 0.
%
%    The explicit form of a plane in 3D is
%
%      the plane through P1, P2 and P3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, D, the implicit plane parameters.
%
%    Output, real P1(3,1), P2(3,1), P3(3,1), three points on the plane.
%
  [ pp, normal ] = plane_imp2normal_3d ( a, b, c, d );

  [ p1, p2, p3 ] = plane_normal2exp_3d ( pp, normal );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plot_mls1p(vPar, vMlsprob, nType)
    % Plot maximum likelihod score (MLS) versus parameter variation
    %  plot_mls1p(vPar, vMlsprob, nType)
    % -----------------------------------------------
    % Plot maximum likelihod score (MLS) versus parameter variation
    % vPar is e.g. a shift, a stretch or a rate change
    %
    % Incoming variable
    % vPar          : Vector of parameter used in search
    % vMlsprob      : Maximum likelihood score for parameter
    % nType    : Parameter identification
    %            1 = dM (Simple shift)
    %            2 = dS (Simple stretch)
    %            3 = Rf (Rate factor)
    %
    % See also plot_mls2p, plot_mls3p
    %
    % Author: J. Woessner, woessner@seismo.ifg,.ethz.ch
    % updated: 28.10.02

    switch nType
        case 1
            sTitle = 'Simple magnitude shift';
            sX = 'dM';
        case 2
            sTitle = 'Simple stretch';
            sX = 'dS';
        case 3
            sTitle = 'Rate factor';
            sX = 'R_f';
        otherwise
            return;
    end
    if exist('mls2_fig','var') &  ishandle(mls2_fig)
        set(0,'Currentfigure',mls2_fig);
    else
        mls2_fig=figure_w_normalized_uicontrolunits('tag','mls2','Name','Max. Likelikehood score','Units','normalized','Nextplot','add',...
            'Numbertitle','off');
        mls2_axs=axes('tag','ax_mls2','Nextplot','add','box','off');
    end

    set(gcf,'tag','mls2');
    set(gca,'tag','ax_mls2','Nextplot','replace','box','off','visible','off');
    plot(vPar, log10(vMlsprob), 'r-^');
    xlabel(sX);
    ylabel('log10(Likelihood score)');
    title(sTitle);
end

author: CelsoReyes


text:
function nd = lagrange_interp_nd_size2 ( m, ind )

%*****************************************************************************80
%
%% LAGRANGE_INTERP_ND_SIZE2 sizes an M-dimensional Lagrange interpolant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer IND(M), the index or level of the 1D rule
%    to be used in each dimension.
%
%    Output, integer ND, the number of points in the product grid.
%
  nd = 1;
  for i = 1 : m
    n = order_from_level_135 ( ind(i) );
    nd = nd * n;
  end

  return
end


text:
function endPointsRoi = dtiCreateRoiFromFiberEndPoints(fg,saveFlag)
% function dtiCreateRoiFromFiberEndPoints(fiberGroup,[saveFlag])
%
% Ths fucntion creates an roi from the endpoints of a given fiber
% group passed in by the user. It borrows code from dtiFiberUI.
% You must first load the fiber group with fg = dtiReadFibers('your.mat'),
% then pass in fg to this fucntion.
%
% If saveFlag == 1 the endPointsRoi will be saved to the current directory
% with the name of the fiber group with _endPointsROI added to the name.
%
% Example:
% fg = dtiReadFibers('your.mat');
% newROI = dtiCreateRoiFromFiberEndPoints(fg);
% dtiWriteRoi(newROI,newROI.name);
%
%
% History:
%  2009.03.25 LMP Wrote it.
%

if notDefined('saveFlag'), saveFlag = 0; end

if(~exist('fg','var') || isempty(fg))
    error 'You must pass in fiber group (fg) struct!';
end

nfibers = length(fg.fibers);
fc = zeros(nfibers*2, 3);
for jj=1:nfibers
    fc((jj-1)*2+1,:) = [fg.fibers{jj}(:,1)'];
    fc((jj-1)*2+2,:) = [fg.fibers{jj}(:,end)'];
end
endPointsRoi = dtiNewRoi([fg.name '_endPointsROI'], fg.colorRgb./255, fc);

if saveFlag == 1
    dtiWriteRoi(endPointsRoi,endPointsRoi.name);
    disp([endPointsRoi.name,' saved to ' pwd]);
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% GAINSANDREFERENCES compute gains matrices, references and other control
%                    related quantities for each state of the state machine.

%% --- Initialization ---

% CoM gains
Gain.KP_CoM =  [50    100  5  % state ==  1  TWO FEET BALANCING
                50    100  5  % state ==  2  COM TRANSITION TO LEFT
                50    100  5  % state ==  3  LEFT FOOT BALANCING
                50    100  5  % state ==  4  YOGA LEFT FOOT
                50    100  5  % state ==  5  PREPARING FOR SWITCHING
                50    100  5  % state ==  6  LOOKING FOR CONTACT
                50    100  5  % state ==  7  TRANSITION TO INITIAL POSITION
                50    100  5  % state ==  8  COM TRANSITION TO RIGHT FOOT
                50    100  5  % state ==  9  RIGHT FOOT BALANCING
                50    100  5  % state == 10  YOGA RIGHT FOOT
                50    100  5  % state == 11  PREPARING FOR SWITCHING
                50    100  5  % state == 12  LOOKING FOR CONTACT
                50    100  5];% state == 13  TRANSITION TO INITIAL POSITION

Gain.KD_CoM = 2*sqrt(Gain.KP_CoM)/20;

% Angular momentum gains
Gain.KI_AngularMomentum = 3;
Gain.KP_AngularMomentum = 2*sqrt(Gain.KI_AngularMomentum)/5;

% Postural task gains
%                   %  TORSO  %%        LEFT ARM   %%       RIGHT ARM   %%         LEFT LEG           %%         RIGHT LEG           %%
Gain.KP_postural = [10   30   20, 10   10    10    8, 10   10    10    8, 30   30   20    20    100 100, 30   50   30    60    100 100  % state ==  1  TWO FEET BALANCING
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   30   20    20    100 100, 30   50   30    60    100 100  % state ==  2  COM TRANSITION TO LEFT
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   50   30    60    100 100, 30   30   20    20    100 100  % state ==  3  LEFT FOOT BALANCING
                    30   30   30, 10   10    10   10, 10   10    10   10,100  200  100   400    100 100,100   50   30   100    100 100  % state ==  4  YOGA LEFT FOOT
                    30   30   30,  5    5    10   10, 10   10    20   10,200  250   20    20     10  10,220  550  220   200     65 300  % state ==  5  PREPARING FOR SWITCHING
                    30   30   30, 10   10    20   10, 10   10    20   10,100  350   20   200     10 100,220  550  220   200     65 300  % state ==  6  LOOKING FOR CONTACT
                    30   30   30, 10   10    10    8, 10   10    10    8, 30   50   60    30      5   5, 30   30   30    20      5   5  % state ==  7  TRANSITION TO INITIAL POSITION
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   50   60    30    100 100, 30   30   30    20    100 100  % state ==  8  COM TRANSITION TO RIGHT FOOT
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   50   30    60    100 100, 30   30   20    20    100 100  % state ==  9  RIGHT FOOT BALANCING
                    30   30   30, 10   10    10   10, 10   10    10   10,100   50   30   100    100 100,100  200  100   100     10  10  % state == 10  YOGA RIGHT FOOT
                    30   30   30, 10   10    10   10, 10   10    10   10,220  550  220   200     65 300,200  250   20    20     10  10  % state == 11  PREPARING FOR SWITCHING
                    30   30   30, 10   10    10   10, 10   10    10   10,220  550  220   200     65 300,100  350   20   200     10 100  % state == 12  LOOKING FOR CONTACT
                    30   30   30, 10   10    10   10, 10   10    10   10,220  550  220   200     65 300,100  350   20   200     10 100];% state == 13  TRANSITION TO INITIAL POSITION

Gain.KD_postural = 2*sqrt(Gain.KP_postural(1,:))/20;

% symmetric gains
Gain.KP_postural(10,18:23) = Gain.KP_postural(10,18:23)*1.5;
Gain.KP_postural(10,1:3)   = Gain.KP_postural(10,1:3)*2;
Gain.KP_postural(4,12:17)  = Gain.KP_postural(4,12:17)*1.5;
Gain.KP_postural(4,1:3)    = Gain.KP_postural(4,1:3)*2;
Gain.KP_postural(4,4:11)   = Gain.KP_postural(4,4:11)*1.5;
Gain.KP_postural(7,2:3)    = Gain.KP_postural(7,2:3)*1.5;

%% Smoothing times

% Smoothing time gain scheduling
Config.SmoothingTimeGainScheduling = 2;

% Smoothing time CoM references
StateMachine.CoMSmoothingTime   = [1;    %% state ==  1  TWO FEET BALANCING
                                   1;    %% state ==  2  COM TRANSITION TO LEFT FOOT
                                   1;    %% state ==  3  LEFT FOOT BALANCING
                                   0.65; %% state ==  4  YOGA LEFT FOOT
                                   2;    %% state ==  5  PREPARING FOR SWITCHING
                                   2;    %% state ==  6  LOOKING FOR CONTACT
                                   1;    %% state ==  7  TRANSITION INIT POSITION
                                   1;    %% state ==  8  COM TRANSITION TO RIGHT FOOT
                                   1;    %% state ==  9  RIGHT FOOT BALANCING
                                   0.65; %% state == 10  YOGA RIGHT FOOT
                                   2;    %% state == 11  PREPARING FOR SWITCHING
                                   5;    %% state == 12  LOOKING FOR CONTACT
                                   10];  %% state == 13  TRANSITION INIT POSITION

% Smoothing time for joints references
StateMachine.jointsSmoothingTime = [1;    %% state ==  1  TWO FEET BALANCING
                                    1;    %% state ==  2  COM TRANSITION TO LEFT FOOT
                                    1;    %% state ==  3  LEFT FOOT BALANCING
                                    0.65; %% state ==  4  YOGA LEFT FOOT
                                    2;    %% state ==  5  PREPARING FOR SWITCHING
                                    2;    %% state ==  6  LOOKING FOR CONTACT
                                    1;    %% state ==  7  TRANSITION INIT POSITION
                                    1;    %% state ==  8  COM TRANSITION TO RIGHT FOOT
                                    1;    %% state ==  9  RIGHT FOOT BALANCING
                                    0.65; %% state == 10  YOGA RIGHT FOOT
                                    2;    %% state == 11  PREPARING FOR SWITCHING
                                    5;    %% state == 12  LOOKING FOR CONTACT
                                    10];  %% state == 13  TRANSITION INIT POSITION

% scale factor smoothing time multiplies the smoothing factor during the
% Yoga (state 4 and 10). The purpose is to reduce the time necessary for
% the reference to converge to the next position, but without changing also
% the valuse stored in Sm.joints_leftYogaRef/Sm.joints_rightYogaRef
StateMachine.scaleFactorSmoothingTime = 0.9;

%% CoM delta

% To be summed to the reference CoM position
StateMachine.CoM_delta  = [% THIS REFERENCE IS USED AS A DELTA W.R.T. THE POSITION OF THE LEFT FOOT
                           0.0,  0.00,  0.0;   %% NOT USED
                           0.0,  0.00,  0.0;   %% state ==  2  COM TRANSITION TO LEFT FOOT
                           0.0,  0.005, 0.0;   %% state ==  3  LEFT FOOT BALANCING
                           0.0,  0.005, 0.0;   %% state ==  4  YOGA LEFT FOOT
                           0.0,  0.00,  0.0;   %% state ==  5  PREPARING FOR SWITCHING
                           0.02, -0.08, 0.0;   %% state ==  6  LOOKING FOR CONTACT
                           0.0,  0.00,  0.0;   %% state ==  7  TWO FEET BALANCING
                           % THIS REFERENCE IS USED AS A DELTA W.R.T. THE POSITION OF THE RIGHT FOOT
                           0.0,  0.00,  0.0;   %% state ==  8  COM TRANSITION TO RIGHT FOOT
                           0.0,  0.00,  0.0;   %% state ==  9  RIGHT FOOT BALANCING
                           0.0, -0.01,  0.0;   %% state == 10  YOGA RIGHT FOOT
                           0.0,  0.00,  0.0;   %% state == 11  PREPARING FOR SWITCHING
                           0.0,  0.09,  0.0;   %% state == 12  LOOKING FOR CONTACT
                           0.0,  0.00,  0.0];  %% NOT USED

%% Joint references
StateMachine.joints_references = [ zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [-0.0348,0.0779,0.0429, ...                          %% state == 2  COM TRANSITION TO LEFT
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630, ...   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151];     %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 3  LEFT FOOT BALANCING
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630, ...   %
                                   0.0005,0.0793,-0.0014,-0.0051,-0.1060,-0.1151];    %
                                   zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [-0.0348,0.0779,0.0429, ...                          %% state == 5  PREPARING FOR SWITCHING
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630, ...   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151];     %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 6  LOOKING FOR CONTACT
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                   0.0107,-0.0741,-0.0001,-0.0120,0.0252,0.1369, ...  %
                                  -0.0026,0.0225,0.0093,-0.0020,0.0027,-0.0277];      %
                                   zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 8  COM TRANSITION TO RIGHT FOOT
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                   0.0107,-0.0741,-0.0001,-0.0120,0.0252,0.1369, ...  %
                                  -0.0026,0.0225,0.0093,-0.0020,0.0027,-0.0277];      %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 9  RIGHT FOOT BALANCING
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151, ...  %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630];      %
                                   zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [-0.0348,0.0779,0.0429, ...                          %% state == 11  PREPARING FOR SWITCHING
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151, ...  %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630];      %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 12  LOOKING FOR CONTACT
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                  -0.0026,0.0225,0.0093,-0.0020,0.0027,-0.0277, ...   %
                                   0.0107,-0.0741,-0.0001,-0.0120,0.0252,0.1369];     %
                                   zeros(1,ROBOT_DOF)];                               %% THIS REFERENCE IS IGNORED

% YOGA MOVESET (joint references during state 4 and 10)
q1 =        [-0.0790,0.2279, 0.4519, ...
             -1.1621,0.6663, 0.4919, 0.9947, ...
             -1.0717,1.2904,-0.2447, 1.0948, ...
              0.2092,0.2060, 0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3484,0.4008,-0.0004,-0.3672,-0.1060,-0.0875];

q2 =        [-0.0790,0.1279, 0.4519, ...
             -1.1621,0.6663, 0.4965, 0.9947, ...
             -1.0717,1.2904,-0.2493, 1.0948, ...
              0.2092,0.2060, 0.0006,-0.1741,-0.1044,0.0700, ...
              0.3714,0.9599, 1.3253,-1.6594,-0.1060,-0.0614];

q3 =        [-0.0852,-0.3273,0.0821,...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092,0.2060, 0.0006,-0.1741,-0.1044,0.0700, ...
              0.3714,0.9599, 1.3253,-1.6594, 0.2443,-0.0614];

q4 =        [-0.0852,-0.4273,0.0821,...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.3473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.2443,-0.0614];

q5 =        [-0.0790,-0.2273, 0.4519, ...
             -1.1621,0.6663, 0.4965, 0.9947, ...
             -1.0717,1.2904,-0.2493, 1.0948, ...
              0.2092, 0.4473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.2443,-0.0614];

q6 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.2443,-0.0614];

q7 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253, -1.6217, 0.2443,-0.0614];

q8 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.2443,-0.0614];

q9 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 0.0107,1.3253,-0.0189, 0.2443,-0.0614];

q10 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.2443,-0.0614];

q11 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 1.3107,1.3253,-0.0189, 0.2443,-0.0614];

q12 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 0.3107,1.3253,-0.0189, 0.2443,-0.0614];

q13 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 1.3107,1.3253,-0.0189, 0.2443,-0.0614];

q14 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 0.0107,1.3253,-0.0189, 0.2443,-0.0614];

q15 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              1.5514, 0.3107,1.3253,-0.0189, 0.2443,-0.0614];

q16 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.2514, 0.0107,1.3253,-0.0189, 0.2443,-0.0614];

q17 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
             -0.3514, 0.3107,1.3253,-0.0189, 0.2443,-0.0614];

StateMachine.joints_leftYogaRef  = [ 0,                                    q1;
                                     1*StateMachine.jointsSmoothingTime(4),q2;
                                     2*StateMachine.jointsSmoothingTime(4),q3;
                                     3*StateMachine.jointsSmoothingTime(4),q4;
                                     4*StateMachine.jointsSmoothingTime(4),q5;
                                     5*StateMachine.jointsSmoothingTime(4),q6;
                                     6*StateMachine.jointsSmoothingTime(4),q7;
                                     7*StateMachine.jointsSmoothingTime(4),q8;
                                     8*StateMachine.jointsSmoothingTime(4),q9;
                                     9*StateMachine.jointsSmoothingTime(4),q10;
                                    10*StateMachine.jointsSmoothingTime(4),q11;
                                    11*StateMachine.jointsSmoothingTime(4),q12;
                                    12*StateMachine.jointsSmoothingTime(4),q13;
                                    13*StateMachine.jointsSmoothingTime(4),q14;
                                    14*StateMachine.jointsSmoothingTime(4),q15;
                                    15*StateMachine.jointsSmoothingTime(4),q16;
                                    16*StateMachine.jointsSmoothingTime(4),q17;
                                    17*StateMachine.jointsSmoothingTime(4),q10;
                                    18*StateMachine.jointsSmoothingTime(4),q11;
                                    19*StateMachine.jointsSmoothingTime(4),q12;
                                    20*StateMachine.jointsSmoothingTime(4),q13;
                                    21*StateMachine.jointsSmoothingTime(4),q14;
                                    22*StateMachine.jointsSmoothingTime(4),q15;
                                    23*StateMachine.jointsSmoothingTime(4),q16;
                                    24*StateMachine.jointsSmoothingTime(4),q17;
                                    25*StateMachine.jointsSmoothingTime(4),q8];

StateMachine.joints_rightYogaRef      = StateMachine.joints_leftYogaRef;
StateMachine.joints_rightYogaRef(:,1) = [0;
                                         1*StateMachine.jointsSmoothingTime(10);
                                         2*StateMachine.jointsSmoothingTime(10);
                                         3*StateMachine.jointsSmoothingTime(10);
                                         4*StateMachine.jointsSmoothingTime(10);
                                         5*StateMachine.jointsSmoothingTime(10);
                                         6*StateMachine.jointsSmoothingTime(10);
                                         7*StateMachine.jointsSmoothingTime(10);
                                         8*StateMachine.jointsSmoothingTime(10);
                                         9*StateMachine.jointsSmoothingTime(10);
                                        10*StateMachine.jointsSmoothingTime(10);
                                        11*StateMachine.jointsSmoothingTime(10);
                                        12*StateMachine.jointsSmoothingTime(10);
                                        13*StateMachine.jointsSmoothingTime(10);
                                        14*StateMachine.jointsSmoothingTime(10);
                                        15*StateMachine.jointsSmoothingTime(10);
                                        16*StateMachine.jointsSmoothingTime(10);
                                        17*StateMachine.jointsSmoothingTime(10);
                                        18*StateMachine.jointsSmoothingTime(10);
                                        19*StateMachine.jointsSmoothingTime(10);
                                        20*StateMachine.jointsSmoothingTime(10);
                                        21*StateMachine.jointsSmoothingTime(10);
                                        22*StateMachine.jointsSmoothingTime(10);
                                        23*StateMachine.jointsSmoothingTime(10);
                                        24*StateMachine.jointsSmoothingTime(10);
                                        25*StateMachine.jointsSmoothingTime(10)];

% if the demo is not "yogaExtended", stop at the 8th move
if ~StateMachine.yogaExtended

    StateMachine.joints_leftYogaRef       = StateMachine.joints_leftYogaRef(1:8,:);
    StateMachine.joints_rightYogaRef      = StateMachine.joints_rightYogaRef(1:8,:);
    StateMachine.joints_leftYogaRef(8,1)  = 15*StateMachine.jointsSmoothingTime(4);
    StateMachine.joints_rightYogaRef(8,1) = 15*StateMachine.jointsSmoothingTime(10);
end

% MIRROR YOGA LEFT MOVESET FOR RIGHT YOGA
for i = 1:size(StateMachine.joints_rightYogaRef,1)

	StateMachine.joints_rightYogaRef(i,2:4)           = [StateMachine.joints_rightYogaRef(i,2) -StateMachine.joints_rightYogaRef(i,3) -StateMachine.joints_rightYogaRef(i,4)];
	rightArm                                          =  StateMachine.joints_rightYogaRef(i,end-15:end-12);
	StateMachine.joints_rightYogaRef(i,end-15:end-12) =  StateMachine.joints_rightYogaRef(i,end-19:end-16);
	StateMachine.joints_rightYogaRef(i,end-19:end-16) =  rightArm;
	rightLeg                                          =  StateMachine.joints_rightYogaRef(i,end-5:end);
	StateMachine.joints_rightYogaRef(i,end-5:end)     =  StateMachine.joints_rightYogaRef(i,end-11:end-6);
	StateMachine.joints_rightYogaRef(i,end-11:end-6)  =  rightLeg;
end

%% References for CoM trajectory (COORDINATOR DEMO ONLY)

% that the robot waits before starting the left-and-right
Config.noOscillationTime       = 0;
Config.directionOfOscillation  = [0;1;0];
Config.amplitudeOfOscillation  = 0.02; % [m]
Config.frequencyOfOscillation  = 0.2;  % [Hz]

author: robotology

repo: whole-body-controllers


text:
% Test file for ADCHEBFUN VOLT

function pass = test_volt

K = @(s, t)  exp(-(s-t).^2);

% List of trigonometric functions to test.
func = @(u) volt(K, u);

% Tolerance for Taylor testing
tolOrder = 1e-2;
tolDiff = 1e-12;
% Initialise vector with pass information
pass = zeros(1, numel(func));

% Call the valueTesting method, which also returns linearity information
[err, lin] = adchebfun.valueTesting(func);

% First, check that the computed function values match what we expect
pass(1) = ( err == 0 );

% Call the taylorTesting method
[order1, order2, nDiff2] = adchebfun.taylorTesting(func,2);

% We expect all elements of ORDER1 to be close to 1. Since the methods being
% tested in this case are all linear, ORDER2 will be noise. However, since
% the methods are indeed linear, we should expect nDiff2 to have values all
% close to machine epsilon, which we can use to check for the correctness of
% the derivative computed.
pass(2) = ( (max(abs(order1 - 1)) < tolOrder) && ...
    (max(abs(nDiff2)) < tolDiff) );

% Check that we received the correct linearity information
pass(3) = ( lin == 1 );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% DEMSWISSROLLLLE2 Demonstrate LLE on the oil data.

% MLTOOLS

[Y, lbls] = lvmLoadData('swissRoll');

options = lleOptions(8);
model = lleCreate(2, size(Y, 2), Y, options);
model = lleOptimise(model, 2);

lvmScatterPlotColor(model, model.Y(:, 2));

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, model.Y(:, 2), 'SwissRoll', 2, true);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%IKINE_SYM  Symbolic inverse kinematics
%
% Q = R.IKINE_SYM(K, OPTIONS) is a cell array (Cx1) of inverse kinematic
% solutions of the SerialLink object ROBOT.  The cells of Q represent the
% solutions for each joint, ie. Q{1} is the solution for joint 1.  A
% cell may contain an array of solutions. The solution is expressed in terms
% of other joint angles and elements of the desired end-point pose which is
% represented by the symbolic matrix (3x4) with elements
%      nx ox ax tx
%      ny oy ay ty
%      nz oz az tz
% where the first three columns specify orientation and the last column
% specifies translation.
%
% K <= N is the number of joint angles solved for.
%
% Options::
%
% 'file',F    Write the solution to an m-file named F
% 'Tpost',T   Add a symbolic 4x4 matrix T to the end of the chain
%
% Example::
%
%         mdl_planar2
%         sol = p2.ikine_sym(2);
%         length(sol)
%
%         q1 = sol{1}   % are the solution for joint 1
%         q2 = sol{2}   % is the solution for joint 2
%         length(q1)
%         ans =
%               2     % there are 2 solutions for this joint
%         q1(1)       % one solution for q1
%         q1(2);      % the other solution for q1
%
% Notes::
% - ignores tool and base transforms.
%
% References::
% - Robot manipulators: mathematics, programming and control
%   Richard Paul, MIT Press, 1981.
% - The kinematics of manipulators under computer control,
%   D.L. Pieper, Stanford report AI 72, October 1968.
%
% Notes::
% - Requires the MATLAB Symbolic Math Toolbox.
% - This code is experimental and has a lot of diagnostic prints.
% - Based on the classical approach using Pieper's method.

% 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

function out = ikine_sym(robot, N, varargin)

    %
    % Given a robot model the following steps are performed:
    % 1. Convert model to symbolic form
    % 2. Find relevant trig equations and solve them for joint angles
    % 3. Write an M-file to implement the solution
    %      xikine(T)
    %      xikine(T, S) where S is a 3 vector with elements 1 or 2 to select
    %       the first or second solution for the corresponding joint.
    %
    % TODO:
    %  - handle the wrist joints, only first 3 joints so far
    %  - handle base and tool transforms
    %  - allow 3DOF solution which is rotation, ie. wrist

    opt.file = [];
    opt.Tpost = [];
    opt.all = false;
    opt = tb_optparse(opt, varargin);

    % make a symbolic representation of the passed robot
    srobot = SerialLink(robot);  % make a deep copy
    srobot = sym(srobot);  % convert to symbolic
    q = srobot.gencoords();

    % define symbolic elements of the homogeneous transform
    syms nx ox ax tx real
    syms ny oy ay ty real
    syms nz oz az tz real
    syms d4 real

    % inits
    Q = {};
    trigsubOld = [];
    trigsubNew = [];

    % loop over each joint variable
    for j=1:N
        fprintf('----- solving for joint %d\n', j);

        % create some equations to sift through
        [left,right] = pieper(srobot, j, 'left', opt.Tpost);

        % decide which equations to look at
        if j <= 3
            % for first three joints only focus on translational part
            left = left(1:3, 4); left = left(:);
            right = right(1:3, 4); right = right(:);
        else
            % for last three joints only focus on rotational part
            left = left(1:3, 1:3); left = left(:);
            right = right(1:3, 1:3); right = right(:);
        end

        % substitute sin/cos for preceding joint as Sj/Cj, essentially removes
        % the joint variables from the equations and treats them as constants.
        if ~isempty(trigsubOld)
            left = rsubs(left, trigsubOld, trigsubNew);
            right = rsubs(right, trigsubOld, trigsubNew);
        end

        % then simplify the LHS
        %   do it after the substitution to prevent sum of angle terms being introduced
        left = simplify(left);

        % search for a solveable equation:
        %    function of current joint variable on the LHS
        %    constant element on the RHS
        k = NaN;
        for i=1:length(left)
            fprintf('%d: %s\n', i, char(left(i) == right(i)));
            if hasonly(left(i), j) && isconstant(right(i))
                k = i;
            end
        end

        eq = [];

        if ~isnan(k)
            % create the equation to solve: LHS-RHS == 0
            fprintf('choosing equation %d\n', k);

            eq = left(k) - right(k);
        else
            % ok, we weren't lucky, try another strategy

            % find all equations:
            %    function of current joint variable on the LHS

            k = [];
            for i=1:length(left)
                % has qj on the left
                if hasonly(left(i), j)
                    k = [k i];
                end
            end

            % hopefully we found at least two of them
            if length(k) < 2
                continue;
            end

            % we did, lets see if the sum square RHS is constant
            for kk = nchoosek(k, 2)'
                fprintf('No simple equation, let''s square and add RHS %d %d\n', kk);
                rhs = simplify(right(kk(1))^2 + right(kk(2))^2); % was simple
                if isconstant( rhs )
                    eq = simplify( expand( left(kk(1))^2 + left(kk(2))^2 ) ) - rhs;
                    break
                end
            end
            if isempty(eq)
                fprintf('** can''t solve this equation, out of options');
                k
                left(k)==right(k)
                error('can''t solve this equation');
            end
            % ensure that S^2+C^2=1 substitutions are made
            eq = subs(eq, trigsubOld, trigsubNew);
            disp(eq)
        end

        % expand the list of joint variable subsitutions
        fprintf('subs sin/cos q%d for S/C\n', j);
        trigsubOld = [trigsubOld mvar('sin(q%d)', j)^2+mvar('cos(q%d)', j)^2 mvar('sin(q%d)', j) mvar('cos(q%d)', j) ];
        trigsubNew = [trigsubNew 1 mvar('S%d', j) mvar('C%d', j) ];

        % now solve the equation
        if srobot.links(j).isrevolute()
            % for revolute joint it will be a trig equation, do we know how to solve it?
            Q{j} = solve_joint(eq, j );
            if isempty(Q)
                warning('RTB:ikine_sym', 'can''t solve this kind of equation');
            end
        else
            fprintf('prismatic case\n')
            q = mvar('q%d', j);
            Q{j} = solve( eq == 0, q);
        end
    end

    % final simplification
    %  get rid of C^2+S^2 and C^4, S^4 terms
    fprintf('**final simplification pass\n')

    trigsubOld = [];
    trigsubNew = [];
    for j=1:N
        trigsubOld = [trigsubOld mvar('S%d', j) mvar('C%d', j) ];
        trigsubNew = [trigsubNew mvar('sin(q%d)', j) mvar('cos(q%d)', j) ];
    end

    Q = simplify_powers(Q, N, trigsubOld, trigsubNew);

    % Q is a cell array of equations for joint variables
    if nargout > 0
        out = Q;
    end

    if ~isempty(opt.file)
        fprintf('**generate MATLAB code\n')
        gencode(Q);
    end
end

%PIEPER Return a set of equations using Pieper's method
%
% [L,R] = pieper(robot, n, which)
%
% If robot has link matrix A1 A2 A3 A4 then returns 12 equations from equating the coefficients of
%
%  A1' T = A2 A3 A4     n=1, which='left'
%  A2' A1' T = A3 A4    n=2, which='left'
%  A3' A2' A1' T = A4   n=3, which='left'
%
%  T A4' = A1 A2 A3     n=1, which='right'
%  T A4' A3' = A1 A2    n=2, which='right'
%  T A4' A3' A2' = A1   n=3, which='right'
%
% A' denotes inversion not transposition
%
% Judicious choice of the equations can lead to joint solutions

function [L,R] = pieper(robot, n, which, Tpost)

    if nargin < 3
        which = 'left';
    end
    if nargin < 4 || isempty(Tpost)
        Tpost = transl(zeros(robot.n, 3));
    end
    assert(n <= robot.n, 'RTB:ikine_sym:badarg', 'N is greater than number of joints');

    syms nx ox ax tx real
    syms ny oy ay ty real
    syms nz oz az tz real

    T = [nx ox ax tx
        ny oy ay ty
        nx oz az tz
        0  0  0  1 ];

    T = inv(robot.base.T) * T * inv(robot.tool.T);

    q = robot.gencoords();

    % Create the symbolic A matrices
    for j=1:robot.n
        A{j} = robot.links(j).A(q(j)).T * Tpost(:,:,j);
    end

    switch which
        case 'left'
            left = T;
            for j=1:n
                left = inv(A{j}) * left ;
            end

            right = eye(4,4);
            for j=n+1:robot.n
                right = right * A{j};
            end

        case 'right'
            left = T;
            for j=1:n
                left = left * inv(A{robot.n-j+1});
            end

            right = eye(4,4);
            for j=1:(robot.n-n)
                right = right * A{j};
            end
    end

    %     left = simple(left);
    %     right = simple(right);

    if nargout == 0
        left == right
    elseif nargout == 1
        L = left;
    elseif nargout == 2
        L = left;
        R = right;
    end
end

%SOLVE_JOINT Solve a trigonometric equation
%
% S = SOLVE_JOINT(EQ, J) solves the equation EQ=0 for the joint variable qJ.
% The result is a vector of symbolic solutions.
%
% The equations must be of the form:
%  A cos(qJ) + B sin(qJ) = 0
%  A cos(qJ) + B sin(qJ) = C
%
% where A, B, C are arbitrarily complex expressions.  qJ can be the only
% joint variable in the expression.
%
% Notes::
% - In general there are two solutions, but if A^2+B^2-C^2 = 0, then only one
%   solution is returned.
% - The one solution case may not be detected symbolically, in which case
%   the two returned solutions will have the same value after numerical
%   substitution
% - The symbolic solution may not be evaluteable numerically with
%   certain parameter values, ie. if A^2+B^2-C^2 < 0

function s = solve_joint(eq, j)

    % see http://petercorke.com/wordpress/solving-trigonometric-equations

    sinj = mvar('sin(q%d)', j);
    cosj = mvar('cos(q%d)', j);

    A = getcoef(eq, cosj);
    B = getcoef(eq, sinj);

    if isempty(A) || isempty(B)
        warning('don''t know how to solve this kind of equation');
    end

    C = -simplify(eq - A*cosj - B*sinj);

    A = simplify_sumsq(A, j);
    B = simplify_sumsq(B, j);
    C = simplify_sumsq(C, j);

    fprintf('A = %s\n', char(A));
    fprintf('B = %s\n', char(B));
    fprintf('C = %s\n', char(C));

    if C == 0
        fprintf('Solve for C == 0\n');
        % A cos(q) + B sin(q) = 0
        s(1) = atan2(A, -B);
        s(2) = atan2(-A, B);
    else
        fprintf('Solve for C != 0\n');
        % A cos(q) + B sin(q) = C
%         r = sqrt(A^2 + B^2 - C^2);
% %         phi = atan2(A, B);
% %
% %         s(2) = atan2(C, r) - phi;
% %         s(1) = atan2(C, -r) - phi;
%         if r == 0
%             s = atan2(B*C, A*C)
%         else
%             s(1) = atan2(B*C + A*r, A*C - B*r);
%             s(2) = atan2(B*C - A*r, A*C + B*r);
%         end
        d = sqrt( simplify(B^2-C^2+A^2) );
        fprintf('d = %s\n', char(d));

        s(1) = atan2(B*C + A*d, A*C - B*d);
        s(2) = atan2(B*C - A*d, A*C + B*d);
    end

%     simplify( A*cos(s(1)) + B*sin(s(1)) - C )
%     simplify( A*cos(s(2)) + B*sin(s(2)) - C )

    if nargout == 0
        try
            eval(s)
        catch
            s
        end
    end
end

function Qout = simplify_sumsq(Q, N)

    tsubOld = [];
    tsubNew = [];
    for j=1:N
        tsubOld = [tsubOld mvar('S%d', j)^2+mvar('C%d', j)^2];
        tsubNew = [tsubNew         1                        ];
    end
    Q = simplify( rsubs(Q, tsubOld, tsubNew) );
    Qout = simplify( rsubs(Q, tsubOld, tsubNew) );
end

function Qout = simplify_powers(Q, N, trigsubOld, trigsubNew)

        % create a list of simplifications
    %  substitute S^2 = 1-C^2, S^4=(1-C^2)^2
    fprintf('power simplification\n');
    tsubOld = [];
    tsubNew = [];
    for j=1:N
        tsubOld = [tsubOld mvar('S%d', j)^2+mvar('C%d', j)^2 mvar('S%d', j)^4];
        tsubNew = [tsubNew         1                        (1-mvar('C%d', j)^2)^2];
    end

    for j=1:length(Q)
        for k=1:5
            % seem to need to iterate this, not quite sure why
            Q{j} = simplify( expand( subs(Q{j}, tsubOld, tsubNew) ) );
        end
        % subs Sx, Cx to sin(qx), cos(qx)
        Qout{j} = simplify( subs(Q{j}, trigsubOld, trigsubNew) );
    end
end

function out = rsubs(eq, old, new)
    for i = 1:length(old)
        eq = subs(eq, old(i), new(i));
    end
    out = eq;
end

%MVAR Create a symbolic variable
%
% V = MVAR(FMT, ARGS) is a symbolic variable created using SPRINTF
%
% eg. mvar('q%d', j)
%
% The symbolic is explicitly declared to be real.

function v = mvar(fmt, varargin)

    if isempty(strfind(fmt, '('))
        % not a function
        v = sym( sprintf(fmt, varargin{:}), 'real' );
    else
        v = str2sym( sprintf(fmt, varargin{:}) );

    end
end

%HASONLY Determine if an expression contains only certain joint variables
%
% S = HASONLY(E L) is true if the joint variables (q1, q2 etc.) in the expression E
% are listed in the vector L.
%
% Eg. hasonly('sin(q1)*cos(q2)*cos(q4)', [1 2 3]) -> true
% Eg. hasonly('sin(q1)*cos(q2)*cos(q4)', [1]) -> false

function s = hasonly(eq, j)

    q = findq(eq);
    if isempty(q)
        s = false;
    else
        s = all(ismember(j, findq(eq)));
    end
end

%ISCONSTANT Determine if an expression is free of joint variables
%
% S = ISCONSTANT(E) is true if the expression E contains no joint variables such
% q1, q2 etc.

function s = isconstant(eq)
    s = isempty(findq(eq));
end

%FINDQ Find the joint variables in expression
%
% Q = FINDQ(E) returns a list of integers indicating the joint variables found
% in the expression E.  For instance an instance of 'q1' would cause a 1 to be
% returned and so on.
%
% Eg. findq('sin(q1)*cos(q2)+S3') -> [1 2]

function q = findq(s)

    q = [];

    for var=symvar(s)
        if isempty(var)
            break
        end
        varname = char(var);
        if varname(1) == 'q'
            q = [q str2num(varname(2:end))];
        end
    end
end

function coef = getcoef(eq, trig)
    z = children( collect(eq, trig) );
    z = children( z(1) );
    coef = z(1);
end

% Output a joint expression to a file

function s = gencode(Q, filename)

    function s = G(s, fmt, varargin)
        s = strvcat(s, sprintf(fmt, varargin{:}));
    end

    s = 'function q = xikine(T, sol)';
    s = G(s, '  if nargin < 2; sol = ones(1, %d); end', length(Q));
    s = G(s, '  px = T(1,4); py = T(2,4); pz = T(3,4);');

    for j=1:3
        Qj = Q{j};   % cast it to subclass
        if length(Qj) == 1
            s = G(s, '  q(%d) = %s', j, matgen2(Qj));
        elseif length(Qj) == 2
            s = G(s, '  if sol(%d) == 1', j);
            s = G(s, '    q(%d) = %s', j, matgen2(Qj(1)));
            s = G(s, '  else');
            s = G(s, '    q(%d) = %s', j, matgen2(Qj(2)));
            s = G(s, '  end');

        end

        s = G(s, '  S%d = sin(q(%d));', j, j);
        s = G(s, '  C%d = cos(q(%d));', j, j);
        s = G(s, ' ');

    end
    s = G(s, 'end');

    fp = fopen(filename, 'w');
    for i=1:numrows(s)
        fprintf(fp, '%s\n', deblank(s(i,:)));
    end
    fclose(fp);

end

% Generate MATLAB code from an expression
%
% Requires a bit of a hack, a subclass of sym (sym2) to do this

function s = matgen2(e)

    s = matgen(sym2(e));

    k = strfind(s, '=');
    s = deblank( s(k+2:end) );
end

author: petercorke

repo: robotics-toolbox-matlab


text:
%detects if we're running Octave
function result = isOctave()
    result = exist('OCTAVE_VERSION') ~= 0;
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function H = mnorm_entropy(U)
warning('This function is deprecated')

% H = mnorm_entropy(U)
%
% U is the Cholesky factor of a covariance matrix C, that is, a lower or
% upper triangylar matrix such that C=U*U' or C=U'*U.

d = size(U,1);

H = d/2*log(2*pi) + 0.5*logdet_chol(U) + d/2;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function bnet = mk_hhmm3(varargin)
% MK_HHMM3 Make a 3 level Hierarchical HMM
% bnet = mk_hhmm3(...)
%
% 3-layer hierarchical HMM where level 1 only connects to level 2, not 3 or obs.
% This enforces sub-models (which differ only in their Q1 index) to be shared.
% Also, we enforce the fact that each model always starts in its initial state
% and only finishes in its final state. However, the prob. of finishing (as opposed to
% self-transitioning to the final state) can be learned.
% The fact that we always finish from the same state means we do not need to condition
% F(i) on Q(i-1), since finishing prob is indep of calling context.
%
% The DBN is the same as Fig 10 in my tech report.
%
%   Q1 ---------->  Q1
%   |              / |
%   |             /  |
%   |  F2 -------    |
%   |  ^         \   |
%   | /|          \  |
%   v  |           v v
%   Q2-| -------->   Q2
%  /|  |             ^
% / |  |            /|
% | |  F3 ---------/ |
% | |  ^           \ |
% | v /              v
% | Q3 ----------->  Q3
% |  |
% \  |
%  v v
%   O
%
%
% Optional arguments in name/value format [default]
%
% Qsizes      - sizes at each level [ none ]
% Osize       - size of O node [ none ]
% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
% Oargs       - cell array of args to pass to the O CPD  [ {} ]
% transprob1  - transprob1(i,j) = P(Q1(t)=j|Q1(t-1)=i)  ['ergodic']
% startprob1  - startprob1(j) = P(Q1(t)=j)  ['leftstart']
% transprob2  - transprob2(i,k,j) = P(Q2(t)=j|Q2(t-1)=i,Q1(t)=k)  ['leftright']
% startprob2  - startprob2(k,j) = P(Q2(t)=j|Q1(t)=k)  ['leftstart']
% termprob2   - termprob2(j,f) = P(F2(t)=f|Q2(t)=j)  ['rightstop']
% transprob3  - transprob3(i,k,j) = P(Q3(t)=j|Q3(t-1)=i,Q2(t)=k)  ['leftright']
% startprob3  - startprob3(k,j) = P(Q3(t)=j|Q2(t)=k)  ['leftstart']
% termprob3   - termprob3(j,f) = P(F3(t)=f|Q3(t)=j)  ['rightstop']
%
% leftstart means the model always starts in state 1.
% rightstop means the model always finished in its last state (Qsize(d)).
%
% Q1:Q3 in slice 1 are of type tabular_CPD
% Q1:Q3 in slice 2 are of type hhmmQ_CPD.
% F2 is of type hhmmF_CPD, F3 is of type tabular_CPD.

ss = 6; D = 3;
Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6;
Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
names = {'Q1', 'Q2', 'Q3', 'F3', 'F2', 'obs'};

intra = zeros(ss);
intra(Q1, Q2) = 1;
intra(Q2, [F2 Q3 obs]) = 1;
intra(Q3, [F3 obs]) = 1;
intra(F3, F2) = 1;

inter = zeros(ss);
inter(Q1,Q1) = 1;
inter(Q2,Q2) = 1;
inter(Q3,Q3) = 1;
inter(F2,[Q1 Q2]) = 1;
inter(F3,[Q2 Q3]) = 1;

% get sizes of nodes
args = varargin;
nargs = length(args);
Qsizes = [];
Osize = 0;
for i=1:2:nargs
  switch args{i},
   case 'Qsizes', Qsizes = args{i+1};
   case 'Osize', Osize = args{i+1};
  end
end
if isempty(Qsizes), error('must specify Qsizes'); end
if Osize==0, error('must specify Osize'); end

% set default params
discrete_obs = 0;
Oargs = {};
startprob1 = 'ergodic';
startprob2 = 'leftstart';
startprob3 = 'leftstart';
transprob1 = 'ergodic';
transprob2 = 'leftright';
transprob3 = 'leftright';
termprob2 = 'rightstop';
termprob3 = 'rightstop';

for i=1:2:nargs
  switch args{i},
   case 'discrete_obs', discrete_obs = args{i+1};
   case 'Oargs',        Oargs = args{i+1};
   case 'Q1args',       Q1args = args{i+1};
   case 'Q2args',       Q2args = args{i+1};
   case 'Q3args',       Q3args = args{i+1};
   case 'F2args',       F2args = args{i+1};
   case 'F3args',       F3args = args{i+1};
  end
end

ns = zeros(1,ss);
ns(Qnodes) = Qsizes;
ns(obs) = Osize;
ns(Fnodes) = 2;

dnodes = [Qnodes Fnodes];
if discrete_obs
  dnodes = [dnodes obs];
end
onodes = [obs];

bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
eclass = bnet.equiv_class;

if strcmp(startprob1, 'ergodic')
  startprob1 = normalise(ones(1,ns(Q1)));
end
if strcmp(startprob2, 'leftstart')
  startprob2 = zeros(ns(Q1), ns(Q2));
  starpbrob2(:, 1) = 1.0;
end
if strcmp(startprob3, 'leftstart')
  startprob3 = zeros(ns(Q2), ns(Q3));
  starpbrob3(:, 1) = 1.0;
end

if strcmp(termprob2, 'rightstop')
  p = 0.9;
  termprob2 = zeros(Qsize(2),2);
  termprob2(:, 2) = p;
  termprob2(:, 1) = 1-p;
  termprob2(1:(Qsize(2)-1), 1) = 1;
end
if strcmp(termprob3, 'rightstop')
  p = 0.9;
  termprob3 = zeros(Qsize(3),2);
  termprob3(:, 2) = p;
  termprob3(:, 1) = 1-p;
  termprob3(1:(Qsize(3)-1), 1) = 1;
end

% SLICE 1

% We clamp untied nodes in the first slice, since their params can't be estimated
% from just one sequence

bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', startprob1, 'adjustable', 0);
bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', startprob2, 'adjustable', 0);
bnet.CPD{eclass(Q3,1)} = tabular_CPD(bnet, Q3, 'CPT', startprob3, 'adjustable', 0);

bnet.CPD{eclass(F2,1)}  = hhmmF_CPD(bnet, F2, Qnodes, 2, D, 'termprob', termprob2);
bnet.CPD{eclass(F3,1)}  = tabular_CPD(bnet, F3, 'CPT', termprob3);

if discrete_obs
  bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, Oargs{:});
else
  bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, Oargs{:});
end

% SLICE 2

bnet.CPD{eclass(Q1,2)} = hhmmQ_CPD(bnet, Q1+ss, Qnodes, 1, D, 'transprob', transprob1, 'startprob', startprob1);
bnet.CPD{eclass(Q2,2)} = hhmmQ_CPD(bnet, Q2+ss, Qnodes, 2, D, 'transprob', transprob2, 'startprob', startprob2);
bnet.CPD{eclass(Q3,2)} = hhmmQ_CPD(bnet, Q3+ss, Qnodes, 3, D, 'transprob', transprob3, 'startprob', startprob3);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function x = r8to_sl ( n, a, b, job )

%*****************************************************************************80
%
%% R8TO_SL solves a R8TO system.
%
%  Discussion:
%
%    The R8TO storage format is used for a Toeplitz matrix, which is constant
%    along diagonals.  Thus, in an N by N Toeplitz matrix, there are at most
%    2*N-1 distinct entries.  The format stores the N elements of the first
%    row, followed by the N-1 elements of the first column (skipping the
%    entry in the first row).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(2*N-1), the R8TO matrix.
%
%    Input, real B(N) the right hand side vector.
%
%    Input, integer JOB,
%    0 to solve A*X=B,
%    nonzero to solve A'*X=B.
%
%    Output, real X(N), the solution vector.
%
  if ( n < 1 )
    x = [];
    return
  end
%
%  Solve the system with the principal minor of order 1.
%
  r1 = a(1);
  x(1) = b(1) / r1;

  if ( n == 1 )
    return
  end
%
%  Recurrent process for solving the system with the Toeplitz matrix.
%
  for nsub = 2 : n
%
%  Compute multiples of the first and last columns of the inverse of
%  the principal minor of order NSUB.
%
    if ( job == 0 )
      r5 = a(n+nsub-1);
      r6 = a(nsub);
    else
      r5 = a(nsub);
      r6 = a(n+nsub-1);
    end

    if ( 2 < nsub )

      c1(nsub-1) = r2;

      for i = 1 : nsub-2
        if ( job == 0 )
          r5 = r5 + a(n+i) * c1(nsub-i);
          r6 = r6 + a(i+1) * c2(i);
        else
          r5 = r5 + a(i+1) * c1(nsub-i);
          r6 = r6 + a(n+i) * c2(i);
        end
      end

    end

    r2 = -r5 / r1;
    r3 = -r6 / r1;
    r1 = r1 + r5 * r3;

    if ( 2 < nsub )

      r6 = c2(1);
      c2(nsub-1) = 0.0E+00;

      for i = 2 : nsub-1
        r5 = c2(i);
        c2(i) = c1(i) * r3 + r6;
        c1(i) = c1(i) + r6 * r2;
        r6 = r5;
      end

    end

    c2(1) = r3;
%
%  Compute the solution of the system with the principal minor of order NSUB.
%
    if ( job == 0 )
      r5 = a(n+1:n+nsub-1) * x(nsub-1:-1:1)';
    else
      r5 = a(2:nsub) * x(nsub-1:-1:1)';
    end

    r6 = ( b(nsub) - r5 ) / r1;

    x(1:nsub-1) = x(1:nsub-1) + c2(1:nsub-1) * r6;
    x(nsub) = r6;

  end

  return
end


text:
function test_failed=test_gabmuleigs
%TEST_GABMULEIGS  Test GABMULEIGS
%
%   Test GABMULEIGS by comparing the output from the iterative and full algorithm.

disp(' ===============  TEST_GABMULEIGS ================');

test_failed=0;

a=20;
M=30;

L=a*M;
N=L/a;

c=randn(M,N);

g=gabtight(a,M,L);

% [V1,D1]=gabmuleigs(10,c,g,a,'iter');
% [V2,D2]=gabmuleigs(10,c,g,a,'full');
F = frame('dgt',g,a,M);
c = framenative2coef(F,c);
[V1,D1]=framemuleigs(F,F,c,10,'iter');
[V2,D2]=framemuleigs(F,F,c,10,'full');

res=norm(D1-D2);

[test_failed,fail]=ltfatdiditfail(res,test_failed);
s=sprintf('GABMULEIGS   L:%3i a:%3i M:%3i %0.5g %s',L,a,M,res,fail);
disp(s);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function d=subdir(rootpath)
% SUBDIR    Returns the paths of all subdirectories of the rootpath
%   d=subdir(rootpath)
%
%   Input:
%       rootpath,   directory to be searched. All subdirectories in the
%                   rootpath are returned
%
%   Output:
%       d,          cell array containing all subdirectories of rootpath
%
%   author:   Tim Myers [myers@metronaviation.com]
%   version:  1.0
%   date:     20-May-2005
%
%   Requires dirdir from P. Macey on Matlab Central File Exchange:
%   http://www.mathworks.com/matlabcentral/files/1570/dirdir.m
%

d=[];
dnew=dirdir(rootpath);

%append subdirectories of root path
for i=1:numel(dnew)
    d=[d; {[rootpath '\' dnew(i).name]}];
end

%initialize search array. 1=directory has been searched.
s=zeros(length(d),1);
while 1
    %find index of subdirectories yet to be searched
    index=find(s==0);
    if isempty(index)
        %all subdirectories have been searched
        return
    else
        %path of subdirectory to be searched
        dtemp=char(d(index(1)));
        %new subdirectories to add to the list
        dnew=dirdir(dtemp);
        %current directory has been searched
        s(index(1))=1;
    end
    %append new directories
    for i=1:numel(dnew)
        d=[d; {[dtemp '\' dnew(i).name]}];
    end
    s=[s; zeros(numel(dnew),1)];
end


text:
 function varargout = num2list(x)
%function varargout = num2list(x)
%
% convert a numeric vector input into a comma separated list of output values
% DOES NOT WORK!

if nargin < 1, num2list_test, ir_usage, end
%nargout
%nargout = numel(x)
%varargout = cell(size(x));

for ii=1:numel(x)
	varargout{ii} = x(ii);
end

function num2list_test
zeros(num2list([2 3 4 5]))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function pdf = discrete_pdf ( x, a, b )

%*****************************************************************************80
%
%% DISCRETE_PDF evaluates the Discrete PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) = B(X) if 1 <= X <= A
%                = 0    otherwise
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the item whose probability is desired.
%
%    Input, integer A, the number of probabilities assigned.
%
%    Input, real B(A), the relative probabilities of
%    outcomes 1 through A.  Each entry must be nonnegative.
%
%    Output, real PDF, the value of the PDF.
%
  b_sum = sum ( b(1:a) );

  if ( 1 <= x & x <= a )
    pdf = b(x) / b_sum;
  else
    pdf = 0.0;
  end

  return
end


text:
function packFonts()
% run extra/genfont.sh first

fonts = dir(fullfile('extra','fonts')) ;
fonts = {fonts([fonts.isdir]).name} ;
fonts(ismember(fonts,{'.','..'})) = [] ;
chars = 'a':'z' ;

im = cell(numel(fonts), numel(chars)) ;
labels = cell(numel(fonts), numel(chars)) ;

for i = 1:numel(fonts)
  for j = 1:numel(chars)
    [p,cmap] = imread(fullfile('extra', 'fonts', fonts{i}, [chars(j) '.png'])) ;
    if isempty(cmap)
      im{i,j} = im2single(p) ;
    else
      im{i,j} = im2single(ind2gray(p,cmap)) ;
    end
    labels{i,j} = j ;
  end
end

imdb.meta.classes = chars ;
imdb.meta.sets = {'train', 'val'} ;
imdb.meta.fonts = fonts ;
imdb.images.id = 1:numel(im) ;
imdb.images.data = cat(3, im{:}) ;
imdb.images.label = cat(2, labels{:}) ;

% Create training and validation sets
sets = ones(numel(fonts), numel(chars)) ;
sets(end-349:end,:) = 2 ;
imdb.images.set = sets(:)' ;

save('data/charsdb.mat', '-struct', 'imdb') ;

author: vedaldi

repo: practical-cnn


text:
function PL=PL_Hata(fc,d,htx,hrx,Etype)
% Hata Model
% Input
%       fc    : carrier frequency[Hz]
%       d     : between base station and mobile station[m]
%       htx   : height of transmitter[m]
%       hrx   : height of receiver[m]
%       Etype : Environment Type('urban','suburban','open')
% output
%       PL    : path loss[dB]
if nargin<5, Etype = 'URBAN'; end
fc=fc/(1e6);
if fc>=150&&fc<=200, C_Rx = 8.29*(log10(1.54*hrx))^2 - 1.1;
 elseif fc>200, C_Rx = 3.2*(log10(11.75*hrx))^2 - 4.97;
 else   C_Rx = 0.8+(1.1*log10(fc)-0.7)*hrx-1.56*log10(fc);
end
PL = 69.55 +26.16*log10(fc) -13.82*log10(htx) -C_Rx ...
     +(44.9-6.55*log10(htx))*log10(d/1000);
EType = upper(Etype);
if EType(1)=='S',  PL = PL -2*(log10(fc/28))^2 -5.4;
  elseif EType(1)=='O'
    PL=PL+(18.33-4.78*log10(fc))*log10(fc)-40.97;
end

author: LyricYang


text:
classdef calc_Omoricross < ZmapVGridFunction
    % CALC_OMORICROSS calculate omori parameters (p, c, k) along a cross section

    properties
        mc_method       McMethods   = McMethods.FixedMc % this function might only know [' Fixed Mc (Mc = Mmin) | Automatic Mc (max curvature) | EMR-method'];
        bootloops       double      = 50
        learningPeriod  duration    = days(100)
        Nmin            double      = 50
        MainShock       ZmapCatalog
        MainShockSelection char {mustBeMember(MainShockSelection,{'Largest','FirstInGlobal','LargestInXsection'})} = 'Largest'
    end

    properties(Constant)
        PlotTag         = 'calc_Omoricross'
        ReturnDetails   = cell2table({ ... VariableNames, VariableDescriptions, VariableUnits
            'p_value',      'p-Value','';...
            'p_value std',  'p-value standard deviation','';...
            'c_value',      'c-value','';...
            'c_value std',  'c-value standard deviation','';...
            'k_value',      'k-value','';...
            'k_value std',  'k-value standard deviation','';...
            'model',        'Chosen fitting model','';...
            'p_value2',     'p-Value2 UNUSED(?)','';...
            'p_value2 std', 'p-value2 standard deviation UNUSED(?)','';...
            'c_value2',     'c-value2  UNUSED(?)','';...
            'c_value2 std', 'c-value2 standard deviation UNUSED(?)','';...
            'k_value2',     'k-value UNUSED(?)','';...
            'k_value2 std', 'k-value standard deviation UNUSED(?)','';...
            'KS_Test H',    'KS-Test (H-value) binary rejection criterion at 95% confidence level','';...
            'KS_Test stat', 'KS-Test statistic for goodness of fit','';...
            'KS_Test P_value','KS-Test p-value','';...
            'RMS',          'RMS value for goodness of fit','';...
            ...'Mc_value','Mc value',''...
            }, 'VariableNames', {'Names','Descriptions','Units'})

        CalcFields      = {'p_value','p_value std', 'c_value','c_value std','k_value','k_value std',...
            'model','p_value2','p_value2 std', 'c_value2','c_value2 std', 'k_value2', 'k_value2 std',...
            'KS_Test H',' KS_Test stat','KS_Test P_value', 'RMS'} % cell array of charstrings, matching into ReturnDetails.Names

        ParameterableProperties = ['nBstSample','Nmin','learningPeriod', 'MainShock']; % array of strings matching into obj.Properties
        References="";
    end

    methods
        function obj=calc_Omoricross(zap, varargin)
            % CALC_OMORICROSS
            % obj = CALC_OMORICROSS() takes catalog, grid, and eventselection from ZmapGlobal.Data
            %
            % obj = CALC_OMORICROSS(ZAP) where ZAP is a ZmapAnalysisPkg

            obj@ZmapVGridFunction(zap, 'p_value');
            obj.MainShock = ZG.maepi.subset(1);
            report_this_filefun();
            unimplemented_error()
            obj.parseParameters(varargin);
            obj.StartProcess();
        end

        function InteractiveSetup(obj)
            % create a dialog that allows user to select parameters neccessary for the calculation

            %% make the interface
            zdlg = ZmapDialog();

            magtype=  ZG.maepi.MagnitudeType(1);
            if isundefined(magtype)
                magtype = "mag";
            end

            zdlg.AddHeader(sprintf('starting %s with an %s %g event', ZG.maepi.Date(1),magtype, ZG.maepi.Magnitude(1)));

            zdlg.AddMcMethodDropdown('mc_method');
            % original list:  [' Fixed Mc (Mc = Mmin) | Automatic Mc (max curvature) | EMR-method']

            zdlg.AddEventSelector('evsel', ZG.GridSelector);

            zdlg.AddGridSpacing('gridOpts',1.00,'km',[],'',1.0,'km');
            zdlg.AddDurationEdit('learningPeriod','Learning Period', obj.learningPeriod, '', @days);

            zdlg.AddCheckbox('useBootstrap',   'Use Bootstrapping',        true,  {'nBstSample'},...
                're takes longer, but provides more accurate results');
            zdlg.AddEdit('nBstSample',         'Number of bootstraps',     obj.bootloops,...
                'Number of bootstraps to determine Mc');
            zdlg.AddEdit('Nmin','Minimum number of events', obj.Nmin);
            zdlg.AddCheckbox('useEventInXsection','Use first largest event in within cross section as mainshock',false,{},'');

            zdlg.Create('Name', 'Omori Parameters [xsec]','WriteToObj', obj,'OkFcn', @obj.doIt);

        end

        function results=Calculate(obj)
            % once the properties have been set, either by the constructor or by interactive_setup
            % get the grid-size interactively and calculate the values in the grid by sorting the
            % seismicity and selecting the appropriate neighbors to each grid point

            if obj.MainShock.Count ~= 1
                error('There must be a single mainshock provided to this function');
            end
            % if fixed magnitude of completeness, request from user
            if obj.mc_method == McMethods.FixedMc
                [~,~,fMcFix] = smart_inputdlg('Fixed Mc input',...
                    struct('prompt','Enter Mc:', 'value', 1.5));
            else
                fMcFix=0;
            end

           [~, mcCalculator] = calc_mc([], obj.mc_method, fBinning, fMcFix);

            obj.gridCalculations(@calculation_function); % Workhorse that calls calculation_function

            if nargout
                results=obj.Result.values;
            end

            % catsave3('calc_Omoricross_orig')

            % View the map
            % view_Omoricross(myvalues, mygrid, 'p-value');

            %% -----

            function out=calculation_function(catalog)
                % calulate values at a single point
                % Grid coordinates

                fMc = mcCalculator(catalog);

                % for some reason this was only associated with radius
                if ~isnan(fMc)
                    catalog = catalog.subset(catalog.Magnitude >= fMc);
                end

                % Number of events per gridnode
                nY=catalog.Count;

                % Calculate the relative rate change, p, c, k, resolution
                if nY >= Nmin  % enough events?
                    nMod = OmoriModel.pck; % Single Omori law
                    [mResult] = calc_Omoriparams(catalog,obj.learningPeriod, timef, obj.bootloops, ZG.maepi, nMod);

                    % Result matrix
                    out = [...
                        mResult.pval1 mResult.pmeanStd1 ...
                        mResult.cval1 mResult.cmeanStd1...
                        mResult.kval1 mResult.kmeanStd1 ...
                        mResult.nMod ... nY fMaxDist...
                        mResult.pval2 mResult.pmeanStd2 ...
                        mResult.cval2 mResult.cmeanStd2...
                        mResult.kval2 mResult.kmeanStd2...
                        mResult.H...
                        mResult.KSSTAT mResult.P mResult.fRMS ...fMc
                        ];
                else
                    out = [NaN NaN NaN NaN NaN NaN NaN ...
                        ...nY fMaxDist
                        NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... fMc
                        ];
                end
            end
        end

        function ModifyGlobals(obj)
        end
    end

    methods(Static)
        function h = AddMenuItem(parent, zapFcn, varargin)
            % create a menu item
            label = 'omori parameters (p-, k-,c-) [xsec]';
            h = uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)XZfun.calc_Omoricross(zapFcn()),...
                varargin{:});
        end

    end
end


text:
function f=comp_iufilterbank_td(c,g,a,Ls,skip,ext)
%COMP_IUFILTERBANK_TD   Synthesis Uniform filterbank by conv2
%   Usage:  f=comp_iufilterbank_td(c,g,a,Ls,skip,ext);
%
%   Input parameters:
%         c    : N*M*W array of coefficients.
%         g    : Filterbank filters - filtLen*M array.
%         a    : Upsampling factor - scalar.
%         Ls   : Output length.
%         skip : Delay of the filters - scalar or array of length M.
%         ext  : Border exension technique.
%
%   Output parameters:
%         f  : Output Ls*W array.
%

%input channel number
W=size(c,3);
%filter number
M=size(g,2);
%length of filters
filtLen = size(g,1);
% Allow filter delay only in the filter support range
if(all(skip>=filtLen) || all(skip<0))
  error('%s: The filter zero index position outside of the filter support.', upper(mfilename));
end

if(numel(skip)==1)
    skip = skip*ones(M,1);
end

% Output memory allocation
f=zeros(Ls,W,assert_classname(c,g));

if(~strcmp(ext,'per'))
    ext = 'zero';
end

skipOut = a*(filtLen-1)+skip;

% W channels are done simultaneously
for m=1:M
   cext = comp_extBoundary(squeeze(c(:,m,:)),filtLen-1,ext,'dim',1);
   ftmp = conv2(g(:,m),comp_ups(cext,a));
   f = f + ftmp(1+skipOut(m):Ls+skipOut(m),:);
end


text: any(cellfun(@exist,{'lsl_create_outlet', 'lsl_loadlib_'})==3)

author: PatternRecognition

repo: OpenBMI


text:
%%*****************************************************************
%% NTscaling: Compute NT scaling matrix
%%
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'.
%% compute W satisfying W*Z*W = X.
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%     important to keep W symmertic.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

 function [W,G,sv,gamx,gamz,dd,ee,ff] = ...
          NTscaling(blk,X,Z,Zchol,invZchol);

    numblk = size(blk,1);
    W = cell(numblk,1); G = cell(numblk,1); sv = cell(numblk,1);
    gamx = cell(numblk,1); gamz = cell(numblk,1);
    dd = cell(numblk,1); ee = cell(numblk,1); ff = cell(numblk,1);
%%
    for p = 1:size(blk,1)
        pblk = blk(p,:);
        numblk = length(pblk{2});
        n = sum(pblk{2});
        if strcmp(pblk{1},'l')
           dd{p} = X{p}./Z{p}; %% do not add perturbation, it badly affects cre-a
        elseif strcmp(pblk{1},'q');
           gamx{p} = sqrt(qops(pblk,X{p},X{p},2));
           gamz{p} = sqrt(qops(pblk,Z{p},Z{p},2));
           w2 = gamz{p}./gamx{p};  w = sqrt(w2);
           dd{p} = qops(pblk,1./w2,ones(n,1),4);
           tt = qops(pblk,1./w,Z{p},3) - qops(pblk,w,X{p},4);
           gamtt = sqrt(qops(pblk,tt,tt,2));
           ff{p} = qops(pblk,1./gamtt,tt,3);
           ee{p} = qops(pblk,sqrt(2)./w,ff{p},4);
        elseif strcmp(pblk{1},'s')
           tmp = Prod2(pblk,Zchol{p},X{p},0);
           tmp = Prod2(pblk,tmp,Zchol{p}',1);
           [sv2,V] = blkeig(pblk,tmp);
           sv2 = max(1e-20,sv2);
           sv{p} = sqrt(sv2);
           tmp  = Prod2(pblk,invZchol{p},V);
           G{p} = Prod2(pblk,spdiags(sqrt(sv{p}),0,n,n),tmp');
           W{p} = Prod2(pblk,G{p}',G{p},1);
        end
    end
%%**********************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function aRatio = computePerClassAspectRatios(posImgDir, posGtDir, pLoad, aRatioType, plot_orientation_dist)

if nargin<5
  plot_orientation_dist = 0;
end

% JMBUENA!!
% Compute the average window size for every subclass of the
% positive metaclass. We are going to keep this average window size
% within the detector classifier. We will train the detector with
% fixed big window (selecting the features in this window) but we
% finally return the "per class" best fit window.
tid  = ticStatus('Per class aspect ratio',1,30);
fs={posImgDir,posGtDir};
fs=bbGt('getFiles',fs); nImg=size(fs,2); assert(nImg>0);
gt   = cell(nImg,1);
lbls = cell(nImg,1);
pLoad2=pLoad;
%if iscell(pLoad)
%  pLoad2 = pLoad{:};
%end

use_mean_aratio = 0;
if iscell(pLoad)
  index = find(strcmp(pLoad2, 'squarify'));
  if ~isempty(index)
    pLoad2{index+1} = {};
  end
elseif isstruct(pLoad2)
  if isfield(pLoad2, 'squarify')
    pLoad2.squarify = {};
  end
end

for i=1:nImg
  [objs_,gt_] = bbGt('bbLoad',fs{2,i},pLoad2);
  indices = gt_(:,5)==0;
  if (sum(indices)>0)
    lbls{i} = [objs_(indices).subclass]';
    gt{i}   = gt_(indices,:);
  end
  tocStatus(tid,i/nImg);
end
gt   = cell2mat(gt);
lbls = cell2mat(lbls);
num_pos_classes = length(unique(lbls));
ratios = gt(:,3)./gt(:,4);
aRatio = zeros(num_pos_classes, 1);
for i=1:num_pos_classes
  if (strcmp(aRatioType, 'mean') == 1)
    aRatio(i) = mean(ratios(lbls==i));
  else %(strcmp(aRatioType, 'median') == 1)
    aRatio(i) = median(ratios(lbls==i));
  end
%  aRatio(i) = prctile(ratios(lbls==i),60);
%   a_ratio_class_std(i) = std(ratios(lbls==i));
end

if plot_orientation_dist
  % Plot orientation distributios
  for i=1:num_pos_classes
    figure;
    haxes = axes;
    data_i = ratios(lbls == i);
    [h, x] = hist(data_i,50);
  %  hf = h ./ length(data_i);
    stairs(x, h, 'g-');
    axis([0 4 0 max(h)]);
    mean_i = mean(data_i);
    median_i = median(data_i);
    hold on;
    line([mean_i, mean_i], get(haxes, 'YLim'), 'Color', [1 0 0]);
    line([median_i, median_i], get(haxes, 'YLim'), 'Color', [1 0 1]);
    hold off;
    title(sprintf('orientation %d', i));
  end;
end;

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


function [ status ] = testGraphical( varargin )
%TESTGRAPHICAL Runs the M2T test suite to produce graphical output
%
% This is quite a thin wrapper around testMatlab2tikz to run the test suite to
% produce a PDF side-by-side report.
%
% Its allowed arguments are the same as those of testMatlab2tikz.
%
% Usage:
%
%     status = testGraphical(...) % gives programmatical access to the data
%
%     testGraphical(...); % automatically invokes makeLatexReport afterwards
%
% See also: testMatlab2tikz, testHeadless, makeLatexReport

    [state] = initializeGlobalState();
    finally_restore_state = onCleanup(@() restoreGlobalState(state));

    [status, args] = testMatlab2tikz('actionsToExecute', @actionsToExecute, ...
                                     varargin{:});

    if nargout == 0
        makeLatexReport(status, args.output);
    end
end
% ==============================================================================
function status = actionsToExecute(status, ipp)
    status = execute_plot_stage(status, ipp);

    if status.skip
        return
    end

    status = execute_save_stage(status, ipp);
    status = execute_tikz_stage(status, ipp);
    %status = execute_hash_stage(status, ipp); %cannot work with files in
    %standalone mode!
    status = execute_type_stage(status, ipp);

    if ~status.closeall && ~isempty(status.plotStage.fig_handle)
        close(status.plotStage.fig_handle);
    else
        close all;
    end
end
% ==============================================================================


text:
function [im_sort,OA]=pixel_order(im)
% Assign strict ordering to image pixels.
%
%   INPUT:
%       - IM is a monochromatic or multispectral image.
%
%   OUPUT:
%       - IM_SORT is an array that has the same dimensions as IM. Its
%       element entries correspond to the order of the grey level pixel in
%       that position. Note that the channels of multispectral images are
%       processed speparately.
%       - OA is a number in the range [0,1] and idicates the fraction of
%       unique filter response combinations (ie. order accuracy). If OA=1
%       the ordering is strict. For multichannel images OA is a vector
%       whose elements are the OAs of the corresponding image channels.
%
%   REFERENCES:
%       1. Coltuc D. and Bolon P., 1999, "Strict ordering on discrete images
%       and applications"
%       2. Coltuc D., Bolon P. and Chassery J-M., 2006, "Exact histogram
%       specification", IEEE Transcations on Image Processing
%       15(5):1143-1152
%
%   AUTHOR: Anton Semechko (asemechk@uoguelph.ca)
%   DATE:   Dec.2009

if nargin<1 || nargin>1
    err='Incorrect number of input arguments';
    error(err)
end

% Image dimensions
M=size(im,1);
N=size(im,2);
P=size(im,3);

% Filters
F2=(1/5)*[0 1 0; 1 1 1; 0 1 0];
F3=(1/9)*ones(3,3);
F4=(1/13)*ones(5,5); F4([1 2 4 5 6 10 16 20 21 22 24 25])=0;
F5=(1/21)*ones(5,5); F5([1,5,21,25])=0;
F6=(1/25)*ones(5,5);
F={F2 F3 F4 F5 F6};

% Convolve filters with the image (one channel at a time) and order
%--------------------------------------------------------------------------
im_sort=[];
OA=[];

for i=1:P
    FR=double(im(:,:,1));
    for j=1:5
        FR_j=imfilter(double(im(:,:,1)),F{j});
        FR=cat(3,FR,FR_j);
        clear FR_j
    end
    im(:,:,1)=[]; % free up some memory

    % Rearange the filter responses
    FR=reshape(FR,M*N,6);

    % Number of unique filter responses and ordering accuracy
    n=size(unique(FR,'rows'),1);
    OA=[OA,n/(M*N)];

    % Sort responses lexicographically
    [FR,idx_pos]=sortrows(FR);
    clear FR

    % idx_pos is pixel position sorted according to ascending pixel order
    % now sort ordered pixels according to pixel position (linear index)
    [idx_pos,idx_o]=sort(idx_pos,'ascend');
    clear idx_pos
    idx_o=reshape(idx_o,M,N);

    im_sort=cat(3,im_sort,idx_o);
    clear idx_o
end

author: Sable

repo: mcbench-benchmarks


text:
function kern = gaussianwhiteKernParamInit(kern, isArd, nInd)

% GAUSSIANWHITEKERNPARAMINIT Gaussian white kernel parameter initialisation.
% The gaussian white kernel used here corresponds to the covariance of an
% output function which has been obtained like the output of the
% convolution operation between a smoothing kernel function T, which follows a
% a Gaussian form and white noise with variance s^2_r. It is given as
%
%    k(x_i, x_j) = s^2_r \int T_i(x_i - z)T_i(x_j - z) dz
%	             = s^2_r N(x_i;x_j, L_{i,r}^{-1} + L_{j,r}^{-1})
%
% where L_{i,r} corresponds to the inverse width in each direction of each
% pseudo point.
%
% FORMAT
% DESC  initialises the gaussian white kernel structure with some default
%       parameters.
% RETURN kern : the kernel structure with the default parameters placed in.
% ARG kern  : the kernel structure which requires initialisation.
% ARG isArd : specifies if the kernel is ARD
% ARG nInd  : number of inducing functions
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009.

% KERN

% By default it assumes the kernel is ARD and only have one inducing kernel
switch nargin
    case 1
        kern.isArd = false;
        kern.nIndFunct = 1;
    case 2
        kern.isArd = isArd;
        kern.nIndFunct = 1;
    case 3
        kern.isArd = isArd;
        kern.nIndFunct = nInd;
    otherwise
        error('Number of inputs is incorrect')
end

kern.sigma2Noise = 1;

if kern.isArd
    if kern.nIndFunct == 1
        kern.precisionT = ones(kern.inputDimension,1);
        kern.nParams = kern.inputDimension + 1;
    else
        kern.precisionT = ones(kern.inputDimension,kern.nIndFunct);
        kern.nParams = numel(kern.precisionT) + 1;
    end
else
    if kern.nIndFunct == 1
        kern.precisionT = 1;  % Number of rows should equal 1 if not ARD
        kern.nParams = 2;
    else
        kern.precisionT = ones(1,kern.nIndFunct);  % Number of rows should equal 1 if not ARD
        kern.nParams = kern.nIndFunct + 1;
    end
end
kern.transforms.index =1:kern.nParams;
kern.transforms.type = optimiDefaultConstraint('positive');
kern.isStationary = true;