text:
function spm_cli(varargin)
% Command line interface for SPM
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_cli.m 6964 2016-12-07 16:50:41Z guillaume $

%-Input arguments
%--------------------------------------------------------------------------
if isempty(varargin), return; end
if nargin == 1 && ismember(varargin{1},{'help','--help'})
    cmd = lower(spm('Ver'));
    fprintf([...
        'Usage: %s [NODE] [arg...]\n',...
        '       %s [NODE] help\n',...
        '       %s [ help | --help ]\n',...
        '       %s help --verbose\n',...
        '\n',...
        '[NODE]         Batch node name\n',...
        '--verbose      List all available batch nodes\n',...
        '--help, help   Print usage statement\n'],...
        cmd, cmd, cmd, cmd);
        return
end
if nargin == 2 ...
        && all(ismember(strtok(varargin(1:2),'-'),{'help','verbose'}))
    c0 = cfg_util('getcfg');
    nodes = find_all_modules(c0.values{1})';
    for i=1:numel(nodes)
        [n, c] = find_module(strsplit(nodes{i},'.'), c0);
        h = normalise_help_text(c.help);
        if ~c.hidden, fprintf('  %s    %s\n',nodes{i},h{1}); end
    end
    return
end
node = strsplit(varargin{1},'.');
opts = varargin(2:end);

%-Get the configuration tree
%--------------------------------------------------------------------------
try
    c0 = cfg_util('getcfg');
catch
    spm_jobman('initcfg');
    c0 = cfg_util('getcfg');
end

%-Search module in configuration tree and allow for node shortcut
% (use first matching module if ambiguity)
%--------------------------------------------------------------------------
[n, c0, found] = find_module(node, c0);

if ~found
    error('Cannot find module %s.\n',strjoin(node,'.'));
else
    n = n(2:end); % first is matlabbatch
    if ~isequal(n,node)
        node = n;
        fprintf('Matching module %s.\n',strjoin(node,'.'));
    end
end

%-Help for a specific node
%--------------------------------------------------------------------------
if numel(opts) == 1 && ismember(opts{1},{'help','--help'})
    cmd = lower(spm('Ver'));
    fprintf([...
        'Usage: %s %s [OPTIONS]\n',...
        '\n',...
        'Options:\n'],...
        cmd,strjoin(node,'.'));
    [flags, h1s, ms] = get_items(c0);
    n = cellfun(@numel,flags);
    for i=1:numel(flags)
        if ms(i), m = '(*) '; else m = ''; end
        fprintf('    --%s%s%s%s\n',flags{i},repmat(' ',1,3+max(n)-n(i)),m,h1s{i});
    end
    fprintf([...
        '\n',...
        '    --help%sPrint usage statement\n'],repmat(' ',1,3+max(n)-4));
    return;
end

%-Get options
%--------------------------------------------------------------------------
i = 1;
params = {};
values = {};
while i <= numel(opts)
    flag = strtok(opts{i},'-');
    opt = {};
    i = i + 1;
    for j=i:numel(opts)
        if opts{j}(1) == '-', break; end
        opt{end+1} = opts{j};
        i = i + 1;
    end
    params{end+1} = flag;
    values{end+1} = opt;
end

%-Create batch job structure
%--------------------------------------------------------------------------
job = struct();
for i=1:numel(params)
    prm = params{i};
    if ~any(prm == '.')
        subs = substruct('.',prm);
    else
        clear subs
        p = strsplit(prm,'.');
        for j=1:numel(p)
            subs(j) = substruct('.',p{j});
        end
        prm = p{end};
    end
    job = subsasgn(job, subs, get_val(values{i},prm,c0));
end

%-Execute batch job
%--------------------------------------------------------------------------
subs = reshape([repmat({'.'},1,numel(node));node],1,[]);
subs = substruct(subs{:});
job  = subsasgn(struct(),subs,job);
cfg_message('none','destination','^matlabbatch:run:jobstart$');
cfg_message('none','destination','^matlabbatch:run:jobdone$');
cfg_message('none','destination','^matlabbatch:run:jobfailed$');
cfg_message('none','destination','^matlabbatch:run:modstart$');
cfg_message('none','destination','^matlabbatch:run:moddone$');
%cfg_message('none','destination','^matlabbatch:run:modfailed$');
try
    [out, job] = spm_jobman('run', {job});
catch
    err = lasterror;
    switch err.identifier
        case 'spm:spm_jobman:jobNotFilled'
            error('Incomplete job.');
        case 'matlabbatch:run:jobfailederr'
            error('Job failed.');
        otherwise
            error('Unknown error during job execution.');
    end
end

%==========================================================================
function [node, c0, found] = find_module(node, c0, found)

if nargin < 3, found = false; end
if found, return; end
if strcmp(node{1}, c0.tag)
    if numel(node) == 1
        if isa(c0,'cfg_exbranch'), found = true; end
    else
        if ~isa(c0,'cfg_exbranch')
            for i=1:numel(c0.values)
                [n, c1 , found] = find_module(node(2:end), c0.values{i});
                if found, node = [c0.tag n]; c0 = c1; return; end
            end
        end
    end
else
    if ~isa(c0,'cfg_exbranch')
        for i=1:numel(c0.values)
            [n, c1, found] = find_module(node, c0.values{i});
            if found, node = [c0.tag n]; c0 = c1; return; end
        end
    end
end

%==========================================================================
function [nodes,n] = find_all_modules(c0,nodes,n)

if nargin < 2, nodes = {}; end
if nargin < 3, n = {}; end
if isa(c0,'cfg_exbranch')
    nodes = [nodes strjoin([n c0.tag],'.')];
else
    m = n;
    for i=1:numel(c0.values)
        [nodes, n] = find_all_modules(c0.values{i},nodes,[m c0.tag]);
    end
end

%==========================================================================
function param = get_val(param,flag,c0)

for i=1:numel(c0.val)
    if strcmp(c0.val{i}.tag, flag)
        switch class(c0.val{i})
            case 'cfg_files'
                param = param(:);
            case 'cfg_entry'
                switch c0.val{i}.strtype
                    case {'n','w','r'}
                        param = str2num(char(param));
                    case {'s'}
                        param = char(param);
                end
            case 'cfg_menu'
                if numel(param) > 1
                    error('Only one value accepted for "%s".',flag);
                end
                param = char(param);
                labels = c0.val{i}.labels;
                j = find(ismember(labels,param));
                if nnz(j)
                    param = c0.val{i}.values{j};
                else
                    opts = sprintf('  "%s"\n',labels{:});
                    error(['Cannot translate option "%s": %s.\n',...
                        'Options are:\n%s'],flag, param, opts);
                end
            otherwise
        end
        return;
    end
end
for i=1:numel(c0.val)
    switch class(c0.val{i})
        case 'cfg_branch'
            param = get_val(param,flag,c0.val{i});
        case 'cfg_repeat'
        case 'cfg_choice'
    end
end

%==========================================================================
function [flags, h1s, ms] = get_items(c0)
flags = {}; h1s = {}; ms = [];
for i=1:numel(c0.val) % go recursive instead
    flag = {}; h1 = {}; m = [];
    if isa(c0.val{i},'cfg_repeat')
        % if more than one, they are not stored yet
        for j=1:numel(c0.val{i}.values)
            [flag, h1, m] = get_items(c0.val{i}.values{j});
            flag = [{c0.val{i}.values{j}.tag} flag];
            h1   = [{char(c0.val{i}.values{j}.help{:})} h1];
            m    = [is_mandatory(c0.val{i}.values{j}) m];
        end
    elseif isa(c0.val{i},'cfg_branch')
        [flag, h1, m] = get_items(c0.val{i});
        flag = [{c0.val{i}.tag} flag];
        h1   = [{char(c0.val{i}.help{:})} h1];
        m    = [is_mandatory(c0.val{i}) m];
    else
        flag = c0.val{i}.tag;
        h1   = char(c0.val{i}.help{:});
        m    = is_mandatory(c0.val{i});
    end
    h1    = normalise_help_text(h1);
    flags = [flags flag];
    h1s   = [h1s h1];
    ms    = [ms m];
end

%==========================================================================
function m = is_mandatory(item)
m = isempty(item.val) && isempty(item.def);

%==========================================================================
function h1 = normalise_help_text(h1)
h1 = cellstr(h1);
for i=1:numel(h1)
    if ~isempty(h1{i}) && h1{i}(end) == '.', h1{i} = h1{i}(1:end-1); end
    if numel(h1{i}) > 70 % use "stty size" to get the number of columns
        h1{i} = [h1{i}(1:70-3) '...'];
    end
end

author: spm


text:
function board = lines(board)

% Add lines to your board by clicking the begin- and endcoordinates. Exit
% by right clicking.

ui = questdlg('Continuing will destroy all line-information that is already in your current board.', ...
              'Continue?', 'Continue', 'Stop');

if ~strcmp(ui, 'Continue')
    return;
end

for i = 1:numel(board)
    p = patch(board(i).xy(:,1), board(i).xy(:,2), [1 1 1]);
    set(p, 'EdgeColor', [0 0 0]);
end

set(gca, 'XTick', [], 'YTick', [], 'box', 'on')
xlim([0, 800])
ylim([0, 600])

eps = 5;

while true
    [x1,y1,button] = ginput(1);
    if button == 3
        break
    end

    [x2,y2,button] = ginput(1);
    if button == 3
        continue
    end

    % Try to make straight lines
    if abs(x1 - x2) < eps
        x1 = 1/2 * (x1 + x2);
        x2 = x1;
    end
    if abs(y1 - y2) < eps
        y1 = 1/2 * (y1 + y2);
        y2 = y1;
    end

    % No lines outside the borders (800x600)
    x = [x1, x2];
    x(x < 0) = 0;
    x(x > 800) = 800;
    y = [y1, y2];
    y(y < 0) = 0;
    y(y > 600) = 600;

    l = line(x, y, 'LineStyle', '--');
    ui = questdlg('OK?', 'OK?', 'Yes', 'No', 'Yes');
    switch ui
        case 'Yes'
            L(end + 1, :) = [x, y]; %#ok<AGROW>
        otherwise
            delete(l)
            continue
    end
end

close(gcf);

author: Sable

repo: mcbench-benchmarks


text:
function [ar,rc,ASAcontrol] = rc2arset(rc,req_order,last)
%RC2ARSET AR reflectioncoefficients to AR models
%   AR = RC2ARSET(RC) converts AR-reflectioncoefficients RC into an AR-
%   model of order LENGTH(RC)-1, with parameter vector AR. The procedure
%   implements the parameter relations of the Levinson-Durbin recursion.
%
%   [SET_AR,SET_RC] = RC2ARSET(RC,REQ_ORDER) returns intermediate AR
%   parameter vectors in the cell array SET_AR and an array SET_RC of
%   reflectioncoefficients, both corresponding to orders requested by
%   REQ_ORDER. REQ_ORDER must be either a row of ascending AR-orders, or
%   a single AR-order.
%
%   RC2ARSET is an ARMASA main function.
%
%   See also: COV2ARSET, AR2ARSET.

%   References: P. Stoica and R.L. Moses, Introduction to Spectral
%               Analysis, Prentice-Hall, Inc., New Jersey, 1997,
%               Chapter 3.

%Header
%=====================================================================

%Declaration of variables
%------------------------

%Declare and assign values to local variables
%according to the input argument pattern
switch nargin
case 1
   if isa(rc,'struct'), ASAcontrol=rc; rc=[];
   else, ASAcontrol=[];
   end
   req_order=[];
case 2
   if isa(req_order,'struct'), ASAcontrol=req_order; req_order=[];
   else, ASAcontrol=[];
   end
case 3
   if isa(last,'struct'), ASAcontrol=last;
   else, error(ASAerr(39))
   end
otherwise
   error(ASAerr(1,mfilename))
end

if isequal(nargin,1) & ~isempty(ASAcontrol)
      %ASAcontrol is the only input argument
   ASAcontrol.error_chk = 0;
   ASAcontrol.run = 0;
end

%ARMASA-function version information
%-----------------------------------

%This ARMASA-function is characterized by
%its current version,
ASAcontrol.is_version = [2000 12 30 20 0 0];
%and its compatability with versions down to,
ASAcontrol.comp_version = [2000 12 30 20 0 0];

%Checks
%------

if ~any(strcmp(fieldnames(ASAcontrol),'error_chk')) | ASAcontrol.error_chk
   %Perform standard error checks
   %Input argument format checks
   ASAcontrol.error_chk = 1;
   if ~isnum(rc)
      error(ASAerr(11,'rc'))
   end
   if ~isavector(rc)
      error(ASAerr(15,'rc'))
   elseif size(rc,1)>1
      rc = rc';
      warning(ASAwarn(25,{'column';'rc';'row'},ASAcontrol))
   end
   if ~isempty(req_order)
      if ~isnum(req_order) | ~isintvector(req_order) |...
            req_order(1)<0 | ~isascending(req_order)
         error(ASAerr(12,{'requested';'req_order'}))
      elseif size(req_order,1)>1
         req_order = req_order';
         warning(ASAwarn(25,{'column';'req_order';'row'},ASAcontrol))
      end
   end

   %Input argument value checks
   if ~isreal(rc)
      error(ASAerr(13))
   end
   if rc(1)~=1
      error(ASAerr(23,{'rc','reflectioncoefficient'}))
   end
   if ~isempty(req_order) & req_order(end) > length(rc)-1
      error(ASAerr(24,'reflectioncoefficients'))
   end
end

if ~any(strcmp(fieldnames(ASAcontrol),'version_chk')) | ASAcontrol.version_chk
      %Perform version check
   ASAcontrol.version_chk = 1;

   %Make sure the requested version of this function
   %complies with its actual version
   ASAversionchk(ASAcontrol);
