text:
function script_rfcn_VOC0712_ResNet101_OHEM_rpn()
% script_rfcn_VOC0712_ResNet101_OHEM_rpn()
% RFCN training and testing with OHEM using ResNet101 model and RPN
% proposals
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

clc;
clear mex;
clear is_valid_handle; % to clear init_key
run(fullfile(fileparts(fileparts(mfilename('fullpath'))), 'startup'));
%% -------------------- CONFIG --------------------
opts.caffe_version          = 'caffe_rfcn';
opts.gpu_id                 = auto_select_gpu;
active_caffe_mex(opts.gpu_id, opts.caffe_version);

% model
model                       = Model.ResNet101_for_RFCN_VOC0712_OHEM();
% cache name
opts.cache_name             = 'rfcn_VOC0712_ResNet101_OHEM_rpn_resnet101';
% config
conf                        = rfcn_config_ohem('image_means', model.mean_image);
% train/test data
fprintf('Loading dataset...')
dataset                     = [];
dataset                     = Dataset.voc0712_trainval_sp(dataset, 'train', conf.use_flipped, 'resnet101');
dataset                     = Dataset.voc2007_test_sp(dataset, 'test', false, 'resnet101');
fprintf('Done.\n');

% do validation, or not
opts.do_val                 = true;

%% -------------------- TRAINING --------------------

opts.rfcn_model        = rfcn_train(conf, dataset.imdb_train, dataset.roidb_train, ...
                                'do_val',           opts.do_val, ...
                                'imdb_val',         dataset.imdb_test, ...
                                'roidb_val',        dataset.roidb_test, ...
                                'solver_def_file',  model.solver_def_file, ...
                                'net_file',         model.net_file, ...
                                'cache_name',       opts.cache_name, ...
                                'caffe_version',    opts.caffe_version);
assert(exist(opts.rfcn_model, 'file') ~= 0, 'not found trained model');

%% -------------------- TESTING --------------------
                          rfcn_test(conf, dataset.imdb_test, dataset.roidb_test, ...
                                'net_def_file',     model.test_net_def_file, ...
                                'net_file',         opts.rfcn_model, ...
                                'cache_name',       opts.cache_name,...
                                'ignore_cache',     true);

end

repo: R-FCN


text:
function test_ft_appendtimelock

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_appendtimelock append_common ft_datatype_timelock

% (C) Johanna Zumer

% make some dummy timelock structures
tlock.label = {'1';'2'};
tlock.time = 1:5;
tlock.dimord = 'rpt_chan_time';
tlock.trial = randn(10, 2, 5);

% concat over trials
tlock2 = tlock;
tlock2.trial = randn(8, 2, 5);

cfg = [];
tlockapp = ft_appendtimelock(cfg, tlock, tlock2);
assert(all( size(tlockapp.trial) == [18 2 5]))

cfg = [];
cfg.channel = {'1'};
tlockapp = ft_appendtimelock(cfg, tlock, tlock2);
assert(all( size(tlockapp.trial) == [18 1 5]))

% concat over channels
clear tlock2
tlock1.label = tlock.label(1);
tlock2.label = tlock.label(2);
tlock1.trial = tlock.trial(:, 1, :);
tlock2.trial = tlock.trial(:, 2, :);
tlock1.time = tlock.time;
tlock2.time = tlock.time;
tlock1.dimord = tlock.dimord;
tlock2.dimord = tlock.dimord;

cfg = [];
tlockapp = ft_appendtimelock(cfg, tlock, tlock2);
assert(all( size(tlockapp.trial) == [20 1 5]))

cfg = [];
tlockapp = ft_appendtimelock(cfg, tlock1, tlock2);
assert(all( size(tlockapp.trial) == [10 2 5]))

%% Test using other combinations of tlock from real data
clear tlock*
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/eventrelatedaveraging/dataFC_LP.mat'));

cfg = [];
tlock1 = ft_timelockanalysis(cfg, dataFC_LP);

cfg = [];
cfg.channel = {'MLT*'};
tlockl = ft_timelockanalysis(cfg, dataFC_LP);
cfg = [];
cfg.channel = {'MRT*'};
tlockr = ft_timelockanalysis(cfg, dataFC_LP);

% should concat over channels
cfg = [];
tlockapp = ft_appendtimelock(cfg, tlockl, tlockr);

try
  cfg = [];
  cfg.appenddim = 'rpt';
  tlockapp = ft_appendtimelock(cfg, tlockl, tlockr);
  catchflag = 0;
catch
  catchflag = 1;
end
assert(logical(catchflag))

% should use only subset of 'l' channels and concat over avg
cfg = [];
tlockapp = ft_appendtimelock(cfg, tlockl, tlock1);

try
  cfg = [];
  cfg.appenddim = 'chan';
  tlockapp = ft_appendtimelock(cfg, tlockl, tlock1);
  catchflag = 0;
catch
  catchflag = 1;
end
assert(logical(catchflag))

% new tlock with keeptrials
cfg = [];
cfg.keeptrials = 'yes';
tlock1 = ft_timelockanalysis(cfg, dataFC_LP);
cfg.channel = {'MLT*'};
tlockl = ft_timelockanalysis(cfg, dataFC_LP);
cfg.channel = {'MRT*'};
tlockr = ft_timelockanalysis(cfg, dataFC_LP);

% should concat over channels
cfg = [];
tlockapp = ft_appendtimelock(cfg, tlockl, tlockr);

try
  cfg = [];
  cfg.appenddim = 'rpt';
  tlockapp = ft_appendtimelock(cfg, tlockl, tlockr);
  catchflag = 0;
catch
  catchflag = 1;
end
assert(logical(catchflag))

% should use only subset of 'l' channels and concat over avg
cfg = [];
tlockapp = ft_appendtimelock(cfg, tlockl, tlock1);

% now test for numerical inaccurracies, should concatenate across 'rpt'
tlock2 = tlock1;
tlock2.time = tlock1.time+0.0000001;
cfg = [];
tlockapp = ft_appendtimelock(cfg, tlock1, tlock2);

%% test for data with labels shuffled around

tlock2 = tlock1;
tlock2.label = tlock1.label(randperm(length(tlock1.label)));

tlock2.trial = randn(size(tlock1.trial));

cfg = [];
tlockshuffled = ft_appendtimelock(cfg, tlock1, tlock2);

% now check whether channels were correctly appended
[a, b] = match_str(tlock1.label, tlockshuffled.label);
x1 = tlock1.trial(:, a(1), :);
y1 = tlockshuffled.trial(1:size(x1, 1), b(1), :);

[a, b] = match_str(tlock2.label, tlockshuffled.label);
x2 = tlock2.trial(:, a(1), :);
y2 = tlockshuffled.trial([size(x1, 1)+1]:end, b(1), :);

if ~all(x1(:) == y1(:)) || ~all(x2(:) == y2(:))
  error('data was wrongly appended when channel labels are differently ordered in input arguments');
end

%% test with sampleinfo

% make some dummy timelock structures without sampleinfo
tlock1 = [];
tlock1.label = {'1';'2'};
tlock1.time = 1:5;
tlock1.dimord = 'rpt_chan_time';
tlock1.trial = randn(10, 2, 5);

tlock2 = [];
tlock2.label = {'1';'2'};
tlock2.time = 1:5;
tlock2.dimord = 'rpt_chan_time';
tlock2.trial = randn(8, 2, 5);

cfg = [];
tlockapp = ft_appendtimelock(cfg, tlock1, tlock2);
assert(~isfield(tlockapp, 'sampleinfo'));
assert(~isfield(tlockapp, 'fsample'));

cfg.keepsampleinfo = 'no';
tlockapp = ft_appendtimelock(cfg, tlock1, tlock2);
assert(~isfield(tlockapp, 'sampleinfo'));
assert(~isfield(tlockapp, 'fsample'));

cfg.keepsampleinfo = 'ifmakessense';
tlockapp = ft_appendtimelock(cfg, tlock1, tlock2);
assert(~isfield(tlockapp, 'sampleinfo'));
assert(~isfield(tlockapp, 'fsample'));

cfg.keepsampleinfo = 'yes';
tlockapp = ft_appendtimelock(cfg, tlock1, tlock2);
assert(isfield(tlockapp, 'sampleinfo')); % here it should exist
assert(~isfield(tlockapp, 'fsample'));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cplx_str,cplx_inc,cplx_lib,cplx_libname,cplx_ver] = opti_FindCplex()
%Finds IBM ILOG CPLEX Libraries and returns path and version

