text:
classdef jsonParseTest < TestCase
    methods
        function this = jsonParseTest(name)
            this = this@TestCase(name);
        end

        function testChar(~)
            assertEqual('', parse('""'));
            assertEqual('a', parse('"a"'));
            assertEqual(sprintf('-\t-\n'), parse('"-\t-\n"'));
            assertEqual('-\t-\n', parse('"-\\t-\\n"'));
            assertEqual(reshape('',[0 1]), parse('[[0,1],""]'));
            assertEqual(reshape('',[3 1 0]), parse('[[3,1,0],""]'));
            assertEqual(['a1';'b2'], parse('[[2,2],"ab12"]'));
        end

        function testLogical(~)
            assertEqual(true, parse('true'));
            assertEqual(false, parse('false'));
            assertEqual(logical.empty, parse('[true,[0,0]]'));
            assertEqual(true(3,1,0), parse('[true,[3,1,0]]'));
            assertEqual(logical([1 0; 1 0]), parse('[true,[2,2],1,1,0,0]'));
        end

        function testDouble(~)
            assertEqual(0, parse('0.0'));
            assertEqual(1, parse('1.0'));
            assertEqual(pi, parse('3.141592653589793'));
            assertEqual(double.empty, parse('[]'));
            assertEqual([3 1 0], parse('[3.0,1.0,0.0]'));
            assertEqual(ones(3,1,0), parse('[[3,1,0]]'));
            assertEqual([-1 -realmax; 0 realmax],...
                parse('[[2,2],-1.0,0.0,-1.7976931348623157E308,1.7976931348623157E308]'));
            assertExceptionThrown(@() parse('NaN'), 'jsonParse:InvalidJSON');
            assertExceptionThrown(@() parse('[[2,2],0,0,NaN,0]'),...
                'jsonParse:InvalidJSON');
        end

        function testSingle(~)
            assertEqual(single(0), parse('[1,[1,1],0.0]'));
            assertEqual(single(pi), parse('[1,[1,1],3.1415927]'));
            assertEqual(single.empty, parse('[1,[0,0]]'));
            assertEqual(ones(3,1,0,'single'), parse('[1,[3,1,0]]'));
            assertEqual(single([-1 -realmax('single'); 0 realmax('single')]),...
                parse('[1,[2,2],-1.0,0.0,-3.4028235E38,3.4028235E38]'));
            assertExceptionThrown(@() parse('[1,[2,2],0,0,NaN,0]'),...
                'jsonParse:InvalidJSON');
        end

        function testInt8(~)
            assertEqual(intmin('int8'), parse('[7,[1,1],-128]'));
            assertEqual(intmax('int8'), parse('[7,[1,1],127]'));
            assertEqual(int8.empty, parse('[7,[0,0]]'));
            assertEqual(ones(3,1,0,'int8'), parse('[7,[3,1,0]]'));
            assertEqual(int8([2 23; -123 0]), parse('[7,[2,2],2,-123,23,0]'));
        end

        function testInt16(~)
            assertEqual(intmin('int16'), parse('[15,[1,1],-32768]'));
            assertEqual(intmax('int16'), parse('[15,[1,1],32767]'));
            assertEqual(int16.empty, parse('[15,[0,0]]'));
            assertEqual(zeros(3,1,0,'int16'), parse('[15,[3,1,0]]'));
            assertEqual(int16([2 23; -234 0]), parse('[15,[2,2],2,-234,23,0]'));
        end

        function testInt32(~)
            assertEqual(intmin('int32'), parse('[31,[1,1],-2147483648]'));
            assertEqual(intmax('int32'), parse('[31,[1,1],2147483647]'));
            assertEqual(int32.empty, parse('[31,[0,0]]'));
            assertEqual(zeros(3,1,0,'int32'), parse('[31,[3,1,0]]'));
            assertEqual(int32([2 23; -234 0]), parse('[31,[2,2],2,-234,23,0]'));
        end

        function testInt64(~)
            assertEqual(intmin('int64'), parse('[63,[1,1],-9223372036854775808]'));
            assertEqual(intmax('int64'), parse('[63,[1,1],9223372036854775807]'));
            assertEqual(int64.empty, parse('[63,[0,0]]'));
            assertEqual(zeros(3,1,0,'int64'), parse('[63,[3,1,0]]'));
            assertEqual(int64([2 23; -234 0]), parse('[63,[2,2],2,-234,23,0]'));
        end

        function testUint8(~)
            assertEqual(intmin('uint8'), parse('[8,[1,1],0]'));
            assertEqual(intmax('uint8'), parse('[8,[1,1],255]'));
            assertEqual(uint8.empty, parse('[8,[0,0]]'));
            assertEqual(zeros(3,1,0,'uint8'), parse('[8,[3,1,0]]'));
            assertEqual(uint8([2 23; 234 0]), parse('[8,[2,2],2,234,23,0]'));
        end

        function testUint16(~)
            assertEqual(intmin('uint16'), parse('[16,[1,1],0]'));
            assertEqual(intmax('uint16'), parse('[16,[1,1],65535]'));
            assertEqual(uint16.empty, parse('[16,[0,0]]'));
            assertEqual(zeros(3,1,0,'uint16'), parse('[16,[3,1,0]]'));
            assertEqual(uint16([2 23; 234 0]), parse('[16,[2,2],2,234,23,0]'));
        end

        function testUint32(~)
            assertEqual(intmin('uint32'), parse('[32,[1,1],0]'));
            assertEqual(intmax('uint32'), parse('[32,[1,1],4294967295]'));
            assertEqual(uint32.empty, parse('[32,[0,0]]'));
            assertEqual(zeros(3,1,0,'uint32'), parse('[32,[3,1,0]]'));
            assertEqual(uint32([2 23; 234 0]), parse('[32,[2,2],2,234,23,0]'));
        end

        function testUint64(~)
            assertEqual(intmin('uint64'), parse('[64,[1,1],0]'));
            assertEqual(intmax('uint64'), parse('[64,[1,1],18446744073709551615]'));
            assertEqual(uint64.empty, parse('[64,[0,0]]'));
            assertEqual(zeros(3,1,0,'uint64'), parse('[64,[3,1,0]]'));
            assertEqual(uint64([2 23; 234 0]), parse('[64,[2,2],2,234,23,0]'));
        end

        function testCellArray(~)
            assertEqual(cell.empty, parse('[[]]'));
            assertEqual(reshape({'jared',0;'jacobs',1}, [2 1 2]),...
                parse('[[2,1,2],["jared","jacobs",0.0,1.0]]'));
            assertEqual({'jane' [] struct('a',0); 432 false -1},...
                parse('[[2,3],["jane",432.0,[],false,{"a":0.0},-1.0]]'));
        end

        function testStruct(~)
            assertEqual(struct.empty, parse('[{},[0,0],[]]'));
            assertEqual(struct('jane', [], 'b', 432),...
                parse('{"b":432.0,"jane":[]}'));
            s = struct;
            s.p = cell(0);
            assertEqual([struct('p', 0) struct('p',true); s struct('p','')],...
                parse('[{},[2,2],["p"],0.0,[[]],true,""]'));
        end

        function testClass(~)
            import edu.stanford.covert.util.SparseMat;
            assertEqual(SparseMat(), parse([
                '["edu.stanford.covert.util.SparseMat",'...
                 '{"siz":[0.0,2.0],"subs":[[0,2]],"vals":[[0,1]]}]']));
            assertEqual(SparseMat([0 0 0; 0 0 0; 0 1 0; 0 0 -1]), parse([
                '["edu.stanford.covert.util.SparseMat",'...
                 '{"siz":[4.0,3.0],"subs":[[2,2],3.0,4.0,2.0,3.0],"vals":[[2,1],1.0,-1.0]}]']));
        end
    end
end

function result = parse(jsonString)
    result = edu.stanford.covert.io.jsonParse(jsonString);
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function test_example_conditional_granger

% MEM 4gb
% WALLTIME 00:10:00

%
%% Conditional Granger causality in the frequency domain
%
% Conditional Granger causality is a derivative of spectral Granger causality that is computed over a triplet of channels (or blocks of channels). It provides the advantage that for this triplet, it allows to differentiate between a delayed parallel drive from sources A to be B and C and a sequential drive from A to B to C.
%
% This example illustrates the simulation and base analysis of the paper
%
% See also: [the connectivity tutorial]/tutorial/connectivity/).
%
%% # Setup and simulating the data sets
%
% First, define parameters under which samples should be simulated.
%
simcfg             = [];
simcfg.ntrials     = 500;
simcfg.triallength = 1;
simcfg.fsample     = 200;
simcfg.nsignal     = 3;
simcfg.method      = 'ar';

% We want to simulate a system with three signals. Their noise is modeled as white noise processes with zero mean and standard deviations
%
%
% They will have the covariances ζ, η and ε. We also require a paramters μ=0.5.
%
% parameters of the model itself
mu                 = 0.5;
absnoise           = [ 1.0   0.2   0.3 ];

% First, we generate the sample for the case of sequential driving.
% We want to incorporate the system
%
%
% which we can do like this:
%
% params(i,j,k): j -> i at t=k
simcfg.params(:,:,1) = [   0      0      0;
                         1.0      0      0;
                           0    1.0     mu];

% Note that the matrix representation for the covariance reads from columns to row, other than the MVAR-model is read intuitively. But we still need to hand the parameters of the noise to the model:
%
% paper defines stds, not cov:
simcfg.noisecov      = diag(absnoise.^2);

data2           = ft_connectivitysimulation(simcfg);

% Now create sample data for the case of differentially delayed driving,
%
%
% which we can write as
%
simcfg.params(:,:,1) = [   0      0      0;
                         1.0      0      0;
                           0      0     mu];
simcfg.params(:,:,2) = [   0      0      0;
                           0      0      0;
                         1.0      0      0];

% We build the actual MVAR-representation...
%
data1           = ft_connectivitysimulation(simcfg);

% #
%
figure
plot(data1.time{1}, data1.trial{1})
legend(data1.label)
xlabel('time (s)')

% Don't be confused that we started with data2 and conclude with data1. This is just to maintain the order the systems have in the paper.
%
%% # MVAR model frequency analysis
% We generate spectral representations from the MVAR representations we defined with data1 and data2. After all, we want to compute spectral Granger causality. Fast Fourier is a good starting point.
%
freq                   = [];
freq.freqcfg           = [];
freq.freqcfg.method    = 'mtmfft';
freq.freqcfg.output    = 'fourier';
freq.freqcfg.tapsmofrq = 2;
freqdata1           = ft_freqanalysis(freq.freqcfg, data1);
freqdata2           = ft_freqanalysis(freq.freqcfg, data2);

%% # "Regular" Granger causality
% Let first compute regular bivariate Granger causality, as this makes the difference clear to what we want.
%
grangercfg = [];
grangercfg.method  = 'granger';
grangercfg.granger.conditional = 'no';
grangercfg.granger.sfmethod = 'bivariate';

gdata = [];
gdata.g1_bivar_reg      = ft_connectivityanalysis(grangercfg, freqdata1);
gdata.g2_bivar_reg      = ft_connectivityanalysis(grangercfg, freqdata2);

%% # Multivariate conditional Granger causality
% However, we clearly want a multivariate approach. Also, we need to define channel combinations, as we now require triplets of inputs.
%
grangercfg.granger.conditional = 'yes';
grangercfg.channelcmb  = {'signal001', 'signal002', 'signal003'};
grangercfg.granger.sfmethod = 'multivariate';
grangercfg.granger.conditional = 'yes';

% block-wise causality
grangercfg.granger.block(1).name   = freqdata1.label{1};
grangercfg.granger.block(1).label  = freqdata1.label(1);
grangercfg.granger.block(2).name   = freqdata1.label{2};
grangercfg.granger.block(2).label  = freqdata1.label(2);
grangercfg.granger.block(3).name   = freqdata1.label{3};
grangercfg.granger.block(3).label  = freqdata1.label(3);

gdata.g1_multi_reg_conditional = ft_connectivityanalysis(grangercfg, freqdata1);
gdata.g2_multi_reg_conditional = ft_connectivityanalysis(grangercfg, freqdata2);