end

if ~any(strcmp(fieldnames(ASAcontrol),'run')) | ASAcontrol.run
      %Run the computational kernel
   ASAcontrol.run = 1;

%Main
%=====================================================

%Recursion initialization
%------------------------

l_rc = length(rc);
max_k = l_rc;
ar = zeros(1,l_rc);
ar(1) = 1;
store = ~isempty(req_order);
if store
   counter = 1;
   max_counter = length(req_order);
   max_k = req_order(max_counter)+1;
   ar_stack = cell(max_counter,1);
   rc_stack = zeros(1,max_counter);
   if req_order(1)==0
      ar_stack{1} = 1;
      rc_stack(1) = 1;
      counter = counter+1;
   end
end

%Levinson Durbin parameter recursion
%-----------------------------------

for k = 2:max_k
   order = k-1;
   rc_temp = rc(k);
   ar(2:order) = ar(2:order)+rc_temp*ar(order:-1:2);
   ar(k) = rc_temp;
   if store & k==req_order(counter)+1
      ar_stack{counter} = ar(1:k);
      rc_stack(counter) = rc_temp;
      counter = counter+1;
   end
end

%Output argument arrangement
%---------------------------

if store
   ar = ar_stack;
   rc = rc_stack;
end

%Footer
%=====================================================

else %Skip the computational kernel
   %Return ASAcontrol as the first output argument
   if nargout>1
      warning(ASAwarn(9,mfilename,ASAcontrol))
   end
   ar = ASAcontrol;
   ASAcontrol = [];
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% former versions        P.M.T. Broersen        p.m.t.broersen@tudelft.nl
% [2000 12 30 20 0 0]    W. Wunderink           wwunderink01@freeler.nl


text:
function h = plotOrbit(color, sma, ecc, inc, raan, arg, minTru, maxTru, gmu, varargin)
%plotOrbit Summary of this function goes here
%   Detailed explanation goes here
    if(minTru >= maxTru)
        h=plotOrbit(color, sma, ecc, inc, raan, arg, minTru, 2*pi, gmu, varargin{:});
        plotOrbit(color, sma, ecc, inc, raan, arg, 0, maxTru, gmu, varargin{:});
        return;
    end

    if(~isempty(varargin) && ~isempty(varargin{1}))
        axisPlotTo = varargin{1};
    else
        axisPlotTo = gca;
    end

    try
        if(~isempty(varargin) && length(varargin)>=2)
            if(length(varargin{2}) == 3)
                rOffset = reshape(varargin{2},3,1);
            else
                rOffset = [0;0;0];
            end
        else
            rOffset = [0;0;0];
        end
    catch
        rOffset = [0;0;0];
    end

    try
        if(~isempty(varargin) && length(varargin)>=3 && length(varargin{3}) == 1)
            lineWidth = varargin{3};
        else
            lineWidth = 1.5;
        end
    catch
        lineWidth = 1.5;
    end

    if(~isempty(varargin) && length(varargin)>=4)
        lineStyle = varargin{4};
    else
        lineStyle = '-';
    end

    tru = linspace(minTru,maxTru,250);
    zeroArr = zeros(size(tru));
    sma = zeroArr + sma;
    ecc = zeroArr + ecc;
    inc = zeroArr + inc;
    raan = zeroArr + raan;
    arg = zeroArr + arg;

    [rVect,~] = vect_getStatefromKepler(sma, ecc, inc, raan, arg, tru, gmu, false);
    x = rVect(1,:) + rOffset(1);
    y = rVect(2,:) + rOffset(2);
    z = rVect(3,:) + rOffset(3);

    hold(axisPlotTo,'on');
    if(ischar(color))
        colorNew = color;
        try
            h = plot3(axisPlotTo, x,y,z, 'Color', colorNew, 'LineWidth', lineWidth, 'LineStyle', lineStyle);
        catch
            h = plot3(axisPlotTo, x,y,z, 'Color', color, 'LineWidth', lineWidth, 'LineStyle', lineStyle);
        end
    else
        h = plot3(axisPlotTo, x,y,z, lineStyle, 'Color', color, 'LineWidth', lineWidth);
    end

    hold(axisPlotTo,'off');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function X = diff(varargin)
% DIFF (overloaded)

Y = varargin{1};
X = Y;
X.basis = [];
for i = 1:size(Y.basis,2)
    base = reshape(full(Y.basis(:,i)),X.dim);
    base = diff(base,varargin{2:end});
    X.basis = [X.basis sparse(base(:))];
end
X.dim = size(base);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [accepted, rejected] = fun_metropolis_hastings(param_init, iterations, data)

x = param_init;

accepted = [];
rejected = [];

for i = 1:iterations
    x_new = fun_transition_model(x);
    x_lik = fun_log_likelihood(x, data);
    x_new_lik = fun_log_likelihood(x_new, data);

    if fun_acceptance(x_lik + log(fun_prior(x)), x_new_lik + log(fun_prior(x_new)))
        x = x_new;
        accepted = [accepted x_new(1)];
    else
        rejected = [rejected x_new(1)];
    end
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function prob_test0276 ( )