%Known CPLEX path locations (Modify to suit your system by adding to cell arrays, or create a new structure for other versions)
CPLX1260.x64 = {'C:\Program Files\IBM\ILOG\CPLEX_Studio126\cplex\'};
CPLX1260.x86 = {'C:\Program Files (x86)\IBM\ILOG\CPLEX_Studio126\cplex\'};
CPLX1260.ver = '12.6';
CPLX1260.lib = 'cplex1260';
CPLX1260.vs  = 'VS2012';

CPLX1251.x64 = {'C:\Program Files\IBM\ILOG\CPLEX_Enterprise_Server1251\CPLEX_Studio\cplex\'};
CPLX1251.x86 = {'C:\Program Files (x86)\IBM\ILOG\CPLEX_Enterprise_Server1251\CPLEX_Studio\cplex\'};
CPLX1251.ver = '12.5.1';
CPLX1251.lib = 'cplex1251';
CPLX1251.vs  = 'VS2012';

CPLX125.x64 = {'C:\Program Files\IBM\ILOG\CPLEX_Studio125\cplex\'};
CPLX125.x86 = {'C:\Program Files (x86)\IBM\ILOG\CPLEX_Studio125\cplex\'};
CPLX125.ver = '12.5';
CPLX125.lib = 'cplex125';
CPLX125.vs  = 'VS2010';

CPLX124.x64 = {'C:\Program Files\IBM\ILOG\CPLEX_Studio_Academic124\cplex\'};
CPLX124.x86 = {'C:\Program Files (x86)\IBM\ILOG\CPLEX_Studio_Academic124\cplex\'};
CPLX124.ver = '12.4';
CPLX124.lib = 'cplex124';
CPLX124.vs  = 'VS2010';

%Add any new structures to below
CPLXLIB = {CPLX1260,CPLX1251,CPLX125,CPLX124};

%DO NOT MODIFY BELOW HERE
for i = 1:length(CPLXLIB)
    %Check for CPLEX
    [ok,cplx_str,cplx_inc,cplx_lib,cplx_libname,cplx_ver] = checkCPLXVer(CPLXLIB{i});
    if(ok), return; end
end

%Still no Luck - user must locate
error('Could not find the IBM ILOG CPLEX Location on your computer. Please modify this file to locate it');

function [ok,cplx_str,cplx_inc,cplx_lib,cplx_libname,cplx_ver] = checkCPLXVer(cplxstr)
%Local check function
ok = false;
%Defaults
cplx_str = []; cplx_inc = []; cplx_lib = []; cplx_libname = []; cplx_ver = [];
%Find local config paths
switch(computer)
    case 'PCWIN'
        if(~isfield(cplxstr,'x86') || isempty(cplxstr.x86))
            error('No x86 CPLEX directories specified to check');
        end
        paths = cplxstr.x86;
    case 'PCWIN64'
        if(~isfield(cplxstr,'x64') || isempty(cplxstr.x86))
            error('No x64 CPLEX directories specified to check');
        end
        paths = cplxstr.x64;
    otherwise
        error('This function is only setup for Windows PCs');
end
for i = 1:length(paths)
    if(exist(paths{i},'dir'))
        dir_cplx = paths{i};
        %Get Include Directory
        cplx_inc = [dir_cplx 'include\ilcplex'];
        if(~exist(cplx_inc,'dir'))
            error('Could not find Cplex include directory. Checked:\n %s',cplx_inc);
        end
        %Get Library Directory
        switch(computer)
            case 'PCWIN'
                switch(cplxstr.vs)
                    case 'VS2010'
                        cplx_lib = [dir_cplx 'lib\x86_windows_vs2010\stat_mda\'];
                    case 'VS2012'
                        cplx_lib = [dir_cplx 'lib\x86_windows_vs2012\stat_mda\'];
                end
            case 'PCWIN64'
                switch(cplxstr.vs)
                    case 'VS2010'
                        cplx_lib = [dir_cplx 'lib\x64_windows_vs2010\stat_mda\'];
                    case 'VS2012'
                        cplx_lib = [dir_cplx 'lib\x64_windows_vs2012\stat_mda\'];
                end
        end
        if(~exist(cplx_lib,'dir'))
            error('Could not find Cplex lib directory. Checked:\n %s\n',cplx_lib);
        end
        %Build complete linker line
        cplx_str = [' -I"' cplx_inc '" -L"' cplx_lib '" -l' cplxstr.lib ' '];
        cplx_libname = cplxstr.lib;
        cplx_ver = cplxstr.ver;
        ok = true;
        return;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ f, g, H ] = opt02_fgh ( x, flag )

%*****************************************************************************80
%
%% OPT02_FGH evaluates F, G and H for test case #2.
%
%  Discussion:
%
%    This example is discussed in Dennis and Schnabel, pages 120 and 138.
%
%    A suggested initial value for X is
%
%      X init = ( 1, 1 ).
%
%    The optimizing value is
%
%      X* = ( 0, 0 )
%
%    and the optimal function value is
%
%      F(X*) = 0.
%
%  Modified:
%
%    02 January 2008
%
%  Author:
%
%    Jeff Borggaard,
%    Gene Cliff,
%    Virginia Tech.
%
%  Reference:
%
%    John Dennis, Robert Schnabel,
%    Numerical Methods for Unconstrained Optimization
%    and Nonlinear Equations,
%    SIAM, 1996,
%    ISBN13: 978-0-898713-64-0,
%    LC: QA402.5.D44.
%
%  Parameters:
%
%    Input, real X(2), the evaluation point.
%
%    Input, string FLAG, indicates what must be computed.
%    'f' means only the value of F is needed,
%    'g' means only the value of G is needed,
%    'all' means F, G and H (if appropriate) are needed.
%    It is acceptable to behave as though FLAG was 'all'
%    on every call.
%
%    Output, real F, the optimization function.
%
%    Output, real G(2,1), the gradient column vector.
%
%    Output, real H(2,2), the Hessian matrix.
%
  n = length ( x );

  if ( n ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'OPT02_FGH - Fatal error!\n' );
    fprintf ( '  The input vector X should have length 2.\n'),
    fprintf ( '  Instead, it has length = %d.\n', n );
    keyboard
  end

  f    = x(1)^4 + x(1)^2 + x(2)^2;

  g    = [ 4*x(1)^3 + 2*x(1); 2*x(2) ];

  H    = [ 12*x(1) + 2 ,   0    ;
                0      ,   2    ];

  return

author: johannesgerer

repo: jburkardt-m


text:
% KATAMA_SSH
% grab time series of water levels from Chen's GOM3 NECOFS forecast model archive
% for points north and south of Katama Bay, using NCTOOLBOX
%
% Quick NCTOOLBOX Install:
% 1. Grab and unzip: https://github.com/acrosby/nctoolbox_recent/zipball/master
% 2. Run "setup_nctoolbox"

% OPeNDAP URL for Chen's Archive of NECOFS GOM3 Forecast
url='http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/archives/necofs_gom3';
nc=ncgeodataset(url);
var='zeta';
zvar=nc.geovariable(var);
start=[2011 6 1 0 0 0];
stop=[2011 9 15 0 0 0];

tdat=zvar.timewindowij(start,stop);
% picked these two spots (north, and south of Katama Bay)
loni=[ -70.4893  -70.4873];
lati=[  41.3973   41.3339];
disp(['Reading time series data from ' url '...'])
dat=nj_tseries(nc,'zeta',loni,lati,'method','nearest',...
      'itime',tdat.index,'ele',1);

save katama_ssh.mat dat
% plot whole time series of water levels
figure(1);
set(gcf,'pos',[40 100 900 400]);
plot(dat.time,dat.vals);datetick
legend('north','south');
title('Water levels from NECOFS GOM3')
grid;set(gcf,'color','white');

% zoom into Irene
figure(2)
set(gcf,'pos',[40 150 900 400]);
ii=date_index(dat.time,[2011 8 23 0 0 0],[2011 9 3 0 0 0]);
plot(dat.time(ii),dat.vals(ii,:));datetick
legend('north','south');
title('Water levels from NECOFS GOM3')
grid;set(gcf,'color','white');

%% plot the model bathymetry and locations
figure(3);
lon = nc.data(zvar.getlonname);
lat = nc.data(zvar.getlatname);
depth=nc{'h'}(:);
gvar=zvar.attribute('mesh'); % find mesh variable
grd=nc{gvar}(:); % get the mesh (connectivity array)
[m,n]=size(grd);
if m==3,grd=grd.';elseif n~=3;disp('Error:triangles only');return;end
%zeta=nc{var}(tdat.index(1),:);
%%
trisurf(grd,lon,lat,zeros(size(depth)),-depth);view(2);...
  shading interp;colorbar;dasp(lat(1));
axis([-70.55 -70.40 41.30  41.44]);
caxis([-20 0])
line(loni,lati,'marker','x','linestyle','none','color','black','markersize',14)
title('Katama Bay area bathymetry in NECOFS GOM3 forecast model (m)')
set(gcf,'color','white');
set(gca,'tickdir','out');

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function [y,deriv] = vectorized_function(w,f,m,direction)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% This template vectorizes the given function F: R^m -> R as follows:
%   k = length(w)/m;
%   If direction=1, X = reshape(w,m,k), y(j) = F(X(:,j)), or
%   if direction=2, X = reshape(w,k,m), y(i) = F(X(i,:)),
%   so that length(y) = k.
%
% Input parameters:
%   w: As with every MV2DF, w can be [], a vector, or a function handle to
%      another MV2DF.
%   f: is a function handle to an m-file that represents the function
%      F: R^m -> R, as well as its first and second derivatives.
%
%   m: The input dimension to F.
%      (optional, default m = 1)
%
%   direction: is used as explained above to determine whether columns,
%              or rows of X are processed by F.
%              (optional, default  direction = 2)
%
%   Function f works as follows:
%     (Note that f, f1 and f2 have to know the required direction, it is
%      not passed to them.)
%   [y,f1] = f(X), where X and y are as defined above.
%
%   Function f1 works as follows:
%   [J,f2] = f1(), where size(J) = size(X).
%                  Column/row i of J is the gradient of y(i) w.r.t.
%                  column/row i of W.
%                  f2 is a function handle to 2nd order derivatives.
%                     If 2nd order derivatives are 0, then f2 should be [].
%
%   Function f2 works as follows:
%   H = f2(dX), where size(dX) = size(X).
%               If direction=1, H(:,j) = H_i * dX(:,j), or
%               if direction=2, H(i,:) = dX(i,:)* H_i, where
%               H_i is Hessian of y(i), w.r.t. colum/row i of X.
%
%

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

if ~exist('m','var')
    m = 1;
    direction = 2;
end

if isempty(w)
    y = @(w)vectorized_function(w,f,m,direction);
    return;
end

if isa(w,'function_handle')
    outer = vectorized_function([],f,m,direction);
    y = compose_mv(outer,w,[]);
    return;
end

if direction==1
    W = reshape(w,m,[]);
elseif direction==2
    W = reshape(w,[],m);
else
    error('illegal direction %i',direction);
end

if nargout==1
    y = f(W);
else
    [y,f1] = f(W);
    deriv = @(dy) deriv_this(dy,f1,direction);
end
y = y(:);

end

function [g,hess,linear] = deriv_this(dy,f1,direction)
if direction==1
    dy = dy(:).';
else
    dy = dy(:);
end
if nargout==1
  J = f1();
  g = reshape(bsxfun(@times,J,dy),[],1);
else
  [J,f2] = f1();
  linear = isempty(f2);
  g = reshape(bsxfun(@times,J,dy),[],1);
  hess = @(d) hess_this(d,f2,J,dy,direction);
end

end

function [h,Jv] = hess_this(dx,f2,J,dy,direction)

dX = reshape(dx,size(J));
if isempty(f2)
    h = [];
else
    h = reshape(bsxfun(@times,dy,f2(dX)),[],1);
end
if nargout>1
    Jv = sum(dX.*J,direction);
    Jv = Jv(:);
end

end

%%%%%%%%%%%%%%%%%%%%  Example function: z = x^2 + y^3 %%%%%%%%%%%%%%%%%%%%

% example function: z = x^2 + y^3
function [z,f1] = x2y3(X,direction)
    if direction==1
        x = X(1,:);
        y = X(2,:);
    else
        x = X(:,1);
        y = X(:,2);
    end
    z = x.^2+y.^3;
    f1 = @() f1_x2y3(x,y,direction);
end

% example function 1st derivative: z = x^2 + y^2
function [J,f2] = f1_x2y3(x,y,direction)
if direction==1
    J = [2*x;3*y.^2];
else
    J = [2*x,3*y.^2];
end
f2 = @(dxy) f2_x2y3(dxy,y,direction);
end

% example function 2nd derivative: z = x^2 + y^2
function H = f2_x2y3(dxy,y,direction)
if direction==1
    H = dxy.*[2*ones(size(y));6*y];
else
    H = dxy.*[2*ones(size(y)),6*y];
end
end

%%%%%%%%%%%%%%%%%%%%  Example function: z = x*y^2 %%%%%%%%%%%%%%%%%%%%

% example function: z = x*y^2
function [z,f1] = xy2(X,direction)
    if direction==1
        x = X(1,:);
        y = X(2,:);
    else
        x = X(:,1);
        y = X(:,2);
    end
    y2 = y.^2;
    z = x.*+y2;
    f1 = @() f1_xy2(x,y,y2,direction);
end

% example function 1st derivative: z = x*y^2
function [J,f2] = f1_xy2(x,y,y2,direction)
if direction==1
    J = [y2;2*x.*y];
else
    J = [y2,2*x.*y];
end
f2 = @(dxy) f2_xy2(dxy,x,y,direction);
end

% example function 2nd derivative: z = x*y^2
function H = f2_xy2(dxy,x,y,direction)
if direction==1
  dx = dxy(1,:);
  dy = dxy(2,:);
  H = [2*y.*dy;2*y.*dx+2*x.*dy];
else
  dx = dxy(:,1);
  dy = dxy(:,2);
  H = [2*y.*dy,2*y.*dx+2*x.*dy];
end
end

function test_this()

k = 5;
m = 2;

dr = 1;
fprintf('Testing x^2+y^2 in direction %i:\n\n',dr);
f = vectorized_function([],@(X)x2y3(X,dr),2,dr);
test_MV2DF(f,randn(k*m,1));

dr = 2;
fprintf('\n\n\n\nTesting x*y^2 in direction %i:\n\n',dr);
f = vectorized_function([],@(X)xy2(X,dr),2,dr);
test_MV2DF(f,randn(k*m,1));

end

author: nesl


text:
function [gETU]= get_gETU(ROIonlyPET,volume,a)

nVoxels = sum(~isnan(ROIonlyPET))
end

author: mvallieres

repo: radiomics


text:
function value = bst_prctile(vector, percentile)
% BST_PRCTILE: Returns the percentile value in vector
%
% USAGE: value = bst_prctile(vector, percentile)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c)2000-2020 University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Martin Cousineau, 2020

% Try to use toolbox function
try
    value = prctile(vector, percentile);
    return;
catch
end

if ~isvector(vector)
    error('Only vectors supported.');
end

% Custom implementation
vector = sort(vector);
rank   = percentile / 100 * (length(vector) + 1);
lowerRank = floor(rank);
upperRank = ceil(rank);
fraction  = rank - lowerRank;

if fraction == 0
    value = vector(rank);
else
    value = fraction * (vector(upperRank) - vector(lowerRank)) + vector(lowerRank);
end


text:
function [ result ] = export_df3( filename, chg )
%EXPORT_DF3 Export volumetric data as a DF3 file.
%   status = export_df3(filename,chg) exports the volumetric data contained
%   the 3 dimensional array chg as a DF3 file. This format is compatible
%   with the POV-Ray ray-tracing package. The data is scaled up by a factor
%   of 1,000,000.

    fid = fopen(filename, 'w');
    if fid==-1
        error(['Error opening ' filename]);
    end
    fwrite(fid, size(chg), 'ushort',0,'ieee-be');
    fwrite(fid, 1e6*chg, 'ulong',0,'ieee-be');
    result = fclose(fid);

end

author: Sable

repo: mcbench-benchmarks


text:
function edges = meshEdges(faces, varargin)
%MESHEDGES Computes array of edge vertex indices from face array.
%
%   EDGES = meshEdges(FACES);
%
%   Example
%     meshEdges
%
%   See also
%     meshes3d, meshEdgeFaces, meshFaceEdges

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-06-28, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

%% Process input arguments

if isstruct(faces) && isfield(faces, 'faces')
    % if input is a mesh structure, extract the 'faces' field
    faces = faces.faces;
elseif nargin > 1
    % if two arguments are given, keep the second one
    faces = varargin{1};
end

if ~iscell(faces)
    %% Process faces given as numeric array
    % all faces have same number of vertices
    nVF = size(faces,2);
    e = nchoosek(1:nVF,2);
    A = sparse(faces(:,e(:,1)),faces(:,e(:,2)),1,max(faces(:)),max(faces(:)));
    [EI,EJ] = find(tril(A+A'));
    edges = [EJ EI];

else
    %% faces are given as a cell array
    % faces may have different number of vertices

    % number of faces
    nFaces  = length(faces);

    % compute the number of edges
    nEdges = 0;
    for i = nFaces
        nEdges = nEdges + length(faces{i});
    end

    % allocate memory
    edges = zeros(nEdges, 2);
    ind = 0;

    % fillup edge array
    for i = 1:nFaces
        % get vertex indices, ensuring horizontal array
        f = faces{i}(:)';
        nVF = length(f);
        edges(ind+1:ind+nVF, :) = [f' f([2:end 1])'];
        ind = ind + nVF;
    end

    % keep only unique edges, and return sorted result
    edges = sortrows(unique(sort(edges, 2), 'rows'));
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function varargout = hexagonalGrid(bounds, origin, size, varargin)
%HEXAGONALGRID Generate hexagonal grid of points in the plane.
%
%   usage
%   PTS = hexagonalGrid(BOUNDS, ORIGIN, SIZE)
%   generate points, lying in the window defined by BOUNDS (=[xmin ymin
%   xmax ymax]), starting from origin with a constant step equal to size.
%   SIZE is constant and is equals to the length of the sides of each
%   hexagon.
%
%   TODO: add possibility to use rotated grid
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/08/2005.
%
size = size(1);
dx = 3*size;
dy = size*sqrt(3);

% consider two square grids with different centers
pts1 = squareGrid(bounds, origin + [0 0],        [dx dy], varargin{:});
pts2 = squareGrid(bounds, origin + [dx/3 0],     [dx dy], varargin{:});
pts3 = squareGrid(bounds, origin + [dx/2 dy/2],  [dx dy], varargin{:});
pts4 = squareGrid(bounds, origin + [-dx/6 dy/2], [dx dy], varargin{:});

% gather points
pts = [pts1;pts2;pts3;pts4];

% eventually compute also edges, clipped by bounds
% TODO : manage generation of edges
if nargout>1
    edges = zeros([0 4]);
    x0 = origin(1);
    y0 = origin(2);

    % find all x coordinate
    x1 = bounds(1) + mod(x0-bounds(1), dx);
    x2 = bounds(3) - mod(bounds(3)-x0, dx);
    lx = (x1:dx:x2)';

    % horizontal edges : first find y's
    y1 = bounds(2) + mod(y0-bounds(2), dy);
    y2 = bounds(4) - mod(bounds(4)-y0, dy);
    ly = (y1:dy:y2)';

    % number of points in each coord, and total number of points
    ny = length(ly);
    nx = length(lx);

    if bounds(1)-x1+dx<size
        disp('intersect bounding box');
    end

    if bounds(3)-x2<size
        disp('intersect 2');
        edges = [edges;repmat(x2, [ny 1]) ly repmat(bounds(3), [ny 1]) ly];
        x2 = x2-dx;
        lx = (x1:dx:x2)';
        nx = length(lx);
    end

    for i=1:length(ly)
        ind = (1:nx)';
        tmpEdges(ind, 1) = lx;
        tmpEdges(ind, 2) = ly(i);
        tmpEdges(ind, 3) = lx+size;
        tmpEdges(ind, 4) = ly(i);
        edges = [edges; tmpEdges];
    end

end

% process output arguments
if nargout>0
    varargout{1} = pts;

    if nargout>1
        varargout{2} = edges;
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% DEMVOWELSFGPLVM3 Model the vowels data with a 2-D FGPLVM using RBF kernel and back constraints, but without PCA initialisation.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'vowels';
experimentNo = 3;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.back = 'mlp';
options.backOptions = mlpOptions;
options.optimiseInitBack = 0;
options.numActive = 200;
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(dataSetName, experimentNo, 'vector')

errors = lvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function geometry_test2069 ( )

%*****************************************************************************80
%
%% TEST2069 tests TRIANGLE_CONTAINS_LINE_PAR_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2007
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  p0 = [ 3.0, 0.0, -7.0 ];
  pd = [ 2.0, 1.0, 5.0 ];
  t = [ ...
    8.0, 4.0, 2.0; ...
    9.0, 0.0, 5.0; ...
    2.0, 1.0, 2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST2069\n' );
  fprintf ( 1, '  TRIANGLE_CONTAINS_LINE_PAR_3D determines whether\n' );
  fprintf ( 1, '    a triangle "contains" a parametric line in 3D.\n' );

  r8mat_transpose_print ( dim_num, 3, t, '  Triangle vertices:' );

  norm = sqrt ( sum ( pd(1:dim_num).^2 ) );
  pd(1:dim_num) = pd(1:dim_num) / norm;

  r8vec_print ( dim_num, p0, '  Parametric base point P0:' );
  r8vec_print ( dim_num, pd, '  Parametric direction PD:' );

  [ inside, pint ] = triangle_contains_line_par_3d ( t, p0, pd );

  if ( inside )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The triangle contains the line.\n' );
    r8vec_print ( dim_num, pint, '  Intersection point:' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The triangle does not contain the line.\n' );
    r8vec_print ( dim_num, pint, '  The intersection point:' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected answer:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    The triangle contains the line, and\n' );
  fprintf ( 1, '    the intersection point is at:\n' );
  fprintf ( 1, '      ( 7, 2, 3 ).\n' );

  return
end


text:
function mpe = find_mpe(engine, evidence)
% FIND_MPE Find the most probable explanation (Viterbi)
% mpe = enter_evidence(engine, evidence, ...)
%
% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
%

obslik = mk_hmm_obs_lik_matrix(engine, evidence);
path = viterbi_path(engine.startprob, engine.transprob, obslik);
bnet = bnet_from_engine(engine);
ns = bnet.node_sizes_slice;
ns(bnet.observed) = 1;
ass = ind2subv(ns, path);
mpe = num2cell(ass');
mpe(bnet.observed,:) = evidence(bnet.observed,:);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function Tim = createTime(Time)

% CREATETIME  Create Tim structure.
%   CREATETIME(Time) returns a structure with the fields
%       .dt             sampling time
%       .firstFrame     first frame to be processed
%       .lastFrame      last frame to be processed
%   from the fields in input structure Time.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

Tim = Time;

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [] = printLabeledList(fid, list, labels, perLine, indent)
% Output the values in list with labels for a specified number of values per line
count = 0;
fprintf(fid, '%s[ ', indent);
for k = 1:length(list)
    fprintf(fid, '%g(%s) ', list(k), labels{k});
    count = count + 1;
    if mod(count, perLine) == 0
        fprintf(fid, '\n%s', indent);
        count = 0;
    end
end
fprintf(']\n');

author: VisLab

repo: EEG-Clean-Tools


text:
% Studies the recovery probability of different algorithms
% at a fixed sparsity level and varying number of signals.

close all;
clear all;
clc;
rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');
% Signal space
N = 256;
% Number of measurements
M = 48;
% Number of signals
Ss = 1:32;
% Sparsity level
K = 16;
num_trials = 100;
num_ss = length(Ss);
success_with_s.somp = zeros(num_ss, 1);
success_with_s.ra_omp = zeros(num_ss, 1);
success_with_s.ra_ormp = zeros(num_ss, 1);
success_with_s.gomp_2_mmv = zeros(num_ss, 1);
success_with_s.gomp_4_mmv = zeros(num_ss, 1);

snr_threshold = 100;

for ns=1:num_ss
    % Current number of signals
    S = Ss(ns);
    num_successes.somp = 0;
    num_successes.ra_omp = 0;
    num_successes.ra_ormp = 0;
    num_successes.gomp_2_mmv = 0;
    num_successes.gomp_4_mmv = 0;
    for nt=1:num_trials
        % Sensing matrix
        Phi = spx.dict.simple.gaussian_dict(M, N);
        % Sparse signal generator
        gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
        % Gaussian distributed non-zero samples
        X = gen.gaussian;
        % Measurement vectors
        Y = Phi * X;

        % Create the solver for simultaneous orthogonal matching pursuit
        solver = spx.pursuit.joint.OrthogonalMatchingPursuit(Phi, K, 2);
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.somp = num_successes.somp + success;

        % Create the solver for rank aware orthogonal matching pursuit
        solver = spx.pursuit.joint.RankAwareOMP(Phi, K);
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.ra_omp = num_successes.ra_omp + success;

        % Create the solver for rank aware order recursive matching pursuit
        solver = spx.pursuit.joint.RankAwareORMP(Phi, K);
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.ra_ormp = num_successes.ra_ormp + success;

        % Create the solver for Generalized OMP MMV with L=2
        solver = spx.pursuit.joint.GOMP(Phi, K, 2);
        solver.L = 2;
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.gomp_2_mmv = num_successes.gomp_2_mmv + success;

        % Create the solver for Generalized OMP MMV with L=4
        solver = spx.pursuit.joint.GOMP(Phi, K, 2);
        solver.L = 4;
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.gomp_4_mmv = num_successes.gomp_4_mmv + success;

        fprintf('S: %d, K=%d, trial=%d\n', S, K, nt);
    end
    success_with_s.somp(ns) = num_successes.somp / num_trials;
    success_with_s.ra_omp(ns) = num_successes.ra_omp / num_trials;
    success_with_s.ra_ormp(ns) = num_successes.ra_ormp / num_trials;
    success_with_s.gomp_2_mmv(ns) = num_successes.gomp_2_mmv / num_trials;
    success_with_s.gomp_4_mmv(ns) = num_successes.gomp_4_mmv / num_trials;
end

save ('bin/success_with_s_comparison.mat');

author: indigits

repo: sparse-plex


text:
function test_bug2265

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_convert_units ft_prepare_sourcemodel

clear all
close all

load(dccnpath('/home/common/matlab/fieldtrip/template/headmodel/standard_bem.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/template/headmodel/standard_mri.mat'));
elecs = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/template/electrode/standard_1020.elc'));

%%
CMMM = 'cm'; scale=1;
CMMM = 'mm'; scale=10; % prepare_leadfield, or another method does not work with mm (a bugg)

vol = ft_convert_units(vol, CMMM); % Convert the vol to cm, or mm

% construct the dipole grid in the template brain coordinates
% the source units are in cm
% the negative inwardshift means an outward shift of the brain surface for inside/outside detection
cfg = [];
cfg.xgrid = scale*[-20:1:20];
cfg.ygrid = scale*[-20:1:20];
cfg.zgrid = scale*[-20:1:20];
cfg.unit = CMMM;
cfg.tight = 'yes';
cfg.inwardshift = -scale*0;
cfg.headmodel = vol;
template_grid = ft_prepare_sourcemodel(cfg);

assert(strcmp(vol.unit, CMMM));
assert(strcmp(template_grid.unit, CMMM));

%% make a figure with the template head model and dipole grid
figure
hold on
ft_plot_headmodel(vol, 'facecolor', 'cortex', 'edgecolor', 'none');alpha 0.5;
camlight;
ft_plot_mesh(template_grid.pos(template_grid.inside,:));

elecs = ft_convert_units(elecs, CMMM);

[vol, elecs] = ft_prepare_vol_sens(vol, elecs);

ft_plot_mesh(elecs.elecpos,'vertexcolor',[1 .3 .3]);

data = [];
cfg.sourcemodel = template_grid;
cfg.headmodel = vol;
cfg.elec = elecs;
cfg.unit = CMMM;
[grid] = ft_prepare_leadfield(cfg, data);

grid = ft_convert_units(grid, CMMM); %% does not really functional!!

%% Dipole simulation
% note that beamformer scanning will be done with a 1cm grid, so you should
% not put the dipole on a position that will not be covered by a grid
% location later
cfg = [];
cfg.headmodel = vol;
cfg.elec = elecs;
cfg.dip.pos = scale*[
  3 -3 6 % dipole 1
  ];
cfg.dip.mom = [ 1 0 0 ]';

cfg.dip.frequency = [10]*1;
cfg.dip.phase=pi/4*[1];
cfg.dip.amplitude=[10];
cfg.relnoise = .1;
cfg.ntrials = 20;
data = ft_dipolesimulation(cfg);

%% compute the data covariance matrix, which will capture the activity of
% the simulated dipole
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = [data.time{1}(1) data.time{1}(end)];
timelock = ft_timelockanalysis(cfg, data);

%% source analysis
cfg = [];
METHOD = 'lcmv'; % 'lcmv' 'rv' 'mne'
cfg.method = METHOD;
cfg.sourcemodel = grid;

cfg.headmodel = vol;
cfg.elec = elecs;
cfg.snr = 10;
source = ft_sourceanalysis(cfg, timelock);

%% interpolate onto the anatomical MRI

cfg = [];
cfg.downsample = 2;
if strcmp(METHOD,'lcmv')
  cfg.parameter = 'avg.pow';
elseif strcmp(METHOD,'rv')
  cfg.parameter = 'avg.rv';
else
  cfg.parameter = 'avg.pow';
end
source = ft_sourceinterpolate(cfg, source , mri); % Note: this function crashes if mne method is used

%% plot the source reconstruction

cfg = [];
if strcmp(METHOD,'lcmv')
  try, source=rmfield(source,'time'); end
  cfg.funparameter = 'avg.pow';
elseif strcmp(METHOD,'rv')
  try, source=rmfield(source,'time'); end
  cfg.funparameter = 'avg.rv';
else
  cfg.funparameter = 'avg.pow';
end

cfg.method = 'ortho';
ft_sourceplot(cfg, source);


text:
function output = calloptidsdp(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

model = yalmip2opticsdp(interfacedata);
if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end

if options.savedebug
    save csdpdebug model
end

solvertime = tic;
model.ops = csdpset;
[y,fvals,exitflag,stats,X] = csdp(model.f,model.A,model.b,model.lb,model.ub,model.sdcone,model.y0,model.ops);
solvertime = toc(solvertime);

% Create dual variable in internal format
if options.saveduals
    if K.l > 0
        top = 1;
        D_struc = X{1};
    else
        top = 0;
        D_struc = [];
    end
    if K.s(1) > 0
        for j = 1:length(K.s)
            D_struc = [D_struc;X{j+top}(:)];
        end
    end
else
    D_struc = [];
end

x = y;  % Our notation do not coincide ...
switch exitflag
    case 0
       problem = 0;
    case 1
        problem = 2;
    case 2
        problem = 1;
    case {3,5,6,7,8,9,10}
        problem = 4;
    case {4,-27}
        problem = 3;
    case -50
        problem = 16;
    case 11
        problem = 7;
    otherwise
        problem = -1;
end
infostr = yalmiperror(problem,interfacedata.solver.tag);

if options.savesolveroutput
	solveroutput.y = y;
    solveroutput.fvals = fvals;
    solveroutput.exitflag = exitflag;
    solveroutput.stats = stats;
    solveroutput.X = X;
else
	solveroutput = [];
end

if options.savesolverinput
	solverinput.model = model;
else
	solverinput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = projPointOnPolygon(point, poly, varargin)
%PROJPOINTONPOLYGON  Compute position of a point projected on a polygon.
%
%   POS = projPointOnPolygon(POINT, POLYGON)
%   Compute the position of the orthogonal projection of a point on a
%   polygon.
%   POINT is a 1-by-2 row vector containing point coordinates
%   POLYGON is a N-by-2 array containing coordinates of polygon vertices
%
%   When POINT is an array of points, returns a column vector with as many
%   rows as the number of points.
%
%   [POS, DIST] = projPointOnPolygon(...)
%   Also returns the distance between POINT and POLYGON. The distance is
%   negative if the point is located inside of the polygon.
%
%   Example
%     poly = [10 10; 20 10;20 20;10 20];
%     projPointOnPolygon([15 0], poly)
%     ans =
%         0.5000
%     projPointOnPolygon([0 16], poly)
%     ans =
%         3.4000
%
%   See also
%   points2d, polygons2d, polygonPoint, projPointOnPolyline
%   distancePointpolygon
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2009-04-30, using Matlab 7.7.0.471 (R2008b)
% Copyright 2009-2022 INRA - Cepia Software Platform

% eventually copy first point at the end to ensure closed polygon
if sum(poly(end, :) == poly(1,:)) ~= 2
    poly = [poly; poly(1,:)];
end

% compute position wrt outline
[pos, minDist] = projPointOnPolyline(point, poly);

% process output arguments
if nargout <= 1
    varargout{1} = pos;
elseif nargout == 2
    varargout{1} = pos;
    if inpolygon(point(:,1), point(:,2), poly(:,1), poly(:,2))
        minDist = -minDist;
    end
    varargout{2} = minDist;
end


text:
clear all;
close all;
clc;
spx.cluster.ssc.util.simulate_subspace_preservation(@ssc_omp, 'ssc_omp');
spx.cluster.ssc.util.print_subspace_preservation_results('ssc_omp');


text:
function M = symmetricfactory(n, k)
% Returns a manifold struct to optimize over k symmetric matrices of size n
%
% function M = symmetricfactory(n)
% function M = symmetricfactory(n, k)
%
% Returns M, a structure describing the Euclidean space of n-by-n symmetric
% matrices equipped with the standard Frobenius distance and associated
% trace inner product, as a manifold for Manopt.
% By default, k = 1. If k > 1, points and vectors are stored in 3D matrices
% X of size nxnxk such that each slice X(:, :, i), for i = 1:k, is
% symmetric.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Jan. 22, 2014.
% Contributors:
% Change log:

    if ~exist('k', 'var') || isempty(k)
        k = 1;
    end

    M.name = @() sprintf('(Symmetric matrices of size %d)^%d', n, k);

    M.dim = @() k*n*(n+1)/2;

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d(:), 'fro');

    M.dist = @(x, y) norm(x(:)-y(:), 'fro');

    M.typicaldist = @() sqrt(k)*n;

    M.proj = @(x, d) multisym(d);

    M.egrad2rgrad = M.proj;

    %M.ehess2rhess = @(x, eg, eh, d) eh;

    M.tangent = @(x, d) d;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

	M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.rand = @() multisym(randn(n, n, k));

    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = multisym(randn(n, n, k));
        u = u / norm(u(:), 'fro');
    end

    M.lincomb = @lincomb;
    function v = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
        if nargin == 3
            v = a1*d1;
        elseif nargin == 5
            v = a1*d1 + a2*d2;
        else
            error('Bad usage of euclidean.lincomb');
        end
    end

    M.zerovec = @(x) zeros(n, n, k);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [m, n]);
    M.vecmatareisometries = @() true;

end

author: andrewssobral

repo: lrslibrary


text:
% This file create a simple regression based network for speech
% dereverberation or enhancement
%
% Created by Xiong Xiao, Temasek Laboratories, NTU, Singapore.
% Last Modified: 08 Feb 2017
%
function layer = genNetworkSeparation_Masking(para)
para.freqBin = (0:1/para.fft_len:0.5)*2*pi; % w = 2*pi*f, where f is the normalized frequency k/N is from 0 to 0.5.
nFreqBin = length(para.freqBin);

% Part 1: generate the BF weight predicting subnet

if para.useWav  % input is row waveform
    layer{1}.name = 'Input';
    layer{end}.inputIdx = 1;
    layer{end}.dim = [1 1]*para.nCh;

    layer{end+1}.name = 'stft';
    layer{end}.prev = -length(layer)+1;
    layer{end}.fft_len = para.fft_len;
    layer{end}.frame_len = para.frame_len;
    layer{end}.frame_shift = para.frame_shift;
    layer{end}.removeDC = para.removeDC;
    layer{end}.win_type = para.win_type;
    layer{end}.dim = [(para.fft_len/2+1)*para.nCh layer{length(layer)+layer{end}.prev}.dim(1)];
    layer{end}.skipBP = 1;  % skip backpropagation

    layer{end+1}.name = 'Affine';       % scaling the Fourier transform
    layer{end}.prev = -1;
    layer{end}.W = [];
    layer{end}.b = [];
    layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
    layer{end}.update = 0;
    layer{end}.skipBP = 1;  % skip backpropagation

    % extract a subset of dimensions for prediction. Not used currently, so we
    % will use all the dimensions.
    layer{end+1}.name = 'ExtractDims';
    layer{end}.prev = -1;
    layer{end}.dimIndex = 1:nFreqBin*para.nCh;
    layer{end}.dim = [length(layer{end}.dimIndex) layer{length(layer)+layer{end}.prev}.dim(1)];
    layer{end}.skipBP = 1;

    % get the log power spectrum and perform CMN
    layer{end+1}.name = 'Power';
    layer{end}.prev = -1;
    layer{end}.dim = [1 1] * layer{end-1}.dim(1);
    layer{end}.skipBP = 1;

    layer{end+1}.name = 'Log';
    layer{end}.const = 1e-2;
    layer{end}.prev = -1;
    layer{end}.dim = [1 1]*layer{end-1}.dim(1);
    layer{end}.skipBP = 1;
else    % input is log spectrogram
    layer{1}.name = 'Input';
    layer{end}.inputIdx = 1;
    layer{end}.dim = [1 1]*(para.fft_len/2+1)*para.nCh;
end

if para.useCMN
    layer{end+1}.name = 'CMN';
    layer{end}.prev = -1;
    layer{end}.dim = [1 1]*layer{end-1}.dim(1);
    layer{end}.skipBP = 1;
end

if strcmpi(para.DeltaGenerationType, 'DeltaByEqn')
    outputDim = nFreqBin;
else
    outputDim = nFreqBin*3;
end

switch para.NetType
    case 'DNN'
        layer{end+1}.name = 'Splice';
        layer{end}.prev = -1;
        layer{end}.context = para.contextSize;
        layer{end}.dim = [layer{end}.context 1]*layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.skipBP = 1;

        layer{end+1}.name = 'Affine';       % scaling the Fourier transform
        layer{end}.prev = -1;
        layer{end}.W = [];
        layer{end}.b = [];
        layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.update = 0;

        layerRegression = genNetworkFeedForward_v2(layer{end}.dim(1), para.hiddenLayerSize, outputDim, 'mse', 'sigmoid');
    case 'LSTM'
        layer{end+1}.name = 'Delta';
        layer{end}.prev = -1;
        layer{end}.delta_order = 2;
        layer{end}.dim = [layer{end}.delta_order+1 1]*layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.skipBP = 1;

        layer{end+1}.name = 'Affine';       % scaling the Fourier transform
        layer{end}.prev = -1;
        layer{end}.W = [];
        layer{end}.b = [];
        layer{end}.dim = [1 1] * layer{length(layer)+layer{end}.prev}.dim(1);
        layer{end}.update = 0;

        tmpTopology.inputDim = layer{end}.dim(1);
        tmpTopology.hiddenLayerSizeLSTM = para.hiddenLayerSize;
        tmpTopology.usePastState = zeros(1,length(para.hiddenLayerSize)); % do not use peeping hole
        tmpTopology.hiddenLayerSizeFF = [];
        tmpTopology.outputDim = outputDim;
        tmpTopology.costFn = 'mse';
        tmpTopology.LastActivation4MSE = 'sigmoid';
        layerRegression = genNetworkLSTM(tmpTopology);
end
layer = [layer layerRegression(2:end-2)];

if ~isempty(para.hiddenLayerSizeFF)
    layerFF = genNetworkFeedForward_v2(para.hiddenLayerSize(end), para.hiddenLayerSizeFF, outputDim, 'mse', 'sigmoid');
    layer = [layer(1:end-1) layerFF(2:end-2)];
end

layer = [layer layer(end-1:end)];
layer{end-1}.prev = -3;
mask_idx = [-2 0] + length(layer);

% apply the masking
if para.useWav  % input is row waveform
    power_idx = ReturnLayerIdxByName(layer, 'power');
else
    input_idx = ReturnLayerIdxByName(layer, 'input');
    layer{end+1}.name = 'exp';
    layer{end}.prev = input_idx(1)-length(layer);
    layer{end}.dim = [1 1]*layer{input_idx(1)}.dim(1);
    power_idx = length(layer);
end

MaskingLayer{1}.name = 'hadamard';
MaskingLayer{end}.dim = [1 1]*layer{end}.dim(1);
MaskingLayer{end}.prev = [mask_idx(1) power_idx(1)] - length(layer) -1;

MaskingLayer{end+1}.name = 'Log';
MaskingLayer{end}.const = 1e-2;
MaskingLayer{end}.prev = -1;
MaskingLayer{end}.dim = [1 1]*MaskingLayer{end-1}.dim(1);

if para.useCMN
    MaskingLayer{end+1}.name = 'CMN';
    MaskingLayer{end}.prev = -1;
    MaskingLayer{end}.dim = [1 1]*MaskingLayer{end-1}.dim(1);
end

if strcmpi(para.DeltaGenerationType, 'DeltaByEqn')
    MaskingLayer{end+1}.name = 'Delta';
    MaskingLayer{end}.prev = -1;
    MaskingLayer{end}.delta_order = 2;
    MaskingLayer{end}.dim = [MaskingLayer{end}.delta_order+1 1]*MaskingLayer{end-1}.dim(1);

    MaskingLayer{end+1}.name = 'Affine';       % scaling the contribution of output dimensions to MSE cost function.
    MaskingLayer{end}.prev = -1;
    MaskingLayer{end}.W = [];
    MaskingLayer{end}.b = [];
    MaskingLayer{end}.dim = [1 1] * MaskingLayer{end-1}.dim(1);
    MaskingLayer{end}.update = 0;
end

layer = [layer MaskingLayer];
output_idx(1) = length(layer);

MaskingLayer2 = MaskingLayer;
MaskingLayer2{1}.prev = [mask_idx(2) power_idx(1)] - length(layer) -1;

layer = [layer MaskingLayer2];
output_idx(2) = length(layer);

if para.useWav
    log_idx = ReturnLayerIdxByName(layer, 'log');
    if para.useCMN
        CleanLayer = layer(1:log_idx(1)+1);
    else
        CleanLayer = layer(1:log_idx(1));
    end
    CleanLayer{1}.inputIdx = 2;
    CleanLayer{1}.dim(:) = 1;
    CleanLayer{2}.dim(:) = CleanLayer{2}.dim(:)/para.nCh;
    CleanLayer{3}.dim(:) = CleanLayer{3}.dim(:)/para.nCh;
    CleanLayer(4) = [];
else
    CleanLayer{1}.name = 'Input';
    CleanLayer{end}.inputIdx = 2;
    CleanLayer{end}.dim = [1 1]*(para.fft_len/2+1);
end

CleanLayer{end+1}.name = 'Delta';
CleanLayer{end}.prev = -1;
CleanLayer{end}.delta_order = 2;
CleanLayer{end}.dim = [CleanLayer{end}.delta_order+1 1]*CleanLayer{length(CleanLayer)+CleanLayer{end}.prev}.dim(1);

CleanLayer{end+1}.name = 'Affine';       % scaling the contribution of output dimensions to MSE cost function.
CleanLayer{end}.prev = -1;
CleanLayer{end}.W = [];
CleanLayer{end}.b = [];
CleanLayer{end}.dim = [1 1] * CleanLayer{length(CleanLayer)+CleanLayer{end}.prev}.dim(1);
CleanLayer{end}.update = 0;
for i=1:length(CleanLayer)
    CleanLayer{i}.skipBP = 1;
end

layer = [layer CleanLayer];
clean1_idx = length(layer);
CleanLayer{1}.inputIdx = 3;
layer = [layer CleanLayer];
clean2_idx = length(layer);

layer{end+1}.name = 'mixture_MSE';
layer{end}.prev = [output_idx clean1_idx clean2_idx] - length(layer);
layer{end}.dim = [1 layer{length(layer)+layer{end}.prev(1)}.dim(1)];

layer = FinishLayer(layer);
end

author: singaxiong

repo: SignalGraph


text:
function y = vpa(varargin);
% (1) vpa in Octave produces a lot of unwanted displayed text. It's directly displayed, not part of the
% function  output. I didn't found a way to remove it (I tried changing
% PAGER). But Dennis uses a clever trick (with "tail") to remove that output in the
% online compiler
% (2) (i) After converting to char, the result is different from Matlab. For
% example, char(vpa(pi^pi)) produces the string
%     'Float(''36.462159607207901501624291995540261'', prec=32)''
% instead of just the string
%     '36.46215960720790150162429199554' as in Matlab.
% (ii) Also, precision is lost by using `char`: char(vpa('-3.4')) gives the
% string
%     'Float(''-3.3999999999999999999999999999999988'', prec=32)'
% Both (i) and (ii) seem to be resolved by `pretty` (which actually calls `disp`):
% `pretty(vpa('-3.4'))` gives the string
%     '-3.4000000000000000000000000000000'
% However, `pretty` formats the output string with spaces and newlines. These
% characters are removed if we are sure they are unwanted: real scalars
y = pretty(builtin('vpa', varargin{:}));
if isscalar(varargin{1}) && isreal(double(varargin{1}))
    y = y(y>32); % remove spaces and newlines produced by `pretty`
end
% In Octave `pretty(vpa('765.0908',20))` gives '765.09080000000000000'. In Matlab it gives '765.0908', and 765 gives '765.0'.
% To remove surplus zeros in Octave, the following seems to work. It
% doesn't work for complex values, but the complex case seems to have
% precision loss anyway
if ~ismember('I',y) % complex case not covered
    y = regexprep(y,'(?<=\.\d*)0+(?=(e|$))','0');
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function eps = eps_RCUs( i_s, n, r )
eps = mean( exp( - max(0, i_s - log(2^(n*r)-1) ) ) ); % avg pr err
% eps = mean( exp( - max(0, i_s - log(2^(n*r)-1)-log(2)) ) ); % max pr err
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
% fig2i HungarianCV matrix plot

i_fish = 6;
[cIX1,gIX1] = LoadCluster_Direct(i_fish,4,2);
[cIX2,gIX2] = LoadCluster_Direct(i_fish,5,2);

% [score,im1] = HungarianCV(cIX1,cIX2,gIX1,gIX2,isPlotFig);
[score,im1] = HungarianCV(cIX2,cIX1,gIX2,gIX1);

%%
figure('Position',[500,300,250,250]);
% subplot(1,2,1)
% imagesc(-im1)
% colormap(bluewhitered)
% axis equal; axis tight;axis xy
%
% subplot(1,2,2)
imagesc(-log(im1))
% colormap(bluewhitered)
axis equal; axis tight;%axis xy
colormap('gray')
ylabel('clusters (1st half)')
xlabel('clusters (2nd half)')
title('Cross-Val: # of cells overlap')
text(size(im1,1)/4,size(im1,2)/20,'score = 0.69');

author: xiuyechen

repo: FishExplorer


text:
function value = r8_aint ( x )

%*****************************************************************************80
%
%% R8_AINT truncates an R8 argument to an integer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the truncated version of X.
%
  if ( x < 0.0 )
    value = - floor ( abs ( x ) );
  else
    value =   floor ( abs ( x ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [post nlZ dnlZ] = infFITC_EP(hyp, mean, cov, lik, x, y)

% FITC-EP approximation to the posterior Gaussian process. The function is
% equivalent to infEP with the covariance function:
%   Kt = Q + G; G = diag(g); g = diag(K-Q);  Q = Ku'*inv(Kuu + snu2*eye(nu))*Ku;
% where Ku and Kuu are covariances w.r.t. to inducing inputs xu and
% snu2 = sn2/1e6 is the noise of the inducing inputs. We fixed the standard
% deviation of the inducing inputs snu to be a one per mil of the measurement
% noise's standard deviation sn. In case of a likelihood without noise
% parameter sn2, we simply use snu2 = 1e-6.
% For details, see The Generalized FITC Approximation, Andrew Naish-Guzman and
%                  Sean Holden, NIPS, 2007.
%
% The implementation exploits the Woodbury matrix identity
%   inv(Kt) = inv(G) - inv(G)*Ku'*inv(Kuu+Ku*inv(G)*Ku')*Ku*inv(G)
% in order to be applicable to large datasets. The computational complexity
% is O(n nu^2) where n is the number of data points x and nu the number of
% inducing inputs in xu.
% The posterior N(f|h,Sigma) is given by h = m+mu with mu = nn + P'*gg and
% Sigma = inv(inv(K)+diag(W)) = diag(d) + P'*R0'*R'*R*R0*P. Here, we use the
% site parameters: b,w=$b,\pi$=tnu,ttau, P=$P'$, nn=$\nu$, gg=$\gamma$
%
% The function takes a specified covariance function (see covFunctions.m) and
% likelihood function (see likFunctions.m), and is designed to be used with
% gp.m and in conjunction with covFITC.
%
% The inducing points can be specified through 1) the 2nd covFITC parameter or
% by 2) providing a hyp.xu hyperparameters. Note that 2) has priority over 1).
% In case 2) is provided and derivatives dnlZ are requested, there will also be
% a dnlZ.xu field allowing to optimise w.r.t. to the inducing points xu. However
% the derivatives dnlZ.xu can only be computed for one of the following eight
% covariance functions: cov{Matern|PP|RQ|SE}{iso|ard}.
%
% Copyright (c) by Hannes Nickisch, 2013-10-29.
%
% See also INFMETHODS.M, COVFITC.M.

persistent last_ttau last_tnu              % keep tilde parameters between calls
tol = 1e-4; max_sweep = 20; min_sweep = 2;     % tolerance to stop EP iterations
inf = 'infEP';
cov1 = cov{1}; if isa(cov1, 'function_handle'), cov1 = func2str(cov1); end
if ~strcmp(cov1,'covFITC'); error('Only covFITC supported.'), end    % check cov
if isfield(hyp,'xu'), cov{3} = hyp.xu; end  % hyp.xu is provided, replace cov{3}

[diagK,Kuu,Ku] = feval(cov{:}, hyp.cov, x);         % evaluate covariance matrix
if ~isempty(hyp.lik)                          % hard coded inducing inputs noise
  sn2 = exp(2*hyp.lik(end)); snu2 = 1e-6*sn2;               % similar to infFITC
else
  snu2 = 1e-6;
end
[n, D] = size(x); nu = size(Kuu,1);
m = feval(mean{:}, hyp.mean, x);                      % evaluate the mean vector

rot180   = @(A)   rot90(rot90(A));                     % little helper functions
chol_inv = @(A) rot180(chol(rot180(A))')\eye(nu);                 % chol(inv(A))
R0 = chol_inv(Kuu+snu2*eye(nu));           % initial R, used for refresh O(nu^3)
V = R0*Ku; d0 = diagK-sum(V.*V,1)';    % initial d, needed for refresh O(n*nu^2)

% A note on naming: variables are given short but descriptive names in
% accordance with Rasmussen & Williams "GPs for Machine Learning" (2006): mu
% and s2 are mean and variance, nu and tau are natural parameters. A leading t
% means tilde, a subscript _ni means "not i" (for cavity parameters), or _n
% for a vector of cavity parameters.

% marginal likelihood for ttau = tnu = zeros(n,1); equals n*log(2) for likCum*
nlZ0 = -sum(feval(lik{:}, hyp.lik, y, m, diagK, inf));
if any(size(last_ttau) ~= [n 1])      % find starting point for tilde parameters
  ttau = zeros(n,1);             % initialize to zero if we have no better guess
  tnu  = zeros(n,1);
  [d,P,R,nn,gg] = epfitcRefresh(d0,Ku,R0,V, ttau,tnu); % compute initial repres.
  nlZ = nlZ0;
else
  ttau = last_ttau;                    % try the tilde values from previous call
  tnu  = last_tnu;
  [d,P,R,nn,gg] = epfitcRefresh(d0,Ku,R0,V, ttau,tnu); % compute initial repres.
  nlZ = epfitcZ(d,P,R,nn,gg,ttau,tnu,d0,R0,Ku,y,lik,hyp,m,inf);
  if nlZ > nlZ0                                           % if zero is better ..
    ttau = zeros(n,1);                    % .. then initialize with zero instead
    tnu  = zeros(n,1);
    [d,P,R,nn,gg] = epfitcRefresh(d0,Ku,R0,V, ttau,tnu);       % initial repres.
    nlZ = nlZ0;
  end
end

nlZ_old = Inf; sweep = 0;               % converged, max. sweeps or min. sweeps?
while (abs(nlZ-nlZ_old) > tol && sweep < max_sweep) || sweep<min_sweep
  nlZ_old = nlZ; sweep = sweep+1;
  for i = randperm(n)       % iterate EP updates (in random order) over examples
    pi = P(:,i); t = R*(R0*pi);                            % temporary variables
    sigmai = d(i) + t'*t; mui = nn(i) + pi'*gg;           % post moments O(nu^2)

    tau_ni = 1/sigmai-ttau(i);          %  first find the cavity distribution ..
    nu_ni = mui/sigmai+m(i)*tau_ni-tnu(i);          % .. params tau_ni and nu_ni

    % compute the desired derivatives of the indivdual log partition function
    [lZ, dlZ, d2lZ] = feval(lik{:}, hyp.lik, y(i), nu_ni/tau_ni, 1/tau_ni, inf);
    ttaui =                            -d2lZ  /(1+d2lZ/tau_ni);
    ttaui = max(ttaui,0);     % enforce positivity i.e. lower bound ttau by zero
    tnui  = ( dlZ + (m(i)-nu_ni/tau_ni)*d2lZ )/(1+d2lZ/tau_ni);
    [d,P(:,i),R,nn,gg,ttau,tnu] = ...                    % update representation
              epfitcUpdate(d,P(:,i),R,nn,gg, ttau,tnu,i,ttaui,tnui, m,d0,Ku,R0);
  end
  % recompute since repeated rank-one updates can destroy numerical precision
  [d,P,R,nn,gg] = epfitcRefresh(d0,Ku,R0,V, ttau,tnu);
  [nlZ,nu_n,tau_n] = epfitcZ(d,P,R,nn,gg,ttau,tnu,d0,R0,Ku,y,lik,hyp,m,inf);
end

if sweep == max_sweep
  warning('maximum number of sweeps reached in function infEP')
end
last_ttau = ttau; last_tnu = tnu;                       % remember for next call

post.sW = sqrt(ttau);                  % unused for FITC_EP prediction with gp.m
dd = 1./(d0+1./ttau);
alpha = tnu./(1+d0.*ttau);
RV = R*V; R0tV = R0'*V;
alpha = alpha - (RV'*(RV*alpha)).*dd;     % long alpha vector for ordinary infEP
post.alpha = R0tV*alpha;       % alpha = R0'*V*inv(Kt+diag(1./ttau))*(tnu./ttau)
B = R0tV.*repmat(dd',nu,1); L = B*R0tV'; B = B*RV';
post.L = B*B' - L;                      % L = -R0'*V*inv(Kt+diag(1./ttau))*V'*R0

if nargout>2                                           % do we want derivatives?
  dnlZ = hyp;                                   % allocate space for derivatives
  RVdd = RV.*repmat(dd',nu,1);
  for i=1:length(hyp.cov)
    [ddiagK,dKuu,dKu] = feval(cov{:}, hyp.cov, x, [], i); % eval cov derivatives
    dA = 2*dKu'-R0tV'*dKuu;                                       % dQ = dA*R0tV
    w = sum(dA.*R0tV',2); v = ddiagK-w;   % w = diag(dQ); v = diag(dK)-diag(dQ);
    z = dd'*(v+w) - sum(RVdd.*RVdd,1)*v - sum(sum( (RVdd*dA)'.*(R0tV*RVdd') ));
    dnlZ.cov(i) = (z - alpha'*(alpha.*v) - (alpha'*dA)*(R0tV*alpha))/2;
  end
  for i = 1:numel(hyp.lik)                                   % likelihood hypers
    dlik = feval(lik{:}, hyp.lik, y, nu_n./tau_n, 1./tau_n, inf, i);
    dnlZ.lik(i) = -sum(dlik);
    if i==numel(hyp.lik)
      % since snu2 is a fixed fraction of sn2, there is a covariance-like term
      % in the derivative as well
      v = sum(R0tV.*R0tV,1)';
      z = sum(sum( (RVdd*R0tV').^2 )) - sum(RVdd.*RVdd,1)*v;
      z = z + post.alpha'*post.alpha - alpha'*(v.*alpha);
      dnlZ.lik(i) = dnlZ.lik(i) + snu2*z;
    end
  end
  [junk,dlZ] = feval(lik{:}, hyp.lik, y, nu_n./tau_n, 1./tau_n, inf);% mean hyps
  for i = 1:numel(hyp.mean)
    dm = feval(mean{:}, hyp.mean, x, i);
    dnlZ.mean(i) = -dlZ'*dm;
  end
  if isfield(hyp,'xu')                   % derivatives w.r.t. inducing points xu
    xu = cov{3};
    cov = cov{2};             % get the non FITC part of the covariance function
    Kpu  = cov_deriv_sq_dist(cov,hyp.cov,xu,x);             % d K(xu,x ) / d D^2
    Kpuu = cov_deriv_sq_dist(cov,hyp.cov,xu);               % d K(xu,xu) / d D^2
    if iscell(cov), covstr = cov{1}; else covstr = cov; end
    if ~ischar(covstr), covstr = func2str(covstr); end
    if numel(strfind(covstr,'iso'))>0              % characteristic length scale
      e = 2*exp(-2*hyp.cov(1));
    else
      e = 2*exp(-2*hyp.cov(1:D));
    end
    B = (R0'*R0)*Ku;

    W = ttau;
    t = W./(1+W.*d0);
    diag_dK = alpha.*alpha + sum(RVdd.*RVdd,1)' - t;
    v = diag_dK+t;                 % BdK = B * ( dnlZ/dK - diag(diag(dnlZ/dK)) )
    BdK = (B*alpha)*alpha' - B.*repmat(v',nu,1);
    BdK = BdK + (B*RVdd')*RVdd;
    A = Kpu.*BdK; C = Kpuu.*(BdK*B'); C = diag(sum(C,2)-sum(A,2)) - C;
    dnlZ.xu = A*x*diag(e) + C*xu*diag(e);    % bring in data and inducing points
  end
end

% refresh the representation of the posterior from initial and site parameters
% to prevent possible loss of numerical precision after many epfitcUpdates
% effort is O(n*nu^2) provided that nu<n
function [d,P,R,nn,gg] = epfitcRefresh(d0,P0,R0,R0P0, w,b)
  nu = size(R0,1);                                   % number of inducing points
  rot180   = @(A)   rot90(rot90(A));                   % little helper functions
  chol_inv = @(A) rot180(chol(rot180(A))')\eye(nu);               % chol(inv(A))
  t  = 1./(1+d0.*w);                                   % temporary variable O(n)
  d  = d0.*t;                                                             % O(n)
  P  = repmat(t',nu,1).*P0;                                            % O(n*nu)
  T  = repmat((w.*t)',nu,1).*R0P0;                % temporary variable O(n*nu^2)
  R  = chol_inv(eye(nu)+R0P0*T');                                    % O(n*nu^3)
  nn = d.*b;                                                              % O(n)
  gg = R0'*(R'*(R*(R0P0*(t.*b))));                                     % O(n*nu)

% compute the marginal likelihood approximation
% effort is O(n*nu^2) provided that nu<n
function [nlZ,nu_n,tau_n] = ...
                        epfitcZ(d,P,R,nn,gg,ttau,tnu, d0,R0,P0, y,lik,hyp,m,inf)
  T = (R*R0)*P;                                             % temporary variable
  diag_sigma = d + sum(T.*T,1)'; mu = nn + P'*gg;       % post moments O(n*nu^2)
  tau_n = 1./diag_sigma-ttau;              % compute the log marginal likelihood
  nu_n  = mu./diag_sigma-tnu+m.*tau_n;            % vectors of cavity parameters
  lZ = feval(lik{:}, hyp.lik, y, nu_n./tau_n, 1./tau_n, inf);
  nu = size(gg,1);
  U = (R0*P0)'.*repmat(1./sqrt(d0+1./ttau),1,nu);
  L = chol(eye(nu)+U'*U);
  ld = 2*sum(log(diag(L))) + sum(log(1+d0.*ttau));
  t = T*tnu; tnu_Sigma_tnu = tnu'*(d.*tnu) + t'*t;
  nlZ = ld/2 -sum(lZ) -tnu_Sigma_tnu/2  ...
    -(nu_n-m.*tau_n)'*((ttau./tau_n.*(nu_n-m.*tau_n)-2*tnu)./(ttau+tau_n))/2 ...
      +sum(tnu.^2./(tau_n+ttau))/2-sum(log(1+ttau./tau_n))/2;

% update the representation of the posterior to reflect modification of the site
% parameters by w(i) <- wi and b(i) <- bi
% effort is O(nu^2)
% Pi = P(:,i) is passed instead of P to prevent allocation of a new array
function [d,Pi,R,nn,gg,w,b] = ...
                           epfitcUpdate(d,Pi,R,nn,gg, w,b, i,wi,bi, m,d0,P0,R0)
  dwi = wi-w(i); dbi = bi-b(i);
  hi = nn(i) + m(i) + Pi'*gg;                   % posterior mean of site i O(nu)
  t = 1+dwi*d(i);
  d(i) = d(i)/t;                                                          % O(1)
  nn(i) = d(i)*bi;                                                        % O(1)
  r = 1+d0(i)*w(i);
  r = r*r/dwi + r*d0(i);
  v = R*(R0*P0(:,i));
  r = 1/(r+v'*v);
  if r>0
    R = cholupdate(R,sqrt( r)*R'*v,'-');
  else
    R = cholupdate(R,sqrt(-r)*R'*v,'+');
  end
  gg = gg + ((dbi-dwi*(hi-m(i)))/t)*(R0'*(R'*(R*(R0*Pi))));            % O(nu^2)
  w(i) = wi; b(i) = bi;                            % update site parameters O(1)
  Pi = Pi/t;                                                             % O(nu)

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function et = start(t)
% TIMETIC/START Start a timer.
%
% et = start(t) starts timer t and returns the current elapsed time.  This
% function does not reset the current elapsed time like tic.  The start
% function is designed to be used with pause to continue accumulating time.
%
% Example:
%   t = timertic;
%   start(t);
%   sphere;
%   axis square;
%   et = pause(t)
%
% See also TIMETIC/PAUSE

%
% David Gleich
% Stanford University
% 27 September 2006
%

t.start();
et = t.elapsed();


text:
function determ = pascal2_determinant ( n )

%*****************************************************************************80
%
%% PASCAL2_DETERMINANT returns the determinant of the PASCAL2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  return
end


text:
function tess_transform( SurfaceFiles, transf )
% TESS_TRANSFORM: Apply a transformation (rotation + translation) to a set of surfaces.
%
% USAGE:  tess_transform( SurfaceFiles, transf )
%
% DESCRIPTION:
%     The transformation structure should contain at least 2 fields: transf.R and transf.T
%     They are the rotation and the translation that will be applied to the coordinates
%     of the vertices of each of the surfaces listed in SurfaceFiles.
%
% INPUT:
%     - SurfaceFiles   : cell array of strings, full path to all the surfaces files to process
%     - Transformation : structure, rotation (.R) + translation (.T) to apply

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009

% Parse inputs
if ischar(SurfaceFiles)
    SurfaceFiles = {SurfaceFiles};
end
% Process all the files
for iFile = 1:length(SurfaceFiles)
    % If relative path name, add the protocol path
    if ~file_exist(SurfaceFiles{iFile})
        SurfaceFiles{iFile} = file_fullpath(SurfaceFiles{iFile});
    end
    % Load vertices
    SurfaceMat = in_tess_bst(SurfaceFiles{iFile});
    % Apply transform
    SurfaceMat.Vertices = ([transf.R, transf.T] * [SurfaceMat.Vertices'; ones(1,size(SurfaceMat.Vertices,1))])';
    % Update comment
    SurfaceMat.Comment = [SurfaceMat.Comment, '_SCS'];
    % History: Apply rotation + translation
    SurfaceMat = bst_history('add', SurfaceMat, 'transform', sprintf('Rotation: [%1.3f,%1.3f,%1.3f; %1.3f,%1.3f,%1.3f; %1.3f,%1.3f,%1.3f]', transf.R'));
    SurfaceMat = bst_history('add', SurfaceMat, 'transform', sprintf('Translation: [%1.3f,%1.3f,%1.3f]', transf.T));
    % Save new vertices positions
    bst_save(SurfaceFiles{iFile}, SurfaceMat, 'v7');
    % Unload surface file
    bst_memory('UnloadSurface', SurfaceFiles{iFile});
end

author: brainstorm-tools


text:
function test_tutorial_clusterpermutationfreq(dataset, datadir)

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqanalysis ft_multiplotTFR ft_singleplotTFR ft_freqstatistics ft_topoplotTFR ft_clustTFRplot ft_megplanar ft_combineplanar

if nargin==0
  dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.ds');
  datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/cluster_permutation_freq');
end

%% PREprocessing
% find the interesting segments of data
cfg = [];                                    % empty configuration
cfg.dataset                 = dataset;       % name of CTF dataset
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 2;
cfg.trialdef.eventvalue     = 3;                    % event value of FIC
cfg = ft_definetrial(cfg);

% remove the trials that have artifacts from the trl
cfg.trl([15, 36, 39, 42, 43, 49, 50, 81, 82, 84],:) = [];

% preprocess the data
cfg.channel   = {'MEG', '-MLP31', '-MLO12'};        % read all MEG channels except MLP31 and MLO12
cfg.demean    = 'yes';                              % do baseline correction with the complete trial

dataFIC = ft_preprocessing(cfg);

% find the interesting segments of data
cfg = [];                                           % empty configuration
cfg.dataset                 = dataset;       % name of CTF dataset
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 2;
cfg.trialdef.eventvalue     = 9;                    % event value of FC
cfg = ft_definetrial(cfg);

% remove the trials that have artifacts from the trl
cfg.trl([2, 3, 4, 30, 39, 40, 41, 45, 46, 47, 51, 53, 59, 77, 85],:) = [];

% preprocess the data
cfg.channel   = {'MEG', '-MLP31', '-MLO12'};        % read all MEG channels except MLP31 and MLO12
cfg.demean    = 'yes';                              % do baseline correction with the complete trial

dataFC = ft_preprocessing(cfg);

%% planar gradient and TFA
cfg = [];
cfg.planarmethod = 'sincos';
% prepare_neighbours determines with what sensors the planar gradient is computed
cfg_neighb.method    = 'distance';
cfg.neighbours       = ft_prepare_neighbours(cfg_neighb, dataFC);

dataFIC_planar = ft_megplanar(cfg, dataFIC);
dataFC_planar  = ft_megplanar(cfg, dataFC);

cfg = [];
cfg.output     = 'pow';
cfg.channel    = 'MEG';
cfg.method     = 'mtmconvol';
cfg.taper      = 'hanning';
cfg.foi        = 20;
cfg.toi        = [-1:0.05:2.0];
cfg.t_ftimwin  = 7./cfg.foi; %7 cycles
cfg.keeptrials = 'yes';

freqFC_planar  = ft_freqanalysis(cfg, dataFC_planar);
freqFIC_planar = ft_freqanalysis(cfg, dataFIC_planar);

cfg = [];
freqFIC_planar_cmb = ft_combineplanar(cfg, freqFIC_planar);
freqFC_planar_cmb = ft_combineplanar(cfg, freqFC_planar);

freqFIC_planar_cmb.grad = dataFIC.grad;
freqFC_planar_cmb.grad = dataFC.grad;

%% Permutation test
cfg = [];
cfg.channel          = {'MEG'};
cfg.latency          = 'all';
cfg.frequency        = 20;
cfg.method           = 'montecarlo';
cfg.statistic        = 'ft_statfun_indepsamplesT';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan        = 2;
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.025;
cfg.numrandomization = 500;
% prepare_neighbours determines what sensors may form clusters
cfg_neighb.method    = 'distance';
cfg.neighbours       = ft_prepare_neighbours(cfg_neighb, dataFC);

design = zeros(1,size(freqFIC_planar_cmb.powspctrm,1) + size(freqFC_planar_cmb.powspctrm,1));
design(1,1:size(freqFIC_planar_cmb.powspctrm,1)) = 1;
design(1,(size(freqFIC_planar_cmb.powspctrm,1)+1):(size(freqFIC_planar_cmb.powspctrm,1)+...
  size(freqFC_planar_cmb.powspctrm,1))) = 2;

cfg.design           = design;
cfg.ivar             = 1;

[stat] = ft_freqstatistics(cfg, freqFIC_planar_cmb, freqFC_planar_cmb);

%% Plotting the results
cfg = [];
cfg.jackknife = 'no';
freqFIC_planar_cmb = ft_freqdescriptives(cfg, freqFIC_planar_cmb);
freqFC_planar_cmb  = ft_freqdescriptives(cfg, freqFC_planar_cmb);

stat.raweffect = freqFIC_planar_cmb.powspctrm - freqFC_planar_cmb.powspctrm;
cfg = [];
cfg.alpha  = 0.025;
cfg.parameter = 'raweffect';
cfg.zlim   = [-1e-27 1e-27];
cfg.layout = 'CTF151.lay';
ft_clusterplot(cfg, stat);

%% Prepropcessing and freqanalysis on planar data
cfg = [];
cfg.toilim = [-1.0 0];
dataFIC_baseline = ft_redefinetrial(cfg, dataFIC);

cfg = [];
cfg.toilim = [0.6 1.6];
dataFIC_activation = ft_redefinetrial(cfg, dataFIC);

dataFIC_baseline.time = dataFIC_activation.time;

cfg = [];
cfg.planarmethod = 'sincos';
% prepare_neighbours determines with what sensors the planar gradient is computed
cfg_neighb.method    = 'distance';
cfg.neighbours       = ft_prepare_neighbours(cfg_neighb, dataFIC_activation);
dataFIC_activation_planar = ft_megplanar(cfg, dataFIC_activation);
dataFIC_baseline_planar   = ft_megplanar(cfg, dataFIC_baseline);

cfg = [];
cfg.output = 'pow';
cfg.channel = 'MEG';
cfg.method = 'mtmconvol';
cfg.taper = 'hanning';
cfg.foi = 20;
cfg.toi = [0.6:0.05:1.6];
cfg.t_ftimwin = 7./cfg.foi; %7 cycles
cfg.keeptrials = 'yes';

freqFIC_activation_planar = ft_freqanalysis(cfg, dataFIC_activation_planar);
freqFIC_baseline_planar   = ft_freqanalysis(cfg, dataFIC_baseline_planar);

cfg = [];
freqFIC_baseline_planar_cmb   = ft_combineplanar(cfg,freqFIC_baseline_planar);
freqFIC_activation_planar_cmb = ft_combineplanar(cfg,freqFIC_activation_planar);

freqFIC_baseline_planar_cmb.grad   = dataFIC.grad;
freqFIC_activation_planar_cmb.grad = dataFIC.grad;

%% Permutation test
cfg = [];
cfg.channel          = {'MEG', '-MLP31', '-MLO12'};
cfg.channel          = {'MEG'};
cfg.latency          = [0.8 1.4];
cfg.method           = 'montecarlo';
cfg.frequency        = 20;
cfg.statistic        = 'ft_statfun_actvsblT';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan        = 2;
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.025;
cfg.numrandomization = 500;
% prepare_neighbours determines what sensors may form clusters
cfg_neighb.method    = 'distance';
cfg.neighbours       = ft_prepare_neighbours(cfg_neighb, freqFIC_activation_planar_cmb);

ntrials = size(freqFIC_activation_planar_cmb.powspctrm,1);
design  = zeros(2,2*ntrials);
design(1,1:ntrials) = 1;
design(1,ntrials+1:2*ntrials) = 2;
design(2,1:ntrials) = [1:ntrials];
design(2,ntrials+1:2*ntrials) = [1:ntrials];

cfg.design   = design;
cfg.ivar     = 1;
cfg.uvar     = 2;

[stat] = ft_freqstatistics(cfg, freqFIC_activation_planar_cmb, freqFIC_baseline_planar_cmb);

%% Plotting the results
cfg = [];
cfg.alpha  = 0.05;
cfg.parameter = 'stat';
cfg.zlim   = [-4 4];
cfg.layout = 'CTF151.lay';
ft_clusterplot(cfg, stat);

%% Within subjects experiment

load(fullfile(datadir, 'GA_TFR_orig.mat'));

cfg = [];
cfg.channel          = {'MEG'};
cfg.latency          = [0 1.8];
cfg.frequency        = 20;
cfg.method           = 'montecarlo';
cfg.statistic        = 'ft_statfun_depsamplesT';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan        = 2;
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.025;
cfg.numrandomization = 500;
% specifies with which sensors other sensors can form clusters
cfg_neighb.method    = 'distance';
cfg.neighbours       = ft_prepare_neighbours(cfg_neighb, GA_TFRFC);

subj = 10;
design = zeros(2,2*subj);
for i = 1:subj
  design(1,i) = i;
end
for i = 1:subj
  design(1,subj+i) = i;
end
design(2,1:subj)        = 1;
design(2,subj+1:2*subj) = 2;

cfg.design   = design;
cfg.uvar     = 1;
cfg.ivar     = 2;

failed = true;
for i=1:100 % this can fail sometime, like every second iteration or so, 100 is a *really* conservative number here
  [stat] = ft_freqstatistics(cfg, GA_TFRFIC, GA_TFRFC);

  % Plotting the results
  pcfg = [];
  pcfg.alpha  = 0.025;
  pcfg.parameter = 'stat';
  pcfg.zlim   = [-4 4];
  pcfg.layout = 'CTF151.lay';
  try
    ft_clusterplot(pcfg, stat);
    failed = false;
    break;
  end
end

if failed
  error('ft_clusterplot fails cause p was too high');
end

author: fieldtrip

repo: fieldtrip


text:
classdef tPanelSettings < matlab.unittest.TestCase

    methods (Test, TestTags = {'Git'})
        function tPanelIRFt56(tc)
            opts = BEARsettings(4, 'ExcelFile', 'data.xlsx');

            fcn = @() setProp(opts, 'IRFt',5);
            tc.verifyError(fcn, 'bear:settings:PANELsettings:UnusedIRFt')

            fcn = @() setProp(opts, 'IRFt',6);
            tc.verifyError(fcn, 'bear:settings:PANELsettings:UnusedIRFt')
        end

        function tPanel1IRFt1(tc)
            opts = BEARsettings(4, 'ExcelFile', 'data.xlsx');
            opts.panel = 1;
            tc.verifyEqual(bear.PANELtype(1), opts.panel);

            fcn = @() setProp(opts, 'IRFt',2);
            tc.verifyError(fcn, 'bear:settings:PANELsettings:WrongIRFt')

            fcn = @() setProp(opts, 'IRFt',3);
            tc.verifyError(fcn, 'bear:settings:PANELsettings:WrongIRFt')

            fcn = @() setProp(opts, 'IRFt',4);
            tc.verifyError(fcn, 'bear:settings:PANELsettings:WrongIRFt')

        end

        function tPanel234IRFt234(tc)
            opts = BEARsettings(4, 'ExcelFile', 'data.xlsx');
            opts.panel = 2;
            tc.verifyEqual(bear.PANELtype(2), opts.panel);

            opts.IRFt = 1;
            tc.verifyEqual(bear.IRFtype(1), opts.IRFt);
            opts.IRFt = 2;
            tc.verifyEqual(bear.IRFtype(2), opts.IRFt);
            opts.IRFt = 3;
            tc.verifyEqual(bear.IRFtype(3), opts.IRFt);
            opts.IRFt = 4;
            tc.verifyEqual(bear.IRFtype(4), opts.IRFt);

            opts.panel = 3;
            tc.verifyEqual(bear.PANELtype(3), opts.panel);

            opts.IRFt = 1;
            tc.verifyEqual(bear.IRFtype(1), opts.IRFt);
            opts.IRFt = 2;
            tc.verifyEqual(bear.IRFtype(2), opts.IRFt);
            opts.IRFt = 3;
            tc.verifyEqual(bear.IRFtype(3), opts.IRFt);
            opts.IRFt = 4;
            tc.verifyEqual(bear.IRFtype(4), opts.IRFt);

            opts.panel = 2;
            tc.verifyEqual(bear.PANELtype(2), opts.panel);

            opts.IRFt = 1;
            tc.verifyEqual(bear.IRFtype(1), opts.IRFt);
            opts.IRFt = 2;
            tc.verifyEqual(bear.IRFtype(2), opts.IRFt);
            opts.IRFt = 3;
            tc.verifyEqual(bear.IRFtype(3), opts.IRFt);
            opts.IRFt = 4;
            tc.verifyEqual(bear.IRFtype(4), opts.IRFt);

        end

        function tPanel56(tc)
            opts = BEARsettings(4, 'ExcelFile', 'data.xlsx');
            opts.panel = 5;
            tc.verifyEqual(bear.PANELtype(5), opts.panel);
        end

        function tSwitchPanels(tc)
            opts = BEARsettings(4, 'ExcelFile', 'data.xlsx');
            opts.panel = 1;
            tc.verifyEqual(bear.PANELtype(1), opts.panel);
            tc.verifyEqual(bear.IRFtype(1), opts.IRFt);

            opts.panel = 3;
            tc.verifyEqual(bear.PANELtype(3), opts.panel);
            tc.verifyEqual(bear.IRFtype(1), opts.IRFt);
            opts.IRFt = 2;
            tc.verifyEqual(bear.IRFtype(2), opts.IRFt);

            opts.panel = 5;
            tc.verifyEqual(bear.PANELtype(5), opts.panel);
            tc.verifyEqual(bear.IRFtype(1), opts.IRFt);
        end
    end

end

function s = setProp(s, prop, value)
s.(prop) = value;
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [c, v, n] = ft_connectivity_psi(input, varargin)

% FT_CONNECTIVITY_PSI computes the phase slope index from a data-matrix
% containing the cross-spectral density. It implements the method described
% in: Nolte et al., Robustly estimating the flow direction of information
% in complex physical systems. Physical Review Letters, 2008; 100; 234101.
%
% Use as
%   [c, v, n] = ft_connectivity_psi(input, ...)
%
% The input data input should be organized as
%   Repetitions x Channel x Channel (x Frequency) (x Time)
% or
%   Repetitions x Channelcombination (x Frequency) (x Time)
%
% The first dimension should be singleton if the input already contains an
% average.
%
% Additional optional input arguments come as key-value pairs:
%   nbin			=	scalar, half-bandwidth parameter: the number of frequency bins
%								across which to integrate
%   hasjack		= 0 or 1, specifying whether the repetitions represent
%               leave-one-out samples (allowing for a variance estimate)
%   feedback	= 'none', 'text', 'textbar' type of feedback showing progress of
%               computation
%   dimord		= string, specifying how the input matrix should be interpreted
%   powindx   =
%   normalize =
%
% The output p contains the phase slope index, v is a variance estimate
% which only can be computed if the data contains leave-one-out samples,
% and n is the number of repetitions in the input data. If the phase slope
% index is positive, then the first chan (1st dim) becomes more lagged (or
% less leading) with higher frequency, indicating that it is causally
% driven by the second channel (2nd dim)
%
% See also FT_CONNECTIVITYANALYSIS

% Copyright (C) 2009-2010 Donders Institute, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% FIXME: interpretation of the slope

hasjack   = ft_getopt(varargin, 'hasjack', 0);
feedback  = ft_getopt(varargin, 'feedback', 'none');
dimord    = ft_getopt(varargin, 'dimord');
powindx   = ft_getopt(varargin, 'powindx');
normalize = ft_getopt(varargin, 'normalize', 'no');
nbin      = ft_getopt(varargin, 'nbin');

if isempty(dimord)
  ft_error('input parameters should contain a dimord');
end

if (length(strfind(dimord, 'chan'))~=2 || contains(dimord, 'pos')>0) && ~isempty(powindx)
  %crossterms are not described with chan_chan_therest, but are linearly indexed

  siz = size(input);

  outsum = zeros(siz(2:end));
  outssq = zeros(siz(2:end));
  pvec   = [2 setdiff(1:numel(siz),2)];

  ft_progress('init', feedback, 'computing metric...');
  %first compute coherency and then phaseslopeindex
  for j = 1:siz(1)
    ft_progress(j/siz(1), 'computing metric for replicate %d from %d\n', j, siz(1));
    c      = reshape(input(j,:,:,:,:), siz(2:end));
    p1     = abs(reshape(input(j,powindx(:,1),:,:,:), siz(2:end)));
    p2     = abs(reshape(input(j,powindx(:,2),:,:,:), siz(2:end)));

    p      = ipermute(phaseslope(permute(c./sqrt(p1.*p2), pvec), nbin, normalize), pvec);

    outsum = outsum + p;
    outssq = outssq + p.^2;
  end
  ft_progress('close');

elseif length(strfind(dimord, 'chan'))==2 || length(strfind(dimord, 'pos'))==2
  %crossterms are described by chan_chan_therest

  siz = size(input);

  outsum = zeros(siz(2:end));
  outssq = zeros(siz(2:end));
  pvec   = [3 setdiff(1:numel(siz),3)];

  ft_progress('init', feedback, 'computing metric...');
  for j = 1:siz(1)
    ft_progress(j/siz(1), 'computing metric for replicate %d from %d\n', j, siz(1));
    p1  = zeros([siz(2) 1 siz(4:end)]);
    p2  = zeros([1 siz(3) siz(4:end)]);
    for k = 1:siz(2)
      p1(k,1,:,:,:,:) = input(j,k,k,:,:,:,:);
      p2(1,k,:,:,:,:) = input(j,k,k,:,:,:,:);
    end
    c      = reshape(input(j,:,:,:,:,:,:), siz(2:end));
    p1     = p1(:,ones(1,siz(3)),:,:,:,:);
    p2     = p2(ones(1,siz(2)),:,:,:,:,:);
    p      = ipermute(phaseslope(permute(c./sqrt(p1.*p2), pvec), nbin, normalize), pvec);
    p(isnan(p)) = 0;
    outsum = outsum + p;
    outssq = outssq + p.^2;
  end
  ft_progress('close');

end

n = siz(1);
c = outsum./n;

if n>1
  n = shiftdim(sum(~isnan(input),1),1);
  if hasjack
    bias = (n-1).^2;
  else
    bias = 1;
  end
  v = bias.*(outssq - (outsum.^2)./n)./(n - 1);
else
  v = [];
end

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

function [y] = phaseslope(x, n, norm)

m   = size(x, 1); %total number of frequency bins
y   = zeros(size(x));
x(1:end-1,:,:,:,:) = conj(x(1:end-1,:,:,:,:)).*x(2:end,:,:,:,:);

if strcmp(norm, 'yes')
  coh = zeros(size(x));
  coh(1:end-1,:,:,:,:) = (abs(x(1:end-1,:,:,:,:)) .* abs(x(2:end,:,:,:,:))) + 1;
  %FIXME why the +1? get the coherence
  for k = 1:m
    begindx = max(1,k-n);
    endindx = min(m,k+n);
    y(k,:,:,:,:) = imag(nansum(x(begindx:endindx,:,:,:,:)./coh(begindx:endindx,:,:,:,:),1));
  end
else
  for k = 1:m
    begindx = max(1,k-n);
    endindx = min(m,k+n);
    y(k,:,:,:,:) = imag(nansum(x(begindx:endindx,:,:,:,:),1));
  end
end

author: spm


text:
function s = zSegment(x,y,z1,z2)

% ZSEGMENT   Vertical segment
%   ZSEGMENT(X,Y,Z1,Z2) makes a vertical segment at the plane
%   location (X,Y) from altitude Z1 to altitude Z2.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

s = [x;y;z1;x;y;z2];

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE


text:
function [nodes edges] = createTestGraph02(varargin)
%CREATETESTGRAPH02  Test graph for geodesic functions
%
%   [nodes edges] = createTestGraph02;
%
%   Example
%     [nodes edges] = createTestGraph02;
%     figure;
%     axis([0 100 10 90]);
%     axis equal;
%     drawGraph(nodes, edges);
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

nodes = [ ...
   10 40; ...
   10 60; ...
   20 50; ...
   40 50; ...
   50 20; ...
   50 40; ...
   50 60; ...
   50 80; ...
   60 50; ...
   80 50; ...
   90 40; ...
   90 60];

edges = [...
    1  3; ...
    2  3; ...
    3  4; ...
    4  6; ...
    4  7; ...
    5  6; ...
    6  9; ...
    7  8; ...
    7  9; ...
    9 10; ...
   10 11; ...
   10 12];

author: mattools

repo: matGeom