%% # Evaluation
% The label combinations are 6x2 cell arrays, containing all 2-permutations
% tuplets from the channels. How to interpret this?
% Is the combination a, b representing Fa→b|c?
% Let's check this. In scenario 2, we should clearly see a higher causality
% from 1→3 | 2 than in scenario 1 of the differentially delayed
% drive. This corresponds to row 4 in the
% gdata.g1_multi_reg_conditional.labelcmb.
% So, let's compare the labelcmb 1, 3 in both scenarios:
%
scenario1_mean = mean(gdata.g1_multi_reg_conditional.grangerspctrm(4, :));
scenario2_mean = mean(gdata.g2_multi_reg_conditional.grangerspctrm(4, :));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%script v0730
while 1
  if or(x(k)==si,y(k)==1)|(v==2)
    break
  end
  xk1=x(k)+1;  yk1=y(k)-1;
  if I(xk1,yk1)
    break
  end
  xk2=xk1-1;yk2=yk1;
  tr=(I(xk2,yk2)>0)|I(xk1,y(k));
  if tr
    per=per+1;
    pr=pr+sqrt(2);
    if xk1==si
      per=per+1;
    else
      per=per+I(xk1+1,yk1);
    end
    if yk1==1
      per=per+1;
    else
      per=per+I(xk1,yk1-1);
    end
    k=k+1;v=6;
    x=[x xk1];  y=[y yk1];
  else
    break
  end
end

author: Sable

repo: mcbench-benchmarks


text:
function [MPa] = Pa2MPa(Pa)
% Convert pressure from pascals to megapascals.
% Chad Greene 2012
MPa = Pa*0.000001;


text:
function[r]=ellrad(varargin)
%ELLRAD  Average and instantaneous ellipse radius.
%
%   R=ELLRAD(KAPPA,LAMBDA,PHI,STR) where KAPPA and LAMBDA are the amplitude
%   and linearity of a time-varying ellise, and PHI is its time-varying
%   phase, returns various measures of the ellipse 'radius'.
%
%   STR determines the radius quantity to be output.
%
%       STR       Symbol     Description
%       ----------------------------------------------------------------
%      'geo'      RM         Geometric mean radius = SQRT(A*B)
%      'ave'      RBAR       Period-averaged distance from the origin
%      'ins'      RI         Instantaneous distance from the origin
%
%   See Lilly and Gascard (2006) for details.
%
%   If STR is omitted, the geometric mean radius is returned by default.
%
%   RM=ELLRAD(KAPPA,LAMBDA) with PHI omitted also works for computing RM.
%   ____________________________________________________________________
%
%   Cell array input/output
%
%   If ELLRAD is given cell array input, it returns cell array output.
%
%   Thus KAPPA, LAMBDA, and PHI may each be cell arrays of the same size,
%   where each element in the cell array is a numerical array.
%
%   The output variables will be also cell arrays of this size.
%   ____________________________________________________________________
%
%   See also ELLVEL, ELLPARAMS, ELLDIFF.
%
%   Usage:  rm=ellrad(kappa,lambda);
%           rbar=ellrad(kappa,lambda,phi,'average');
%           ri=ellrad(kappa,lambda,phi,'instantaneous');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2005--2014 J.M. Lilly --- type 'help jlab_license' for details

if ischar(varargin{end})
    str=varargin{end}(1:3);
    varargin=varargin(1:end-1);
else
    str='geo';
end

kappa=varargin{1};
lambda=varargin{2};

if length(varargin)==2
    if ~aresame(str,'geo')
           error('Only the geometric mean radius can be computed without PHI input.')
    else
        phi=kappa;  %It doesn't matter what we set this to as it will not be used.
    end
else
    phi=varargin{3};
end

if ~isempty(kappa)
    if ~iscell(kappa)
        r=ellrad_one(kappa,lambda,phi,str);
    else
        r=kappa;
        for i=1:length(kappa)
            r{i}=ellrad_one(kappa{i},lambda{i},phi{i},str);
        end
    end
else
    r=kappa;
end

function[r]=ellrad_one(kappa,lambda,phi,str)
if aresame(str,'geo')
     [a,b]=kl2ab(kappa,lambda);
     r=sqrt(a.*abs(b));   %r=kappa^2*sqrt(1-lambda^2);
elseif aresame(str,'ave')
     ecc=ecconv(lambda,'lin2ecc');
     [K,E]=ellipke(ecc.^2);
     r=frac(2*kappa,pi).*sqrt(1+abs(lambda)).*E;
elseif aresame(str,'ins')
     r=kappa.*sqrt(1+abs(lambda).*cos(2*phi));
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function G = edge_detect(im)
% Here we directly use canny

    [BW,~,G]=edge_canny(rgb2gray(im),'canny');
    BW = imdilate(BW,strel('diamond',5));
    BW = single(BW);
    BW(BW==0) = 0.5;
    BW(1:10,:) = 0;
    BW(end-10:end,:) = 0;
    BW(:,1:10) = 0;
    BW(:,end-10:end) = 0;
    G = G.*BW;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function dataTYPES = CreateNewDataTypes(mrSESSION)
% Create top-level dataTYPES strycture.
%
%  dataTYPES = CreateNewDataTypes(mrSESSION);
%
% We should use the dtGet/Set/Create suite.
%
% We shouldn't set such specific parameters for the default; we should set
% the values empty.  Also, it is awkward to have the same values
% represented in dataTYPES and mrSESSION because the two may lose
% synchronization. (BW)
%
% See also:
%
% Programming Notes
% djh 9/26/2001
% ras, 03/04: added more eventAnalysisParams, since it's assumed they'll
% be created anyway down the line (e.g., if detrend option is 2, mrLoadRet
% looks for the detrendFrames in eventAnalysisParams).
% ras, 02/06: also added the parfile and scanGroup fields for event-
% related analyses, so they're there even if they're not used.

disp('CreateNewDataTypes is obsolete:  Use dtCreate');
evalin('caller','mfilename')

dataTYPES.name = 'Original';

for iScan = 1:length(mrSESSION.functionals)
    dataTYPES.scanParams(iScan).annotation = '';
    dataTYPES.scanParams(iScan).nFrames = mrSESSION.functionals(iScan).nFrames;
    dataTYPES.scanParams(iScan).framePeriod = mrSESSION.functionals(iScan).framePeriod;
    dataTYPES.scanParams(iScan).slices = mrSESSION.functionals(iScan).slices;
    dataTYPES.scanParams(iScan).cropSize = mrSESSION.functionals(iScan).cropSize;
    dataTYPES.scanParams(iScan).PfileName = mrSESSION.functionals(iScan).PfileName;
    dataTYPES.scanParams(iScan).parfile = '';
    dataTYPES.scanParams(iScan).scanGroup = '';
    dataTYPES.blockedAnalysisParams(iScan).blockedAnalysis = 1;
    dataTYPES.blockedAnalysisParams(iScan).detrend = 1;
    dataTYPES.blockedAnalysisParams(iScan).inhomoCorrect = 1;
    dataTYPES.blockedAnalysisParams(iScan).temporalNormalization = 1;
    dataTYPES.blockedAnalysisParams(iScan).nCycles = 8;
    dataTYPES.eventAnalysisParams(iScan) = er_defaultParams;
end

dataTYPES.retinotopyModelParams = [];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [string, file_name] = user_string(string_name, string)
%USER_STRING  Get/set a user specific string
%
% Examples:
%   string  = user_string(string_name)
%   isSaved = user_string(string_name, new_string)
%   [value, file_name] = user_string(...)
%
% Function to get and set a string in a system or user specific file. This
% enables, for example, system specific paths to binaries to be saved.
%
% The specified string will be saved in a file named <string_name>.txt,
% either in a subfolder named .ignore under this file's folder, or in the
% user's prefdir folder (in case this file's folder is non-writable).
%
% IN:
%   string_name - String containing the name of the string required, which
%                 sets the filename storing the string: <string_name>.txt
%   new_string  - The new string to be saved in the <string_name>.txt file
%
% OUT:
%   string  - The currently saved string. Default: ''
%   isSaved - Boolean indicating whether the save was succesful
%   file_name - path of the .txt file used to contain the requested string

% Copyright (C) Oliver Woodford 2011-2014, Yair Altman 2015-

% This method of saving paths avoids changing .m files which might be in a
% version control system. Instead it saves the user dependent paths in
% separate files with a .txt extension, which need not be checked in to
% the version control system. Thank you to Jonas Dorn for suggesting this
% approach.