%*****************************************************************************80
%
%% TEST0276 tests CARDIOID_MEAN, CARDIOID_SAMPLE, CARDIOID_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  sample_num = 1000;

  a = 0.0;
  b = 0.25;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0276\n' );
  fprintf ( 1, '  For the Cardioid PDF:\n' );
  fprintf ( 1, '  CARDIOID_MEAN computes the mean;\n' );
  fprintf ( 1, '  CARDIOID_SAMPLE samples;\n' );
  fprintf ( 1, '  CARDIOID_VARIANCE computes the variance.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %f\n', a );
  fprintf ( 1, '  PDF parameter B = %f\n', b );

  if ( ~cardioid_check ( a, b ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return;
  end

  mean = cardioid_mean ( a, b );
  variance = cardioid_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF mean =                    %f\n', mean );
  fprintf ( 1, '  PDF variance =                %f\n', variance );

  for i = 1 : sample_num
    [ x(i), seed ] = cardioid_sample ( a, b, seed );
  end

  mean = r8vec_mean ( sample_num, x );
  variance = r8vec_variance ( sample_num, x );
  xmax = max ( x(1:sample_num) );
  xmin = min ( x(1:sample_num) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', sample_num );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gpuids = getGpuIds(gpuname)
    deviceCount = getGpuCount_mex();
    gpuids = int32(0:(deviceCount-1));
    for idx = 1:deviceCount
        name = getGpuName_mex(gpuids(idx));
        if ~strcmp(name, gpuname)
            gpuids(idx) = -1;
        end
    end
    gpuids = gpuids(gpuids>=0);
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function [niiC,doseFileNameC] = dose2nii(planC, doseNumV, scanNum, tmpDirPath,reorientFlag)

if ~exist('tmpDirPath','var') || isempty(tmpDirPath) || ~exist(tmpDirPath,'dir')
    tmpDirPath = fullfile(getCERRPath, 'ImageRegistration', 'tmpFiles');
end

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

[affineMat,~, voxel_size, ~, dose3MC] = getPlanCAffineMat(planC, scanNum, reorientFlag, [], doseNumV);
qOffset = affineMat(1:3,end)';

for i = 1:numel(doseNumV)
    dose3M = dose3MC{i};
    [scanUniqName, ~] = genScanUniqName(planC,scanNum);
    doseFileName = fullfile(tmpDirPath, ['dose_' num2str(doseNumV(i)) '_' scanUniqName '.nii']);
    niiC{i} = vol2nii(dose3M,affineMat,qOffset,voxel_size,doseFileName);
    doseFileNameC{i} = doseFileName;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function x = hermitian_semidefinite( n )

%HERMITIAN_SEMIDEFINITE   Complex Hermitian positive semidefinite matrices.
%    HERMITIAN_SEMIDEFINITE(N), where N is an integer, creates a complex
%    Hermitian matrix variable of size [N,N] and constrains it to be
%    positive semidefinite. Therefore, given the declaration
%       variable x(n,n) Hermitian
%    the constraint
%       x == hermitian_semidefinite(n)
%    is equivalent to
%       lambda_min(x) >= 0;
%    In fact, lambda_min is implemented in CVX using HERMITIAN_SEMIDEFINITE
%    for complex matrices.
%
%    HERMITIAN_SEMIDEFINITE(SX), where SX is a valid size vector, creates
%    an array variable of size SX and constrains each subarray along the
%    leading two dimensions to be positive semidefinite. SX(1) and SX(2)
%    must be equal. Therefore, given the declaration
%       variable x(sx) Hermitian
%    the constraint
%       x == hermitian_semidefinite(sx)
%    is equivalent to
%       for k = 1:prod(sx(3:end)),
%          lambda_min(x(:,:,k)) >= 0;
%       end
%
%   Disciplined convex programming information:
%       SEMIDEFINITE is a cvx set specification. See the user guide for
%       details on how to use sets.

x = semidefinite( n, true );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
function f_x = ParFor5(in1)
%PARFOR5
%    F_X = PARFOR5(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.1.
%    16-Jul-2019 11:57:25

u = in1(:,1);
uxx = in1(:,5);
f_x = (u.*9.981631376631412e-1+uxx.*9.965561702440027+u.*uxx.*1.006065455272619)./(u+9.978905414318433);

author: dynamicslab

repo: SINDy-PI


text:
function pick = nms_IoMin(boxes, overlap)
% top = nms(boxes, overlap)
% Non-maximum suppression. (FAST VERSION)
% Greedily select high-scoring detections and skip detections
% that are significantly covered by a previously selected
% detection.
%
% NOTE: This is adapted from Pedro Felzenszwalb's version (nms.m),
% but an inner loop has been eliminated to significantly speed it
% up in the case of a large number of boxes

% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if isempty(boxes)
  pick = [];
  return;
end

x1 = boxes(:,1);
y1 = boxes(:,2);
x2 = boxes(:,3);
y2 = boxes(:,4);
s = boxes(:,end);

area = (x2-x1+1) .* (y2-y1+1);
[vals, I] = sort(s);

pick = s*0;
counter = 1;
while ~isempty(I)
  last = length(I);
  i = I(last);
  pick(counter) = i;
  counter = counter + 1;

  xx1 = max(x1(i), x1(I(1:last-1)));
  yy1 = max(y1(i), y1(I(1:last-1)));
  xx2 = min(x2(i), x2(I(1:last-1)));
  yy2 = min(y2(i), y2(I(1:last-1)));

  w = max(0.0, xx2-xx1+1);
  h = max(0.0, yy2-yy1+1);

  inter = w.*h;
%   o = inter ./ (area(i) + area(I(1:last-1)) - inter);
  o = inter ./ min(area(i),area(I(1:last-1)));

  I = I(find(o<=overlap));
end

pick = pick(1:(counter-1));

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [state] = CliffFrontLeftSensorRoomba(serPort);
%[state] = CliffFrontLeftSensorRoomba(serPort)
% Specifies the state of the Cliff Front Left sensor
% Either triggered or not triggered

% By; Joel Esposito, US Naval Academy, 2011

%Initialize preliminary return values
state = nan;
try

%Flush Buffer
N = serPort.BytesAvailable();
while(N~=0)
fread(serPort,N);
N = serPort.BytesAvailable();
end

warning off
global td

fwrite(serPort, [142]);  fwrite(serPort,10);
CliffFrntLft = dec2bin(fread(serPort, 1));
state = bin2dec(CliffFrntLft(end));

pause(td)
catch
    disp('WARNING:  function did not terminate correctly.  Output may be unreliable.')
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function r8mat_indicator_test ( )

%*****************************************************************************80
%
%% R8MAT_INDICATOR_TEST tests R8MAT_INDICATOR_TEST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_INDICATOR_TEST\n' );
  fprintf ( 1, '  R8MAT_INDICATOR returns an indicator matrix.\n' );

  a = r8mat_indicator ( m, n );

  r8mat_print ( m, n, a, '  Indicator matrix:' );

  return
end


text:
function sendToSamiam(pedigree, factorList, alleleList, phenotypeList, positions, outputFileName)
% This function creates a file for Samiam using a list of factors and a
% pedigree.
% Input:
%   pedigree: Data structure that includes names and parent-child
%   relationships
%   factorList: List of factors
%	alleleList: List of alleles
%   phenotypeList: List of phenotypes
%   pedigreePositionsFileName: List with positions for the nodes for
%   each person in the family -- each row is for a person, and the contents
%   of each row is (genotype horizontal position, genotype vertical
%   position, phenotype horizontal position, phenotype vertical position)
%   outputFileName: Name of the file that will go to Samiam, should end in
%   .net

numPeople = length(pedigree.names);

outputFileName = strcat(outputFileName, '.net');
fid = fopen(outputFileName, 'w+');
fprintf(fid, '%s', 'net');
fprintf(fid, '\n', '');
fprintf(fid, '%s', '{');
fprintf(fid, '\n', '');
fprintf(fid, '\t', '');
fprintf(fid, '%s', 'node_size = (90 36);');
fprintf(fid, '\n', '');
fprintf(fid, '%s', '}');
fprintf(fid, '\n', '');
fprintf(fid, '\n', '');

for i = 1:numPeople
    % Make a genotype node and a phenotype node for each person
    fprintf(fid, '%s', 'node ');
    name = pedigree.names{i};
    fprintf(fid, '%s', name);
    fprintf(fid, '%s', 'Genotype');
    fprintf(fid, '\n', '');
    fprintf(fid, '%s', '{');
    fprintf(fid, '\n', '');
    fprintf(fid, '\t', '');
    fprintf(fid, '%s', 'label = "');
    fprintf(fid, '%s', name);
    fprintf(fid, '%s', 'Genoytpe";');
    fprintf(fid, '\n', '');
    fprintf(fid, '\t', '');
    fprintf(fid, '%s', 'position = (');
    fprintf(fid, '%d', positions(i,1));
    fprintf(fid, '%s', ' ');
    fprintf(fid, '%d', positions(i,2));
    fprintf(fid, '%s', ');');
    fprintf(fid, '\n', '');
    fprintf(fid, '\t', '');
    fprintf(fid, '%s', 'states = (');
    for j = 1:length(alleleList)
        % Iterate through alleles
        for k = j:length(alleleList)
            % Iterate through pairs of alleles
            fprintf(fid, '%s', '"');
            fprintf(fid, '%s', alleleList{j});
            fprintf(fid, '%s', alleleList{k});
            fprintf(fid, '%s', '"');
            if (j ~= length(alleleList)) || (k ~= length(alleleList))
                % Put in a space because there will be another genotype
                fprintf(fid, '%s', ' ');
            end
        end
    end
    fprintf(fid, '%s', ');');
    fprintf(fid, '\n', '');
    fprintf(fid, '%s', '}');
    fprintf(fid, '\n', '');
    fprintf(fid, '\n', '');

    fprintf(fid, '%s', 'node ');
    fprintf(fid, '%s', name);
    fprintf(fid, '%s', 'Phenotype');
    fprintf(fid, '\n', '');
    fprintf(fid, '%s', '{');
    fprintf(fid, '\n', '');
    fprintf(fid, '\t', '');
    fprintf(fid, '%s', 'label = "');
    fprintf(fid, '%s', name);
    fprintf(fid, '%s', 'Phenotype";');
    fprintf(fid, '\n', '');
    fprintf(fid, '\t', '');
    fprintf(fid, '%s', 'position = (');
    fprintf(fid, '%d', positions(i,3));
    fprintf(fid, '%s', ' ');
    fprintf(fid, '%d', positions(i,4));
    fprintf(fid, '%s', ');');
    fprintf(fid, '\n', '');
    fprintf(fid, '\t', '');
    fprintf(fid, '%s', 'states = (');
    for j = 1:length(phenotypeList)
        % Iterate through phenotypes
        fprintf(fid, '%s', '"');
        fprintf(fid, '%s', phenotypeList{j});
        fprintf(fid, '%s', '"');
        if j ~= length(phenotypeList)
            fprintf(fid, '%s', ' ');
        end
    end
    fprintf(fid, '%s', ');');
    fprintf(fid, '\n', '');
    fprintf(fid, '%s', '}');
    fprintf(fid, '\n', '');
    fprintf(fid, '\n', '');
end

for i = 1:2*numPeople
    % Iterate through factors and make a potential for each
    factor = factorList(i);
    fprintf(fid, '%s', 'potential (');
    firstNode = factor.var(1);
    if firstNode <= numPeople
        % Genotype node
        firstName = pedigree.names{firstNode};
        fprintf(fid, '%s', firstName);
        fprintf(fid, '%s', 'Genotype |');
        if pedigree.parents(firstNode, 1) ~= 0
            % Person has parents
            parentOneName = pedigree.names{factor.var(2)};
            parentTwoName = pedigree.names{factor.var(3)};
            fprintf(fid, '%s', ' ');
            fprintf(fid, '%s', parentOneName);
            fprintf(fid, '%s', 'Genotype ');
            fprintf(fid, '%s', parentTwoName);
            fprintf(fid, '%s', 'Genotype');
        end
        fprintf(fid, '%s', ')');
        fprintf(fid, '\n', '');
        fprintf(fid, '%s', '{');
        fprintf(fid, '\n', '');
        fprintf(fid, '\t', '');
        fprintf(fid, '%s', 'data = ');
        fprintf(fid, '%s', '(');
        for k = 1:length(factor.val)
            % Iterate through the values of the factor
            fprintf(fid, '%f', factor.val(k));
            if k == length(factor.val)
                % At end of factor
                fprintf(fid, '%s', ')');
                fprintf(fid, '%s', ';');
                fprintf(fid, '\n', '');
                fprintf(fid, '%s', '}');
                fprintf(fid, '\n', '');
                fprintf(fid, '\n', '');
            else
                fprintf(fid, '%s', ' ');
            end
        end
    else
        nameNode = factor.var(2);
        name = pedigree.names{nameNode};
        fprintf(fid, '%s', name);
        fprintf(fid, '%s', 'Phenotype | ');
        fprintf(fid, '%s', name);
        fprintf(fid, '%s', 'Genotype)');
        fprintf(fid, '\n', '');
        fprintf(fid, '%s', '{');
        fprintf(fid, '\n', '');
        fprintf(fid, '\t', '');
        fprintf(fid, '%s', 'data = ');
        for k = 1:length(factor.var)
            fprintf(fid, '%s', '(');
        end
        varIndexes = ones(1, length(factor.var));
        for k = 1:length(factor.val)
            % Iterate through the values of the factor
            fprintf(fid, '%f', factor.val(k));
            closedParen = 0;
            j = length(factor.var);
            while j > 0
                % Iterate through variables to see what closed parentheses are
                % necessary
                if varIndexes(j) == factor.card(j)
                    % Need to close parentheses
                    closedParen = closedParen + 1;
                    fprintf(fid, '%s', ')');
                    varIndexes(j) = 1;
                else
                    varIndexes(j) = varIndexes(j) + 1;
                    break
                end
                j = j - 1;
            end
            if (closedParen > 0) && (closedParen < length(factor.var))
                % Need a new line
                fprintf(fid, '\n', '');
                fprintf(fid, '\t', '');
                fprintf(fid, '\t', '');
                for j = 1:closedParen
                    % Make new open parentheses
                    fprintf(fid, '%s', '(');
                end
            elseif closedParen == length(factor.var)
                % At end of factor
                fprintf(fid, '%s', ';');
                fprintf(fid, '\n', '');
                fprintf(fid, '%s', '}');
                fprintf(fid, '\n', '');
                fprintf(fid, '\n', '');
            else
                fprintf(fid, '%s', ' ');
            end
        end
    end
end

fclose(fid);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [cMap]=fire(varargin)

% function [cMap]=fire(n)
% ------------------------------------------------------------------------
% Creates the colormap data for n levels for the fire colormap.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2015/01/01
%------------------------------------------------------------------------

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

% cMap=[0 0 0; 1 0 0; 1 1 0; 1 1 1]; %Simple version
cMap=[0,0,0;0.112,0.0130,0;0.213,0.0240,0.00100;0.323,0.0370,0.00100;...
       0.421,0.0520,0.00100;0.526,0.0720,0.00200;0.620,0.100,0.00300;...
       0.693,0.131,0.0110;0.770,0.172,0.0190;0.837,0.217,0.0320;0.894,0.261,0.0430;...
       0.954,0.315,0.0520;0.993,0.379,0.0630;0.998,0.484,0.0780;1,0.574,0.110;...
       1,0.644,0.152;1,0.711,0.211;1,0.773,0.274;1,0.837,0.322;1,0.914,0.369;...
       1,0.953,0.443;1,0.982,0.539;0.999,1,0.646;0.985,1,0.772;1,1,1];

[cMap]=resampleColormap(cMap,n);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [err, time] = spincomp(S, N, dt, pref)
%SPINCOMP  Compare time-stepping schemes in 1D/2D/3D.
%   SPINCOMP(S, N, DT, PREF) compares the time-stepping schemes in PREF.SCHEME
%   applied to the PDE specified by the SPINOP/SPINOP2/SPINOP3 S, using a fixed
%   number N of grid points and different time-steps DT. The results are shown
%   in a pair of log-log plots. PREF is a SPINPREF/SPINPREF2/SPINPREF3 object.
%
%   [ERR, TIME] = SPINCOMP(S, N, DT, PREF) returns the errors and the computer
%   times.
%
% Example 1: Compare the ETDRK schemes for 1D Kuramoto-Sivashinsky equation
%
%   dom = [0 32*pi]; tspan = [0 30];
%   u0 = chebfun('cos(x/16).*(1 + sin(x/16))',dom,'trig');
%   S = spinop('KS');
%   S.init = u0; S.tspan = tspan;
%   N = 256; dt = [1, 5e-1, 2e-1, 1e-1, 5e-2, 2e-2];
%   pref = spinpref;
%   pref.scheme = {'etdrk4', 'exprk5s8', 'krogstad', 'friedli', ...
%      'hochbruck-ostermann', 'minchev', 'strehmel-weiner'};
%   spincomp(S, N, dt, pref);
%
% Example 2: Compare the PREDICTOR-CORRECTOR schemes for 2D Gray-Scott equations
%
%   G = 3; dom = G*[0 1 0 1]; tspan = [0 10];
%   u01 = chebfun2(@(x,y) 1-exp(-150*((x-G/2).^2 + (y-G/2).^2)), dom, 'trig');
%   u02 = chebfun2(@(x,y) exp(-150*((x-G/2).^2 + 2*(y-G/2).^2)), dom, 'trig');
%   u0 = chebmatrix(u01);
%   u0(2,1) = u02;
%   S = spinop2('GS'); S.init = u0; S.tspan = tspan;
%   N = 128; dt = [1, 5e-1, 2e-1, 5e-2, 2e-2];
%   pref = spinpref2;
%   pref.scheme = {'pec423', 'pecec433', 'pec524', 'pecec534', 'pec625', ...
%      'pecec635', 'pec726', 'pecec736'};
%   spincomp(S, N, dt, pref);
%
% See also SPINOP, SPINOP2, SPINOP3, SPINSCHEME.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Get the schemes:
schemes = pref.scheme;
nSchemes = length(schemes);

% Get the time-steps:
timesteps = dt;
nTimesteps = length(timesteps);

% Convert the initial conditiona to a CHEBMATRIX:
u0 = S.init;
if ( isa(u0, 'chebfun') == 1 || isa(u0, 'chebfun2') == 1 || ...
    isa(u0, 'chebfun3') == 1)
    u0 = chebmatrix(u0);
elseif ( isa(u0, 'chebfun2v') == 1 || isa(u0, 'chebfun3v') == 1 )
    temp = chebmatrix(u0(1));
    for k = 2:size(u0, 1)
        temp(k,1) = u0(k);
    end
    u0 = temp;
end
nVars = length(u0);

% Create a grid to compare exact and computed solutions:
dom = u0{1}.domain;
xx = trigpts(32, dom);
if ( isa(u0{1}, 'chebfun') == 1 )
    dim = 1;
elseif ( isa(u0{1}, 'chebfun2') == 1 )
    dim = 2;
    [xx, yy] = meshgrid(xx);
elseif ( isa(u0{1}, 'chebfun3') == 1 )
    dim = 3;
    [xx, yy, zz] = meshgrid(xx);
end

% First, estimate the exact solution using a very small time step (half the
% smallest time-step) and with PECEC736 (7th-order multistep scheme):
if ( isa(pref, 'spinpref') == 1 )
    prefu = spinpref();
elseif ( isa(pref, 'spinpref2') == 1 )
    prefu = spinpref2();
elseif ( isa(pref, 'spinpref3') == 1 )
    prefu = spinpref3();
end
prefu.scheme = 'pecec736';
prefu.plot = 'off';
prefu.M = pref.M;
uexact = spinoperator.solvepde(S, N, min(timesteps)/2, prefu);
if ( isa(uexact, 'chebmatrix') == 0 )
    uexact = chebmatrix(uexact);
end

% Scale (i.e., maximum amplitude) of the exact solution:
scale = 0;
for i = 1:nVars
    if ( dim == 1 )
        scale = max(scale, max(abs(uexact{i}(xx))));
    elseif( dim == 2 )
        scale = max(scale, max(max(abs(uexact{i}(xx,yy)))));
    elseif ( dim == 3 )
        scale = max(scale, max(max(max(abs(uexact{i}(xx,yy,zz))))));
    end
end

% Second, compute solutions for different schemes and time-steps:
err = zeros(nTimesteps, nSchemes);
time = zeros(nTimesteps, nSchemes);
for k = 1:nTimesteps
    for l = 1:nSchemes
        prefu.scheme = schemes{l};
        [u, ~, t] = spinoperator.solvepde(S, N, timesteps(k), prefu);
        time(k,l) = t;
        if ( isa(u, 'chebmatrix') == 0 )
            u = chebmatrix(u);
        end
        isNan = isNanTest(u, dim);
        if ( isNan == 1 )
            err(k,l) = NaN;
        else
            for i = 1:nVars
                if ( dim == 1 )
                    temp = abs(u{i}(xx) - uexact{i}(xx));
                elseif( dim == 2 )
                    temp = abs(u{i}(xx,yy) - uexact{i}(xx,yy));
                elseif ( dim == 3 )
                    temp = abs(u{i}(xx,yy,zz) - uexact{i}(xx,yy,zz));
                end
                temp = max(temp(:));
                err(k,l) = max(err(k,l), temp);
            end
        end
    end
end
err = err/scale; % Relative error
tspan = S.tspan;
TF = tspan(end);
timesteps = timesteps/TF; % Relative time-steps

% Plot Accuarcy vs time-step:
figure, subplot(1, 2, 1)
labels = cell(nSchemes, 1);
for l = 1:nSchemes
    labels{l} = schemes{l};
    if ( l > 1 )
        hold on
    end
    if ( l < 8 )
        loglog(timesteps, err(:, l), '.-', 'linewidth', 2, 'markersize', 30)
    elseif ( l < 15 )
        loglog(timesteps, err(:, l), 'x-', 'linewidth', 2, 'markersize', 10)
    else
        loglog(timesteps, err(:, l), 'o-', 'linewidth', 2, 'markersize', 10)
    end
end
left = 10^floor(log10(min(timesteps)));
right = 10^ceil(log10(max(timesteps)));
down = max(min(1e-10, 10^floor(log10(min(err(:))))), 1e-12);
up = min(max(1e0, 10^ceil(log10(max(err(:))))), 1e2);
set(gca, 'fontsize', 16), axis([left right down up])
xlabel('Relative time-step'), ylabel(sprintf('Relative error at t = %.3f', TF))
legend(labels, 'Location', 'NorthWest')

% Plot Accuarcy vs computer time:
subplot(1, 2, 2)
for l = 1:nSchemes
    if ( l > 1 )
        hold on
    end
    if ( l < 8 )
        loglog(time(:,l), err(:, l), '.-', 'linewidth', 2, 'markersize', 30)
   elseif ( l < 15 )
        loglog(time(:,l), err(:, l), 'x-', 'linewidth', 2, 'markersize', 10)
    else
        loglog(time(:,l), err(:, l), 'o-', 'linewidth', 2, 'markersize', 10)
    end
end
left = 10^(floor(log10(min(min(time)))));
right = 10^(ceil(log10(max(max(time)))));
set(gca, 'fontsize', 16), axis([left right down up])
xlabel('Computer time (s)'), ylabel(sprintf('Relative error at t = %.3f', TF))
legend(labels, 'Location', 'NorthEast')

end

function out = isNanTest(u, dim)

B = u.blocks;
if ( dim == 1 )
    out = any(cell2mat(cellfun(@(C) isnan(C), B, 'UniformOutput', 0)));
elseif ( dim == 2 )
    out = any(cell2mat(cellfun(@(C) isnan(C.cols), B, 'UniformOutput', 0)));
    out = out || ...
        any(cell2mat(cellfun(@(C) isnan(C.rows), B, 'UniformOutput', 0)));
elseif ( dim == 3 )
    out = any(cell2mat(cellfun(@(C) isnan(C.cols), B, 'UniformOutput', 0)));
    out = out || ...
        any(cell2mat(cellfun(@(C) isnan(C.rows), B, 'UniformOutput', 0)));
    out = out || ...
        any(cell2mat(cellfun(@(C) isnan(C.tubes), B, 'UniformOutput', 0)));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%                               431-400 Year Long Project
%                           LA1 - Medical Image Processing 2003
%  Supervisor     :  Dr Lachlan Andrew
%  Group Members  :  Alister Fong    78629   a.fong1@ugrad.unimelb.edu.au
%                    Lee Siew Teng   102519  s.lee1@ugrad.unimelb.edu.au
%                    Loh Jien Mei    103650  j.loh1@ugrad.unimelb.edu.au
%
%  File and function name : test_sort_coord.m
%  Version                : 2.0
%  Date of completion     : 18 December 2003
%  Written by    :   Alister Fong    78629   a.fong1@ugrad.unimelb.edu.au
%
% Inputs :
%           None
%
% Outputs :
%           None
%
% Description :
%   Test the clockwise and anti-clockwise sorting function on a circle
%
% Usage >> test_sort_coord

function test_sort_coord
clc
outlined_image = zeros(100);
outlined_image_coord = plot_circle(50,50,20,'bresenham');
assigned_values = outlined_image_coord(:,1);
assigned_values(:) = 1;
outlined_image = place_values(outlined_image_coord,assigned_values,outlined_image);

[Ycoord,Xcoord] = find(outlined_image>0);

disp(' ');
disp('--------------');
disp(strcat('Number of pixels :',num2str(length(Ycoord))));
disp('--------------');
disp(' ');

temp_image = zeros(size(outlined_image));

% -------------------------------------------------------------------
disp('Display the image via the "find" function');
f = figure;
for n = 1:1:length(Ycoord)
    temp_image(Ycoord(n),Xcoord(n)) = 1;
    imshow(temp_image);
    title('"find" function');
    pause(0.01);
end

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

% Clockwise rotation
disp('Press any key to continue....');
pause;

disp('sort_coord_pixels clockwise');
tic
output_coord = sort_coord_pixel([Xcoord,Ycoord],'clockwise');
toc

disp('Press any key to continue....');
pause;
disp('Display the image via "sort_coord_pixel" function');
temp_image(:) = 0;
for n = 1:1:length(output_coord(:,1))
    temp_image(output_coord(n,2),output_coord(n,1)) = 1;
    imshow(temp_image);
    title('"sort-coord-pixel" function - clockwise');
    pause(0.01);
end

% -------------------------------------------------------------------
% Anti-clockwise rotation
disp('Press any key to continue....');
pause;

disp('sort_coord_pixels anti-clockwise');
tic
output_coord = sort_coord_pixel([Xcoord,Ycoord],'anti-clockwise');
toc

disp('Press any key to continue....');
pause;
disp('Display the image via "sort_coord_pixel" function');
temp_image(:) = 0;
for n = 1:1:length(output_coord(:,1))
    temp_image(output_coord(n,2),output_coord(n,1)) = 1;
    imshow(temp_image);
    title('"sort-coord-pixel" function - anti-clockwise');
    pause(0.01);
end

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

% -------------------------------------------------------------------
disp('Press any key to continue....');
pause;
disp('Get image "test_sort1.tif"');
temp_image = imread('test_sort1.tif');
[Ycoord,Xcoord] = find(temp_image~=0);
disp('sort_coord_pixels clockwise');
input_coord = [Xcoord,Ycoord];
disp('------------------------');
tic
output_coord = sort_coord_pixel(input_coord,'clockwise','progress track');
toc
disp('------------------------');
disp('Press any key to continue....');
pause;
disp('Display the image via "sort_coord_pixel" function');
temp_image(:) = 0;
for n = 1:1:length(output_coord(:,1))
    temp_image(output_coord(n,2),output_coord(n,1)) = 1;
    imshow(temp_image);
    title('"sort-coord-pixel" function - clockwise');
    pause(0.01);
end

% -------------------------------------------------------------------
disp('Press any key to continue....');
pause;
disp('Get image "test_sort1.tif"');
temp_image = imread('test_sort1.tif');
[Ycoord,Xcoord] = find(temp_image~=0);
disp('sort_coord_pixels anti-clockwise');
input_coord = [Xcoord,Ycoord];
disp('------------------------');
tic
output_coord = sort_coord_pixel(input_coord,'anti-clockwise','progress track');
toc
disp('------------------------');
disp('Press any key to continue....');
pause;
disp('Display the image via "sort_coord_pixel" function');
temp_image(:) = 0;
for n = 1:1:length(output_coord(:,1))
    temp_image(output_coord(n,2),output_coord(n,1)) = 1;
    imshow(temp_image);
    title('"sort-coord-pixel" function - anti-clockwise');
    pause(0.01);
end

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

% -------------------------------------------------------------------
disp('Press any key to continue....');
pause;
disp('Get image "test_sort2.tif"');
temp_image = imread('test_sort2.tif');
[Ycoord,Xcoord] = find(temp_image~=0);
disp('sort_coord_pixels clockwise, 4 connectivity');
input_coord = [Xcoord,Ycoord];
disp('------------------------');
tic
output_coord = sort_coord_pixel(input_coord,'clockwise',4);
toc
disp('------------------------');
disp('Press any key to continue....');
pause;
disp('Display the image via "sort_coord_pixel" function');
temp_image(:) = 0;
for n = 1:1:length(output_coord(:,1))
    temp_image(output_coord(n,2),output_coord(n,1)) = 1;
    imshow(temp_image);
    title('"sort-coord-pixel" function - clockwise');
    pause(0.01);
end

% -------------------------------------------------------------------
disp('Press any key to continue....');
pause;
disp('Get image "test_sort2.tif"');
temp_image = imread('test_sort2.tif');
[Ycoord,Xcoord] = find(temp_image~=0);
disp('sort_coord_pixels anti-clockwise, 4 connectivity');
input_coord = [Xcoord,Ycoord];
disp('------------------------');
tic
output_coord = sort_coord_pixel(input_coord,'anti-clockwise',4);
toc
disp('------------------------');
disp('Press any key to continue....');
pause;
disp('Display the image via "sort_coord_pixel" function');
temp_image(:) = 0;
for n = 1:1:length(output_coord(:,1))
    temp_image(output_coord(n,2),output_coord(n,1)) = 1;
    imshow(temp_image);
    title('"sort-coord-pixel" function - anti-clockwise');
    pause(0.01);
end

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

disp('Press any key to Close figure....');
pause;
close(f)


text:
%% Recording algorithm behavior with MatlabBGL
% In this example, we will write a simple visitor that outputs an
% algorithm's behavior.  The algorithm we will examine is dijkstra_sp.  To
% examine the runtime behavior we will use a visitor which outputs a string
% every time a function is called.

%% Setup
% To begin, we load a graph.

load ../graphs/clr-25-2.mat

%%
% Next, let's check the documentation to see which functions to implement
% for the visitor

help dijkstra_sp

%%
% The help states that dijkstra_sp allows visitors functions for
% initialize_vertex, discover_vertex, examine_vertex, examine_edge,
% edge_relaxed, edge_not_relaxed, and finish_vertex.
%
% Rather than implementing 7 functions ourselves, we define two helper
% functions.  These helper functions return functions themselves.  There is
% one helper that returns a vertex visitor function and one helper than
% returns an edge visitor function.

vertex_vis_print_func = @(str) @(u) ...
    fprintf('%s called on %s\n', str, char(labels{u}));
edge_vis_print_func = @(str) @(ei,u,v) ...
    fprintf('%s called on (%s,%s)\n', str, char(labels{u}), char(labels{v}));

%%
% These anonymous functions return functions themselves.

ev_func = vertex_vis_print_func('examine_vertex');
ev_func(1)

%%
% I hope you see how these functions are useful in saving quite a bit of
% typing.

%% Calling dijkstra_sp
% We are almost done.  Now, we just have to setup the visitor structure to
% pass to the dijkstra_sp call.

vis = struct();
vis.initialize_vertex = vertex_vis_print_func('initialize_vertex');
vis.discover_vertex = vertex_vis_print_func('discover_vertex');
vis.examine_vertex = vertex_vis_print_func('examine_vertex');
vis.finish_vertex = vertex_vis_print_func('finish_vertex');
vis.examine_edge = edge_vis_print_func('examine_edge');
vis.edge_relaxed = edge_vis_print_func('edge_relaxed');
vis.edge_not_relaxed = edge_vis_print_func('edge_not_relaxed');

%%
% With the visitor setup, there is hardly any work left.

dijkstra_sp(A,1,struct('visitor', vis));

%% Understanding the output
% To understand the output, we find it helpful to have a copy of
% Introduction to Algorithms by Cormen, Leiserson, and Rivest.  The source
% for the graph is Figure 25-2 in that book and the authors use the graph
% to illustrate how Dijkstra's algorithm runs.  In particular, Figure 25-5
% shows a sample run of Dijkstra's algorithm.
%
% Perhaps the first thing to notice is that the initialize vertex visitor
% is never called.  This results from an error in the MatlabBGL and Boost
% documentation.  Once it is resolved, we will update the MatlabBGL
% documentation to match the Boost graph library.
%
% The results: discover_vertex is called before examine_vertex.  For the
% edges, examine_edge is always called before either edge_relaxed
% or edge_not_relaxed.  The edges that are relaxed are the shaded edges in
% Figure 25-5.
%
% Finally, finish vertex is called on a vertex after all of its edges have
% been examined and possibly relaxed.

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function ex3bvp
%EX3BVP  Example 3 of the BVP tutorial.
%   This is the example for D02KAF from the NAG library.  D02KAF is a code
%   for Sturm-Liouville problems that takes advantage of special features of
%   the problem.  The task is to find the fourth  eigenvalue of Mathieu's
%   equation
%
%      y'' + (lambda -2*q*cos(2*x))*y = 0
%
%   on the interval [0, pi] with boundary conditions y'(0) = 0, y'(pi) = 0.
%   The parameter q = 5.
%
%   A code that exploits fully the special nature of the Sturm-Liouville
%   problem can compute a specific eigenvalue.  Of course using BVP4C we can
%   only compute an eigenvalue near to a guessed value.  We can make it much
%   more likely that we compute the desired eigenfunction by supplying a
%   guess that has the correct qualitative behavior.  We use here the same
%   guess lambda = 15 as the example.  The eigenfunction is determined only
%   to a constant multiple, so the normalizing condition y(0) = 1 is used to
%   specify a particular solution.
%
%   Plotting the solution on the mesh found by BVP4C does not result in a
%   smooth graph.  The solution S(x) is continuous and has a continuous
%   derivative.  It can be evaluated inexpensively using DEVAL at as many
%   points as necessary to get a smooth graph.

% Copyright 2002, The MathWorks, Inc.

% BVPINT is used to form an initial guess for a mesh of 10 equally
% spaced points.  The guess cos(4x) for y(x) and its derivative as guess
% for y'(x) are evaluated in EX3INIT.  The desired eigenvalue is the one
% nearest the guess lambda = 15.  A guess for unknown parameters is the
% (optional) last argument of BVPINIT.
solinit = bvpinit(linspace(0,pi,10),@ex3init,15);

options = bvpset('stats','on');

sol = bvp4c(@ex3ode,@ex3bc,solinit,options);

% BVP4C returns the solution as the structure 'sol'. The computed eigenvalue
% is returned in the field sol.parameters.
fprintf('\n');
fprintf('D02KAF computed lambda = 17.097.\n')
fprintf('bvp4c  computed lambda =%7.3f.\n',sol.parameters)

figure
plot(sol.x,sol.y(1,:),sol.x,sol.y(1,:),'*')
axis([0 pi -1 1])
title('Eigenfunction for Mathieu''s equation.')
xlabel('Solution at mesh points only.')

% Plotting the solution just at the mesh points does not result in a
% smooth graph near the ends of the interval.  The approximate solution
% S(x) is continuous and has a continuous derivative.  DEVAL is used to
% evaluate it at enough points to get a smooth graph.
figure
xint = linspace(0,pi);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
axis([0 pi -1 1])
title('Eigenfunction for Mathieu''s equation.')
xlabel('Solution evaluated on a finer mesh with DEVAL.')

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

function dydx = ex3ode(x,y,lambda)
%EX3ODE  ODE function for Example 3 of the BVP tutorial.
q = 5;
dydx = [              y(2)
         -(lambda - 2*q*cos(2*x))*y(1) ];

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

function res = ex3bc(ya,yb,lambda)
%EX3BC  Boundary conditions for Example 3 of the BVP tutorial.
res = [ ya(2)
        yb(2)
        ya(1) - 1 ];

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

function v = ex3init(x)
%EX3INIT  Guess for the solution of Example 3 of the BVP tutorial.
v = [   cos(4*x)
      -4*sin(4*x) ];


text:
function [varargout]=cunique(A)

% function [A_uni,ind1,ind2,Ac]=cunique(A)
%-------------------------------------------------------------------------
%This function is similar to MATLAB's unique function. There are three
%differences: 1) An additional 4th optional output is available providing
%the count, or number of occurances, for each element in the input array.
%2) The 2nd output mathces the size of the first input, 3) The 3rd output
%is reshaped to be the size of the input variable.
%
% See also: unique
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2018/03/21: Created
%-------------------------------------------------------------------------

%%

[A_uni,ind1,ind2]=unique(A);

varargout{1}=A_uni;
varargout{2}=reshape(ind1,size(A_uni));
varargout{3}=reshape(ind2,size(A));

if nargout==4
    [subInd] = ind2subn(size(A_uni),ind2);
    Ac=accumarray(subInd,ones(numel(ind2),1),size(A_uni));
    Ac=reshape(Ac(ind2),size(A));
    varargout{4}=Ac;
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function jac = p30_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P30_JAC evaluates the jacobian for problem p30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  csum = 0.0;
  for i = 1 : 19
    csum = csum + i.^(4.0/3.0);
  end

  pprime = 0.0;
  for i = 1 : 19
    pprime = pprime + ( 4.0 / 3.0 ) * r8_cube_root ( t - i );
  end

  jac(1,1) = pprime / csum;

  return
end


text:
%% Copyright 2011-2012 The MathWorks, Inc.
function ReplaceCard()
% This is a simple demo to show a match & replace of the video
% data.
%
% To try this,
% 1. Get a deck of playing cards with "MATLAB & Simulink" logo or print
% the file "MWqueen_crop_small.bmp" on a color printer.
% 2. Run this script in MATLAB
% 3. Show printed queen of diamond in front of the usb web camera.
% 4. Press either 1 or 2 of your keyboard to see the image changing
% 5. Press 0 of your keyboard to reset it back to queen of diamonds
%
% Original version created by Takuya Otani
% Senior Application Engineer, MathWorks, Japan
%

clear all; close all; clc; imaqreset;
global process_flag;
process_flag = false;

%% Define Geometric Transformation Objects
gte = vision.GeometricTransformEstimator;
gte.Method = 'Random Sample Consensus (RANSAC)';
gte.Transform = 'Nonreflective similarity';
gte.NumRandomSamplingsMethod = 'Desired confidence';
gte.MaximumRandomSamples = 500;
gte.DesiredConfidence = 99.8;

agt = vision.GeometricTransformer;
agt.OutputImagePositionSource = 'Property';
agt.OutputImagePosition = [1 1 480 640];

halphablender = vision.AlphaBlender( ...
    'Operation', 'Binary mask', 'MaskSource', 'Input port');

%% Setup Image Acquisition
hCamera = videoinput('winvideo', 1, 'RGB24_640x480');
hCamera.ReturnedColorspace = 'rgb';
hCamera.FramesPerTrigger =  1;

% Initialization of the camera
triggerconfig(hCamera, 'manual');
start(hCamera);

%% Load reference image data (queen of diamond) Pad the image data
%% to match 480x640 dimension
ref_img = imread('MWqueen_crop_small.bmp');
ref_img_pad = padarray(ref_img, [140 0], 'pre'); %340x240
ref_img_pad = padarray(ref_img_pad, [0 400], 'pre');
ref_img_gray = rgb2gray(ref_img_pad);
ref_keys = detectSURFFeatures(ref_img_gray, 'MetricThreshold', 1000);
% Extract SURF Features
[ref_features  ref_validPts] = extractFeatures(ref_img_gray,  ref_keys, 'SURFSize', 64);

% Open a Figure window which accepts keyboard inputs
figure('KeyPressFcn', @keycontrol);
himg = imshow(zeros(480, 640, 3, 'uint8'));

%% image data to replace
global rep_img_pad;

while ishandle(himg)

  %% Capture data from the camera, and extrace SURF features, and descriptors
  vid_img = getsnapshot(hCamera);
  vid_img = im2single(vid_img);  vid_img_gray = rgb2gray(vid_img);
  vid_keys = detectSURFFeatures(vid_img_gray, 'MetricThreshold', 1000);
  [vid_features vid_validPts] = extractFeatures(vid_img_gray, vid_keys, 'SURFSize', 64);

  %Here, I am using a MATLAB Coder compiled version of matchFeatures
  %To speed things up.
  %index_pairs = matchFeatures(ref_features, vid_features, 'MatchThreshold', 2);
  index_pairs = matching_fcn_mex(ref_features, vid_features, 2);

  if process_flag == true & size(index_pairs,1) > 60

    ref_matched_points = ref_validPts(index_pairs(:,1));
    vid_matched_points = vid_validPts(index_pairs(:,2));

    % Estimate the transformation matrix
    [tform_matrix inlierIdx] = step(gte, ...
                                    ref_matched_points.Location, ...
                                    vid_matched_points.Location);

    % peform a geometric transformation
    rep_img_trans = step(agt, rep_img_pad, tform_matrix);

    % Replace the card
    mask = rep_img_trans(:,:,1) > 0;
    mosaic = step(halphablender, vid_img, rep_img_trans, mask);

    set(himg, 'CData', mosaic);
  else
    set(himg, 'CData', vid_img);
  end

  drawnow;
end

% When keyboard is pressed, replace rep_img data.
function keycontrol(src,evnt)
global rep_img_pad;
global process_flag;

if evnt.Character == '0'
    process_flag = false;
elseif evnt.Character == '1'
    process_flag = true;
    rep_img = imread('MWJoker_crop_small.bmp');
    rep_img = im2single(rep_img);
    rep_img_pad = padarray(rep_img, [140 0], 'pre');
    rep_img_pad = padarray(rep_img_pad, [0 400], 'pre');
elseif evnt.Character == '2'
    process_flag = true;
    rep_img = imread('MWspade_crop_small.bmp');
    rep_img = im2single(rep_img);
    rep_img_pad = padarray(rep_img, [140 0], 'pre');
    rep_img_pad = padarray(rep_img_pad, [0 400], 'pre');
end


text:
function [y,newp,retcode]=ssfile(obj,y,p,d,id) %#ok<INUSL>

% risessfile --  computes the steady state of ... analytically
%
% ::
%
%
%   [y,newp,retcode]=risessfile(obj,y,p,d,id)
%
% Args:
%
%    - **obj** [rise|dsge]: model object (not always needed)
%
%    - **y** [vector]: endo_nbr x 1 vector of initial steady state
%
%    - **p** [struct]: parameter structure
%
%    - **d** [struct]: definitions
%
%    - **id** [vector]: location of the variables to calculate
%
% Returns:
%    :
%
%    - **y** []: endo_nbr x 1 vector of updated steady state
%
%    - **newp** [struct]: structure containing updated parameters if any
%
%    - **retcode** [0|number]: return 0 if there are no problems, else return
%      any number different from 0
%
% Note:
%
%    - this is new approach has three main advantages relative to the previous
%      one:
%      - The file is valid whether we have many regimes or not
%      - The user does not need to know what regime is being computed
%      - It is in sync with the steady state model
%
% Example:
%
%    See also:

% Conversion of Dynare file [Dynare_edo_steadystate.m] into RISE file [risessfile.m]
%
%  Done 23-Jul-2017 00:31:15.

retcode=0;

if nargin==1
    % list of endogenous variables to be calculated
    %----------------------------------------------
    % y=get(obj,'endo_list');
    y={'RC','RK','WC','WK','YC','YK','MCC','MCK','KC','KK','PKB','R','L','QK',...
        'HC','HSC','HK','HSK','UHC','UHSC','UHK','UHSK','empC','HrC',...
        'empK','HrK','empSC','HrSC','empSK','HrSK','unemp','EIK','EC','INFWC',...
        'INFWK','INFC','INFK','DIFFNORMGDP','NORMINFGDP','DIFFREALGDP',...
        'DIFFREALEC','DIFFREALEIK','DIFFREALW','AH','INFGDP','INFCNA','INFCOR',...
        'GAP','PFGAP','INFC10','ECD','KD','RCD','QCD','KCH','RCH','ECH','QCH',...
        'LAGKD','LAGKCH','UK','UC','DIFFREALECH','DIFFREALECD','betas','XiL',...
        'Lpref','EFFK','MUZK','MUZM','HG','MUC','MUK','EFFECD','EFFECH','STAR',...
        'RL1','RL2','RL3','RL4','RL5','RL6','RL7','RT2','DIFFREALGDP_obs',...
        'DIFFREALEC_obs','DIFFREALEIK_obs','DIFFREALECD_obs','DIFFREALECH_obs',...
        'DIFFREALW_obs','AH_obs','INFCNA_obs','INFCOR_obs','INFK_obs','R_obs',...
        'RT2_obs','unemp_obs'};
    % list of parameters to be computed during steady state calculation
    %-------------------------------------------------------------------
    newp={'AA','AHSS','A_HC','A_HK','DD','DIFFREALECDSS','DIFFREALECDSS_obs',...
        'DIFFREALECHSS','DIFFREALECHSS_obs','DIFFREALECSS','DIFFREALECSS_obs',...
        'DIFFREALEIKSS','DIFFREALEIKSS_obs','DIFFREALGDPSS','DIFFREALGDPSS_obs',...
        'DIFFREALWSS','DIFFREALWSS_obs','ECDSS','ECHSS','ECSS','EIKSS','HCSS',...
        'HKSS','HSCSS','HSKSS','HSS','HrCSS','HrKSS','HrSCSS','HrSKSS',...
        'IMPHSSS','INFC10SS','INFCNASS','INFCNASS_obs','INFCORSS','INFCORSS_obs',...
        'INFCSS','INFGDPSS','INFKSS','INFKSS_obs','INFWCSS','INFWKSS','KCDSS',...
        'KCHSS','KCSS','KKSS','LSS','MCCSS','MCKSS','MUCSS','MUCSShabit','MUKSS',...
        'MUKSShabit','MUZCSS','ONE','PKBSS','PYSS','QCDSS','QCHSS','QKSS',...
        'RCDSS','RCHSS','RCSS','RKSS','RL1SS','RL2SS','RL3SS','RL4SS','RL5SS',...
        'RL6SS','RL7SS','RR','RSS','RSS_obs','RT2SS','RT2SS_obs','Rnr','UCSS',...
        'UHCSS','UHKSS','UHSCSS','UHSKSS','UKSS','USS','WCSS','WKSS','YCSS',...
        'YKSS','YYSS','beta_','beta_0','beta_2','empCSS','empKSS','empSCSS',...
        'empSKSS','eta_cd','eta_cd_eta_cnn','eta_ch','eta_ch_eta_cnn','eta_cnn',...
        'hc_hk','mu_','s_c_ech','s_ecdc','s_k','s_k_ecd','s_k_eik','s_yc',...
        'theta_wc','theta_wk','unempSS','unempSS_obs','xsi_HrC','xsi_HrK',...
        'xsi_NC','xsi_NK','ycbi','ycbi_ykb','ykb'};
else

    newp=struct();

    %start_steady_state;

    newp.beta_0 = p.pbeta;
    newp.beta_2 = p.pbeta*p.rpr; % s.s. funds rate premium
    newp.beta_ = newp.beta_2;
    newp.MUZCSS=1;
    newp.ONE=1;
    newp.USS=1;
    newp.MUKSS=p.MUZKSS*p.MUZMSS;
    newp.MUCSS=p.MUZKSS^p.alpha_*p.MUZMSS;
    newp.MUKSShabit=newp.MUKSS;
    newp.MUCSShabit=newp.MUCSS;
    newp.PKBSS=p.theta_k/(p.theta_k-1)*(p.theta_c-1)/p.theta_c;
    newp.PYSS=1;
    newp.MCCSS=(p.theta_c-1)/p.theta_c;
    newp.MCKSS=(p.theta_k-1)/p.theta_k;
    newp.RKSS=newp.MUKSS/newp.beta_2-(1-p.delta_);
    newp.RCSS=newp.MUKSS/newp.beta_2-(1-p.delta_);
    newp.RCHSS=newp.MUCSS/newp.beta_2-(1-p.delta_ch); % Housing sector
    newp.RCDSS=newp.MUKSS/newp.beta_2-(1-p.delta_cd); % Durable sector
    newp.USS=1;
    newp.mu_=newp.RCSS;
    newp.AA=p.alpha_/newp.RKSS*newp.MCKSS;
    newp.DD = 0.135;
    newp.RR = 0.075;
    newp.eta_cnn=1;
    newp.eta_cd_eta_cnn=newp.DD/((newp.MUKSShabit-newp.beta_2*p.h_cd)/(1-newp.beta_2*p.h/newp.MUCSShabit)*(1-p.h/newp.MUCSShabit)/(1-p.h_cd/newp.MUKSShabit)*(1-(1-p.delta_cd)/newp.MUKSS)/newp.RCDSS);
    newp.eta_ch_eta_cnn=newp.RR/((newp.MUCSShabit-newp.beta_2*p.h_ch)/(1-newp.beta_2*p.h/newp.MUCSShabit)*(1-p.h/newp.MUCSShabit)/(1-p.h_ch/newp.MUCSShabit)*(1-(1-p.delta_ch)/newp.MUCSS)/newp.RCHSS);
    newp.eta_ch=newp.eta_ch_eta_cnn;
    newp.eta_cd=newp.eta_cd_eta_cnn;
    newp.DD=newp.eta_cd_eta_cnn*(newp.MUKSShabit-newp.beta_2*p.h_cd)/(1-newp.beta_2*p.h/newp.MUCSShabit)*(1-p.h/newp.MUCSShabit)/(1-p.h_cd/newp.MUKSShabit)*(1-(1-p.delta_cd)/newp.MUKSS)/newp.RCDSS;
    newp.RR=newp.eta_ch_eta_cnn*(newp.MUCSShabit-newp.beta_2*p.h_ch)/(1-newp.beta_2*p.h/newp.MUCSShabit)*(1-p.h/newp.MUCSShabit)/(1-p.h_ch/newp.MUCSShabit)*(1-(1-p.delta_ch)/newp.MUCSS)/newp.RCHSS;
    newp.Rnr=(1-(1-p.delta_)/newp.MUKSS)*newp.AA*newp.MUKSS;
    newp.ycbi_ykb=((1-p.s_AS)-newp.Rnr)/((newp.DD*(1-p.s_AS)/(1+newp.RR))+newp.Rnr);
    newp.hc_hk=newp.ycbi_ykb*(newp.RCSS*newp.MCKSS/(newp.RKSS*newp.MCCSS))^(p.alpha_/(1-p.alpha_));
    newp.HSS=0.25;
    newp.AHSS=newp.HSS;
    newp.HKSS=newp.HSS/(1+newp.hc_hk);
    newp.HCSS=newp.HSS-newp.HKSS;
    newp.HrCSS=1/3;
    newp.HrKSS=1/3;
    newp.empCSS=newp.HCSS/newp.HrCSS;
    newp.empKSS=newp.HKSS/newp.HrKSS;
    newp.ycbi=newp.HCSS*(newp.AA)^(p.alpha_/(1-p.alpha_));
    newp.ykb=newp.HKSS*(newp.AA)^(p.alpha_/(1-p.alpha_));
    newp.YCSS=newp.ycbi;
    newp.YKSS=newp.ykb;
    newp.KCSS=newp.AA*newp.ycbi*newp.MUKSS;
    newp.KKSS=newp.AA*newp.ykb*newp.MUKSS;
    newp.ECHSS=newp.RR/(1+newp.RR)*newp.ycbi*(1-p.s_AS);
    newp.ECSS=1/(1+newp.RR)*newp.ycbi*(1-p.s_AS);
    newp.ECDSS=newp.DD*newp.PKBSS*newp.ECSS;
    newp.EIKSS=(1-(1-p.delta_)/newp.MUKSS)*(newp.KCSS+newp.KKSS);
    newp.KCDSS=newp.ECDSS/(1-(1-p.delta_cd)/newp.MUKSS);
    newp.KCHSS=newp.ECHSS/(1-(1-p.delta_ch)/newp.MUCSS);
    newp.YYSS=(newp.YCSS+newp.YKSS*newp.PKBSS)/newp.PYSS;
    newp.s_k_ecd=newp.ECDSS/newp.YKSS;
    newp.s_c_ech=newp.ECHSS/newp.YCSS;
    newp.s_k_eik=newp.EIKSS/newp.YKSS;
    newp.s_yc = (newp.YCSS/newp.YYSS);
    newp.s_ecdc=newp.PKBSS*newp.ECDSS/(newp.ECSS+newp.PKBSS*newp.ECDSS+(newp.MUCSS/newp.beta_2-1+p.delta_ch)*newp.KCHSS);
    newp.INFCNASS=exp(.02/4);
    newp.INFCSS = newp.INFCNASS*((newp.MUZCSS/p.MUZKSS)^(1-p.alpha_))^(-newp.s_ecdc);
    newp.INFCORSS=newp.INFCNASS;
    newp.INFKSS=newp.INFCSS*(newp.MUZCSS/p.MUZKSS)^(1-p.alpha_);
    newp.INFWCSS=newp.INFCSS*p.MUZKSS^p.alpha_*p.MUZMSS;
    newp.INFWKSS=newp.INFWCSS;
    newp.RSS=newp.INFCSS/newp.beta_0*newp.MUCSS;
    newp.RT2SS=exp(p.tp2)*newp.RSS;
    newp.INFC10SS = newp.INFCNASS;
    newp.IMPHSSS = newp.RCHSS*newp.KCHSS;
    newp.s_k=newp.PKBSS*newp.YKSS/newp.YYSS;
    newp.INFGDPSS=newp.INFCSS^(newp.YCSS/newp.YYSS)*newp.INFKSS^(newp.YKSS*newp.PKBSS/(newp.YYSS));
    newp.LSS=newp.eta_cnn/(newp.ECSS*(1-p.h/newp.MUCSShabit))-newp.eta_cnn*newp.beta_2*p.h/(newp.ECSS*(newp.MUCSShabit-p.h));
    newp.WCSS=newp.MCCSS*(1-p.alpha_)*newp.YCSS/newp.HCSS;
    newp.WKSS=newp.MCKSS*(1-p.alpha_)*newp.YKSS/newp.HKSS;
    % xsiN_xsiH_C = ((newp.HrCSS/newp.empCSS)^(1+p.sigmah))/(1+1/p.sigmah);
    % xsiN_xsiH_K = ((newp.HrKSS/newp.empKSS)^(1+p.sigmah))/(1+1/p.sigmah);
    % gC = (1/(1+p.sigman) + 1/p.sigmah)*(xsiN_xsiH_C*(1+p.sigmah)/p.sigmah)^(-(1+p.sigman)/(1+p.sigman+p.sigmah));
    % markup_xsiN_C = (newp.HCSS^((1+p.sigmah)*(1+p.sigman)/(1+p.sigmah+p.sigman)-1))*gC/(newp.LSS*newp.WCSS);
    % gK = (1/(1+p.sigman) + 1/p.sigmah)*(xsiN_xsiH_K*(1+p.sigmah)/p.sigmah)^(-(1+p.sigman)/(1+p.sigman+p.sigmah));
    % markup_xsiN_K = (newp.HKSS^((1+p.sigmah)*(1+p.sigman)/(1+p.sigmah+p.sigman)-1))*gK/(newp.LSS*newp.WKSS);
    % markup_w = (1-newp.unempSS)^((1+p.sigmah+p.sigman)/(1+p.sigmah) - 1 - p.sigman);
    markup_w = (1-p.unempSS)^((1+p.sigmah+p.sigman)/(1+p.sigmah) - 1 - p.sigman);
    newp.theta_wc = markup_w/(markup_w -1); newp.theta_wk = newp.theta_wc;
    newp.A_HC=newp.LSS*(newp.theta_wc-1)/newp.theta_wc*newp.WCSS/(((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))*newp.HCSS^(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah))));
    newp.A_HK=newp.LSS*(newp.theta_wk-1)/newp.theta_wk*newp.WKSS/(((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))*newp.HKSS^(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah))));
    newp.xsi_NC=newp.A_HC/((1/(1+p.sigman)+1/p.sigmah)*(newp.HCSS^p.sigman/newp.HrCSS^(1+p.sigman+p.sigmah))^((1+p.sigman)/(1+p.sigman+p.sigmah)));
    newp.xsi_NK=newp.A_HK/((1/(1+p.sigman)+1/p.sigmah)*(newp.HKSS^p.sigman/newp.HrKSS^(1+p.sigman+p.sigmah))^((1+p.sigman)/(1+p.sigman+p.sigmah)));
    newp.xsi_HrC=newp.xsi_NC*(1+p.sigmah)/p.sigmah*(newp.HCSS^p.sigman/newp.HrCSS^(1+p.sigman+p.sigmah));
    newp.xsi_HrK=newp.xsi_NK*(1+p.sigmah)/p.sigmah*(newp.HKSS^p.sigman/newp.HrKSS^(1+p.sigman+p.sigmah));
    newp.UHCSS=newp.A_HC*((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))*newp.HCSS^(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah)))/newp.LSS;
    newp.UHKSS=newp.A_HK*((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))*newp.HKSS^(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah)))/newp.LSS;
    newp.HSCSS=(newp.WCSS*newp.LSS/(newp.A_HC*((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))))^(1/(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah))));
    newp.HSKSS=(newp.WKSS*newp.LSS/(newp.A_HK*((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))))^(1/(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah))));
    newp.empSCSS=((1+p.sigmah)/p.sigmah*newp.xsi_NC/newp.xsi_HrC)^(-1/(1+p.sigmah+p.sigman))*newp.HSCSS^(1/(1+p.sigman/(1+p.sigmah)));
    newp.empSKSS=((1+p.sigmah)/p.sigmah*newp.xsi_NK/newp.xsi_HrK)^(-1/(1+p.sigmah+p.sigman))*newp.HSKSS^(1/(1+p.sigman/(1+p.sigmah)));
    newp.HrSCSS=newp.HSCSS/newp.empSCSS;
    newp.HrSKSS=newp.HSKSS/newp.empSKSS;
    newp.UHSCSS=newp.A_HC*((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))*newp.HSCSS^(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah)))/newp.LSS;
    newp.UHSKSS=newp.A_HK*((1+p.sigman)/(1+p.sigman/(1+p.sigmah)))*newp.HSKSS^(-1+(1+p.sigman)/(1+p.sigman/(1+p.sigmah)))/newp.LSS;
    newp.unempSS=(newp.empSCSS+newp.empSKSS-(newp.empCSS+newp.empKSS))/(newp.empSCSS+newp.empSKSS);
    newp.QKSS=1;
    newp.QCDSS=1;
    newp.QCHSS=1;
    newp.UCSS=1;
    newp.UKSS=1;
    % XiBSS=1;
    % XiDSS=1;
    % XiHSS=1;
    newp.RL1SS=newp.RSS;
    newp.RL2SS=newp.RSS;
    newp.RL3SS=newp.RSS;
    newp.RL4SS=newp.RSS;
    newp.RL5SS=newp.RSS;
    newp.RL6SS=newp.RSS;
    newp.RL7SS=newp.RSS;
    newp.DIFFREALECSS =exp( log(newp.MUCSS));
    newp.DIFFREALEIKSS =exp( log(newp.MUKSS));
    newp.DIFFREALECDSS =exp( log(newp.MUKSS));
    newp.DIFFREALECHSS =exp( log(newp.MUCSS));
    newp.DIFFREALWSS =exp( log(newp.MUCSS) );
    newp.DIFFREALGDPSS =exp( (1-newp.s_k)*log(newp.MUCSS)+(newp.s_k)*log(newp.MUKSS));

    %end_steady_state;

    %trends;

    newp.DIFFREALGDPSS_obs=(1-newp.s_k)*log(newp.MUCSS)*100+(newp.s_k)*log(newp.MUKSS)*100;
    newp.DIFFREALECSS_obs=log(newp.MUCSS)*100;
    newp.DIFFREALEIKSS_obs=log(newp.MUKSS)*100;
    newp.DIFFREALECDSS_obs=log(newp.MUKSS)*100;
    newp.DIFFREALECHSS_obs=log(newp.MUCSS)*100;
    newp.DIFFREALWSS_obs=log(newp.MUCSS)*100;
    newp.INFCNASS_obs=(1-newp.s_ecdc)*log(newp.INFCSS)*100+newp.s_ecdc*log(newp.INFKSS)*100;
    newp.INFCORSS_obs=(1-newp.s_ecdc)*log(newp.INFCSS)*100+newp.s_ecdc*log(newp.INFKSS)*100;
    newp.INFKSS_obs=log(newp.INFCSS)*100-log(newp.MUKSS)*100+log(newp.MUCSS)*100;
    newp.RSS_obs=log(newp.RSS)*100;
    newp.RT2SS_obs=log(newp.RT2SS)*100;
    newp.unempSS_obs=100*log(newp.unempSS);

    %end_trends;

    ys = [
        newp.RCSS
        newp.RKSS
        newp.WCSS
        newp.WKSS
        newp.YCSS
        newp.YKSS
        newp.MCCSS
        newp.MCKSS
        newp.KCSS
        newp.KKSS
        newp.PKBSS
        newp.RSS
        newp.LSS
        newp.QKSS
        newp.HCSS
        newp.HSCSS
        newp.HKSS
        newp.HSKSS
        newp.UHCSS
        newp.UHSCSS
        newp.UHKSS
        newp.UHSKSS
        newp.empCSS
        newp.HrCSS
        newp.empKSS
        newp.HrKSS
        newp.empSCSS
        newp.HrSCSS
        newp.empSKSS
        newp.HrSKSS
        newp.unempSS
        newp.EIKSS
        newp.ECSS
        newp.INFWCSS
        newp.INFWKSS
        newp.INFCSS
        newp.INFKSS
        newp.ONE
        newp.ONE
        newp.DIFFREALGDPSS
        newp.DIFFREALECSS
        newp.DIFFREALEIKSS
        newp.DIFFREALWSS
        newp.AHSS
        newp.INFGDPSS
        newp.INFCNASS
        newp.INFCORSS
        newp.ONE
        newp.ONE
        newp.INFC10SS
        newp.ECDSS
        newp.KCDSS
        newp.RCDSS
        newp.QCDSS
        newp.KCHSS
        newp.RCHSS
        newp.ECHSS
        newp.QCHSS
        newp.KCDSS
        newp.KCHSS
        newp.USS
        newp.USS
        newp.DIFFREALECHSS
        newp.DIFFREALECDSS
        newp.beta_
        newp.ONE
        newp.ONE
        newp.ONE
        p.MUZKSS
        p.MUZMSS
        newp.ONE
        newp.MUCSS
        newp.MUKSS
        newp.ONE
        newp.ONE
        newp.ONE
        newp.RL1SS
        newp.RL2SS
        newp.RL3SS
        newp.RL4SS
        newp.RL5SS
        newp.RL6SS
        newp.RL7SS
        newp.RT2SS
        newp.DIFFREALGDPSS_obs
        newp.DIFFREALECSS_obs
        newp.DIFFREALEIKSS_obs
        newp.DIFFREALECDSS_obs
        newp.DIFFREALECHSS_obs
        newp.DIFFREALWSS_obs
        newp.ONE
        newp.INFCNASS_obs
        newp.INFCORSS_obs
        newp.INFKSS_obs
        newp.RSS_obs
        newp.RT2SS_obs
        newp.unempSS_obs
        ];

    y(id)=ys;

end

end

author: jmaih


text:
function [success] = manageParallelWorkers(varargin)
	% Manages loading and stopping parallel processing workers.
	% Biafra Ahanonu
	% started: 2015.12.01

	% changelog
		% 2022.02.28 [18:36:15] - Added ability to input just the number of workers to open as 1st single input argument that aliases for the "setNumCores" Name-Value input, still support other input arguments as well.

	if length(varargin)==1
		[success] = ciapkg.io.manageParallelWorkers(varargin{1});
	else
		[success] = ciapkg.io.manageParallelWorkers('passArgs', varargin);
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [DEMc,MASK] = crop(DEM,varargin)

%CROP crop an instance of GRIDobj with axis-aligned minimum bounding box
%
% Syntax
%
%     DEMc = crop(DEM)
%     DEMc = crop(DEM,I)
%     DEMc = crop(DEM,I,fillval)
%     DEMc = crop(DEM,ix)
%     DEMc = crop(DEM,x,y)
%     DEMc = crop(DEM,'interactive');
%
% Description
%
%     crop removes outer parts of a grid based on missing values (NaNs),
%     a mask, linear indices or coordinate pairs such that the axis-aligned
%     minimum bounding rectangle remains.
%
% Input arguments
%
%     DEM      instance of GRIDobj
%     I        mask (logical) as instance of GRIDobj
%     fillval  scalar; values in the output grid where I.Z==0 obtain this
%              value (crop and clip by setting fillval to nan)
%     ix       linear index into the DEM
%     x,y      coordinate vectors
%
%
% See also: IND2SUB, GRIDobj/pad, GRIDobj/resample
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 2. September, 2017

narginchk(1,3);

if nargin == 1;
    MASK = isnan(DEM.Z);

    if ~any(MASK(:));
        DEMc = DEM;
        if nargout == 2;
            MASK = DEM;
            MASK.Z = true(DEM.size);
            MASK.zunit = '';
            MASK.name = 'mask for cropping';
        end

        return
    end
    MASK = ~MASK;
    MASK = bwperim(MASK);
    IX  = find(MASK);

elseif nargin >= 2;
    if isa(varargin{1},'GRIDobj') || isa(varargin{1},'logical');
        % GRIDobj
        validatealignment(DEM,varargin{1});
        if isa(varargin{1},'GRIDobj')
            MASK = logical(varargin{1}.Z);
        else
            MASK = varargin{1};
        end

        if nargin == 3
            validateattributes(varargin{2},{'numeric'},{'scalar'},'crop','fillval',3);
            if isnan(varargin{2}) && ~(isa(DEM.Z,'single') || isa(DEM.Z,'double'))
                DEM.Z(~MASK) = 0;
                warning('TopoToolbox:GRIDobj',...
                    ['fillval set to zero since data type of the input grid\newline' ...
                     'does not support nans'])
            else
                DEM.Z(~MASK) = varargin{2};
            end
        end
        MASK = bwperim(MASK);
        IX  = find(MASK);
    elseif nargin == 2;
        if isnumeric(varargin{1})
            % indices are supplied
            IX  = varargin{1};
            if numel(IX) < 2;
                error('TopoToolbox:GRIDobj',...
                    'At least two indices are required to crop the grid.')
            end

            if any(IX<1) || any(IX>prod(DEM.size));
                error('TopoToolbox:GRIDobj',...
                    ['Index must range between 1 and ' num2str(prod(DEM.size)) '.'])
            end
        else
            % interactive part
            try
            imagesc(DEM);
            c = uicontrol('Style','Text','Units','normalized','position',[0 0 1 0.05],...
                'String','Draw rectangle and double click when finished.');
            h = imrect;
            addNewPositionCallback(h,@(pos) set(c,'String',...
                ['LX:' num2str(round(pos(1)),'%d') ', LY:' num2str(round(pos(2)),'%d') ...
                 ', UX:' num2str(round(pos(1)+pos(3)),'%d') ', UY:' num2str(round(pos(2)+pos(4)),'%d')]));
            [~] = wait(h);
            MASK  = createMask(h);
            delete(h)
            close
            catch ME
                error('TopoToolbox:crop','output variable undefined')
%                 warning('empty matrix returned')
%                 DEMc = [];
                return
            end

            MASK  = bwperim(MASK);
            IX    = find(MASK);
        end

    else

        [X,Y] = refmat2XY(DEM.refmat,DEM.size);
        x = varargin{1};
        y = varargin{2};

        x = max(min(X),x);
        x = min(max(X),x);

        y = max(min(Y),y);
        y = min(max(Y),y);

        IX = coord2ind(X,Y,x,y);
    end

end

if nargout == 2
    MASK = DEM;
    MASK.Z = false(DEM.size);
    MASK.Z(IX) = true;
    MASK.name = 'mask for cropping';
    MASK.zunit = '';
end

% nr of dimensions
siz  = DEM.size;
k    = [1 cumprod(siz(1:end-1))];

% preallocate subsref structure
S    = substruct('()',cell(1,2));

% subset size
sizout = zeros(1,2);

% loop through dimensions (see ind2sub)
% and get subscripts of minimum bounding rectangle/box/...
for r = 2:-1:1
    IX2       = rem(IX-1,k(r))+1;
    subdim    = (IX-IX2)/k(r)+1;
    S.subs{r} = min(subdim):max(subdim);
    sizout(r) = numel(S.subs{r});
    IX        = IX2;
end

DEMc      = DEM;
DEMc.Z    = reshape(subsref(DEM.Z,S),sizout);
DEMc.size = sizout;
[x,y]     = sub2coord(DEM,S.subs{1}(1),S.subs{2}(1));
DEMc.refmat(3,:) = [x y] - DEMc.cellsize*[1 -1];

DEMc.name   = [DEM.name ' (cropped)'];
DEMc.zunit  = DEM.zunit;
DEMc.xyunit = DEM.xyunit;

if ~isempty(DEM.georef)
    % Copy all referencing information
    DEMc.georef = DEM.georef;
    DEMc.georef.SpatialRef = refmatToMapRasterReference(DEMc.refmat, DEMc.size);

end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
% Copyright (C) 2008  David Bateman
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, see <http://www.gnu.org/licenses/>.
%
% -*- texinfo -*-
% @deftypefn {Function File} {@var{w} =} window (@var{f}, @var{n}, @var{opts})
% Create a @var{n}-point windowing from the function @var{f}. The
% function @var{f} can be for example @code{@@blackman}. Any additional
% arguments @var{opt} are passed to the windowing function.
% @end deftypefn

function wout = window (f, n, varargin)
  if (nargin == 0)
    error ('window: UI tool not supported');
  elseif (nargin > 1)
    w = feval (f, n, varargin{:});
    if (nargout > 0)
      wout = w;
    end % if
  else
    help(mfilename);
  end % if
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Example
% PTP motion conditional torques, KUKA iiwa 7 R 800

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise matlab toolbox

% First start the server on the KUKA iiwa controller
% Then run this script using Matlab

% Copyright: Mohammad SAFEEA, 8th of April 2018

% Important: Be careful when runnning the script,

close all,clear all;clc;

warning('off')

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  disp('Connection could not be establised, script aborted');
  return;
end
disp('An example on using the interruptible PTP motion functions for Human robot collaboration')
%% Move PTP in joint space to some inital position
jPos={0,0,0,-pi/2,0,pi/2,0};
relVel=0.15;
movePTPJointSpace( t_Kuka , jPos, relVel); % move to initial configuration
%% Move PTP in joint space with torque limits
jPos{1}=pi/2;
joints_indices=[1,4]; % joint index start from one
max_torque=[5,5];
min_torque=[-5,-5];
res=movePTP_ConditionalTorque_JointSpace( t_Kuka...
    , jPos, relVel,joints_indices,max_torque,min_torque);
if(res==1)
    disp('Motion completed successfully');
elseif(res==0)
    disp('Motion interrupted');
elseif(res==-1)
    disp('Error while performing the motion');
end
%% Turn off the server
net_turnOffServer( t_Kuka );
fclose(t_Kuka);

repo: KST-Kuka-Sunrise-Toolbox


text:
function test_pull1377

% MEM 12gb
% WALLTIME 1:30:00
% DEPENDENCY ft_prepare_sourcemodel headsurface ft_prepare_leadfield ft_freqanalysis ft_sourceanalysis

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

% This function creates a set of input-structures to be used for testing source
% analysis steps prior to inverse solution. The structure of this script is more or less
% 1. get the data (MEG and EEG)
% 2. preprocess the data (MEG and EEG)
%   a. timelock with covariance for lcmv and mne
%   b. timelock without keeptrials for dipolefitting
%   c. freq for dics and pcc
% 3. create leadfield
%   a. internally computed leadfield (prepare volume conductor first and leadfield afterwards (MEG localspheres, EEG single sphere))
%   b. load externally computed leadfield and convert to appropriate format
% 4. perform source analysis (MEG and EEG): LCMV, MNE, DICS, PCC, dipolefitting

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In general, variables in one function workspace are not available to other
% functions. However, nested functions can access and modify variables in the
% workspaces of the functions that contain them.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% prevent errors from cfg.mne.keepleadfield, etc
global ft_default
ft_default.checkconfig = 'loose';

datameg = []; % touch these variables to make them shared between the main function and the nexted function
dataeeg = [];

loading_raw_data

headmodel_megsphere = [];
headmodel_eegsphere = [];
sourcemodel_megsphere = [];
sourcemodel_eegsphere = [];
sourcemodel_precomputed_meg = [];
sourcemodel_precomputed_eeg = [];
headmodel_simbio = [];
sourcemodel_simbio = [];
headmodel_openmeeg = [];
sourcemodel_openmeeg = [];

create_leadfields

MEG_tlck = [];
MEG_tlck_df = [];
MEG_freq = [];
EEG_tlck = [];
EEG_tlck_df = [];
EEG_freq = [];

preprocessing

sourcemodel_meg = sourcemodel_megsphere;
headmodel_meg = headmodel_megsphere;
sourcemodel_eeg = sourcemodel_eegsphere;
headmodel_eeg = headmodel_eegsphere;
perform_source_analysis

sourcemodel_meg = sourcemodel_precomputed_meg;
headmodel_meg = [];
sourcemodel_eeg = sourcemodel_precomputed_eeg;
headmodel_eeg = [];
perform_source_analysis

sourcemodel_meg = [];
headmodel_meg = [];
sourcemodel_eeg = sourcemodel_simbio;
headmodel_eeg = headmodel_simbio;
perform_source_analysis

sourcemodel_meg = [];
headmodel_meg = [];
sourcemodel_eeg = sourcemodel_openmeeg;
headmodel_eeg = headmodel_openmeeg;
perform_source_analysis

% !! create worse-case scenario, whereby order and nr of chans don't match across inputs (grad/elec and leadfields)
% Step 1: remove 2-3 random chans from both MEG and EEG raw data
% Step 2: shuffle the order of the channels (note that this breaks the mapping between channels and data)
cfg = [];
cfg.channel = randperm(length(datameg.label)-2);
datameg = ft_selectdata(cfg, datameg);
datameg.label = datameg.label(randperm(length(datameg.label)));

cfg.channel = randperm(length(dataeeg.label)-2);
dataeeg = ft_selectdata(cfg, dataeeg);
dataeeg.label = dataeeg.label(randperm(length(dataeeg.label)));

preprocessing

sourcemodel_meg = sourcemodel_megsphere;
headmodel_meg = headmodel_megsphere;
sourcemodel_eeg = sourcemodel_eegsphere;
headmodel_eeg = headmodel_eegsphere;
perform_source_analysis

sourcemodel_meg = sourcemodel_precomputed_meg;
headmodel_meg = [];
sourcemodel_eeg = sourcemodel_precomputed_eeg;
headmodel_eeg = [];
perform_source_analysis

sourcemodel_meg = [];
headmodel_meg = [];
sourcemodel_eeg = sourcemodel_simbio;
headmodel_eeg = headmodel_simbio;
perform_source_analysis

sourcemodel_meg = [];
headmodel_meg = [];
sourcemodel_eeg = sourcemodel_openmeeg;
headmodel_eeg = headmodel_openmeeg;
perform_source_analysis

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. get the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function loading_raw_data

    clc
    disp('1: Get the raw data')

    % for MEG data + sensor info

    load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf151.mat'), 'data');
    datameg = data;
    clear data

    % for EEG data + channel info

    cfg = [];
    cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/preprocessing_erp/s04.eeg');
    hdr = ft_read_header(cfg.dataset);
    event = ft_read_event(cfg.dataset);
    EVsample = [event.sample]';
    EVvalue = {event.value}';
    Word = find(strcmp('S141', EVvalue)==1);

    begsample = EVsample(Word) - 10;
    endsample = EVsample(Word) + 500;
    offset = -100*ones(size(endsample));
    cfg.trl = [begsample endsample offset];

    dataeeg = ft_preprocessing(cfg);

    % Fix EEG channel names, but note that the actual recording was not done with a Biosemi system
    cfg = [];
    cfg.layout = 'biosemi64.lay';
    lay = ft_prepare_layout(cfg);
    dataeeg.label = lay.label(1:end-2);
    clear lay

    dataeeg.elec = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/template/electrode/standard_1020.elc'), 'senstype', 'eeg');

  end % loding raw data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. preprocess the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function preprocessing

    clc
    disp('2: preprocess the data')

    % create timelock structure with covariance for lcmv and mne
    % and timeloch without keeptrials for dipolefitting

    % for MEG
    cfg = [];
    cfg.covariance = 'yes';
    cfg.keeptrials = 'yes';
    cfg.channel = 'MEG';
    MEG_tlck = ft_timelockanalysis(cfg, datameg);

    cfg = [];
    cfg.covariance = 'yes';
    % cfg.keeptrials = 'yes'; %if this is not commented, the .avg field necessary in dipolefitting is missing
    cfg.channel = 'MEG';
    MEG_tlck_df = ft_timelockanalysis(cfg, datameg);

    % for EEG
    cfg = [];
    cfg.covariance = 'yes';
    cfg.keeptrials = 'yes';
    EEG_tlck = ft_timelockanalysis(cfg, dataeeg);

    % for EEG (dimension problems in dipolefitting!!)
    cfg = [];
    cfg.covariance = 'yes';
    % cfg.keeptrials = 'yes';
    EEG_tlck_df = ft_timelockanalysis(cfg, dataeeg);

    % create freq structure for dics beamforming and pcc beamforming
    % for MEG
    cfg = [];
    cfg.method = 'mtmfft';
    cfg.output = 'fourier';
    cfg.tapsmofrq = 4;
    cfg.foilim = [0 20];
    cfg.channel = 'MEG';
    MEG_freq = ft_freqanalysis(cfg, datameg);

    % for EEG
    cfg = [];
    cfg.method = 'mtmfft';
    cfg.output = 'fourier';
    cfg.tapsmofrq = 4;
    cfg.foilim = [0 20];
    EEG_freq = ft_freqanalysis(cfg, dataeeg);

  end % function preprocessing

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3. create leadfield
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function create_leadfields

    clc
    disp('3: create leadfield')

    %% 3.a internally

    % get volume conductor model
    % for MEG, localsphere
    headmodel_megsphere = ft_read_headmodel(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/vol/Subject01vol_localspheres.mat'));

    % for EEG, singlesphere
    headmodel_eegsphere.o = [0 0 2];
    headmodel_eegsphere.r = 8;
    headmodel_eegsphere.c = 1;
    headmodel_eegsphere.unit = 'cm';

    % create leadfield (internally to ft)
    % for MEG
    cfg = [];
    cfg.grad = datameg.grad;
    cfg.headmodel = headmodel_megsphere;
    cfg.channel = 'MEG';
    cfg.resolution = 1.5;
    sourcemodel_megsphere = ft_prepare_leadfield(cfg);

    cfg = [];
    cfg.elec = dataeeg.elec;
    cfg.headmodel = headmodel_eegsphere;
    cfg.channel = 'EEG';
    cfg.resolution = 1.5;
    sourcemodel_eegsphere = ft_prepare_leadfield(cfg); % inconsistent with dataeeg.labels

    %% 3.b externally (mimic externally created leadfields)
    clc;
    disp('3b: Externally computed leadfields for EEG and MEG')

    Ndip = 1000;

    % for MEG
    % here we are assuming that the externally computed leadfield has the structure Ndip x Nsens
    leadf = randn(3*Ndip, length(datameg.grad.label));

    sourcemodel_precomputed_meg = [];
    sourcemodel_precomputed_meg.leadfield = cell(Ndip,1);
    sourcemodel_precomputed_meg.unit = 'cm';
    sourcemodel_precomputed_meg.leadfielddimord = '{pos}_chan_ori';
    sourcemodel_precomputed_meg.label = datameg.grad.label;

    for d =1:Ndip
      sourcemodel_precomputed_meg.leadfield{d} = [leadf(d,:); leadf(d+Ndip,:); leadf(d+2*Ndip,:)]';
    end

    sourcemodel_precomputed_meg.inside = ones(size(sourcemodel_precomputed_meg.leadfield));
    sourcemodel_precomputed_meg.pos = randn(length(sourcemodel_precomputed_meg.leadfield),3);

    % for EEG
    % here we are assuming that the externally computed leadfield has the structure Ndip x Nsens
    leadf = randn(3*Ndip, length(dataeeg.elec.label));

    sourcemodel_precomputed_eeg = [];
    sourcemodel_precomputed_eeg.leadfield = cell(Ndip,1);
    sourcemodel_precomputed_eeg.unit = 'cm';
    sourcemodel_precomputed_eeg.leadfielddimord = '{pos}_chan_ori';
    sourcemodel_precomputed_eeg.label = dataeeg.elec.label;

    for d =1:Ndip
      sourcemodel_precomputed_eeg.leadfield{d,1} = [leadf(d,:); leadf(d+Ndip,:); leadf(d+2*Ndip,:)]';
    end

    sourcemodel_precomputed_eeg.inside = ones(size(sourcemodel_precomputed_eeg.leadfield));
    sourcemodel_precomputed_eeg.pos = randn(length(sourcemodel_precomputed_eeg.leadfield),3);

    %% 3.c simbio leadfield for EEG

    segm = [];
    segm.brain = logical(floor(3*rand(10,10,10)));
    segm.dim = size(segm.brain);
    segm.unit = 'mm';
    segm.coordsys = 'ctf';
    segm.transform = eye(4);

    cfg = [];
    cfg.method = 'hexahedral';
    cfg.resolution = 1; % this is in mm
    mesh_vol = ft_prepare_mesh(cfg, segm);

    ft_plot_mesh(mesh_vol, 'surfaceonly', false);
    view(120, 30)

    cfg = [];
    cfg.method ='simbio';
    cfg.conductivity = 0.33; % order follows mesh.tissuelabel
    headmodel_simbio = ft_prepare_headmodel(cfg, mesh_vol);

    cfg = [];
    cfg.elec = dataeeg.elec;
    cfg.headmodel = headmodel_simbio;
    cfg.channel = 'EEG';
    cfg.resolution = 1.5;
    sourcemodel_simbio = ft_prepare_leadfield(cfg); % inconsistent with dataeeg.labels

    %% 3.c OpenMEEG leadfield for EEG

    load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/SubjectSEF_vol.mat'), 'vol')
    bnd_openmeeg = vol;

    cfg = [];
    cfg.conductivity = 0.33;
    cfg.method = 'openmeeg';
    headmodel_openmeeg = ft_prepare_headmodel(cfg, bnd_openmeeg);

    cfg = [];
    cfg.elec = dataeeg.elec;
    cfg.headmodel = headmodel_openmeeg;
    cfg.channel = 'EEG';
    cfg.resolution = 1.5;
    sourcemodel_openmeeg = ft_prepare_leadfield(cfg); % inconsistent with dataeeg.labels

  end % function create_leadfields

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4. perform source analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function perform_source_analysis

    clc;
    disp('4: inverse solution, ft_sourceanalysis and ft_dipolefitting')

    %% for MEG
    if ~isempty(sourcemodel_meg) || ~isempty(headmodel_meg)

      % do LCMV beamforming
      cfg = [];
      cfg.method = 'lcmv';
      cfg.keepleadfield = 'yes';
      cfg.lcmv.keepfilter = 'yes';
      cfg.lcmv.keepcov = 'yes';
      cfg.lcmv.lambda = '5%';
      cfg.sourcemodel = sourcemodel_meg;
      cfg.headmodel = headmodel_meg;
      sourcelcmv3d1 = ft_sourceanalysis(cfg, MEG_tlck);
      sourcelcmv3d1 = ft_datatype_source(sourcelcmv3d1, 'version', 'latest');

      % project through computed filter
      cfg.rawtrial = 'yes';
      cfg.sourcemodel = keepfields(sourcelcmv3d1, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'});
      ft_sourceanalysis(cfg, MEG_tlck);

      % do MNE
      cfg = [];
      cfg.method = 'mne';
      cfg.mne.keepleadfield = 'yes';
      cfg.mne.keepfilter = 'yes';
      cfg.mne.lambda = 1e4;
      cfg.sourcemodel = sourcemodel_meg;
      cfg.headmodel = headmodel_meg;
      sourcemne3d1 = ft_sourceanalysis(cfg, MEG_tlck);
      sourcemne3d1 = ft_datatype_source(sourcemne3d1, 'version', 'latest');

      cfg.rawtrial = 'yes';
      cfg.sourcemodel = keepfields(sourcemne3d1, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'});
      ft_sourceanalysis(cfg, MEG_tlck);

      % do DICS
      cfg = [];
      cfg.method = 'dics';
      cfg.dics.keepfilter = 'yes';
      cfg.dics.keepleadfield = 'yes';
      cfg.dics.keepcsd = 'yes';
      cfg.dics.lambda = '5%';
      cfg.frequency = 10;
      cfg.sourcemodel = sourcemodel_meg;
      cfg.headmodel = headmodel_meg;
      sourcedics3d1 = ft_sourceanalysis(cfg, MEG_freq);
      sourcedics3d1 = ft_datatype_source(sourcedics3d1, 'version', 'latest');

      cfg.rawtrial = 'yes';
      cfg.sourcemodel = keepfields(sourcedics3d1, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'});
      ft_sourceanalysis(cfg, MEG_freq);

      % do PCC
      cfg = [];
      cfg.method = 'pcc';
      cfg.pcc.keepfilter = 'yes';
      cfg.pcc.keepleadfield = 'yes';
      cfg.pcc.keepcsd = 'yes';
      cfg.pcc.keepmom = 'yes';
      cfg.pcc.lambda = '5%';
      cfg.frequency = 10;
      cfg.sourcemodel = sourcemodel_meg;
      cfg.headmodel = headmodel_meg;
      ft_sourceanalysis(cfg, MEG_freq);

      % do dipolefit
      cfg = [];
      cfg.numdipoles = 1;
      % cfg.nonlinear = 'no';
      cfg.sourcemodel = sourcemodel_meg;
      cfg.headmodel = headmodel_meg;
      cfg.latency = 0.025;
      cfg.channel = [1 3:151];
      ft_dipolefitting(cfg, MEG_tlck_df);
    end

    %% EEG
    if ~isempty(sourcemodel_eeg) || isempty(headmodel_eeg)

      % do LCMV beamforming
      cfg = [];
      cfg.method = 'lcmv';
      cfg.lcmv.keepleadfield = 'yes';
      cfg.lcmv.keepfilter = 'yes';
      cfg.lcmv.keepcov = 'yes';
      cfg.lcmv.lambda = '5%';
      cfg.sourcemodel = sourcemodel_eeg;
      cfg.headmodel = headmodel_eeg;
      sourcelcmv3d1 = ft_sourceanalysis(cfg, EEG_tlck);
      sourcelcmv3d1 = ft_datatype_source(sourcelcmv3d1, 'version', 'latest');

      % project through computed filter
      cfg.rawtrial = 'yes';
      cfg.sourcemodel = keepfields(sourcelcmv3d1, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'});
      ft_sourceanalysis(cfg, EEG_tlck);

      % do MNE
      cfg = [];
      cfg.method = 'mne';
      cfg.mne.keepleadfield = 'yes';
      cfg.mne.keepfilter = 'yes';
      cfg.mne.lambda = 1e4;
      cfg.sourcemodel = sourcemodel_eeg;
      cfg.headmodel = headmodel_eeg;
      sourcemne3d1 = ft_sourceanalysis(cfg, EEG_tlck);
      sourcemne3d1 = ft_datatype_source(sourcemne3d1, 'version', 'latest');

      cfg.rawtrial = 'yes';
      cfg.sourcemodel = keepfields(sourcemne3d1, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'});
      ft_sourceanalysis(cfg, EEG_tlck);

      % do DICS
      cfg = [];
      cfg.method = 'dics';
      cfg.dics.keepfilter = 'yes';
      cfg.dics.keepleadfield = 'yes';
      cfg.dics.keepcsd = 'yes';
      cfg.dics.lambda = '5%';
      cfg.frequency = 10;
      cfg.sourcemodel = sourcemodel_eeg;
      cfg.headmodel = headmodel_eeg;
      sourcedics3d1 = ft_sourceanalysis(cfg, EEG_freq);

      cfg.rawtrial = 'yes';
      cfg.sourcemodel = keepfields(sourcedics3d1, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'});
      ft_sourceanalysis(cfg, EEG_freq);

      % do PCC
      cfg = [];
      cfg.method = 'pcc';
      cfg.pcc.keepfilter = 'yes';
      cfg.pcc.keepleadfield = 'yes';
      cfg.pcc.keepcsd = 'yes';
      cfg.pcc.keepmom = 'yes';
      cfg.pcc.lambda = '5%';
      cfg.frequency = 10;
      cfg.sourcemodel = sourcemodel_eeg;
      cfg.headmodel = headmodel_eeg;
      ft_sourceanalysis(cfg, EEG_freq);

      % do dipolefit
      cfg = [];
      cfg.numdipoles = 1;
      % cfg.nonlinear = 'no';
      cfg.sourcemodel = sourcemodel_eeg;
      cfg.headmodel = headmodel_eeg;
      cfg.latency = 0.025;
      ft_dipolefitting(cfg, EEG_tlck_df);
    end

  end % function perform_source_analysis

end % main function