% 10/01/2013 - Access files in text, not binary mode, as latter can cause
%              errors. Thanks to Christian for pointing this out.
% 29/05/2015 - Save file in prefdir if current folder is non-writable (issue #74)
% 09/01/2018 - Fix issue #232: if the string looks like a file/folder path, ensure it actually exists
% 15/01/2020 - Added file_name output argument

    if ~ischar(string_name)
        error('string_name must be a string.');
    end
    % Create the full filename
    fname = [string_name '.txt'];
    dname = fullfile(fileparts(mfilename('fullpath')), '.ignore');
    file_name = fullfile(dname, fname);
    default_file_name = file_name;
    if nargin > 1
        % Set string
        if ~ischar(string)
            error('new_string must be a string.');
        end
        % Make sure the save directory exists
        %dname = fileparts(file_name);
        if ~exist(dname, 'dir')
            % Create the directory
            try
                if ~mkdir(dname)
                    string = false;
                    return
                end
            catch
                string = false;
                return
            end
            % Make it hidden
            try
                fileattrib(dname, '+h');
            catch
            end
        end
        % Write the file
        fid = fopen(file_name, 'wt');
        if fid == -1
            % file cannot be created/updated - use prefdir if file does not already exist
            % (if file exists but is simply not writable, don't create a duplicate in prefdir)
            if ~exist(file_name,'file')
                file_name = fullfile(prefdir, fname);
                fid = fopen(file_name, 'wt');
            end
            if fid == -1
                string = false;
                file_name = default_file_name;
                return
            end
        end
        try
            fprintf(fid, '%s', string);
        catch
            fclose(fid);
            string = false;
            return
        end
        fclose(fid);
        string = true;
    else
        % Get string
        fid = fopen(file_name, 'rt');
        if fid == -1
            % file cannot be read, try to read the file in prefdir
            file_name = fullfile(prefdir, fname);
            fid = fopen(file_name, 'rt');
            if fid == -1
                string = '';
                file_name = default_file_name;
                return
            end
        end
        string = fgetl(fid);
        fclose(fid);

        % Fix issue #232: if the string looks like a file/folder path, ensure it actually exists
        if ~isempty(string) && any(string=='\' | string=='/') && ~exist(string) %#ok<EXIST>
            string = '';
        end
    end
end

author: altmany


text:
function callbackMFLDA (t, f, xi, gamma, Phi, auxdata)
  fprintf('%3d [%0.3f] \n',t,-f);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%error('Tracker not configured! Please edit the tracker_MFT.m file.'); % Remove this line after proper configuration

% The human readable label for the tracker, used to identify the tracker in reports
% If not set, it will be set to the same value as the identifier.
% It does not have to be unique, but it is best that it is.
tracker_label = ['MFT'];
tracker_repo_path = 'your MFT path'
% For MATLAB implementations we have created a handy function that generates the appropritate
% command that will run the matlab executable and execute the given script that includes your
% tracker implementation.
%
% Please customize the line below by substituting the first argument with the name of the
% script (not the .m file but just the name of the script as you would use it in within Matlab)
% of your tracker and also provide the path (or multiple paths) where the tracker sources
% are found as the elements of the cell array (second argument).
tracker_command = generate_matlab_command('MFT()', {[tracker_repo_path, 'vot2018_main']});

tracker_interpreter = 'matlab';

% tracker_linkpath = {}; % A cell array of custom library directories used by the tracker executable (optional)

repo: MFT


text:
function startHandGuiding( t )
%% About:
% This function is used to start KUKA's off-the-shelf hand guiding

%% Syntax:
% startHandGuiding( t )

%% Areguments:
% t: is the TCP/IP connection object

% Copyright, Mohammad SAFEEA, 9th of June 2017

theCommand='handGuiding';
fprintf(t, theCommand);

fprintf('Hand gudining functionality is started \n');
fprintf('1- Press white button to activate the hand guiding \n');
fprintf('2- Release the white button so the robot stops in its configuration \n');
fprintf('3- Release the white button so the robot stops in its configuration \n');
fprintf('4- To reactivate the handguding functionality, press the green button once, for less than 1 sec  \n');
fprintf('5- Repeat from 1  \n');
fprintf('5- To terminate the hand guding function, press the green button for more than 1.5 sec \n keep pressing until the red light starts to flicker then release your hand, \n the hand guiding mode will be terminated ');
    message='';
    w = warning ('off','all');
    readingFlag=false;
    while readingFlag==false
        message=fgets(t);

        if checkAcknowledgment(message)
            readingFlag=true;
        end
    end
     w = warning ('on','all');
    fprintf('Hand gudining was terminated by user at robot side \n');
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function value = r8_gami ( a, x )

%*****************************************************************************80
%
%% R8_GAMI evaluates the incomplete gamma function for an R8 argument.
%
%  Discussion:
%
%    GAMI = Integral ( 0 <= T <= X ) exp ( - t ) * t^( a - 1 )  dt
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real A, the parameter.
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the value of the incomplete
%    gamma function.
%
  if ( a <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_GAMI - Fatal error!\n' );
    fprintf ( 1, '  A <= 0.\n' );
    error ( 'R8_GAMI - Fatal error!' )
  end

  if ( x < 0.0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_GAMI - Fatal error!\n' );
    fprintf ( 1, '  X < 0.\n' );
    error ( 'R8_GAMI - Fatal error!' )

  elseif ( x == 0.0 )

    value = 0.0;

  else

    factor = exp ( r8_lngam ( a ) + a * log ( x ) );

    value = factor * r8_gamit ( a, x );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cg = demo()
%DEMO   Return a random BVP CHEBGUI demo.
%   CG = CHEBGUI.DEMO() returns a random BVP CHEBGUI demo, on a CHEBGUI object
% format.

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

% Find the folders which demos are stored in. The chebguiDemos folder lives in
% the trunk folder, find the path of the Chebfun trunk.
trunkPath = chebfunroot();

% Append directory information:
bvppath = fullfile(trunkPath, 'chebguiDemos', 'bvpdemos');

% Setup ODEs:
D = dir(bvppath);

% Find how many demos there are available:
numberOfDemos = length(D)-2; % First two entries are . and ..

% Create a random integer in 1:numberOfDemos. We avoid use of randi which was
% introduced in v2008b.
randnum = floor(numberOfDemos*rand) + 1;
selected = randnum + 2; % Need to shift the random integer by 2

% The full path to the demo file:
demoPath = fullfile(bvppath, D(selected).name);

% Convert the .guifile to a CHEBGUI object:
cg = chebgui.demo2chebgui(demoPath);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% IM = mkSquare(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN, TWIDTH)
%      or
% IM = mkSine(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN, TWIDTH)
%
% Compute a matrix of dimension SIZE (a [Y X] 2-vector, or a scalar)
% containing samples of a 2D square wave, with given PERIOD (in
% pixels), DIRECTION (radians, CW from X-axis, default = 0), AMPLITUDE
% (default = 1), and PHASE (radians, relative to ORIGIN, default = 0).
% ORIGIN defaults to the center of the image.  TWIDTH specifies width
% of raised-cosine edges on the bars of the grating (default =
% min(2,period/3)).
%
% In the second form, FREQ is a 2-vector of frequencies (radians/pixel).

% Eero Simoncelli, 6/96.

% TODO: Add duty cycle.  Allow soft threshold (using rcosFn)

function [res] = mkSquare(sz, per_freq, dir_amp, amp_phase, phase_orig, orig_twidth, twidth)

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (prod(size(per_freq)) == 2)
  frequency = norm(per_freq);
  direction = atan2(per_freq(1),per_freq(2));
  if (exist('dir_amp') == 1)
    amplitude = dir_amp;
  else
    amplitude = 1;
  end
  if (exist('amp_phase') == 1)
    phase = amp_phase;
  else
    phase = 0;
  end
  if (exist('phase_orig') == 1)
    origin = phase_orig;
  end
  if (exist('orig_twidth') == 1)
    transition = orig_twidth;
  else
    transition = min(2,2*pi/(3*frequency));
  end
  if (exist('twidth') == 1)
    error('Too many arguments for (second form) of mkSine');
  end
else
  frequency = 2*pi/per_freq;
  if (exist('dir_amp') == 1)
    direction = dir_amp;
  else
    direction = 0;
  end
  if (exist('amp_phase') == 1)
    amplitude = amp_phase;
  else
    amplitude = 1;
  end
  if (exist('phase_orig') == 1)
    phase = phase_orig;
  else
    phase = 0;
  end
  if (exist('orig_twidth') == 1)
    origin = orig_twidth;
  end
  if (exist('twidth') == 1)
    transition = twidth;
  else
    transition = min(2,2*pi/(3*frequency));
  end

end

%------------------------------------------------------------

if (exist('origin') == 1)
  res = sin(mkRamp(sz, direction, frequency, phase, origin));
else
  res = sin(mkRamp(sz, direction, frequency, phase));
end

[Xtbl,Ytbl] = rcosFn(2*sin(frequency*transition/2),0,[-amplitude amplitude]);

res = pointOp(res,Ytbl,Xtbl(1),Xtbl(2)-Xtbl(1),0);

% OLD threshold version:
%res = amplitude * (mod(res,2*pi) < pi);


text:
%% IMAGE_POOL demonstrates MATLAB's SPMD command for parallel programming.
%
%  Discussion:
%
%    The function reads in an image and adds noise to it.
%    It then divides the image up into R, G and B portions and expects
%    3 workers to filter the portions.
%    The function then returns the original, noisy, and filtered versions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2010
%
%  Author:
%
%    John Burkardt
%    Gene Cliff
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'IMAGE_POOL\n' );
  fprintf ( 1, '  Demonstrate the use of MATLAB''s SPMD command\n' );
  fprintf ( 1, '  to carry out a parallel computation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The client reads a noisy RGB image.\n' );
  fprintf ( 1, '  Each worker filters one color.\n' );
  fprintf ( 1, '  The client reassembles the data.\n' );
  fprintf ( 1, '\n' );

  n = 3;

  matlabpool ( 'open', 'local', n );

  [ x, y, z ] = image_fun ( );

  matlabpool close
%
%  Display the original, noisy, and filtered versions.
%
  figure;

  subplot ( 1, 3, 1 );
  imshow ( x );
  title ( 'Original image' );

  subplot ( 1, 3, 2 );
  imshow( y );
  title( 'Noisy Image' );

  subplot ( 1, 3, 3 );
  imshow ( z );
  title ( 'Median Filtered Image' );


text:
function [ mO ] = ImageFilteringFrequencyDomain( mI, mH, paddingMode )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = ImageFilteringFrequencyDomain( mI, mH, paddingMode )
% Applies Image Filtering in the Frequency Domain.
% Input:
%   - mI                -   Input Image.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - mH                -   Filtering Kernel.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paddingMode       -   Padding Mode.
%                           Sets whether the padding is by Zeros,
%                           Symmetric, Replicate or Circular mode.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, 4}.
% Output:
%   - mI                -   Output Image.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'imfilter()' - https://www.mathworks.com/help/images/ref/imfilter.html.
% Remarks:
%   1.  A
% TODO:
%   1.
%   Release Notes:
%   -   1.0.000     04/04/2019  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

PADDING_MODE_ZEROS      = 1;
PADDING_MODE_SYMMETRIC  = 2;
PADDING_MODE_REPLICATE  = 3;
PADDING_MODE_CIRCULAR   = 4;

numRows     = size(mI, 1);
numCols     = size(mI, 2);

numRowsKernel = size(mH, 1);
numColsKernel = size(mH, 2);

radiusRows = floor(numRowsKernel / 2);
radiusCols = floor(numColsKernel / 2);

vPadRadius = [radiusRows; radiusCols];

switch(paddingMode)
    case(PADDING_MODE_ZEROS)
        % Size of the Linear Convolution Support
        numRowsL = numRows + numRowsKernel - 1;
        numColsL = numCols + numColsKernel - 1;

        % Equivalent of Full Linear Convolution (Zero Padding Built In).
        % This is due the fact padding with zeors at the end with symmetric
        % assumption and shifting yield correct result.
        mO = ifft2(fft2(mI, numRowsL, numColsL) .* fft2(mH, numRowsL, numColsL), 'symmetric');

        firstRowIDx = ceil((numRowsKernel + 1) / 2);
        firstColIDx = ceil((numColsKernel + 1) / 2);

        mO = mO(firstRowIDx:(firstRowIDx + numRows - 1), firstColIDx:(firstColIDx + numCols - 1));
    case({PADDING_MODE_SYMMETRIC, PADDING_MODE_REPLICATE})
        % Padding to apply "Linear Convolution" on the image pixels.
        mI = PadArray2D(mI, vPadRadius, paddingMode);

        numRowsPad = numRows + (2 * radiusRows);
        numColsPad = numCols + (2 * radiusCols);

        mHC = mH;
        mHC(numRowsPad, numColsPad) = 0;
        mHC = circshift(mHC, [-radiusRows, -radiusCols]);

        mO = ifft2(fft2(mI) .* fft2(mHC), 'symmetric');
        mO = mO((radiusRows + 1):(radiusRows + numRows), (radiusCols + 1):(radiusCols + numCols));
    case(PADDING_MODE_CIRCULAR)
        mHC = mH;
        mHC(numRows, numCols) = 0;
        mHC = circshift(mHC, [-radiusRows, -radiusCols]);

        mO = ifft2(fft2(mI) .* fft2(mHC), 'symmetric');
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
clear classes

W = diag([0.1, 1*pi/180].^2);
P0 = diag([0.005, 0.005, 0.001].^2);
V = diag([0.005, 0.5*pi/180].^2);

randinit
map = Map(20, 10);
veh = Vehicle(V);
RandomPath(veh, map.dim);
sensor = RangeBearingSensor(veh, map, W);
sensor.interval = 5;
ekf = EKF(veh, [], P0, sensor, W, []);
ekf.verbose = true;

randinit
ekf.run(1000);

f1
clf
map.visualize()
veh.plot_xy('b');
ekf.plot_map();
grid on
xyzlabel

author: petercorke

repo: robotics-toolbox-matlab


text:
function res = PDD_Profile_res(params, doseV_obj, doseV_array,doseV_arrayFF, doseV_e, ...
    doseProfile1_obj, doseP1p5cm_array, doseP1p5cm_arrayFF, doseP1p5cm_e, ...
    doseProfile2_obj, doseP5cm_array, doseP5cm_arrayFF,  doseP5cm_e, ...
    doseProfile3_obj, doseP10cm_array, doseP10cm_arrayFF, doseP10cm_e, ...
    doseProfile4_obj, doseP20cm_array, doseP20cm_arrayFF, doseP20cm_e, ...
    energy, numBin, extraBin, numBinFF)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Start the function.
%Calculate the residue of the PDD and dose profile difference between DPM and the measurement.

doseV_sum = zeros(size(doseV_obj));

ener = energy/numBin: energy/numBin:energy*(1+extraBin/numBin);
a = Fatigue(params(1:4), ener);
a(find(isnan(a))) = 0;

cutoff = 0.85;
Ef = energy*cutoff;
kt = energy*0.15;  % change kt
f = 1./(1+exp((ener-Ef)/kt));
a = a.*f;

% Sum the contribution of primary and "OnlyHorn" effect
doseV_sum = (doseV_array) * real(a)';
doseProfile1_sum = (doseP1p5cm_array) * real(a)';
doseProfile2_sum = (doseP5cm_array) * real(a)';
doseProfile3_sum = (doseP10cm_array)  * real(a)';
doseProfile4_sum = (doseP20cm_array) * real(a)';

% Get the weights/spectrum for the flattening filter
enerFF = energy/numBin : energy/numBin : numBinFF*energy/numBin;
aFF = Fatigue(params(1:4), params(5)* enerFF);
aFF(find(isnan(aFF))) = 0;

% Get modified Flattening filter spectrum/weights
f = 1./(1+exp((params(5)*enerFF-Ef)/kt));
aFF = aFF.* f;
aFF = aFF * (sum(a)/sum(aFF)*params(6)) ;  % scale FF weight respect to sum(a)

%calculate the dose contribution of the flattening filter
doseV_sumFF = doseV_arrayFF * real(aFF)';
doseProfile1_sumFF = doseP1p5cm_arrayFF * real(aFF)';
doseProfile2_sumFF = doseP5cm_arrayFF * real(aFF)';
doseProfile3_sumFF = doseP10cm_arrayFF * real(aFF)';
doseProfile4_sumFF = doseP20cm_arrayFF * real(aFF)';

% Nov. 1 2006
% Add the asd
filterSize = 7;
filterWindow = [-3 -2 -1 0 1 2 3];
G1 = gauss(filterWindow, params(9));  %p(9) is sigma
G1 = G1/sum(G1);   % normalize, to make sum(G1)==1;

doseV_DPM = doseV_sum + doseV_sumFF + params(8) * sum(a)* doseV_e;
doseProfile1_DPM = doseProfile1_sum + doseProfile1_sumFF + params(8) * sum(a)* doseP1p5cm_e;
doseProfile1_DPM = conv(G1, doseProfile1_DPM);
doseProfile2_DPM = doseProfile2_sum + doseProfile2_sumFF + params(8) * sum(a)* doseP5cm_e;
doseProfile2_DPM = conv(G1, doseProfile2_DPM);
doseProfile3_DPM = doseProfile3_sum + doseProfile3_sumFF + params(8) * sum(a)* doseP10cm_e;
doseProfile3_DPM = conv(G1, doseProfile3_DPM);
doseProfile4_DPM = doseProfile4_sum + doseProfile4_sumFF + params(8) * sum(a)* doseP20cm_e;
doseProfile4_DPM = conv(G1, doseProfile4_DPM);

% Use relative weight for PDD as 5; for lateral dose profile, use weight
% as 1.
% for doseV_obj, dose different absolute dose value make a difference?
% Should it be weighed by the dose values at those points?

diffsqr = 1.0*sum ((doseV_obj - doseV_DPM).^2) ...
      + sum((doseProfile1_obj - doseProfile1_DPM(4:end-3)).^2) ...
    + sum((doseProfile2_obj - doseProfile2_DPM(4:end-3)).^2) ...
    + sum((doseProfile3_obj - doseProfile3_DPM(4:end-3)).^2) ...
      + sum((doseProfile4_obj - doseProfile4_DPM(4:end-3)).^2);
res = sqrt(diffsqr);
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ readerobj ] = open_tiff_reader(filename)
%OPEN_TIFF_READER Intiates a reader object for TIFF imagery
%
% Works for most generic TIFFs, but really intended for SAR datasets stored
% in TIFF format.  For some GeoTIFFs (like RADARSAT-2 and Sentinel-1) this
% will also find and process metadata into the SICD format.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Sensor-specific metadata
% A TIFF file by itself is usually light on any SAR metadata.  Here we look
% for sensor-specific metadata files that might be associated with this
% TIFF file and use them to compute a SICD metadata structure.
tiff_sensor = guess_tiff_sensor(filename);
meta2sicd_tiff_fun = ['meta2sicd_tiff' tiff_sensor];
if ~isempty(tiff_sensor) && exist(meta2sicd_tiff_fun,'file')==2 % Sensor-specific
    [meta, symmetry] = feval(meta2sicd_tiff_fun, filename);
else
    meta = struct();
    symmetry = [0 0 0]; % Default orientation if none can be determined
end

%% Open reader
tiff_meta = imfinfo(filename);
if isfield(tiff_meta,'SampleFormat') && ... % Normal MATLAB TIFF reader does not handle complex data types
        any(strncmpi(tiff_meta.SampleFormat,'complex',7))
    symmetry(3) = ~symmetry(3); % Our TIFF reader is transposed from MATLAB's
    if iscell(meta) % TIFF should be split into multiple parts
        % This case really only works for Sentinel1.  We should make a more
        % generic solution that allows the meta2sicd function to define how
        % TIFF is split.
        last_col = 0;
        base_reader = open_ctiff_reader_noxml(filename,symmetry);
        for i = 1:numel(meta)
            readerobj{i} = subset_reader(base_reader,...
                    last_col + [1 meta{i}.ImageData.NumCols],[1 meta{i}.ImageData.NumRows]);
            last_col = last_col + meta{i}.ImageData.NumCols;
        end
    else
        readerobj = open_ctiff_reader_noxml(filename,symmetry);
    end
else
    readerobj = open_tiff_reader_noxml(filename,symmetry);
end
if iscell(meta)
    for i = 1:numel(meta)
        meta{i} = setstructfields(meta{i},readerobj{i}.get_meta()); % Merge metadata
        readerobj{i}.get_meta=@() meta{i};
    end
else
    meta = setstructfields(meta,readerobj.get_meta()); % Merge metadata
    readerobj.get_meta=@() meta;
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [ cuboid_rect_score, segment_rect_score] = computeCuboidImgFea( rotImg, allhyps, checkrect, checkseg, bufName )
%COMPUTECUBOIDIMGFEA Summary of this function goes here
%   Detailed explanation goes here

%% project image to 6 views and compute pixel-wise hog
rotImg  = im2double(rotImg);
vp = [-1  0  0; ...
       1  0  0; ...
       0 -1  0; ...
       0  1  0; ...
       0  0 -1; ...
       0  0  1];
uv = xyz2uvN(vp);
fov = 160/180*pi;
cutSize = 2000;

hasBuff = exist('bufName','var') && exist(bufName,'file') && ~isempty(bufName);
% hasBuff = false;
if hasBuff
    load(bufName);
else
    all_hog = cell(6,1);
    [sepSceneSix] = separatePano( rotImg, fov, uv(:,1), uv(:,2), cutSize);
end
% load './rectangleDetector/finalModel.mat'; % load in allModel
% rectModel = allModel{4}; % after 2nd negative mining
% load './rectangleDetector/selectModel.mat';

% compute response
detconfig;
load(config.modelfile);
partTemplate = reshape(rectModel.belta, 5, 5, 31, 9);

max_scale = config.max_scale;
min_scale = config.min_scale;
for i = 1:6
    if ~hasBuff
        img = double(sepSceneSix(i).img .* 255);
        feature = feature_pyramid( img, max_scale, min_scale );
        all_hog{i} = feature;
    else
        feature = all_hog{i};
    end
    for pid = 1:9
        for sid = 1:length(feature.feat)
            match{i, pid, sid} = fconvblas(feature.feat{sid}, {partTemplate(:,:,:,pid)}, 1, 1);
            match{i, pid, sid} = match{i, pid, sid}{1};
        end
    end
end

%% checkrect
% hyps = allhyps(checkrect);
ptIDofSurf = [8 5 1 4; 6 7 3 2; 5 6 2 1; 7 8 4 3; 1 2 3 4; 8 7 6 5];
PROJECTVIEW = [2 1 4 3 6 5];
cuboid_rect_score = cell(length(allhyps),1);
for hid = find(checkrect')
%     fprintf('HID: %d\n', hid);
    hyps = allhyps{hid};
    score = repmat(struct('surfscore',[],'surfID',[]), length(hyps), 1);
    for oid = 1:length(hyps)
%         fprintf('OID: %d\n', oid);
        p1 = hyps(oid).align(1,:);
        p7 = hyps(oid).align(7,:);
        if hyps(oid).type == 1
            tmpstd = std(hyps(oid).align, 1);
            [~,id] = min(tmpstd);
            priorvisibility = false(6,1);
            priorvisibility([id*2-1 id*2]) = true;
        else
            priorvisibility = true(6,1);
        end
        visibility = find(sum(vp .* [p1;p7;p1;p7;p1;p7], 2)<0 & priorvisibility);
        surfscore = -100*ones(length(visibility),1);
        for i = 1:length(visibility)
            scalescore = -100*ones(length(feature.feat),1);
            sid = PROJECTVIEW(visibility(i));
            pid = ptIDofSurf(visibility(i),:);
            [out2D, valid, ~] = projectPoint2SeparateView( hyps(oid).align(pid,:), vp(sid,:), fov, cutSize );
            minscalex1 = 2*8/(out2D(1,1)-1);
            minscalex2 = 3*8/(cutSize - out2D(3,1));
            minscaley1 = 2*8/(out2D(1,2)-1);
            minscaley2 = 3*8/(cutSize - out2D(3,2));
            valid_scale_ID = find(feature.scale>max([minscalex1 minscalex2 minscaley1 minscaley2]));
            if min([minscalex1 minscalex2 minscaley1 minscaley2])<0
                valid_scale_ID = [];
            end

            for j = valid_scale_ID%1:length(feature.feat)
                [xIDmin, yIDmin] = getGridID(out2D(1,1), out2D(1,2), feature.scale(j));
                [xIDmax, yIDmax] = getGridID(out2D(3,1)+1, out2D(3,2)+1, feature.scale(j));
                xID1 = xIDmin;  yID1 = yIDmin;
                xID2 = round((xIDmin+xIDmax)/2);    yID2 = yIDmin;
                xID3 = xIDmax;    yID3 = yIDmin;
                xID4 = xIDmin;  yID4 = round((yIDmin+yIDmax)/2);
                xID5 = round((xIDmin+xIDmax)/2);    yID5 = round((yIDmin+yIDmax)/2);
                xID6 = xIDmax;    yID6 = round((yIDmin+yIDmax)/2);
                xID7 = xIDmin;  yID7 = yIDmax;
                xID8 = round((xIDmin+xIDmax)/2);    yID8 = yIDmax;
                xID9 = xIDmax;    yID9 = yIDmax;

                scalescore(j) = match{sid,1,j}(yID1,xID1) + match{sid,2,j}(yID2,xID2) + match{sid,3,j}(yID3,xID3) ...
                              + match{sid,4,j}(yID4,xID4) + match{sid,5,j}(yID5,xID5) + match{sid,6,j}(yID6,xID6) ...
                              + match{sid,7,j}(yID7,xID7) + match{sid,8,j}(yID8,xID8) + match{sid,9,j}(yID9,xID9);
            end
            surfscore(i) = max(scalescore) - rectModel.b;
        end
        score(oid).surfscore = surfscore;
        score(oid).surfID = PROJECTVIEW(visibility);
    end
    cuboid_rect_score{hid} = score;
end

clear rectModel

%% compute multiple layer of segmentation
rotImg  = im2double(rotImg);
thresholds = [200 500 800 1200 2000 3000];
if ~hasBuff
    for i = 1:length(thresholds)
        all_seg(i).labelMap = gbPanoSegment( im2uint8(rotImg), 0.5, thresholds(i), 50);
        all_seg(i).thresh = thresholds(i);
    end
end

%% checkseg
% load('./rectangleDetector/segmentation/uniformvector_lvl6.mat');
[ coor, tri ] = getUniformVector( 6 );
[sH, sW, ~] = size(all_seg(1).labelMap);
coor2D = uv2coords(xyz2uvN(coor), sW, sH);
coor2D_indi = sub2ind([sH sW], coor2D(:,2), coor2D(:,1));

segment_rect_score = cell(length(allhyps),1);
for hid = find(checkseg')
    hyps = allhyps{hid};
    vector_num = size(coor, 1);
    indicator = false(vector_num, length(hyps));
    for i = 2:length(hyps)
        if hyps(i).type==1
            p = hyps(i).out_points_w([1 2 3 4],:);
        else
            p = hyps(i).align;
        end
        p = p./repmat(sqrt(sum(p.^2,2)),1,3);
        vp = sum(p,1); vp = vp./norm(vp); %vp(3) = 0;
        [out2D, valid, ~] = projectPoint2SeparateView( p, vp, pi/3, 100 );
        if any(~valid)
            vp = sum(p,1); vp(1:2) = 0; vp = vp./norm(vp);
            [out2D, valid, ~] = projectPoint2SeparateView( p, vp, pi/3, 100 );
            if any(~valid)
                vp = sum(p,1); vp(3) = 0; vp = vp./norm(vp);
                [out2D, valid, ~] = projectPoint2SeparateView( p, vp, pi/3, 100 );
            end
        end
        K = convhull( out2D(:,1), out2D(:,2));
        [ inside, ~, ~ ] = insideCone( p(K(end:-1:2),:), coor, 0 );
        if ~any(inside)
            fprintf('%d:%d\n', hid, i);
        end
%         assert(any(inside));
        indicator(:,i) = inside;
    end

    seg_score = zeros( length(all_seg), length(hyps));
    for sid = 1:length(all_seg)
        label = all_seg(sid).labelMap(coor2D_indi);
        MINLAB = min(label);
        MAXLAB = max(label);
        N = hist(label, MINLAB:MAXLAB);
        for j = 1:length(hyps)
            lab = label(indicator(:,j));
            if isempty(lab)
                continue;
            end
            M = hist(lab, MINLAB:MAXLAB);
            I = M>0.05*sum(M);
            seg_score(sid,j) = sum(M(I))./sum(N(I));
        end
    end
    max_seg_score = max(seg_score, [], 1);
    segment_rect_score{hid} = max_seg_score;
end

if ~hasBuff
    save(bufName, 'all_hog', 'all_seg');
end

end

function [xID, yID] = getGridID(x, y, s)
% xID = round( (x-1) * s / 8 ) + 1;
% yID = round( (y-1) * s / 8 ) + 1;
xID = ceil( (x-1) * s / 8 ) - 2;
yID = ceil( (y-1) * s / 8 ) - 2;
end

author: zouchuhang

repo: LayoutNet


text:
%CALCOPTICALFLOWDF  DeepFlow optical flow algorithm implementation
%
%     flow = cv.calcOpticalFlowDF(I0, I1)
%
% ## Input
% * __I0__ First 8-bit single-channel grayscale input image.
% * __I1__ Second input image of the same size and type as `I0`.
%
% ## Output
% * __flow__ computed flow image that has the same size as `I0` and type
%   `single` (2-channels). Flow for `(x,y)` is stored in the third dimension.
%
% The class implements the DeepFlow optical flow algorithm described in
% [Weinzaepfel2013]. See also
% [DeepMatching](http://lear.inrialpes.fr/src/deepmatching/).
% Parameters - class fields - are:
%
% * __Alpha__ Smoothness assumption weight. default 1.0
% * __Delta__ Color constancy assumption weight. default 0.5
% * __Gamma__ Gradient constancy weight. default 5.0
% * __Sigma__ Gaussian smoothing parameter. default 0.6
% * __MinSize__ Minimal dimension of an image in the pyramid (next, smaller
%   images in the pyramid are generated until one of the dimensions reaches
%   this size). default 25
% * __DownscaleFactor__ Scaling factor in the image pyramid (must be <1).
%   default 0.95
% * __FixedPointIterations__ How many iterations on each level of the pyramid.
%   default 5
% * __SorIterations__ Iterations of Successive Over-Relaxation (solver).
%   default 25
% * __Omega__ Relaxation factor in SOR. default 1.6
%
% ## References
% [Weinzaepfel2013]:
% > Philippe Weinzaepfel, Jerome Revaud, Zaid Harchaoui, and Cordelia Schmid.
% > "Deepflow: Large displacement optical flow with deep matching".
% > In Computer Vision (ICCV), 2013 IEEE International Conference on,
% > pages 1385-1392. IEEE, 2013.
%
% See also: cv.calcOpticalFlowSF, cv.calcOpticalFlowFarneback
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function result = p00_exp_transform ( problem, order )

%*****************************************************************************80
%
%% P00_EXP_TRANSFORM applies an exponential transform and Gauss-Legendre rule.
%
%  Discussion:
%
%    To approximate:
%
%      Integral ( alpha <= x < Infinity ) f(x) dx
%
%    Transform:
%
%      u = exp ( -x )
%      du = - exp ( -x ) dx
%
%      x = - log ( u )
%      dx = - du / u
%
%      x = alpha    => u = exp ( -alpha )
%      x = Infinity => u = 0
%
%    Transformed integral:
%
%      Integral ( 0 <= u <= exp ( -alpha ) ) f ( -log(u) ) du / u
%
%    We apply a Gauss-Legendre rule here, but we could easily use any rule
%    that avoids evaluation at U = 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer PROBLEM, the index of the problem.
%
%    Input, integer ORDER, the order of the Gauss-Legendre rule
%    to apply.
%
%    Output, real RESULT, the approximate integral.
%
  alpha = p00_alpha ( problem );
%
%  Get the abscissas and weights for Gauss-Legendre quadrature.
%
  [ u, weight ] = legendre_compute ( order );
%
%  Modify the weights from [-1,1] to [0,exp(-alpha)].
%
  weight(1:order) = exp ( -alpha ) * weight(1:order) / 2.0;
%
%  Linear transform of abscissas from [-1,1] to [0,exp(-alpha)].
%
  u(1:order) = ( ( 1.0 + u(1:order) ) * exp ( - alpha ) ...
               + ( 1.0 - u(1:order) ) * 0.0 )           ...
               / ( 2.0              );
%
%  Define U_LOG = - log ( U )
%
  u_log(1:order) = - log ( u(1:order) );
%
%  Evaluate F ( -LOG(U) ).
%
  f_vec = p00_fun ( problem, order, u_log );
%
%  The integrand is F ( -LOG(U) ) / U
%
  f_vec(1:order) = f_vec(1:order) ./ u(1:order);
%
%  Sum.
%
  result = weight(1:order) * f_vec(1:order)';

  return
end


text:
function [vp f] = compute_vp(lines, imgsize, camcalib)
% compute vanishing points 'vp' using 'lines'
% assign line class to 'lines' and 'linesmore'

vpmultitype = 0;

% if camcalib.provided==1     % true camera parameters provided
%     [normal_vec, vp] = vanishline(lines, camcalib);
% else                        % images from the web

    if vpmultitype==0
        [vp f] = vanish_from_minevidence(lines, imgsize);

    elseif vpmultitype==1
        [normal_vec, vp] = vanishline_weakorthoconst(lines, camcalib);

    elseif vpmultitype==2
        vpset = vanishline_multi(lines, camcalib);
        % script_dispmultivp(vpset, lines, linesmore, img);
        % figure;

    elseif vpmultitype==3
        vpset = vanishline_multi3(lines, camcalib);
        % script_dispmultivp(vpset, lines, linesmore, img);
        % figure;

    elseif vpmultitype==4
        vpset = vanishline_multi4(lines, camcalib);
        % script_dispmultivp(vpset, lines, linesmore, img);
    end

% end


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

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

% Last Modified by GUIDE v2.5 13-Mar-2004 17:41:39

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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

% --- Executes on slider movement.
function Slider1_Callback(hObject, eventdata, handles)
% hObject    handle to Slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Find the value of the slider
value = get(handles.Slider1,'Value');

% Place the value in the text field
str = sprintf('%.2f',value);
set (handles.Label1,'String',str);

% --- Executes during object creation, after setting all properties.
function Slider1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


text:
function out = spm_run_fmri_data(job)
% Set up the design matrix and run a design
% SPM job execution function
% takes a harvested job data structure and call SPM functions to perform
% computations on the data.
% Input:
% job    - harvested job data structure (see matlabbatch help)
% Output:
% out    - computation results, usually a struct variable.
%__________________________________________________________________________
% Copyright (C) 2005-2019 Wellcome Trust Centre for Neuroimaging

% $Id: spm_run_fmri_data.m 7739 2019-12-02 14:00:18Z guillaume $

original_dir = pwd;
cd(spm_file(job.spmmat{1},'fpath'));

load(fullfile(pwd,'SPM.mat'));

%-Image filenames
%--------------------------------------------------------------------------
SPM.xY.P = char(job.scans);

%-Let SPM configure the design
%--------------------------------------------------------------------------
SPM = spm_fmri_spm_ui(SPM);

if ~isempty(job.mask{1})
    SPM.xM.VM         = spm_data_hdr_read(job.mask{:});
    SPM.xM.xs.Masking = [SPM.xM.xs.Masking, '+explicit mask'];
end

%-Save SPM.mat
%--------------------------------------------------------------------------
fprintf('%-40s: ','Saving SPM configuration')                           %-#
fmt = spm_get_defaults('mat.format');
s = whos('SPM');
if s.bytes > 2147483647, fmt = '-v7.3'; end
save('SPM.mat','SPM', fmt);
fprintf('%30s\n','...SPM.mat saved')                                    %-#

out.spmmat{1} = fullfile(pwd, 'SPM.mat');

cd(original_dir);

author: spm


text:
%RUNSCRIPT Run an M-file in interactive fashion
%
% RUNSCRIPT(FNAME, OPTIONS) runs the M-file FNAME and pauses after every
% executable line in the file until a key is pressed.  Comment lines are shown
% without any delay between lines.
%
% Options::
% 'delay',D    Don't wait for keypress, just delay of D seconds (default 0)
% 'cdelay',D   Pause of D seconds after each comment line (default 0)
% 'begin'      Start executing the file after the comment line %%begin (default false)
% 'dock'       Cause the figures to be docked when created
% 'path',P     Look for the file FNAME in the folder P (default .)
% 'dock'       Dock figures within GUI
%
% Notes::
% - If no file extension is given in FNAME, .m is assumed.
% - If the executable statement has comments immediately afterward (no blank lines)
%   then the pause occurs after those comments are displayed.
% - A simple '-' prompt indicates when the script is paused, hit enter.
% - If the function cprintf() is in your path, the display is more
%   colorful, you can get this file from MATLAB Central.
% - If the file has a lot of boilerplate, you can skip over and not display
%   it by giving the 'begin' option which searchers for the first line
%   starting with %%begin and commences execution at the line after that.
%
% See also eval.

% Copyright (C) 1993-2014, 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 runscript(fname, varargin)

    opt.path = [];
    opt.delay = [];
    opt.begin = false;
    opt.cdelay = 0;
    opt.dock = false;

    opt = tb_optparse(opt, varargin);

    close all

    prevDockStatus = get(0,'DefaultFigureWindowStyle');
    if opt.dock
        set(0,'DefaultFigureWindowStyle','docked');
    else
        set(0,'DefaultFigureWindowStyle','normal');
    end

    if ~isempty(opt.path)
        fname = fullfile(opt.path, [fname '.m']);
    else
        fname = [fname '.m'];
    end

    fp = fopen(fname, 'r');

    clc
    fprintf('--- runscript <-- %s\n', fname);

    running = false;
    shouldPause = false;
    savedText = [];

    if ~opt.begin
        running = true;
    end

    lineNum = 1;

    % stashMode
    %  0 normal
    %  1 loop
    %  2 continuation
    continMode = false;
    compoundDepth = 0;

    while 1
        % get the next line from the file, bail if EOF
        line = fgetl(fp);
        if line == -1
            break
        end
        lineNum = lineNum+1;

        % logic to skip lines until we see one beginning with %%begin
        if ~running
            if strcmp(line, '%%begin')
                running = true;
            else
                continue;
            end
        end;

        if length(strtrim(line)) == 0
            % blank line
            fprintf('\n');
            if shouldPause
                scriptwait(opt);
                shouldPause = false;
            end
            continue
        elseif startswith(strtrim(line), '%')
            % line was a comment
            disp(line)
            pause(opt.cdelay)  % optional comment delay
            continue;
        else
            if shouldPause
                scriptwait(opt);
                shouldPause = false;
            end
        end

        % if the start of a loop, stash the text for now
        if startswith(line, 'for') || startswith(line, 'while') || startswith(line, 'if')
            % found a compound block, don't eval it until we get to the end
            compoundDepth = compoundDepth + 1;
        end
        % if the statement has a continuation
        if endswith(line, '...')  && compoundDepth == 0
            % found a compound statement, don't eval it until we get to the end
            continMode = true;
        end

        if compoundDepth == 0 && ~continMode
            prompt = '>> ';
        else
            prompt = '';
        end

        % display the line with a pretend MATLAB prompt
        if exist('cprintf')
            cprintf('blue', '%s%s\n', prompt, line)
        else
            fprintf('%s', prompt); disp(line)
        end

        if compoundDepth > 0 || continMode
            % we're in stashing mode
            savedText = strcat(savedText, '\n', line);
        end

        if compoundDepth > 0 && startswith(line, 'end')
            % the compound block is fully unnested

            compoundDepth = compoundDepth - 1;
            if compoundDepth == 0
                evalSavedText(savedText, lineNum);
                savedText = '';
                shouldPause = true;
            end
            continue

        elseif continMode && ~endswith(line, '...')
            % no longer in continuation mode

            evalSavedText(savedText, lineNum);
            savedText = '';
            continMode = false;
            shouldPause = true;
            continue
        end

        if compoundDepth == 0 && ~continMode
            % it's a simple executable statement, execute it
            fprintf(' \n');
            evalSavedText(line, lineNum);
            shouldPause = true;
        end
    end
    fprintf('------ done --------\n');
    % restore the docking mode if we set it
    set(0,'DefaultFigureWindowStyle', prevDockStatus)
end

    function evalSavedText(text, lineNum)
        if length(strtrim(text)) == 0
            return
        end

        text = sprintf(text);

        try
            evalin('base', text);
        catch m
            fprintf('error in script %s at line %d', lineNum);
            m.rethrow();
        end
    end

% delay or prompt according to passed options
function scriptwait(opt)
    if isempty(opt.delay)
        %a = input('-', 's');
        prompt = 'continue?';
        if exist('cprintf')
            cprintf('red', prompt); pause;
            for i=1:length(prompt)
                cprintf('text', '\b');
            end
        else
            fprintf(prompt); pause;
            for i=1:length(prompt)
                fprintf('\b');
            end
        end
    else
        pause(opt.delay);
    end
end

% test if s2 is at the start of s1 (ignoring leading spaces)
function res = startswith(s1, s2)

    s1 = strtrim(s1);  % trim leading white space
    r = strfind(s1, s2);
    res = false;
    if ~isempty(r) && (r(1) == 1)
        res = true;
    end
end

% test if s2 is at the end of s1
function res = endswith(s1, s2)

    if length(s1) < length(s2)
        res = false;
    else
        n2 = length(s2)-1;
        res = strcmp(s1(end-n2:end), s2);
    end

end

author: Allopart

repo: rbpf-gmapping


text:
function [ x, error_norm, iter, flag ]  = gs ( A, x, b, max_it, tol )

%*****************************************************************************80
%
%% GS solves the linear system Ax=b using the Gauss-Seidel Method.
%
%  Modified:
%
%    02 April 2006
%
%  Reference:
%
%    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
%      June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
%      Charles Romine, Henk van der Vorst
%    Templates for the Solution of Linear Systems: Building Blocks for
%      Iterative Methods,
%    SIAM Publications, 1993.
%
%  Parameters:
%
%    Input, real A(N,N), the symmetric positive definite matrix.
%
%    Input, real X(N), the initial guess vector.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, integer MAX_IT, the maximum number of iterations.
%
%    Input, real TOL, an error tolerance.
%
%    Output, real X(N), the solution.
%
%    Output, real ERROR_NORM, the norm of the error.
%
%    Output, integer ITER, the number of iterations performed.
%
%    Output, integer FLAG, the return flag.
%    0 = the solution was found to within the specified tolerance.
%    1 = a satisfactory solution was not found.  The iteration limit
%        was exceeded.
%

%
%  Initialization.
%
  flag = 1;
  iter = 0;

  bnrm2 = norm ( b );
  if ( bnrm2 == 0.0 )
    bnrm2 = 1.0
  end

  r = b - A * x;
  error_norm = norm ( r ) / bnrm2;
  errorhist = [ ];
  errorhist(1) = error_norm;

  if ( error_norm < tol )
    flag = 0;
    return
  end
%
%  Split the matrix.
%
  w = 1.0;
  [ M, N, b ] = split ( A, b, w, 3 );

  for iter = 1 : max_it

    x_1 = x;
%
%  Update the approximation.
%
    x = M \ ( N * x + b );
%
%  Compute the error.
%
    error_norm = norm ( x - x_1 ) / norm ( x );
    errorhist(iter+1) = error_norm;
%
%  Check for convergence.
%
    if ( error_norm <= tol )
      flag = 0;
      break
    end

  end

  error_norm = errorhist;

  return
end


text:
function legendre_zeros_test ( )

%*****************************************************************************80
%
%% LEGENDRE_ZEROS_TEST tests LEGENDRE_ZEROS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_ZEROS_TEST:\n' );
  fprintf ( 1, '  LEGENDRE_ZEROS computes the zeros of the N-th Legendre\n' );
  fprintf ( 1, '  polynomial.\n' );

  for n = 1 : 7
    l = legendre_zeros ( n );
    r8vec_print ( n, l, '  Legendre zeros' );
  end

  return
end


text:
function [Y] = diagonalize(X,samplesize)
% ----------------------- Input ------------------------------
% X: samples of all tasks (each row is a sample)
% samplesize: the i-th entry is the sample size of the i-th task
% ----------------------- Output -----------------------------
% Y: sparse data matrix which is diagonal

tasknum = length(samplesize); % the number of tasks
[totalnum,dim] = size(X);
row = zeros(dim*totalnum,1); col = row; datavec = row;
accumsize = 0; accumind = 0;
for i = 1:tasknum
    accumsize = accumsize + samplesize(i);
    accumind = accumind + dim*samplesize(i);
    indsample = (accumsize-samplesize(i)+1 : accumsize)';
    indnz = (accumind-dim*samplesize(i)+1 : accumind)';
    row(indnz) = repmat(indsample,dim,1);
    col(indnz) = reshape(repmat(dim*(i-1)+1:dim*i,samplesize(i),1),dim*samplesize(i),1);
%     Xi = X(indsample,:)/sqrt(tasknum*samplesize(i));
    Xi = X(indsample,:);
    datavec(indnz) = Xi(:);
end
Y = sparse(row,col,datavec,totalnum,dim*tasknum);

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function status = besa_save2PmgPos(file_path, file_name, ...
    channel_coords, channel_orientations, channel_labels)
% besa_save2PmgPos writes the channel coordinates and the corresponding
% labels into the file 'file_name'.
% Brief *.pos and *.pmg file format description:
% - one sensor per line
% - magnetometers: label (optional), six coordinates per line (location,
%   orientation)
% - gradiometers: label (optional), nine coordinates per line (location of
%   primary sensor, location of secondary sensor, orientation).
%
% Parameters:
%     [file_path]
%         Full path to the folder where the file should be saved.
%
%     [file_name]
%         The name of the file where the output should be written.
%
%     [channel_coords]
%         Cell containing the channel coordinates.
%
%     [channel_orientations]
%         Cell containing the channel coordinates.
%
%     [channel_labels]
%         Optional parameter. Cell (size: number of channels x 1)
%         containing the channel labels.
%
% Return:
%     [status]
%         The status of the writing process. Still unused.
%

% Copyright (C) 2013, BESA GmbH
%
% File name: besa_save2PmgPos.m
%
% This file is part of MATLAB2BESA.
%
%    MATLAB2BESA 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.
%
%    MATLAB2BESA 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 MATLAB2BESA. If not, see <http://www.gnu.org/licenses/>.
%
% Author: Robert Spangler
% Created: 2014-01-11

n_channelCoords = size(channel_coords, 1);
n_channelOrientations = size(channel_orientations, 1);
% Check if number of coordinates and orientations match
if n_channelCoords ~= n_channelOrientations
    error('Mismatch in number of sensor coordinates and orientations!');
end;

% Check number of orientations is 3
if size(channel_orientations, 2) ~= 3
    error('Three orientations per sensor required!');
end;

% Check if number of coordinates and labels match
n_channelLabels = 0;
if exist('channel_labels', 'var')
    n_channelLabels = size(channel_labels, 1);
    % Check if number of channel coordinates matches number of channel
    % labels. If they do not match, do not write any labels to file.
    if n_channelCoords ~= n_channelLabels
        n_channelLabels = 0;
    end;
end;

% Change to the directory where the data should be saved.
cd(file_path)

% Open file for writing.
fid = fopen(file_name, 'w');

% Write the channel coordinates and the labels to the file.
for i=1:n_channelCoords
    % Write label
    if n_channelLabels > 0
        fprintf(fid, '%s\t', channel_labels{i});
    end;
    % Write coordinates
    if size(channel_coords{i}, 2) == 3
        % Write magnetometer info: three coordinates
        fprintf(fid, '%2.6f\t%2.6f\t%2.6f\t', channel_coords{i});
    elseif size(channel_coords{i}, 2) == 6
        % Write gradiometer info: six coordinates
        fprintf(fid, '%2.6f\t%2.6f\t%2.6f\t%2.6f\t%2.6f\t%2.6f\t',...
            channel_coords{i});
    else
        % Number of coordinates and orientations does not mach file format
        % description!
        error('Invalid number of coordinates and/or orientations!');
    end;

     % Write orientations
     fprintf(fid, '%2.6f\t%2.6f\t%2.6f\n', channel_orientations(i, :));
end

fclose(fid);
status = 1;

author: fieldtrip

repo: fieldtrip


text:
function r8mat_exp_story ( test )

%*****************************************************************************80
%
%% R8MAT_EXP_STORY prints explanatory text for each problem.
%
%  Discussion:
%
%     1) Diagonal example
%     2) Symmetric example
%     3) Laub
%     4) Moler and Van Loan
%     5) Moler and Van Loan
%     6) Moler and Van Loan
%     7) Moler and Van Loan
%     8) Wikipedia example
%     9) NAG F01ECF
%    10) Ward #1
%    11) Ward #2
%    12) Ward #3
%    13) Ward #4
%    14) Moler example
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Laub,
%    Review of "Linear System Theory" by Joao Hespanha,
%    SIAM Review,
%    Volume 52, Number 4, December 2010, page 779-781.
%
%    Cleve Moler, Charles VanLoan,
%    Nineteen Dubious Ways to Compute the Exponential of a Matrix,
%    Twenty-Five Years Later,
%    SIAM Review,
%    Volume 45, Number 1, March 2003, pages 3-49.
%
%    Cleve Moler,
%    Cleve's Corner: A Balancing Act for the Matrix Exponential,
%    July 23rd, 2012.
%
%    Robert Ward,
%    Numerical computation of the matrix exponential with accuracy estimate,
%    SIAM Journal on Numerical Analysis,
%    Volume 14, Number 4, September 1977, pages 600-610.
%
%  Parameters:
%
%    Input, integer TEST, the index of the test case.
%
  if ( test == 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This matrix is diagonal.\n' );
    fprintf ( 1, '  The calculation of the matrix exponential is simple.\n' );
  elseif ( test == 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This matrix is symmetric.\n' );
    fprintf ( 1, '  The calculation of the matrix exponential is straightforward.\n' );
  elseif ( test == 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This example is due to Laub.\n' );
    fprintf ( 1, '  This matrix is ill-suited for the Taylor series approach.\n' );
    fprintf ( 1, '  As powers of A are computed, the entries blow up too quickly.\n' );
  elseif ( test == 4 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This example is due to Moler and Van Loan.\n' );
    fprintf ( 1, '  The example will cause problems for the series summation approach,\n' );
    fprintf ( 1, '  as well as for diagonal Pade approximations.\n' );
  elseif ( test == 5 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This example is due to Moler and Van Loan.\n' );
    fprintf ( 1, '  This matrix is strictly upper triangular\n' );
    fprintf ( 1, '  All powers of A are zero beyond some (low) limit.\n' );
    fprintf ( 1, '  This example will cause problems for Pade approximations.\n' );
  elseif ( test == 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This example is due to Moler and Van Loan.\n' );
    fprintf ( 1, '  This matrix does not have a complete set of eigenvectors.\n' );
    fprintf ( 1, '  That means the eigenvector approach will fail.\n' );
  elseif ( test == 7 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This example is due to Moler and Van Loan.\n' );
    fprintf ( 1, '  This matrix is very close to example 5.\n' );
    fprintf ( 1, '  Mathematically, it has a complete set of eigenvectors.\n' );
    fprintf ( 1, '  Numerically, however, the calculation will be suspect.\n' );
  elseif ( test == 8 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This matrix was an example in Wikipedia.\n' );
  elseif ( test == 9 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This example is due to the NAG Library.\n' );
    fprintf ( 1, '  It is an example for function F01ECF.\n' );
  elseif ( test == 10 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This is Ward''s example #1.\n' );
    fprintf ( 1, '  It is defective and nonderogatory.\n' );
    fprintf ( 1, '  The eigenvalues are 3, 3 and 6.\n' );
  elseif ( test == 11 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This is Ward''s example #2.\n' );
    fprintf ( 1, '  It is a symmetric matrix.\n' );
    fprintf ( 1, '  The eigenvalues are 20, 30, 40.\n' );
  elseif ( test == 12 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This is Ward''s example #3.\n' );
    fprintf ( 1, '  Ward''s algorithm has difficulty estimating the accuracy\n' );
    fprintf ( 1, '  of its results.  The eigenvalues are -1, -2, -20.\n' );
  elseif ( test == 13 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This is Ward''s example #4.\n' );
    fprintf ( 1, '  This is a version of the Forsythe matrix.\n' );
    fprintf ( 1, '  The eigenvector problem is badly conditioned.\n' );
    fprintf ( 1, '  Ward''s algorithm has difficulty estimating the accuracy\n' );
    fprintf ( 1, '  of its results for this problem.\n' );
  elseif ( test == 14 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  This is Moler''s example.\n' );
    fprintf ( 1, '  This badly scaled matrix caused problems for MATLAB''s expm().\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_EXP_STORY - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of TEST = %d\n', test );
    error ( 'R8MAT_EXP_STORY - Fatal error!' );
  end

  return
end


text:
function eyeInDeg = eyetrack2deg_roughBars(eyeTrackMatFile)
% eyetrack2deg_roughBars - convert fMRI pixel based eyemovement data
% with the roughBars stimuli to degrees of visual angle
%
% eyeInDeg = eyetrack2deg_roughBars(eyeTrackMatFile);
%
% 2007/08 SOD: wrote it.

if ~exist('eyeTrackMatFile','var') || isempty(eyeTrackMatFile)
    error('need eyeTrackMatFile');
end

% frames per second
framerate = 30;

% calculate number of time points
stimrep = 5;
tr      = 1.5;
%prescan = 1
%ntime = (14.*stimrep.*tr).*framerate;

% copied from makeRetinotopyStimulus_roughBars, but we removed the
% last mean luminance presentation because we clipped the eye movie
% from the first stimulus presentation to the last.
sequence = [5:8 9:12 13:16 21,21 17:20 9:12 3,4,1,2 21,21 15,16 13,14 11,12 9,10 7,8 5,6 21,21 19,20 17,18 11,12 9,10 1:4];
fixseq   = ceil(sequence./4);
fixseq = reshape(fixseq,2,numel(fixseq)./2)';% pair up
fixseq = fixseq(:,1);
fixseq = repmat(fixseq,1,tr.*framerate.*stimrep);
fixseq = fixseq';
fixseq = fixseq(:);

% fill in fixseq with nearest values
while sum(fixseq==6),
    a = diff(fixseq==6);
    ii=find(a==1);
    fixseq(ii+1)= fixseq(ii);
    ii=find(a==-1);
    fixseq(ii) = fixseq(ii+1);
end

% load eyetrackdata
if ischar(eyeTrackMatFile)
    load(eyeTrackMatFile);
else
    et = eyeTrackMatFile;
end

% compute error:
slip = abs(numel(fixseq)-size(et,1));
fprintf(1,'[%s]:Timing error (difference) is about: %.2f seconds.\n',mfilename,slip/framerate);

% simply cut the difference
if numel(fixseq)<size(et,1)
    et = et(1:numel(fixseq),:);
else
    et = [et; zeros(slip,size(et,2))];
end

% interpolate missing data
t  = (0:numel(fixseq)-1)'./framerate;
eyeClosed = et(:,1)==0 | et(:,2)==0 | ~isfinite(sum(et(:,1:2),2));

% interpolate single points (probably due to fit failure)
singlePoints = eyeClosed + ...
               [eyeClosed(2:end); 0]+...
               [0; eyeClosed(1:end-1)];
singlePoints = singlePoints==1;
et(:,1) = interp1(t(~singlePoints),et(~singlePoints,1),t,'linear',mean(et(~singlePoints,1)));
et(:,2) = interp1(t(~singlePoints),et(~singlePoints,2),t,'linear',mean(et(~singlePoints,2)));
eyeClosed = eyeClosed & ~singlePoints;

% ignore 10 seconds around eyeblick
goodData = 1-min(eyeClosed+...
                [eyeClosed(2:end); 0]+...
                [eyeClosed(3:end); 0; 0]+...
                [eyeClosed(4:end); 0; 0; 0]+...
                [eyeClosed(5:end); 0; 0; 0; 0]+...
                [eyeClosed(6:end); 0; 0; 0; 0; 0]+...
                [0; eyeClosed(1:end-1)]+...
                [0; 0; eyeClosed(1:end-2)]+...
                [0; 0; 0; eyeClosed(1:end-3)]+...
                [0; 0; 0; 0; eyeClosed(1:end-4)]+...
                [0; 0; 0; 0; 0; eyeClosed(1:end-5)],1);%sigh...
goodData = logical(goodData);
%et_old = et;
et(:,1) = interp1(t(goodData),et(goodData,1),t,'linear',mean(et(goodData,1)));
et(:,2) = interp1(t(goodData),et(goodData,2),t,'linear',mean(et(goodData,2)));

% smooth estimates
% 1: median filter (jumps of 1/30th of seconds are unlikely)
% 2: gaussian filter (smoothness constraint fwhm=100ms)
met = ones(size(et,1),1)*mean(et(:,1:2));
et(:,1:2) = blurTC(et(:,1:2)-met,3,3,10)+met;

% remove data around fixation changes = 0; remove [-1 2] seconds
% find fixation changes
lim = [-1 2];
fc = find(diff(fixseq)~=0)-framerate*abs(min(lim));
fcii = fc(:)*ones(1,framerate*diff(lim))+...
    ones(size(fc(:)))*(0:framerate*diff(lim)-1);
fcii = fcii';
fcii = fcii(:);
badCalData = false(size(fixseq));
badCalData(fcii) = true;

% find good calibration data
goodCalData = goodData & ~badCalData;

%--- limit to good calibration data
f =fixseq(goodCalData);
% measured fixation in pixels
X = et(goodCalData,1:2);
% time
t = t(goodCalData);

Y = [-1 -1;...
     +1 -1;...
      0  0;...
     -1 +1;...
     +1 +1].*5;
Xsummary = zeros(max(fixseq),2);
fixseq_deg = zeros(size(fixseq,1),2);
for n=1:5,
    ii = fixseq==n;
    fixseq_deg(ii,:) = ones(sum(ii),1)*Y(n,:);
    ii = f==n;
    Xsummary(n,:) = mean(X(ii,:));
end

Xt = et(:,1:2);
Xt2 = Xt(goodCalData,:);
for n=1:5,
    ii = f==n;
    Xsummary(n,:) = mean(Xt2(ii,:));
end

% linear xfm
[Xi xfm] = linear_xfm(Xsummary,Y);
xfm(abs(xfm)<eps) = 0;
X2  = [Xt(:,1:2) ones(size(et,1),1)];
Xt = X2*xfm;
Xt = Xt(:,1:2);
xfm1 = [xfm';0 0 1];
xfm2 = xfm1;

c = logical([ones(2,3); zeros(1,3); 1 1 0]);

id = [eye(3); 0 0 1];

a=zeros(1001,9);
a(1,:)=[xfm1(c)' sqrt(sum(skewness(Xt(goodCalData,1:2)).^2))];

% now carefully add perspective scaling with boosting and monitoring the
% skewness we like it if the skewness goes down but not up.
iii = 1;
while(1),
    Xt2 = Xt(goodCalData,:);
    for n=1:5,
        ii = f==n;
        Xsummary(n,:) = mean(Xt2(ii,:));
    end

    % perspective xfm
    [Xi xfm] = perspective_xfm(Xsummary,Y);

    tmp = xfm(c)-id(c);
    tmp(abs(tmp)<10*eps)=0;
    if ~any(tmp), break; end
    [vM iM]  = max(abs(tmp));
    iM  = iM(1);
    vM  = tmp(iM);
    tmp = zeros(size(tmp));
    tmp(iM) = min(abs(vM),1e-3).*sign(vM);
    tmp = tmp + id(c);

    xfm = id;
    xfm(c) = tmp;

    X2  = [Xt(:,1:2) ones(size(et,1),1)]';
    Xtt = ((xfm(1:3,1:3)*X2) ./ (ones(3,1)*(xfm(4,:)*X2)))';
    Xtt = Xtt(:,1:2);

    a(iii+1,:) = [xfm2(c)' sqrt(sum(skewness(Xtt(goodCalData,1:2)).^2))];

    % monitor the skewness
    if round(a(iii+1,9)*1e4)>round(max(1.05*a(iii,9)*1e4,0.2*1e4))
        %[xfm1 xfm2],
        break
    else
        Xt = Xtt;
        xfm2 = xfm2 + xfm - id;
    end
    iii = iii + 1;
end
% % linear xfm
% [Xi xfm] = linear_xfm(Xsummary,Y);
% X2  = [et(:,1:2) ones(size(et,1),1)];
% Xt = (X2*xfm);
% Xt = Xt(:,1:2);
%
% Xt2 = Xt(goodCalData,:);
% for n=1:5,
%     ii = f==n;
%     Xsummary(n,:) = mean(Xt2(ii,:));
% end
% % perspective xfm
% [Xi xfm] = perspective_xfm(Xsummary,Y);
% xfm{1} =
% xfm{2}
% X2  = [Xt(:,1:2) ones(size(et,1),1)]';
% %Xt = ((xfm{1}*X2) ./ (ones(3,1)*(xfm{2}*X2)))';
% Xt = Xt(:,1:2);
%
%

Xi = Xt(goodCalData,:);

% output
eyeInDeg.blinks    = eyeClosed;
eyeInDeg.eye       = Xt;
eyeInDeg.framerate = framerate;
eyeInDeg.fixation  = fixseq_deg;
eyeInDeg.goodData  = goodCalData;

% compute derivative
der = diff(eyeInDeg.eye)*framerate;
der = ([ones(1,size(der,2));der] + [der; ones(1,size(der,2))])./2;
vel = hypot(der(:,1),der(:,2));

eyeInDeg.velocity = vel;

% plot
if ~nargout,
    colors = {'m.','y.','b.','g.','r.'};

    figure(1);clf;hold on;
    for n=1:5,
        plot(X(f==n,1),X(f==n,2),colors{n},'markersize',5);
    end
    ylabel('y (pixels)');
    xlabel('x (pixels)');

    figure(2);clf;hold on;
    for n=1:5,
        plot(Xi(f==n,1),Xi(f==n,2),colors{n},'markersize',5);
    end
    for n=1:5,
        plot(Y(n,1),Y(n,2),'kx','markersize',15,'linewidth',5);
        plot(mean(Xi(f==n,1)),mean(Xi(f==n,2)),'ko','markersize',15,'linewidth',3);
    end
    ylabel('y (deg)');
    xlabel('x (deg)');
    axis([-20 20 -20 20])

    figure(3);clf;
    subplot(2,1,1);hold on;
    for n=1:5,
        plot(t(f==n),Xi(f==n,1),colors{n},'markersize',5);
    end
    title('x');
    ylabel('position (deg)');
    axis([1 max(t) -20 20]);

    subplot(2,1,2);hold on;
    for n=1:5,
        plot(t(f==n),Xi(f==n,2),colors{n},'markersize',5);
    end
    title('y')
    xlabel('time (sec)');
    ylabel('position (deg)');
    axis([1 max(t) -20 20]);

    figure(4);clf;
    subplot(2,1,1);plot(a(1:iii,1:8)-ones(iii,1)*a(1,1:8));
    subplot(2,1,2);plot(a(1:iii,9),'k');
    h=axis;
    axis([h(1) h(2) 0-.01 max(h(4),0.2)+.01]);
end


text:
function l = list(varargin)

%tstoolbox/@list/list
%   Syntax:
%     * l = list
%       creates empty list
%     * l = list('Hello world')
%       create list with one entry, 'Hello world'
%     * l = list('Hello', 'My' , 'World')
%       create list with three entries
%     * l = list('Hello', 'My' , 'World')
%       create list with three entries
%
%   An object of type list contains a list of strings.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

if nargin == 0 			% create an empty list
	l.data = {};		% implementation as one-dimensional cell array (columnwise)
	l.len = 0;
	l = class(l, 'list');
else					% create list from arguments
	if iscellstr(varargin)
		l = list(varargin);
	elseif ischar(varargin{1})
	    l = list(cellstr(varargin{1}));
	elseif isa(varargin{1}, 'list')
		l = varargin{1};
	elseif iscellstr(varargin{1})
		l.data = varargin{1}(:);
		l.len = length(varargin{1});
		l = class(l, 'list');
	else
		error('Wrong type of argument(s) given');
	end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [DCM] = VBA_spm_dcm_fmri_check(P)
% post-hoc diagnostics for DCM (bilinear or nonlinear) of fMRI data
% FORMAT [DCM] = spm_dcm_fmri_check(DM)
% DCM - DCM structure or its filename
%
% This routine provides some diagnostics to ensure model inversion has
% converged. It plots the predicted and observed responses over all regions
% and provides the coefficient of determination ? or percent variance
% explained. This should normally be above 10%. An abnormally low
% coefficient of determination is highlighted in red. Quantitatively, one
% would normally expect to see one or more extrinsic (between source)
% connections with the strength of 1/8 Hz or greater. If all the extrinsic
% posterior expectations are below this value, then this suggests a failure
% of convergence or that the data are very noisy (possibly due to using
% very small regions of interest to summarise regional responses). Finally,
% the posterior correlations among all parameters are shown in terms of a
% correlation matrix. The number of effective parameters estimated is
% reported in terms of the (KL) divergence between the posterior and
% prior densities over parameters. This is divided by the log of the
% number of observations, by appealing to the Bayesian information
% criterion. The divergence corresponds to complexity or Bayesian
% surprise. Normally, one would expect the posterior and prior to diverge
% in a non-trivial fashion.
%
% Posterior densities are shown as bars with 90% confidence intervals in
% pink. An informed model inversion would normally provide posterior
% densities with confidence intervals that are, for some connections,
% displaced from prior expectations (at or around zero).
%
% This routine is compatible with DCM8, DCM10 and DCM12 files.
%__________________________________________________________________________
% Copyright (C) 20012 Wellcome Trust Centre for Neuroimaging
% Karl Friston

% $Id: spm_dcm_fmri_check.m karl $

%-Load DCM structure
%--------------------------------------------------------------------------
if ~nargin
  uiopen('load');
elseif isstruct(P)
  DCM = P;
else
  load(P)
end

% Assemble diagnostics
%==========================================================================
% coefficient of determination (percent variance explained)
%--------------------------------------------------------------------------
PSS  = sum(sum(DCM.y.^2));
RSS  = sum(sum(DCM.R.^2));
D(1) = 100*PSS/(PSS + RSS);

% largest absolute posterior expectation (extrinsic connections)
%--------------------------------------------------------------------------
try
   A = DCM.Ep.A;
catch
   A = DCM.A;
end
D(2) = max(max(abs(A - diag(diag(A)))));

% complexity and effective number of parameters estimated
%--------------------------------------------------------------------------
qE = VBA_spm_vec(DCM.Ep);
pE = VBA_spm_vec(DCM.M.pE);
qC = DCM.Cp;
pC = DCM.M.pC;
k  = rank(full(pC));
pC = VBA_spm_inv(pC);

D(3) = (trace(pC*qC) + (pE - qE)'*pC*(pE - qE) - VBA_spm_logdet(qC*pC) - k)/2;
D(3) = D(3)/log(DCM.v);

% Plot summary of inversion
%==========================================================================
VBA_spm_figure('GetWin','DCM diagnostics'); clf

% plot predicted and observed regional responses
%--------------------------------------------------------------------------
subplot(2,1,1);
t = (1:DCM.v)*DCM.Y.dt;

plot(t,DCM.y,t,DCM.y + DCM.R,':');
str = sprintf('variance explained %0.0f%%', D(1));
str = {'responses and predictions',str};
if D(1) > 10
  title(str,'FontSize',16);
else
  title(str,'FontSize',16,'Color','r');
end
xlabel('time {seconds}');

% posterior densities over A parameters
%--------------------------------------------------------------------------
try
  i = VBA_spm_fieldindices(DCM.Ep,'A');
catch
  i = 1 + (1:DCM.n^2);
end
qE = VBA_spm_vec(DCM.Ep);
qC = DCM.Cp;
subplot(2,2,3)
VBA_spm_plot_ci(qE(i),qC(i,i)), hold on
str = sprintf('largest connection strength %0.2f', D(2));
str = {'intrinsic and extrinsic connections',str};
if D(2) > 1/8
  title(str,'FontSize',16);
else
  title(str,'FontSize',16,'Color','r');
end
xlabel('parameter}');
axis square

% posterior correlations among all parameters
%--------------------------------------------------------------------------
subplot(2,2,4)
imagesc(VBA_cov2corr(DCM.Cp))
title('posterior correlations','FontSize',16)
str = sprintf('estimable parameters %0.0f', D(3));
str = {'posterior correlations',str};
if D(3) > 1
  title(str,'FontSize',16);
else
  title(str,'FontSize',16,'Color','r');
end
axis square

author: MBB-team

repo: VBA-toolbox


text:
% Compute the error of a pose-pose constraint
% x1 3x1 vector (x,y,theta) of the first robot pose
% x2 3x1 vector (x,y,theta) of the second robot pose
% z 3x1 vector (x,y,theta) of the measurement
%
% You may use the functions v2t() and t2v() to compute
% a Homogeneous matrix out of a (x, y, theta) vector
% for computing the error.
%
% Output
% e 3x1 error of the constraint
% A 3x3 Jacobian wrt x1
% B 3x3 Jacobian wrt x2
function [e, A, B] = linearize_pose_pose_constraint(x1, x2, z)

  % TODO compute the error and the Jacobians of the error
    X1 = v2t(x1);
    X2 = v2t(x2);
    Z = v2t(z);
    e = t2v(Z\(X1\X2));

    Rij = Z(1:3,1:3);
    Ri = X1(1:3,1:3);

    thi = atan2(Ri(2,1),Ri(1,1));
    thij = atan2(Rij(2,1),Rij(1,1));

    xi = x1(1);
    yi = x1(2);
    xj = x2(1);
    yj = x2(2);

%    A = [-Rij'*Ri' 0 0; 0 -Rij'*Ri' 0; 0 0 -1];
%    B = [Rij'*Ri' 0 0; 0 Rij'*Ri' 0; 0 0 1];

    A = [-cos(thi)*cos(thij)+sin(thi)*sin(thij) -sin(thi)*cos(thij)-cos(thi)*sin(thij) 0; cos(thi)*sin(thij)+sin(thi)*cos(thij) sin(thi)*sin(thij)-cos(thi)*cos(thij) 0; 0 0 -1];

    A(1:2,3) = [cos(thij)*(-sin(thi)*(xj-xi)+cos(thi)*(yj-yi))+sin(thij)*(-cos(thi)*(xj-xi)-sin(thi)*(yj-yi)); -sin(thij)*(-sin(thi)*(xj-xi)+cos(thi)*(yj-yi))+cos(thij)*(-cos(thi)*(xj-xi)-sin(thi)*(yj-yi))];

    B = [cos(thi)*cos(thij)-sin(thi)*sin(thij) sin(thi)*cos(thij)+cos(thi)*sin(thij) 0; -cos(thi)*sin(thij)-sin(thi)*cos(thij) -sin(thi)*sin(thij)+cos(thi)*cos(thij) 0; 0 0 1];

end;

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function out=prod(x,dim)

if ~isempty(x)

 precision=x(1).precision;
 % can only handle up to 2 dimensions
 s=size(x);
 if nargin==1
  if any(s==1)
   out=mp(1,precision);
   for ii=1:numel(x)
    out=out*x(ii);
   end
  else
   out=mp(ones(1,s(2)),precision);
   for j=1:s(2)
    for i=1:s(1)
     out(1,j)=out(1,j)*x(i,j);
    end
   end
  end
 elseif nargin==2
  if dim==1
   out=mp(ones(1,s(2)),precision);
   for j=1:s(2)
    for i=1:s(1)
     out(1,j)=out(1,j)*x(i,j);
    end
   end
  elseif dim==2
   out=mp(ones(s(1),1),precision);
   for j=1:s(2)
    for i=1:s(1)
     out(i,1)=out(i,1)*x(i,j);
    end
   end
  else
   error('mp prod can only handle up to 2-D arrays (11/04)');
  end
 end

else
 out=mp(1);
end

author: opencobra

repo: cobratoolbox


text:
function [u,p] = StokesDGSRT0(A,B,u,p,f,g,itStep,auxMat)
%% STOKESDGS
%  Created by Ming Wang. Discussed with Long Chen.
% Reference:

Su = auxMat.trilA;
Bt = auxMat.Bt;
Sp = auxMat.Sp;
invMp = auxMat.invMp;

%% DGS relaxation step
for k = 1: itStep
    % Step 1: relax Momentum eqns twice
    for i = 1:2
        u = u + Su\(f-Bt*p-A*u);
    end
    % Step 2: relax Continuity eqns
    rp = g - B*u;
%     rp = rp - mean(rp);  Whether need this step ?????
    dq = Sp\rp;
    % Step 3: update u and p
    u = u + Bt*dq;
    p = p - invMp*rp;
end

repo: ifem


text:
function o = nifti2struc
% Create a data structure describing NIFTI-2 headers
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: nifti2struc.m 7147 2017-08-03 14:07:01Z spm $

persistent org;
if ~isempty(org)
    o = org;
    return;
end
t = struct(...
    'conv',{@char  , @uint8 , @int16 , @int32 , @int64 , @single , @double },...
    'prec',{'uint8', 'uint8', 'int16', 'int32', 'int64', 'single', 'double'},...
    'size',{   1   ,    1   ,    2   ,    4   ,    8   ,    4    ,    8    });
c = t(1);
b = t(2);
s = t(3);
i = t(4);
l = t(5);
f = t(6);
d = t(7);

table = {...
    i,  1, 'sizeof_hdr', 540
    c,  8, 'magic', ['ni2' char(0) sprintf('\r\n\032\n')]
    s,  1, 'datatype', 2
    s,  1, 'bitpix', 8
    l,  8, 'dim', [3 0 0 0 1 1 1 1]
    d,  1, 'intent_p1', 0
    d,  1, 'intent_p2', 0
    d,  1, 'intent_p3', 0
    d,  8, 'pixdim', [0 1 1 1]
    l,  1, 'vox_offset', 0
    d,  1, 'scl_slope', 1
    d,  1, 'scl_inter', 0
    d,  1, 'cal_max', []
    d,  1, 'cal_min', []
    d,  1, 'slice_duration', []
    d,  1, 'toffset', []
    l,  1, 'slice_start', []
    l,  1, 'slice_end', []
    c, 80, 'descrip', 'NIFTI-2 Image'
    c, 24, 'aux_file', ''
    i,  1, 'qform_code', 0
    i,  1, 'sform_code', 0
    d,  1, 'quatern_b', 0
    d,  1, 'quatern_c', 0
    d,  1, 'quatern_d', 0
    d,  1, 'qoffset_x', 0
    d,  1, 'qoffset_y', 0
    d,  1, 'qoffset_z', 0
    d,  4, 'srow_x', [1 0 0 0]
    d,  4, 'srow_y', [0 1 0 0]
    d,  4, 'srow_z', [0 0 1 0]
    i,  1, 'slice_code', []
    i,  1, 'xyzt_units', 10
    i,  1, 'intent_code', 0
    c, 16, 'intent_name', ''
    b,  1, 'dim_info', []
    c, 15, 'unused_str', ''};

org = struct('label',table(:,3),'dtype',table(:,1),'len',table(:,2),...
    'offset',0,'def',table(:,4));
os  = 0;
for j=1:length(org)
    os  = org(j).dtype.size*ceil(os/org(j).dtype.size);
    fun = org(j).dtype.conv;
    if ischar(org(j).def), z = char(0); else z = 0; end
    def = [org(j).def repmat(z,1,org(j).len-length(org(j).def))];
    org(j).def    = feval(fun,def);
    org(j).offset = os;
    os  = os + org(j).len*org(j).dtype.size;
end
o = org;


text:
function value = r4_gamr ( x )

%*****************************************************************************80
%
%% R4_GAMR evaluates the reciprocal gamma function of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the value of the reciprocal gamma
%    function at X.
%
  if ( x <= 0.0 && r4_aint ( x ) == x )

    value = 0.0;

  elseif ( abs ( x ) <= 10.0 )

    value = 1.0 / r4_gamma ( x );

  else

    [ alngx, sgngx ] = r4_lgams ( x );
    value = sgngx * exp ( - alngx );

  end

  return
end


text:
classdef prtOptionsComputation

    properties

        largestMatrixSize = 9e6; % Total number of elements in the matrix (3000 x 3000)
        verboseStorage = true;   % Save's the training dataset within each
                                 % action. True enables easier plotting
                                 % and data exploration.
        showProgressBar = true;
    end

    methods
        function obj = prtOptionsComputation(varargin)
            obj = prtUtilAssignStringValuePairs(obj,varargin{:});
        end
    end

    methods
        function obj = set.largestMatrixSize(obj,val)
            assert(prtUtilIsPositiveInteger(val),'largestMatrixSize must be a positive integer');

            if val < 1e3
                warning('prt:prtOptionsComputation','Although valid, this value for largestMatrixSize is small. Consider using a larger value.')
            end

            obj.largestMatrixSize = val;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT