text:
function [rec,prec,ap] = wsddnVOCevaldet(VOCopts,cls,res,draw)

% load test set
tic;
VOCopts.annocachepath=[VOCopts.localdir '%s_anno_cache.mat'];
cp=sprintf(VOCopts.annocachepath,VOCopts.testset);
if exist(cp,'file')
  fprintf('%s: pr: loading ground truth\n',cls);
  load(cp,'gtids','recs');
else
  [gtids,t]=textread(sprintf(VOCopts.imgsetpath,VOCopts.testset),'%s %d');
  for i=1:length(gtids)
    % display progress
    if toc>1
      fprintf('%s: pr: load: %d/%d\n',cls,i,length(gtids));
      drawnow;
      tic;
    end

    % read annotation
    recs(i)=PASreadrecord(sprintf(VOCopts.annopath,gtids{i}));
  end
  save(cp,'gtids','recs');
end

fprintf('%s: pr: evaluating detections\n',cls);

% hash image ids
hash=wsddnVOChash_init(gtids);

% extract ground truth objects

npos=0;
gt(length(gtids))=struct('BB',[],'diff',[],'det',[]);
for i=1:length(gtids)
  % extract objects of class
  clsinds=strmatch(cls,{recs(i).objects(:).class},'exact');
  gt(i).BB=cat(1,recs(i).objects(clsinds).bbox)';
  gt(i).diff=[recs(i).objects(clsinds).difficult];
  gt(i).det=false(length(clsinds),1);
  npos=npos+sum(~gt(i).diff);
end

% load results
ids        = res.ids;
confidence = res.confidence;
BB         = res.bbox';

% sort detections by decreasing confidence
[sc,si]=sort(-confidence);
ids=ids(si);
BB=BB(:,si);

% assign detections to ground truth objects
nd=length(confidence);
tp=zeros(nd,1);
fp=zeros(nd,1);
tic;
for d=1:nd
  % display progress
  if toc>1
    fprintf('%s: pr: compute: %d/%d\n',cls,d,nd);
    drawnow;
    tic;
  end

  % find ground truth image
  i=wsddnVOChash_lookup(hash,ids{d});
  if isempty(i)
    error('unrecognized image "%s"',ids{d});
  elseif length(i)>1
    error('multiple image "%s"',ids{d});
  end

  % assign detection to ground truth object if any
  bb=BB(:,d);
  ovmax=-inf;
  for j=1:size(gt(i).BB,2)
    bbgt=gt(i).BB(:,j);
    bi=[max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
    iw=bi(3)-bi(1)+1;
    ih=bi(4)-bi(2)+1;
    if iw>0 & ih>0
      % compute overlap as area of intersection / area of union
      ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
        (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
        iw*ih;
      ov=iw*ih/ua;
      if ov>ovmax
        ovmax=ov;
        jmax=j;
      end
    end
  end
  % assign detection as true positive/don't care/false positive
  if ovmax>=VOCopts.minoverlap
    if ~gt(i).diff(jmax)
      if ~gt(i).det(jmax)
        tp(d)=1;            % true positive
        gt(i).det(jmax)=true;
      else
        fp(d)=1;            % false positive (multiple detection)
      end
    end
  else
    fp(d)=1;                    % false positive
  end
end

% compute precision/recall
fp=cumsum(fp);
tp=cumsum(tp);
rec=tp/npos;
prec=tp./(fp+tp);

ap=wsddnVOCap(rec,prec);

if draw
  % plot precision/recall
  plot(rec,prec,'-');
  grid;
  xlabel 'recall'
  ylabel 'precision'
  title(sprintf('class: %s, subset: %s, AP = %.3f',cls,VOCopts.testset,ap));
end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function out = param_map(GA)

% rundesignsimISI2
% Tor Wager, 11/24/01
% if your input var name is already a variable in the wkspace, uses that...
% otherwise creates one with zeros.

warning off
index = 1;

% -------------------------------------------------
% load GA field names
% -------------------------------------------------
fN = fieldnames(GA);
for i = 1:size(fN,1)
    eval([fN{i} ' = GA.' fN{i}])
end

% -------------------------------------------------
% common user variables
% -------------------------------------------------
tickRes = [.1 .1];
ISIrange = [.1 16]
FCrange = [0 0];
GA.sizegenerations = 100;

myISI = ISIrange(1):tickRes(1):ISIrange(2);
myFC = FCrange(1):tickRes(2):FCrange(2);

ISIsize = length(myISI);
FCsize = length(myFC);

ISItick = 1:max(round(ISIsize./25),1):ISIsize;
ISIticklab = myISI(ISItick);
FCtick = 1:2:FCsize;
FCticklab = myFC(FCtick);

disp(['Starting ' num2str(prod([ISIsize FCsize])) ' surface points.'])

out.ga = GA;
out.myISI = myISI;
out.myFC = myFC;

% ----------------------------------------------------------------
% * HRF GAMMA FUNCTION - spm 99
% ----------------------------------------------------------------

   HRF = spm_hrf(.1);						% SPM HRF sampled at .1 s
   HRF = HRF/ max(HRF);

% ----------------------------------------------------------------
% * make the map
% ----------------------------------------------------------------

for RR = 1:size(myISI,2)

    % set ISI
    GA.ISI = myISI(RR);
    numsamps = ceil(numStim*GA.ISI/TR);

    out.numsamps(find(myISI == GA.ISI)) = numsamps;

    % ----------------------------------------------------------------
    % * get smoothing matrix and autocorrelation matrix
    % ----------------------------------------------------------------
    if ~isfield(GA,'LPsmooth'),GA.LPsmooth = 1;,end
    [S,Vi,svi] = getSmoothing(HPlength,GA.LPsmooth,TR,numsamps,xc);
    clear Vi

	for CC = 1:size(myFC,2)
		disp(['Starting model ' num2str(index)])

		% set rest frequency
		GA.freqConditions(end+1) = myFC(CC);
		GA.freqConditions(1:end-1) = (1 - myFC(CC)) / (size(GA.freqConditions,2)-1);

		disp(['FREQC and ISI are:	' num2str(GA.freqConditions) ' and ' num2str(GA.ISI)])

        % ----------------------------------------------------------------
        % * generate an unsorted list of conditions
        % ----------------------------------------------------------------
        numStim = ceil(scanLength / (GA.ISI));
        if ~isempty(restevery) & ~isempty(restlength),
            numRestStim = (ceil(numStim/(mean(restevery)+restlength)) - 1) * restlength;
        else
	        numRestStim = 0;
        end
        numStimEachCond = ceil((numStim-numRestStim) * freqConditions);			            % row vector of stim in each cond
        numBlanks = numStim - numRestStim - sum(numStimEachCond);

        unsortedList = [];
        for i = 1:size(conditions,2)
            unsortedList = [unsortedList; conditions(i) * ones(numStimEachCond(i),1)];
        end
        unsortedList = [unsortedList; zeros(numBlanks,1)];

        % ----------------------------------------------------------------
        % * randomize order and test each list
        % ----------------------------------------------------------------

	    clear listMatrix
        for z = 1:GA.sizeGenerations

  	        stimList = getRandom(unsortedList);

            % -------------------------------------------------------------------------------------------------
		    % * counterbalancing
		    % -------------------------------------------------------------------------------------------------

            cbal(z) = getCounterBal(stimList, maxOrder,conditions,freqConditions);

            % -------------------------------------------------------------------------------------------------
		    % * efficiency
		    % -------------------------------------------------------------------------------------------------

			model = designvector2model(stimList,ISI,HRF,TR,numsamps,nonlinthreshold,S);
        	xtxitx = pinv(model);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
			eff(z) = calcEfficiency(contrastweights,contrasts,xtxitx,svi);

            % -------------------------------------------------------------------------------------------------
		    % * HRF shape estimation efficiency
		    % -------------------------------------------------------------------------------------------------

            delta = [];
            for i = 1:max(stimList(:,1))
                delta(:,i) = (stimList == i);
            end

			[model] = tor_make_deconv_mtx2(delta,5,TR / ISI);
            if ~isempty(S), model = S * model;,end

        	xtxitx = pinv(model);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
			hrf_eff(z) = calcEfficiency([],[],xtxitx,svi);

            % -------------------------------------------------------------------------------------------------
		    % calculate average time between stimulus repetitions in each trial type
		    % -------------------------------------------------------------------------------------------------

            for i = 1:max(stimList),mytimeBtwn(i) = mean(diff(find(stimList==i))) .* GA.ISI;, end
            timeBtwn(z) = mean(mytimeBtwn);

            % old method:
            %[fitness,convar,xtxi,models] = testlist(listMatrix,GA,HRF,'svi',svi,'S',S);
            %convar = 1./convar;
	        %eval([name '(RR,CC) = mean(convar,2);']);

        end % testing at this parameter combination

        out.cbal(RR,CC) = mean(cbal);
        out.eff(RR,CC) = mean(eff);
        out.hrf_eff(RR,CC) = mean(hrf_eff);

        out.cbal_se(RR,CC) = std(cbal) / sqrt(sizeGenerations);
        out.eff_se(RR,CC) = std(eff) / sqrt(sizeGenerations);
        out.hrf_eff_se(RR,CC) = std(hrf_eff) / sqrt(sizeGenerations);
        out.timeBtwn(RR,CC) = mean(timeBtwn);
        out.timeBtwn_se(RR,CC) = std(timeBtwn) / sqrt(sizeGenerations);

		%eval(['save ' name num2str(index) ' ' name ' SIM'])
		index = index + 1;
	end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function error_string = yalmiperror(errorcode,solver)
%YALMIPERROR Creates YALMIP error message based on error code
%
%   s = YALMIPERROR(ID) gives a textual description of an error
%   code generated by YALMIP (typically in SOLVESDP)
%
%   The complete set of error codes are
%
%   -12 Solver license cannot be located
%   -11 Solver license expired
%   -10 NaN in model data
%    -9 Specified solver name not recognized
%    -8 Problem does not satisfy geometric programming rules
%    -7 Solver does not return error codes
%    -6 Search space not bounded (bound all variables)
%    -5 License problems in solver
%    -4 Solver not applicable
%    -3 Solver not found in MATLAB path
%    -2 Successfully solved
%    -1 Unknown error
%     0 Successfully solved
%     1 Infeasible problem
%     2 Unbounded objective function
%     3 Maximum #iterations or time-limit exceeded
%     4 Numerical problems
%     5 Lack of progress
%     6 Initial solution infeasible
%     7 YALMIP sent incorrect input to solver
%     8 Feasibility cannot be determined
%     9 Unknown problem in solver
%    10 bigM failed (obsolete)
%    11 Other identified error (use savesolveroutput and refer to solver
%    documentation)
%    12 Infeasible or unbounded
%    13 YALMIP cannot determine status in solver
%    14 Model creation failed
%    15 Problem either infeasible or unbounded
%    16 User terminated
%    17 Presolve recovery failed
%    18 Missing non-negativity bounds in GP formulation
%    19 Convexity requirements not met
%    20 Solver complains about bad data
%    21 Failed to initialize bisection space (typically means infeasibility)
%    22 Ill-posed problem according to solver
%
%   See also OPTIMIZE

if nargin ==0
    help yalmiperror
    return
end

if nargin==1
    solver = '';
else
    solver = ['(' strrep(solver,'+','') ')'];
end

switch errorcode
case -12
  error_string = ['Solver license cannot be located ' solver];
case -11
  error_string = ['Solver license expired ' solver];
case -10
  error_string = ['NaN in model data ' solver];
case -9
  error_string = ['Specified solver name not recognized ' solver];
case -8
  error_string = ['Problem does not satisfy geometric programming rules'];
case -7
  error_string = ['Solver does not return error codes ' solver];
case -6
  error_string = ['Search space not bounded ' solver];
case -5
  error_string = ['License problems in solver ' solver];
 case -4
  error_string = ['Solver not applicable ' solver];
case -3
  error_string = ['Solver not found ' solver];
 case -2
  error_string = 'No suitable solver';
 case -1
  error_string = ['Unknown error ' solver];
 case 0
  error_string = ['Successfully solved ' solver ];
 case 1
  error_string = ['Infeasible problem ' solver ];
 case 2
  error_string = ['Unbounded objective function ' solver ];
 case 3
  error_string = ['Maximum iterations or time limit exceeded ' solver ];
 case 4
  error_string = ['Numerical problems ' solver ];
 case 5
  error_string = ['Lack of progress ' solver ];
 case 6
  error_string = ['Initial solution infeasible ' solver ];
 case 7
  error_string = ['YALMIP called solver with incorrect input ' solver ];
 case 8
  error_string = ['Feasibility cannot be determined ' solver ];
 case 9
  error_string = ['Unknown problem in solver (Turn on ''debug'' in sdpsettings) ' solver ];
 case 10
  error_string = ['bigM failed, increase sp.Mfactor ' solver ];
 case 11
  error_string = ['Other identified error ' solver ];
 case 12
  error_string = ['Either infeasible or unbounded ' solver ];
 case 13
  error_string = ['YALMIP cannot determine status in solver ' solver ];
 case 14
  error_string = ['Model creation failed ' solver ];
 case 15
  error_string = ['Infeasible or unbounded problem ' solver ];
 case 16
  error_string = ['User terminated ' solver ];
 case 17
  error_string = ['Presolve recovery failed ' solver ];
 case 18
  error_string = ['Missing non-negativity bounds in GP formulation ' solver ];
 case 19
  error_string = ['Convexity requirements not met ' solver ];
case 20
  error_string = ['Solver complains about bad data ' solver ];
case 21
  error_string = ['Failed to initialize bisection space (probably infeasible)' solver ];
case 22
  error_string = ['Ill-posed problem according to solver ' solver ];

 otherwise
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function tapas_sem_trapp_int_check_input(u, theta, ptheta)
% Check whether the input is compatible with mpdcm. If not an error is
% returned.
%
% Input:
% u         -- Structure. Inputs to DCM in mpdcm format.
% theta     -- Structure. Model parameters in mpdcm format.
% ptheta    -- Structure. Priors of the model in mpdcm format.
%
% Output:
%

% aponteeduardo@gmail.com
%
% Author: Eduardo Aponte, TNU, UZH & ETHZ - 2015
% Copyright 2016 by Eduardo Aponte <aponteeduardo@gmail.com>
%

tapas_mpdcm_check_input_u(u);
tapas_mpdcm_check_input_ptheta(ptheta);

check_theta(theta)

su = size(u);
st = size(theta);

assert(su(1) == st(1), ...
    'mpdcm:fmri:int:input:dmatch', ...
    'Dimensions of u and theta doesn''t match')

su = size(u{1});
du = theta{1}.dim_u;
assert(su(1) == du, ...
    'mpdcm:fmri:int:input:theta:dim_u:match_u', ...
    'theta.dim_u doesn''t match u.');
end

function check_theta(theta)
%% Throws an error if something is wrong with theta

assert(iscell(theta), ...
    'mpdcm:fmri:int:input:theta:not_cell', ...
    'theta should be a cell array')
assert(ndims(theta) == 2, ...
    'mpdcm:fmri:int:input:theta:ndim', ...
    'theta should be two dimensional, number of dimensions is %d', ...
    ndims(theta))

ofD = 0;
nfD = 0;

for i = 1:numel(theta)

    assert(isstruct(theta{i}), ...
        'mpdcm:fmri:int:input:theta:cell:not_struct', ...
        'theta{%d} is not struct', i)

    assert(isscalar(theta{i}), ...
        'mpdcm:fmri:int:input:theta:cell:ndim', ...
        'theta{%d} should be 1 x 1', i)

    tapas_mpdcm_check_input_matrix(theta, [1, 1], 'dim_x', i);
    tapas_mpdcm_check_input_matrix(theta, [1, 1], 'dim_u', i);

    if i == 1
        dx = theta{i}.dim_x;
        du = theta{i}.dim_u;
    end

    assert(dx == theta{i}.dim_x, ...
        'mpdcm:fmri:int:input:theta:cell:dim_x:not_match', ...
        'theta{%d}.dim_x doesn''t match previous values', i);

    assert(du == theta{i}.dim_u, ...
        'mpdcm:fmri:int:input:theta:cell:dim_u:not_match', ...
        'theta{%d}.dim_u doesn''t match previous values', i);

    tapas_mpdcm_check_input_matrix(theta, [dx, dx], 'A', i);
    tapas_mpdcm_check_input_matrix(theta, [dx, dx], 'B', i);
    tapas_mpdcm_check_input_matrix(theta, [dx, du], 'C', i);

    tapas_mpdcm_check_input_matrix(theta, [dx, 1], 'x0', i);
    tapas_mpdcm_check_input_matrix(theta, [dx, 1], 'beta', i);
    tapas_mpdcm_check_input_matrix(theta, [dx, 1], 'theta', i);
    tapas_mpdcm_check_input_matrix(theta, [1, 1], 'tau', i);

end

end

author: translationalneuromodeling

repo: tapas


text:
function plot_createMenu(hFig, firstCall, idxPC)

% Menu for all point clouds ----------------------------------------------------

if firstCall

    hFig.MenuBar = 'none';

    % Create menu for all PCs
    hFig.UserData.hMenuAll = uimenu(hFig, 'Label', 'All PCs');
    hFig.UserData.hMenuAll.UserData.idxPC = idxPC; % save index of PC in UserData of menu
    fillmenu(hFig.UserData.hMenuAll, hFig);

    % Menu for adding a new PC
    hFig.UserData.hMenuAddPC = uimenu(hFig, 'Label', 'Add PC...', 'Callback', @plot_addPC);

else

    % Add only idxPC to menu for all PCs
    hFig.UserData.hMenuAll.UserData.idxPC = [hFig.UserData.hMenuAll.UserData.idxPC idxPC]; % save index of PC in UserData of menu

end

% Menu for actual point cloud --------------------------------------------------

hFig.UserData.PC{idxPC}.U.hMenu = uimenu(hFig, 'Label', hFig.UserData.PC{idxPC}.label);
hFig.UserData.PC{idxPC}.U.hMenu.UserData.idxPC = idxPC; % save index of PC in UserData of menu
hFig.UserData.PC{idxPC}.U.hMenu.Position = hFig.UserData.PC{idxPC}.U.hMenu.Position - 1; % move menu
pause(0.001) % without this line position of menu doesn't change (bug???)

fillmenu(hFig.UserData.PC{idxPC}.U.hMenu, hFig);

end

function fillmenu(parentMenu, hFig)

% Menu for all PCs?
if strcmp(parentMenu.Label, 'All PCs'), menu4allPCs = true; else menu4allPCs = false; end

% Visible ----------------------------------------------------------------------

if ~menu4allPCs, h = uimenu(parentMenu, 'Label', 'Visible', 'Checked', 'On', 'Callback', @plot_setVisible); end

% Color ------------------------------------------------------------------------

h = uimenu(parentMenu, 'Label', 'Color');
if ~menu4allPCs, h.Separator = 'On'; end

% Unicolor
hUnicolor = uimenu(h, 'Label', 'Unicolor');
for i = 1:size(hFig.UserData.colors,1)
    uimenu('Parent', hUnicolor, 'Label', hFig.UserData.colors{i,2}, 'Callback', @plot_setColor, 'ForegroundColor', hFig.UserData.colors{i,1});
end
uimenu('Parent', hUnicolor, 'Label', 'by PC' , 'Callback', @plot_setColor, 'Separator', 'On');
uimenu('Parent', hUnicolor, 'Label', 'custom', 'Callback', @plot_setColor);

% Attributes
hAttributes = uimenu(h, 'Label', 'Attribute');
uimenu('Parent', hAttributes, 'Label', 'x'  , 'Callback', @plot_setColor);
uimenu('Parent', hAttributes, 'Label', 'y'  , 'Callback', @plot_setColor);
uimenu('Parent', hAttributes, 'Label', 'z'  , 'Callback', @plot_setColor);
uimenu('Parent', hAttributes, 'Label', 'RGB', 'Callback', @plot_setColor);

if ~menu4allPCs

    % Add attribute to menu for actual PC
    if isstruct(hFig.UserData.PC{parentMenu.UserData.idxPC}.A)
        att = fields(hFig.UserData.PC{parentMenu.UserData.idxPC}.A);
        if ~isempty(att)
            for i = 1:numel(att)
                h = uimenu('Parent', hAttributes, 'Label', att{i}, 'Callback', @plot_setColor);
                if i == 1, h.Separator = 'On'; end
            end
        end
    end

    % Update attributes in menu for all PCs
    updateAttributesInMenu4allPCs(hFig)

else % if menu4allPCs is true

    hFig.UserData.hMenuAllAttributes = hAttributes; % save handle to attributes menu (for adding attributes later)

end

% MarkerSize -------------------------------------------------------------------

h = uimenu(parentMenu, 'Label', 'MarkerSize');
for i = 1:10
    uimenu(h, 'Label', num2str(i), 'Callback', @plot_setMarkerSize);
end
uimenu(h, 'Label', 'custom', 'Callback', @plot_setMarkerSize, 'Separator', 'On');

% MaxPoi -----------------------------------------------------------------------

h = uimenu(parentMenu, 'Label', 'MaxPoints', 'Callback', @plot_setMaxPoints);

% ColormapName -----------------------------------------------------------------

colormaps = {'parula' 'jet' 'jetinv' 'hsv' 'hot' 'cool' 'spring' 'summer' 'autumn' 'winter' 'gray' 'bone' 'copper' 'pink' 'lines' 'colorcube' 'prism' 'flag' 'difpal' 'classpal'}; % from 'doc colormap'
h = uimenu(parentMenu, 'Label', 'ColormapName');
for i = 1:numel(colormaps)
    uimenu('Parent', h, 'Label', colormaps{i}, 'Callback', @plot_setColormapName);
end
uimenu('Parent', h, 'Label', 'custom', 'Callback', @plot_setColormapName, 'Separator', 'On');

% CAxisLim ---------------------------------------------------------------------

h = uimenu(parentMenu, 'Label', 'CAxisLim', 'Callback', @plot_setCAxisLim);

% Select -----------------------------------------------------------------------

h = uimenu(parentMenu, 'Label', 'Select', 'Separator', 'on');

% Limits
hLimits = uimenu(h, 'Label', 'Limits');

           uimenu('Parent', hLimits, 'Label', 'Set'     , 'Callback', @plot_selectLimits);
hChoose  = uimenu('Parent', hLimits, 'Label', 'Choose'  , 'Callback', @plot_selectLimits);
           uimenu('Parent', hLimits, 'Label', 'Zoom out', 'Callback', @plot_selectLimits);
hShowAll = uimenu('Parent', hLimits, 'Label', 'Show all', 'Callback', @plot_selectLimits);
if menu4allPCs % shortcuts
    hChoose.Accelerator  = 's';
    hShowAll.Accelerator = 'a';
end

% Profile
hProfile = uimenu('Parent', h, 'Label', 'Profile', 'Callback', @plot_selectProfile);

% InPolygon
hInPolygon = uimenu('Parent', h, 'Label', 'InPolygon', 'Callback', @plot_selectInPolygon);

% Modify -----------------------------------------------------------------------

h = uimenu(parentMenu, 'Label', 'Modify');

% transform
hProfile = uimenu('Parent', h, 'Label', 'transform', 'Callback', @plot_transform);

% ecef2mapTrafo
hProfile = uimenu('Parent', h, 'Label', 'ecef2mapTrafo', 'Callback', @plot_ecef2mapTrafo);

% map2ecefTrafo
hProfile = uimenu('Parent', h, 'Label', 'map2ecefTrafo', 'Callback', @plot_map2ecefTrafo);

% Plot normals -----------------------------------------------------------------

if ~menu4allPCs

    if isstruct(hFig.UserData.PC{parentMenu.UserData.idxPC}.A)
        att = fields(hFig.UserData.PC{parentMenu.UserData.idxPC}.A);
        if ~isempty(att)
            if any(strcmpi(att, 'nx'))
                uimenu(parentMenu, 'Label', 'Plot normals', 'Separator', 'on', 'Callback', @plot_plotNormals);
            end
        end
    end

end

end

function updateAttributesInMenu4allPCs(hFig)
% Update attributes in menu for all PCs
% (not in plot_updateMenus, as plot_updateMenus is called on each plot and this is only necessary when a new PC is added)

for i = 1:numel(hFig.UserData.hMenuAllAttributes.Children) % delete all attributes (except x, y, z, RGB)
    if i == 1, idx2del = []; end
    if all(~strcmp(hFig.UserData.hMenuAllAttributes.Children(i).Label, {'x' 'y' 'z' 'RGB'}))
        idx2del = [idx2del; i];
    end
end
delete(hFig.UserData.hMenuAllAttributes.Children(idx2del));

for idxPC = 1:numel(hFig.UserData.PC) % find attributes common to all PCs
    if isstruct(hFig.UserData.PC{idxPC}.A)
        att{idxPC} = fields(hFig.UserData.PC{idxPC}.A);
    else
        att{idxPC} = {};
    end
    if idxPC == 1
        commonAtt = att{idxPC};
    else
        commonAtt = intersect(commonAtt, att{idxPC});
    end
end

for i = 1:numel(commonAtt) % add attributes common to all PCs
    if i == 1, separator = 'On'; else separator = 'Off'; end
    uimenu(hFig.UserData.hMenuAllAttributes, 'Label', commonAtt{i}, 'Callback', @plot_setColor, 'Separator', separator);
end

end

author: pglira


text:
function pass = test_nonlinSys2_C1(pref)
% Test 2x2 system (sin/cos). This is the same problem as test_nonlinSystem1,
% but using CHEBMATRIX {} syntax.
%
% Asgeir Birkisson, April 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

% Smooth domain:
d = [-pi pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% Colloc1
pref.discretization = @chebcolloc1;

A = chebop(@(x,u) [u{1} - diff(u{2}, 2) + u{1}.^2; diff(u{1}) + sin(u{2})], d);
A.lbc = @(u) u{1} - 1;
A.rbc = @(u) [u{2} - 1/2; diff(u{2})];

u = mldivide(A, f, pref);

% Check the residual
err1 = norm(A(u) - f );

% Want to check BCs as well.
bcFunLeft = A.lbc(u);
bcFunRight = chebfun(A.rbc(u));
err2 = [norm(bcFunLeft(d(1))), norm(bcFunRight(d(end)))];

pass(1) = err1 < tol;
pass(2) = all( err2 < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function vt_video_main(visualDataFile, varargin)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

    %% initial some variables
    % define a variable of visual template
    global VT;

    % define the number of visual templates
    global NUM_VT;
    NUM_VT = 1;

    % define a variable of previous vt id
    global PREV_VT_ID;
    PREV_VT_ID = 1;

    % define history of visual template
    global VT_HISTORY;
    global VT_HISTORY_FIRST;
    global VT_HISTORY_OLD;

    VT_HISTORY = [0];

    % for showing the comparing result between current imgs with templates
    global MIN_DIFF_CURR_IMG_VTS;
    MIN_DIFF_CURR_IMG_VTS = [0];

    global DIFFS_ALL_IMGS_VTS;
    DIFFS_ALL_IMGS_VTS = [0];

    % define a variable for showing sub vt imgage
    global SUB_VT_IMG;

    %% Process the parameters
    % parser some pameters

    % define the threshold for vt matching
    % This threshold determines whether a new view template is generated
    global VT_MATCH_THRESHOLD;

    % define the x, y range of visual template img
    global VT_IMG_CROP_Y_RANGE; % Row range
    global VT_IMG_CROP_X_RANGE; % Column range

    % Define resized range of visual template img for reducing resolution
    global VT_IMG_RESIZE_X_RANGE;
    global VT_IMG_RESIZE_Y_RANGE;

    % define x, y shift range for template matching
    % This determines how many +- pixels (therefore rotation) will be tested for match
    global VT_IMG_X_SHIFT;
    global VT_IMG_Y_SHIFT;

    % define decay of vt for updating vts
    % VT_GLOBAL_DECAY is subtracted for all the view templates at each time step
    % VT_ACTIVE_DECAY is added the best matching view template
    global VT_GLOBAL_DECAY;
    global VT_ACTIVE_DECAY;

     % define the patch size in x, y for patch normalisation
    global PATCH_SIZE_Y_K;
    global PATCH_SIZE_X_K;
    global BLOCK_READ;

    for i=1:(nargin-1)
        if ischar(varargin{i})
            switch varargin{i}
                case 'VT_MATCH_THRESHOLD', VT_MATCH_THRESHOLD = varargin{i+1};

                case 'VT_IMG_CROP_Y_RANGE', VT_IMG_CROP_Y_RANGE = varargin{i+1};
                case 'VT_IMG_CROP_X_RANGE', VT_IMG_CROP_X_RANGE = varargin{i+1};

                case 'VT_IMG_RESIZE_X_RANGE', VT_IMG_RESIZE_X_RANGE = varargin{i+1};
                case 'VT_IMG_RESIZE_Y_RANGE', VT_IMG_RESIZE_Y_RANGE = varargin{i+1};

                case 'VT_IMG_X_SHIFT', VT_IMG_X_SHIFT = varargin{i+1};
                case 'VT_IMG_Y_SHIFT', VT_IMG_Y_SHIFT = varargin{i+1};

                case 'VT_GLOBAL_DECAY', VT_GLOBAL_DECAY = varargin{i+1};
                case 'VT_ACTIVE_DECAY', VT_ACTIVE_DECAY = varargin{i+1};

                case 'PATCH_SIZE_Y_K', PATCH_SIZE_Y_K = varargin{i+1};
                case 'PATCH_SIZE_X_K', PATCH_SIZE_X_K = varargin{i+1};

                case 'BLOCK_READ', BLOCK_READ = varargin{i+1};
            end
        end
    end

    % define half offset for reversed matching
    global VT_IMG_HALF_OFFSET;
    VT_IMG_HALF_OFFSET = [0 floor(size(VT_IMG_CROP_X_RANGE) / 2)];

    gc_x = 1;
    gc_y = 1;
    gc_z = 1;

    hdc_yaw = 0;
    hdc_pitch = 0;

    VT(1).id = 1;
    VT(1).template(1) = 1;
    VT(1).decay = 0.7;
    VT(1).gc_x = gc_x;
    VT(1).gc_y = gc_y;
    VT(1).gc_z = gc_z;
    VT(1).hdc_yaw = hdc_yaw;
    VT(1).hdc_pitch = hdc_pitch;
    VT(1).first = 1;           % don't want to inject energy as the vt is been created
    VT(1).numExp = 1;
    VT(1).exps(1).id = 1;

    VT(1).template = zeros(size(VT_IMG_CROP_Y_RANGE, 2), size(VT_IMG_CROP_X_RANGE, 2));

    %% read imgs as VT input

    % specify the movie and the frames to read
    [videoInfo, iStartFrame, totalFramesNum] = get_video_data_info(visualDataFile);

    for curFrame = 1 : totalFramesNum

        % read the visual data in blocks
        if (mod(curFrame, BLOCK_READ) == 1)
            [curVideoBlock] = read_video_block(videoInfo,curFrame, iStartFrame, totalFramesNum, BLOCK_READ);
        end

        % visual templates and visual odometry uses intensity so convert to grayscale
        [curImg] = get_current_frame_image(curVideoBlock, curFrame, BLOCK_READ);

        curGrayImg = rgb2gray(curImg);
        rawImg = im2double(curGrayImg);

        % get the most active view template
        [vt_id] = visual_template(rawImg, gc_x, gc_y,gc_z, hdc_yaw, hdc_pitch);

        % render the raw img
        subplot(3, 3, 1, 'replace');
        imshow(curGrayImg, []);
        title('Raw img');
        axis on

        % render the history of visual templates
        subplot(3, 3, [2 3 5 6], 'replace');
        plot(VT_HISTORY, '.r');
        title('Frame vs View Template');
        axis on

        % render the resized image of visual templates
        subplot(3, 3, 4, 'replace');
        imshow(SUB_VT_IMG, []);
        title('Resized img');
        axis on

        % render the history of visual templates
        subplot(3, 3, 7, 'replace');
        plot(MIN_DIFF_CURR_IMG_VTS, '.');
        title('Min SAD with each template');
        axis on

        % render the history of visual templates
        subplot(3, 3, [8 9], 'replace');
        plot(DIFFS_ALL_IMGS_VTS, '.');
        title('Min SAD of each current img');
        axis on

        DIFFS_ALL_IMGS_VTS = vt_id;

        drawnow;
    end

end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
% Script for testing the string array feature of SOFA

%% Test Strings as application-specific variable
% Load some arbritrary HRTFs
hrtf = SOFAload(fullfile(SOFAdbPath, 'database','ari','dtf_nh2.sofa'));
% Add a string array
str={};
for ii=1:hrtf.API.M
  str{ii,1}=['String' num2str(round(rand(1,1)*10000))];
end
% SOFAaddVariable(Obj,Name,Dim,Value)
hrtf2 = SOFAaddVariable(hrtf,'Test','MS',str);
% Save as SOFA
SOFAsave('stringtest_applicationvar.sofa',hrtf2);
% Reload the file
hrtf = SOFAload('stringtest_applicationvar.sofa');
% compare the strings
if prod(strcmp(hrtf.Test,hrtf2.Test))
    disp('SimpleFreeFieldHRIR: String Load-Reload: OK');
    delete('stringtest_applicationvar.sofa');
else
    error('String comparison showed differences');
end
clear all

%% Test with conventions GeneralString
% Create an empty object
Obj = SOFAgetConventions('GeneralString');
% Create numeric data with M=15, R=2, N=10
Obj.Data.Double=rand(15,2,10);
% Create string arrays
str2={}; str={};
for ii=1:15
  id = num2str(round(rand(1,1)*1000000));
  str{ii,1}=['X' id];
  str2{ii,1}=['Left' id];
  str2{ii,2}=['Right' id];
end
Obj.String2 = str2;      % String1=[MRS]
Obj.Data.String1 = str;  % Data.String1=[MS]
Obj.Data.String2 = str2; % Data.String2=[MRS]
% Add a new string with dimensions [RS]
strn={'left ear'; 'right ear'};
Obj = SOFAaddVariable(Obj, 'Ears', 'RS', strn);
% Update dimensions
Obj = SOFAupdateDimensions(Obj);
% Save as SOFA
SOFAsave('stringtest_generalstring.sofa',Obj);
% Reload the file
Obj2 = SOFAload('stringtest_generalstring.sofa');
% Compare the strings
if ~prod(strcmp(Obj2.Data.String2,Obj.Data.String2))
    error('Data.String2: Comparison showed differences');
end
if ~prod(strcmp(Obj2.String2,Obj.String2))
    error('String2: Comparison showed differences');
end
if ~prod(strcmp(Obj2.Data.String1,Obj.Data.String1))
    error('Data.String1: Comparison showed differences');
end
if ~prod(strcmp(Obj2.Ears,Obj.Ears))
    error('Ears: Comparison showed differences');
end
disp('GeneralString: String1, String2, Data, Ears: Load-Reload: OK');
clear all
delete('stringtest_generalstring.sofa');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% LFBuild4DFreqHyperfan - construct a 4D hyperfan passband filter in the frequency domain
%
% Usage:
%
%     [H, FiltOptions] = LFBuild4DFreqHyperfan( LFSize, Slope1, Slope2, BW, FiltOptions )
%     H = LFBuild4DFreqHyperfan( LFSize, Slope1, Slope2, BW )
%
% This file constructs a real-valued magnitude response in 4D, for which the passband is a hyperfan.
% This is useful for selecting objects over a range of depths from a lightfield, i.e. volumetric
% focus.
%
% Once constructed the filter must be applied to a light field, e.g. using LFFilt4DFFT. The
% LFDemoBasicFilt* files demonstrate how to contruct and apply frequency-domain filters.
%
% A more technical discussion, including the use of filters for denoising and volumetric focus, and
% the inclusion of aliases components, is included in:
%
% [2] D.G. Dansereau, O. Pizarro, and S. B. Williams, "Linear Volumetric Focus for Light Field
% Cameras," to appear in ACM Transactions on Graphics (TOG), vol. 34, no. 2, 2015.
%
% Inputs:
%
%     LFSize : Size of the frequency-domain filter. This should match or exceed the size of the
%     light field to be filtered. If it's larger than the input light field, the input is
%     zero-padded to match the filter's size by LFFilt4DFFT.
%
%     Slope : The slope of the planar passband. If different slopes are desired in s,t and u,v,
%     the optional aspect parameter should be used.
%
%     BW : 3-db Bandwidth of the planar passband.
%
%     [optional] FiltOptions : struct controlling filter construction
%            HyperfanMethod : 'sweep' or 'direct', controls how the hyperfan is constructed. 'sweep'
%                             concatenates frequency planes, while 'direct' evaluates a single 4D
%                             equation; default 'direct'. Warning: sweep can be very slow.
%               HyperconeBW : Sets the underlying hypercone's BW separately from the dual-fan,
%                             default is to use the BW parameter for both; only applies when
%                             HyperfanMethod is 'direct'
%               SlopeMethod : 'Skew' or 'Rotate' default 'skew'
%                 Precision : 'single' or 'double', default 'single'
%                   Rolloff : 'Gaussian' or 'Butter', default 'Gaussian'
%                     Order : controls the order of the filter when Rolloff is 'Butter', default 3
%                  Aspect4D : aspect ratio of the light field, default [1 1 1 1]
%                    Window : Default false. By default the edges of the passband are sharp; this adds
%                             a smooth rolloff at the edges when used in conjunction with Extent4D or
%                             IncludeAliased.
%                  Extent4D : controls where the edge of the passband occurs, the default [1 1 1 1]
%                             is the edge of the Nyquist box. When less than 1, enabling windowing
%                             introduces a rolloff after the edge of the passband. Can be greater
%                             than 1 when using IncludeAliased.
%            IncludeAliased : default false; allows the passband to wrap around off the edge of the
%                             Nyquist box; used in conjunction with Window and/or Extent4D. This can
%                             increase processing time dramatically, e.g. Extent4D = [2,2,2,2]
%                             requires a 2^4 = 16-fold increase in time to construct the filter.
%                             Useful when passband content is aliased, see [2].
%                SweepSteps : For HyperfanMethod 'sweep', how many steps to use in the sweep.
%
% Outputs:
%
%                 H : real-valued frequency magnitude response
%       FiltOptions : The filter options including defaults, with an added PassbandInfo field
%                     detailing the function and time of construction of the filter
%
% User guide: LFToolbox.pdf
% See also:  LFDemoBasicFiltGantry, LFDemoBasicFiltIllum, LFDemoBasicFiltLytroF01,
% LFBuild2DFreqFan, LFBuild2DFreqLine, LFBuild4DFreqDualFan, LFBuild4DFreqHypercone,
% LFBuild4DFreqHyperfan, LFBuild4DFreqPlane, LFFilt2DFFT, LFFilt4DFFT, LFFiltShiftSum

% Copyright (c) 2013-2020 Donald G. Dansereau

function [H, FiltOptions] = LFBuild4DFreqHyperfan( LFSize, Slope1, Slope2, BW, FiltOptions )

FiltOptions = LFDefaultField('FiltOptions', 'HyperfanMethod', 'Direct'); % 'Direct', 'Sweep'

switch( lower(FiltOptions.HyperfanMethod ))
	case 'sweep'
		FiltOptions = LFDefaultField('FiltOptions', 'SweepSteps', 5);
		SweepVec = linspace(Slope1, Slope2, FiltOptions.SweepSteps);
		[H, FiltOptions] = LFBuild4DFreqPlane( LFSize, SweepVec(1), BW, FiltOptions );
		for( CurSlope = SweepVec(2:end) )
			H = max(H, LFBuild4DFreqPlane( LFSize, CurSlope, BW, FiltOptions ));
		end

	case 'direct'
		FiltOptions = LFDefaultField('FiltOptions', 'HyperconeBW', BW);
		[H, FiltOptions] = LFBuild4DFreqHypercone( LFSize, FiltOptions.HyperconeBW, FiltOptions );
		[Hdf, FiltOptions] = LFBuild4DFreqDualFan( LFSize, Slope1, Slope2, BW, FiltOptions );
		H = H .* Hdf;

	otherwise
		error('Unrecognized hyperfan construction method');
end

TimeStamp = datestr(now,'ddmmmyyyy_HHMMSS');
FiltOptions.PassbandInfo = struct('mfilename', mfilename, 'time', TimeStamp, 'VersionStr', LFToolboxVersion);

repo: LFToolbox


text:
% Calculate Kendall coefficient of agreement from a winning matrix
% according to Wei-Sheng La et al. "A Comparative Study for Single Image
% Blind Deblurring", CVPR 2016
function u = coefficientOfAgreement(C)

    % Number of methods described by the winning matrix.
    numMethods = sum(sum(C + C', 2) > 0);
    % Number of observers for the winning matrix.
    numUsers = round(sum(C(:)) / nchoosek(numMethods, 2) );

    C = C(:);
    C = C(C > 1);
    w = 0;
    for i = 1:length(C)
        w = w + nchoosek(C(i), 2);
    end

    u = (2 * w) / ( nchoosek(numMethods, 2) * nchoosek(numUsers, 2) ) - 1;

author: thomas-koehler

repo: SupER


text:
%MDL_UR5 Create model of Universal Robotics UR3 manipulator
%
% MDL_UR5 is a script that creates the workspace variable ur3 which
% describes the kinematic characteristics of a Universal Robotics UR3 manipulator
% using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         arm along +ve x-axis configuration
%
% Reference::
% - https://www.universal-robots.com/how-tos-and-faqs/faq/ur-faq/actual-center-of-mass-for-robot-17264/
%
% Notes::
% - SI units of metres are used.
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also mdl_ur5, mdl_ur10, mdl_puma560, SerialLink.

% MODEL: Universal Robotics, UR3, 6DOF, standard_DH

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

function r = mdl_ur3()

    deg = pi/180;

    % robot length values (metres)
    a = [0, -0.24365, -0.21325, 0, 0, 0]';

    d = [0.1519, 0, 0, 0.11235, 0.08535, 0.0819]';

    alpha = [1.570796327, 0, 0, 1.570796327, -1.570796327, 0]';

    theta = zeros(6,1);

    DH = [theta d a alpha];

    mass = [2, 3.42, 1.26, 0.8, 0.8, 0.35];

    center_of_mass = [
        0,-0.02, 0
        0.13, 0, 0.1157
        0.05, 0, 0.0238
        0, 0, 0.01
        0, 0, 0.01
        0, 0, -0.02 ];

    % and build a serial link manipulator

    % offsets from the table on page 4, "Mico" angles are the passed joint
    % angles.  "DH Algo" are the result after adding the joint angle offset.

    robot = SerialLink(DH, ...
        'name', 'UR3', 'manufacturer', 'Universal Robotics');

    % add the mass data, no inertia available
    links = robot.links;
    for i=1:6
        links(i).m = mass(i);
        links(i).r = center_of_mass(i,:);
    end

    % place the variables into the global workspace
    if nargin == 1
        r = robot;
    elseif nargin == 0
        assignin('caller', 'ur3', robot);
        assignin('caller', 'qz', [0 0 0 0 0 0]); % zero angles
        assignin('caller', 'qr', [180 0 0 0 90 0]*deg); % vertical pose as per Fig 2
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [xUpdate,PUpdate,innov,Pzz,W]=KalmanUpdate(xPred,PPred,z,R,H,condIssue)
%KALMANUPDATE Perform the measurement update step in the standard linear
%             Kalman filter. This assumes that the predicted target state
%             and the measurement are both multivariate Gaussian
%             distributed.
%
%INPUTS: xPred The xDimX1 predicted target state.
%        PPred The xDimXxDim predicted state covariance matrix.
%            z The zDimX1 measurement vector.
%            R The zDimXzDim measurement covariance matrix.
%            H The optional zDimXxDim measurement matrix. The measurement
%              is modeled as z=H*x+noise. If this parameter is omitted or
%              an empty matrix is passed, then H will be taken as a
%              zDimXzDim identity matrix followed by columns of zeros
%              (Assuming that zDim<=xDim. Otherwise, H must be provided).
%    condIssue An optional parameter indicating whether a conditioning
%              issue with the measurement covariance matrix. Singularity
%              can occur if two elements are perfectly correlated. This
%              replaces the matrix inverse with a pseudoinverse.
%
%OUTPUTS: xUpdate The xDim X 1 updated target state vector.
%         PUpdate The updated xDim X xDim state covariance matrix.
%      innov, Pzz The zDimX1 innovation and the zDimXzDim innovation
%                 covariance matrix are returned in case one wishes to
%                 analyze the consistency of the estimator or use those
%                 values in gating or likelihood evaluation.
%               W The xDimXzDim gain used in the update. This can be
%                 useful when gating and using the function
%                 calcMissedGateCov.
%
%Given a prediction of the target state xPred with an associated covariance
%matrix PPred, and assuming that they are the mean and covariance matrix of
%a Gaussian distribution, under the measurement model
%z=H*x+w
%where z is the measurement, x is the true target state, H is a zDimXxDim
%matrix and w is additive Gaussian noise, this function updates the state
%estimate and its associated covariance estimate.
%
%The Joseph-form covariance update is used for improved numerical
%stability. The algorithm is derived in Chapter 5 of [1].
%
%EXAMPLE:
%This example simulates a linear dynamic process that extracts Cartesian
%measurements. it then looks at the measurement error reduction factor
%(MERF), and the normalized estimation error squared (NEES). The MERF
%should be less than 1 or else there is no point in using the filter-
%connecting the dots works better. The NEES should be close to 1 to
%indicate covariance consistency.
% numMonteCarlo=500;
% numSteps=50;
% %Power spectral density of the process noise.
% q0=processNoiseSuggest('PolyKal-ROT',2*9.8,1);
% xDim=4;%2D position and velocity.
% T=1;%prediction interval.
% F=FPolyKal(T,xDim,1);
% Q=QPolyKal(T,xDim,1,q0);
% SQ=chol(Q,'lower');
% R=[50^2,-1000;
%    -1000,50^2];
% SR=chol(R,'lower');
% H=[eye(2,2),zeros(2,2)];
% zDim=size(H,1);
%
% z=zeros(zDim,numSteps,numMonteCarlo);
% xTrue=zeros(xDim,numSteps,numMonteCarlo);
% %The initial state.
% xTrue(:,1,:)=repmat([1e3;0;75;-50],[1,1,numMonteCarlo]);
% maxVel=400;%Used for single-point initiation.
% higherDerivStdDev=maxVel/sqrt(2);%used for single-point initiation.
%
% xEst=zeros(xDim,numSteps,numMonteCarlo);
% PEst=zeros(xDim,xDim,numSteps,numMonteCarlo);
% for curRun=1:numMonteCarlo
%     z(:,1,curRun)=H*xTrue(:,1,curRun)+SR*randn(zDim,1);
%     for curStep=2:numSteps
%         xTrue(:,curStep,curRun)=F*xTrue(:,curStep-1,curRun)+SQ*randn(xDim,1);
%         z(:,curStep,curRun)=H*xTrue(:,curStep,curRun)+SR*randn(zDim,1);
%     end
%
%     %One-point initialization. This works in this instance, because H
%     %extracts the Cartesian components.
%     [xInit,PInit]=onePointCartInit(z(:,1,curRun),R,higherDerivStdDev,1);
%     xEst(:,1,curRun)=xInit;
%     PEst(:,:,1,curRun)=PInit;
%     for curStep=2:numSteps
%         [xPred,PPred]=discKalPred(xEst(:,curStep-1,curRun,1),PEst(:,:,curStep-1,curRun,1),F,Q);
%         [xUpdate,PUpdate]=KalmanUpdate(xPred,PPred,z(:,curStep,curRun),R,H);
%         xEst(:,curStep,curRun)=xUpdate;
%         PEst(:,:,curStep,curRun)=PUpdate;
%     end
% end
%
% posTrue=xTrue(1:2,:,:);
% %The MERF requires the measurement to be in the same coordinates are the
% %corresponding element sof the state.
% MERF=calcMERF(posTrue,z,xEst(1:2,:,:),0,true);
% NEES=calcNEES(xTrue,xEst,PEst);
% figure(1)
% clf
% hold on
% plot(MERF,'-k','linewidth',2)
% h1=xlabel('step');
% h2=ylabel('Measurement Error Reduction Factor');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
% figure(2)
% clf
% hold on
% plot(NEES,'-r','linewidth',2)
% h1=xlabel('step');
% h2=ylabel('NEES');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%
%July 2012 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(xPred,1);

if(nargin<6||isempty(condIssue))
    condIssue=false;
end

if(nargin<5||isempty(H))
	zDim=size(z,1);
    H=[eye(zDim,zDim),zeros(zDim,xDim-zDim)];
end

zPred=H*xPred;
innov=z-zPred;

Pzz=R+H*PPred*H';
%Ensure symmetry
Pzz=(Pzz+Pzz')/2;

if(condIssue==false)
    opts.SYM=true;
    opts.RECT=false;
    opts.TRANSA=false;
    W=linsolve(Pzz,H*PPred,opts)';%W=PPred*H'/Pzz;
else
    %If there is possibly poor conditioning.
    W=PPred*H'*pinv(Pzz);
end

xUpdate=xPred+W*innov;

temp=W*H;
temp=eye(xDim,xDim)-temp;
PUpdate=temp*PPred*temp'+W*R*W';

%Ensure symmetry
PUpdate=(PUpdate+PUpdate')/2;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function MIMICII_data_importer
% MIMICII_data_importer imports data from the MIMIC II database for use
% with RRest, a toolbox of respiratory rate algorithms.
%
%               MIMICII_data_importer
%
%	Inputs:
%       none
%                       just specify the relevant paths in the
%                       "universal_parameters" function below.
%
%	Outputs:
%       a single file containing all the data is written to the path
%       specified by "up.paths.analysis_path".
%
%   Requirements:
%       This requires the WFDB Toolbox, which can be downloaded from:
%           https://physionet.org/physiotools/matlab/wfdb-app-matlab/
%       data download is performed automatically using the script below.
%
%   Further Information:
%       This version of the MIMICII_data_importer is provided to facilitate
%       reproduction of the analysis performed in:
%           Charlton P.H. et al., "Waveform Analysis to Estimate
%           Respiratory Rate" [In Press]
%       Further information on this study can be obtained at:
%           http://peterhcharlton.github.io/RRest/waveform_analysis.html
%       In addition, further information on RRest, including future
%       versions, can be obtained at:
%           http://peterhcharlton.github.io/RRest/index.html
%
%   Comments, Questions, Criticisms, Feedback, Contributions:
%       See: http://peterhcharlton.github.io/RRest/contributions.html
%
%   Version:
%       v.1 - published on 23rd Feb 2016 by Peter Charlton
%
%   Licence:
%       please see the accompanying file named "LICENSE"
%

%% Setup
up = universal_parameters;

%% Download data
download_data(up);

%% Identify patient stays
pt_stays = identify_pt_stays(up);

%% Extract data from files on computer
data = extract_from_mimic_ii(pt_stays, up);

%% Save data in common format
save_data_in_common_format(data, up);

%% Identify sub-groups of neonates and adults
% This isn't necessary to obtain the data, but creates a plot used in the
% book chapter.
results = sub_group_identification(up);

end

function up = universal_parameters

fprintf('\n -- Setting up Universal Parameters')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% PARAMETERS TO BE SPECIFIED %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Specify the root data directory (where the data will be stored)
up.paths.data_root = 'C:\Documents\Data\mimicii\';
up.paths.data_save_folder = 'C:\Documents\Data\';

% Specify the web address of the data to be downloaded
rel_database = 'mimic2wdb';
up.paths.database_dir = ['http://physionet.org/physiobank/database/', rel_database];

% Please note that your system may require the slashes in file paths to be
% of the opposite direction. In which case, change the following:
slash = '\';

% if you want .eps illustrations, then do as follows:
up.eps_figs = 0;  % set this to 1
if up.eps_figs
    % download 'export_fig' from:
    % http://uk.mathworks.com/matlabcentral/fileexchange/23629-export-fig
    % add the path of your download:
    export_fig_dir_path = 'C:\Users\pc13\Documents\GitHub\phd\Tools\Other Scripts\export_fig\';
    addpath(export_fig_dir_path)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%   OTHER PARAMETERS   %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%   (don't need editing)   %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% extraction definitions
up.extraction.period = 600;    % period of waveform data to extract (in s)
up.extraction.rel_sigs = {'II', 'PLETH', 'RESP'};    % required signals
up.extraction.rel_nums = {'HR', 'PULSE', 'RESP'};    % required numerics
up.no_pt_stays = 100;

% database definitions
up.mimic_db.version = 3;
up.mimic_db.part = 0;

% paths
up.paths.file_location = [up.paths.data_root, 'physionet.org\physiobank\database\', rel_database, '\', num2str(up.mimic_db.version), num2str(up.mimic_db.part), '\'];
up.paths.analysis_path = [up.paths.data_root, 'Analysis_files', slash 'Data_for_Analysis', slash ];
up.paths.plots = [up.paths.data_root, 'Analysis_files', slash 'Results', slash 'Figures', slash ];
up.paths.records_filepath = [up.paths.analysis_path, 'RECORDS'];

% url paths
up.paths.url.records = [up.paths.database_dir, '/', num2str(up.mimic_db.version), num2str(up.mimic_db.part), '/', 'RECORDS'];

% Add all functions within the directory
addpath(genpath(fileparts(mfilename('fullpath'))));

% check that all the required folders exist. If not, create them:
req_folders = {up.paths.data_root, up.paths.analysis_path, up.paths.plots};
displayed_msg = 0;
for req_folder_no = 1 : length(req_folders)
    curr_folder = req_folders{req_folder_no};
    if ~exist(curr_folder, 'dir')
        if displayed_msg == 0
            fprintf('\n -- Creating folder in which to store data and analysis')
            displayed_msg = 1;
        end
        mkdir(curr_folder)
    end
end

% Check that the WFDB Matlab Toolbox has been installed
if ~exist('getWfdbClass', 'file')
    error('Couldn''t find the WFDB Matlab Toolbox. Please install as described at the top of the file.')
end

end

function download_data(up)

fprintf('\n -- Downloading Data')

%% set current directory to that of wget:
cd(up.paths.data_root)

%% Find out what records are available:
% download list of records
outfilename = websave(up.paths.records_filepath,up.paths.url.records);

% form list of records from downloaded file
records = load(up.paths.records_filepath);

%% Download each record
rel_records = records(1:up.no_pt_stays);

for rec_no = 1 : length(rel_records)

    % specify record
    record = rel_records(rec_no);

    % specify download location
    down_loc = [up.paths.database_dir, '/', num2str(up.mimic_db.version), num2str(up.mimic_db.part), '/', num2str(record), '/', 'RECORDS'];

    % set save location
    save_folder = [up.paths.file_location, num2str(record), '\'];
    save_loc = [save_folder, 'RECORDS'];
    if ~exist(save_folder, 'dir')
        mkdir(save_folder)
    end

    % download RECORDS file for this record
    outfilename = websave(save_loc,down_loc);

    % find out what individual files are in this record:
    fileID = fopen(save_loc);
    file_names = textscan(fileID, '%s %*[^\n]'); file_names = file_names{1};
    fclose(fileID);

    for file_no = 1 : length(file_names)

        % download header file
        curr_file_name = [file_names{file_no}, '.hea'];
        down_loc = [up.paths.database_dir, '/', num2str(up.mimic_db.version), num2str(up.mimic_db.part), '/', num2str(record), '/', curr_file_name];
        save_loc = [save_folder, curr_file_name];
        outfilename = websave(save_loc,down_loc);

        if ~strcmp(file_names{file_no}, num2str(record))
            % also download data file
            curr_file_name = [file_names{file_no}, '.dat'];
            down_loc = [up.paths.database_dir, '/', num2str(up.mimic_db.version), num2str(up.mimic_db.part), '/', num2str(record), '/', curr_file_name];
            save_loc = [save_folder, curr_file_name];
            outfilename = websave(save_loc,down_loc);
        end
    end

end

end

function pt_stays = identify_pt_stays(up)

fprintf('\n -- Identifying Patient Stays')

% Folders (each one is a patient stay)
dirs = dir(up.paths.file_location);

elim_els = ones(length(dirs),1);
for s = 1:length(dirs)
    if length(dirs(s).name) < 3
        elim_els(s) = 0;
    end
end

dirs = dirs(logical(elim_els'));

pt_stays = extractfield(dirs, 'name');

end

function data = extract_from_mimic_ii(pt_stays, up)
%
% Based on RDSAMP, part of the WFDB Toolbox. Further details are available
% at:
%       http://physionet.org/physiotools/matlab/wfdb-app-matlab/
%
% Reads MIMIC II data and outputs it as a structure of .t and .v values.

fprintf('\n -- Extracting Data from MIMIC II files')

if up.no_pt_stays > length(pt_stays)
    warning('There are fewer ICU stay datasets than requested')
end
no_pt_stays = min([up.no_pt_stays, length(pt_stays)]);

extraction_stats.no_pt_stays = no_pt_stays;
extraction_stats.no_records_searched = 0;
extraction_stats.no_records_without_signals = 0;
extraction_stats.no_records_without_num = 0;
extraction_stats.no_records_with_sigs_and_num = 0;

for pt_stay_no = 1 : no_pt_stays

    % set current dir
    dir_path = [up.paths.file_location, pt_stays{pt_stay_no}, '\'];
    cd(dir_path)

    % check to see if there are any files
    temp = dir([dir_path, pt_stays{pt_stay_no}, '_0*.hea']);
    if ~isempty(temp)
        records = extractfield(temp, 'name'); clear temp
    else
        extraction_stats.no_records_searched = extraction_stats.no_records_searched+1;
        extraction_stats.no_records_without_signals = extraction_stats.no_records_without_signals+1;
        continue
    end
    clear dir_path

    remaining_time_to_extract = up.extraction.period;

    for record_no = 1 : length(records)

        % setup java
        persistent javaWfdbExec config
        if(isempty(javaWfdbExec))
            [javaWfdbExec,config]=getWfdbClass('rdsamp');
        end

        %Remove file extension
        record_str = records{record_no}(1:(end-4));

        % Extract file information (including which signals are present)
        [siginfo,Fs]=extract_file_info(record_str);

        % Identify relevant signals
        sig_inds = 1:siginfo.no_sigs;
        sig_log = zeros(length(sig_inds),1);
        for sig_no = 1 : length(sig_inds)
            if sum(strcmp(up.extraction.rel_sigs, siginfo.sig_names(sig_no)))
                sig_log(sig_no) = 1;
            end
        end

        % skip if this doesn't have both the required signals
        if sum(sig_log) ~= length(up.extraction.rel_sigs)
            clear current_waves extracted_data Fs ListCapacity siginfo signalList javaWfdbExec config wfdb_argument
            continue
        end

        % create wfdb argument
        N0 = 1;
        wfdb_argument={'-r',record_str,'-Ps','-f',['s' num2str(N0-1)]};

        % Make command to only extract the required signals
        signalList = sig_inds(logical(sig_log));
        if(~isempty(signalList))
            wfdb_argument{end+1}='-s ';
            %-1 is necessary because WFDB is 0 based indexed.
            for sInd=1:length(signalList)
                wfdb_argument{end+1}=[num2str(signalList(sInd)-1)];
            end
        end

        % determine how many samples to extract from this file
        remaining_samples_to_extract = remaining_time_to_extract*Fs;
        N=siginfo.no_samps;
        if remaining_samples_to_extract < N
            N = remaining_samples_to_extract;
        end
        remaining_time_to_extract = remaining_time_to_extract - (N/Fs);

        % skip if we don't require any more data
        if N == 0
            clear current_waves extracted_data Fs ListCapacity siginfo signalList javaWfdbExec config wfdb_argument
            break
        end

        % additional bits of wfdb_argument
        wfdb_argument{end+1}='-t';
        wfdb_argument{end+1}=['s' num2str(N)];

        % extract data
        extracted_data=javaWfdbExec.execToDoubleArray(wfdb_argument);
        clear wfdb_argument

        %% Re-format data for output

        waves.fs = siginfo.fs;

        % store this new bit of the wave
        current_waves = fieldnames(waves);
        signalNames = siginfo.sig_names(signalList);
        for wave_no = 1 : length(up.extraction.rel_sigs)
            rel_sig_ind = find(strcmp(signalNames, up.extraction.rel_sigs{wave_no}));
            if ~strcmp(current_waves, up.extraction.rel_sigs{wave_no})
                eval(['waves.' up.extraction.rel_sigs{wave_no} '.t = extracted_data(:, 1);']);
                eval(['waves.' up.extraction.rel_sigs{wave_no} '.v = extracted_data(:, rel_sig_ind+1);']);
            else
                eval(['current_last_time = waves.' up.extraction.rel_sigs{wave_no} '.t(end);']);
                eval(['waves.' up.extraction.rel_sigs{wave_no} '.t = [waves.' up.extraction.rel_sigs{wave_no} '.t; current_last_time + (1/Fs) + extracted_data(:, 1)];']);
                eval(['waves.' up.extraction.rel_sigs{wave_no} '.v = [waves.' up.extraction.rel_sigs{wave_no} '.v; extracted_data(:, rel_sig_ind+1)];']);
            end
        end
        clear current_waves extracted_data Fs ListCapacity siginfo signalList javaWfdbExec config

    end

    %% Extract corresponding numerics

    if exist('waves', 'var')

        % setup java
        persistent javaWfdbExec_num config_num
        if(isempty(javaWfdbExec_num))
            [javaWfdbExec_num,config_num]=getWfdbClass('rdsamp');
        end

        filepath = [pt_stays{pt_stay_no}, 'n'];

        if exist([up.paths.file_location, pt_stays{pt_stay_no}, '\', filepath, '.hea'], 'file')

            % Extract file information (including which signals are present)
            [numinfo,Fs]=extract_file_info(filepath);
            loc = numinfo.loc;

            % Identify relevant numerics
            num_inds = 1:numinfo.no_sigs;
            num_log = zeros(length(num_inds),1);
            for num_no = 1 : length(num_inds)
                if sum(strcmp(up.extraction.rel_nums, numinfo.sig_names(num_no)))
                    num_log(num_no) = 1;
                end
            end

            % skip if this doesn't have both the required signals
            if sum(num_log) == length(up.extraction.rel_nums)

                numerics.fs = numinfo.fs;

                % create wfdb argument
                N0 = 1;
                N = up.extraction.period*numerics.fs;
                wfdb_argument={'-r',filepath,'-Ps','-f',['s' num2str(N0-1)]};

                % Make command to only extract the required signals
                numList = num_inds(logical(num_log));
                if(~isempty(numList))
                    wfdb_argument{end+1}='-s ';
                    %-1 is necessary because WFDB is 0 based indexed.
                    for sInd=1:length(numList)
                        wfdb_argument{end+1}=[num2str(numList(sInd)-1)];
                    end
                end

                % Additional wfdb_argument bits
                wfdb_argument{end+1}='-t';
                wfdb_argument{end+1}=['s' num2str(N)];

                % extract data
                extracted_data=javaWfdbExec_num.execToDoubleArray(wfdb_argument);
                clear wfdb_argument

                %% Re-format data for output

                % store this new bit of the wave
                current_numerics = fieldnames(numerics);
                numNames = numinfo.sig_names(numList);
                for wave_no = 1 : length(up.extraction.rel_nums)
                    rel_num_ind = find(strcmp(numNames, up.extraction.rel_nums{wave_no}));
                    if ~strcmp(current_numerics, up.extraction.rel_nums{wave_no})
                        eval(['numerics.' up.extraction.rel_nums{wave_no} '.t = extracted_data(:, 1);']);
                        eval(['numerics.' up.extraction.rel_nums{wave_no} '.v = extracted_data(:, rel_num_ind+1);']);
                    else
                        eval(['current_last_time = numerics.' up.extraction.rel_nums{wave_no} '.t(end);']);
                        eval(['numerics.' up.extraction.rel_nums{wave_no} '.t = [numerics.' up.extraction.rel_nums{wave_no} '.t; current_last_time + (1/Fs) + extracted_data(:, 1)];']);
                        eval(['numerics.' up.extraction.rel_nums{wave_no} '.v = [numerics.' up.extraction.rel_nums{wave_no} '.v; extracted_data(:, rel_num_ind+1)];']);
                    end
                end

            end

        end

        clear current_numerics extracted_data Fs ListCapacity numinfo numList javaWfdbExec_num config_num

    end

    %% Find extraction statistics
    extraction_stats.no_records_searched = extraction_stats.no_records_searched+1;
    if ~exist('numerics', 'var')
        extraction_stats.no_records_without_num = extraction_stats.no_records_without_num+1;
    end
    if ~exist('waves', 'var')
        extraction_stats.no_records_without_signals = extraction_stats.no_records_without_signals+1;
    end

    %% Store data for this patient stay in the output
    if exist('waves', 'var') && exist('numerics', 'var')
        extraction_stats.no_records_with_sigs_and_num = extraction_stats.no_records_with_sigs_and_num+1;
        data{extraction_stats.no_records_with_sigs_and_num}.loc = loc;
        data{extraction_stats.no_records_with_sigs_and_num}.ID = pt_stays{pt_stay_no};
        data{extraction_stats.no_records_with_sigs_and_num}.waves = waves;
        data{extraction_stats.no_records_with_sigs_and_num}.numerics = numerics;
        clear waves numerics
    end

end

%% Find out how many records had the required duration signals
total_records = length(data);
recs_w_req_duration = zeros(total_records,1);
for s = 1  : length(data)
    if sum(strcmp(fieldnames(data{1,s}), 'waves'))
        if length(data{1,s}.waves.II.t) == up.extraction.period*data{1,s}.waves.fs ...
            && length(data{1,s}.waves.PLETH.t) == up.extraction.period*data{1,s}.waves.fs ...
            && length(data{1,s}.waves.RESP.t) == up.extraction.period*data{1,s}.waves.fs
            recs_w_req_duration(s) = 1;
        end
    end
end

%% Only output those records with the required duration
extraction_stats.no_records_with_req_duration = sum(recs_w_req_duration);
data = data(logical(recs_w_req_duration));

%% Create Pie Chart of extraction stats

% bar chart
ftsize = 14;
figure('Position', [200, 200, 800, 450])
labels = {'Signals or numerics absent', 'Shorter than 10 mins', 'Included in analysis'};
piedata(3) = extraction_stats.no_records_with_req_duration;
piedata(2) = extraction_stats.no_records_with_sigs_and_num - extraction_stats.no_records_with_req_duration;
piedata(1) = extraction_stats.no_records_searched - extraction_stats.no_records_with_sigs_and_num;
p = pie(piedata);
for s = 2:2: length(p)
    p(s).FontSize = ftsize;
end
set(gca, 'FontSize', ftsize)
legend(labels, 'Location', 'southeastoutside')
title('Reasons for exclusion of records', 'FontSize', ftsize)
% save
save_name = 'record_exclusion';
savepath = [up.paths.plots, save_name, '.png'];
set(gcf,'PaperUnits','inches');
set(gcf,'PaperPosition',[1 1 8 4.5]);
print(gcf, savepath, '-dpng')
close

end

function [headers, Fs] = extract_file_info(filepath)

%% Obtain info on signals

fileID = fopen([filepath, '.hea']);
header_lines = textscan(fileID,'%s','Delimiter','\n'); header_lines = header_lines{1,1};
fclose all;

% Extract the relevant header information
temp.top_line = textscan(header_lines{1}, '%s', 'delimiter', sprintf(' ')); temp.top_line = temp.top_line{1,1};  % adapted from http://uk.mathworks.com/matlabcentral/newsreader/view_thread/249016
headers.no_sigs = str2double(temp.top_line{2});
headers.fs = str2double(temp.top_line{3});
headers.no_samps = str2double(temp.top_line{4});

% extract signal names
headers.sig_names = cell(0);
for sig_no = 1 : headers.no_sigs
    temp.curr_line = textscan(header_lines{1+sig_no}, '%s', 'delimiter', sprintf(' ')); temp.curr_line = temp.curr_line{1,1};  % adapted from http://uk.mathworks.com/matlabcentral/newsreader/view_thread/249016
    headers.sig_names = [headers.sig_names; temp.curr_line{end}];
end

% extract location
loc_line = header_lines{end};
headers.loc = loc_line(13:end);

% putting in a format for rdsamp

Fs = headers.fs;

end

function save_data_in_common_format(old_data, up)

fprintf('\n -- Saving data in appropriate format')

for subj_el = 1:length(old_data)

    SID = old_data{1,subj_el}.ID;
    loc = old_data{1,subj_el}.loc;

    % insert group name
    if strcmp(loc, 'nicu')
        data(1,subj_el).group = 'neonate';
    else
        data(1,subj_el).group = 'adult';
    end

    % Insert fixed params

    data(1,subj_el).fix.id = SID;
    data(1,subj_el).fix.loc = loc;
    data(1,subj_el).fix.ventilation = 'unknown';
    data(1,subj_el).fix.recording_conditions = 'critical care';

    % insert PPG signal
    ppg = old_data{1,subj_el}.waves.PLETH;
    data(1,subj_el).ppg.v = ppg.v;
    clear ppg
    data(1,subj_el).ppg.fs = old_data{1,subj_el}.waves.fs;
    % insert EKG signal
    ekg = old_data{1,subj_el}.waves.II;
    data(1,subj_el).ekg.v = ekg.v;
    clear ekg
    data(1,subj_el).ekg.fs = old_data{1,subj_el}.waves.fs;
    % insert RESP signal
    resp = old_data{1,subj_el}.waves.RESP;
    data(1,subj_el).ref.resp_sig.imp.v = resp.v;
    clear resp
    data(1,subj_el).ref.resp_sig.imp.fs = old_data{1,subj_el}.waves.fs;
    data(1,subj_el).ref.resp_sig.imp.method = 'thoracic impedance measured using clinical monitor';

    % insert numerics
    data(1,subj_el).ref.params.rr.t = old_data{1,subj_el}.numerics.RESP.t;
    data(1,subj_el).ref.params.rr.v = old_data{1,subj_el}.numerics.RESP.v;
    data(1,subj_el).ref.params.rr.method = 'derived from thoracic impedance measured using clinical monitor';
    data(1,subj_el).ref.params.rr.units.t = 's';
    data(1,subj_el).ref.params.rr.units.v = 'breaths/min';
    data(1,subj_el).ref.params.hr.t = old_data{1,subj_el}.numerics.HR.t;
    data(1,subj_el).ref.params.hr.v = old_data{1,subj_el}.numerics.HR.v;
    data(1,subj_el).ref.params.hr.method = 'derived from ecg measured using clinical monitor';
    data(1,subj_el).ref.params.hr.units.t = 's';
    data(1,subj_el).ref.params.hr.units.v = 'beats/min';
    data(1,subj_el).ref.params.pr.t = old_data{1,subj_el}.numerics.PULSE.t;
    data(1,subj_el).ref.params.pr.v = old_data{1,subj_el}.numerics.PULSE.v;
    data(1,subj_el).ref.params.pr.method = 'derived from ppg measured using clinical monitor';
    data(1,subj_el).ref.params.pr.units.t = 's';
    data(1,subj_el).ref.params.pr.units.v = 'beats/min';

end

% Save to file
save([up.paths.data_save_folder, 'mimicii_data'], 'data')

end

function results = sub_group_identification(up)

%% Load data
load([up.paths.data_save_folder, 'mimicii_data'], 'data')

%% Split into sub-groups and identify data for each
sub_group = zeros(length(data),1);   % zero for neonates, one for adults
n_rrs = []; a_rrs = [];
for s = 1 : length(data)
    if strcmp(data(1,s).group, 'neonate')
        n_rrs = [n_rrs; data(1,s).ref.params.rr.v];
    else
        a_rrs = [a_rrs; data(1,s).ref.params.rr.v];
        sub_group(s) = 1;
    end
end
results.sub_group = sub_group;

%% Plot histograms

% generate data of prop of rrs in each category:
xbins1 = -1:120;
rrs.a.tot = sum(~isnan(a_rrs));
rrs.a.bins = nan(length(xbins1),1);
rrs.n.tot = sum(~isnan(n_rrs));
rrs.n.bins = nan(length(xbins1),1);
for s = 1:length(xbins1)
    rel_rr_range = [(xbins1(s)-0.5),(xbins1(s)+0.5)];
    rel_a_rrs = a_rrs(a_rrs>rel_rr_range(1) & a_rrs<=rel_rr_range(2));
    rel_n_rrs = n_rrs(n_rrs>rel_rr_range(1) & n_rrs<=rel_rr_range(2));
    rrs.a.bins(s) = length(rel_a_rrs)/rrs.a.tot;
    rrs.n.bins(s) = length(rel_n_rrs)/rrs.n.tot;
end

% bar chart
ftsize = 14;
paper_size = [9, 3.5];
figure('Position', [200, 200, 100*paper_size(1), 100*paper_size(2)])
b = bar(xbins1, [rrs.a.bins, rrs.n.bins], 1.2);
b(1).FaceColor = [0,0,1];
b(1).EdgeColor = [0,0,1];
b(2).FaceColor = [1, 0.6, 0.6];
b(2).EdgeColor = [1, 0.6, 0.6];
xlim([-1 xbins1(end)])
ylim([0 1.1*max([rrs.a.bins; rrs.n.bins])])
xlabel('Impedance RR numeric [bpm]', 'FontSize', ftsize)
ylabel('Proportion of Measurements', 'FontSize', ftsize)
set(gca, 'YTick', [], 'FontSize', ftsize)
legend(b, {'Adults', 'Neonates'})

%% save histograms
save_name = 'adults_vs_neonates';
savepath = [up.paths.plots, save_name];
PrintFigs(gcf, paper_size, savepath, up)

end

function PrintFigs(h, paper_size, savepath, up)
set(h,'PaperUnits','inches');
set(h,'PaperSize', [paper_size(1), paper_size(2)]);
set(h,'PaperPosition',[0 0 paper_size(1) paper_size(2)]);
print(h,'-dpdf',savepath)
print(h,'-dpng',savepath)

if up.eps_figs
    export_fig(savepath, '-eps')
end
close all
end

author: peterhcharlton

repo: RRest


text:
function errMsg = validateNumber(x, numberName, lb, ub, isInt, prevErrMsgs, enteredStr)
%validateNumber Summary of this function goes here
%   Detailed explanation goes here

    errMsg=prevErrMsgs;
    if(not(isnumeric(x)) || isnan(x))
        enteredStr = ['(Entered: ', enteredStr, ')'];
        errMsg{end+1} = [numberName, ' must be numeric.  ', enteredStr];
        return;
    end

    enteredStr = ['(Entered: ', num2str(x), ')'];

    if(isInt==true)
        if(not(ceil(x) == floor(x)))
            errMsg{end+1} = [numberName, ' must be an integer.  ', enteredStr];
        end
    end

    if(x<lb || x>ub)
        rngStr = ['[',num2str(lb), ', ', num2str(ub)   ,']'];
        errMsg{end+1} = [numberName, ' must be within the range ', rngStr, '.  ', enteredStr];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function y=synth(freq,dur,amp,Fs,type)
% y=synth(freq,dur,amp,Fs,type)
%
% Synthesize a single note
%
% Inputs:
%  freq - frequency in Hz
%  dur - duration in seconds
%  amp - Amplitude in range [0,1]
%  Fs -  sampling frequency in Hz
%  type - string to select synthesis type
%         current options: 'fm', 'sine', or 'saw'

% Copyright (c) 2009 Ken Schutte
% more info at: http://www.kenschutte.com/midi

if nargin<5
  error('Five arguments required for synth()');
end

N = floor(dur*Fs);

if N == 0
  warning('Note with zero duration.');
  y = [];
  return;

elseif N < 0
  warning('Note with negative duration. Skipping.');
  y = [];
  return;
end

n=0:N-1;
if (strcmp(type,'sine'))
  y = amp.*sin(2*pi*n*freq/Fs);

elseif (strcmp(type,'saw'))

  T = (1/freq)*Fs;     % period in fractional samples
  ramp = (0:(N-1))/T;
  y = ramp-fix(ramp);
  y = amp.*y;
  y = y - mean(y);

elseif (strcmp(type,'fm'))

  t = 0:(1/Fs):dur;
  envel = interp1([0 dur/6 dur/3 dur/5 dur], [0 1 .75 .6 0], 0:(1/Fs):dur);
  I_env = 5.*envel;
  y = envel.*sin(2.*pi.*freq.*t + I_env.*sin(2.*pi.*freq.*t));

else
  error('Unknown synthesis type');
end

% smooth edges w/ 10ms ramp
if (dur > .02)
  L = 2*fix(.01*Fs)+1;  % L odd
  ramp = bartlett(L)';  % odd length
  L = ceil(L/2);
  y(1:L) = y(1:L) .* ramp(1:L);
  y(end-L+1:end) = y(end-L+1:end) .* ramp(end-L+1:end);
end

author: kts

repo: matlab-midi

save_path: github-repos/MATLAB/kts-matlab-midi


text:
function [normed_traindata,normed_testdata ] = normalization(traindata,testdata)
% Calculating mean and standard deviation of train samples on each band
% Using calculated mean and standard deviation to perform normalization on
% both train samples and test samples

mean_val=mean(traindata,1);
sigma_val=std(traindata,0,1);

traindata=bsxfun(@minus,traindata,mean_val);
normed_traindata=bsxfun(@rdivide,traindata,sigma_val);

testdata=bsxfun(@minus,testdata,mean_val);
normed_testdata=bsxfun(@rdivide,testdata,sigma_val);

%selTrainData=bsxfun(@rdivide,selTrainData,sum(selTrainData,2));
%selTestData=bsxfun(@rdivide,selTestData,sum(selTestData,2));

% mean_val = mean(traindata,1);
% normed_traindata = bsxfun(@minus,traindata,mean_val);
%
% sigma_val = std(normed_traindata,0,1);
% normed_traindata = bsxfun(@rdivide,normed_traindata,sigma_val);
% normed_testdata = bsxfun(@minus,testdata,mean_val);
% normed_testdata = bsxfun(@rdivide,normed_testdata,sigma_val);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [M_right_up, M_right_down] = index_Matrix(N)
n = log2(N);
M_right_up = zeros(N/2, n);
M_right_down = zeros(N/2, n);
for i = 1 : n
    for j = 1 : 2^(i - 1)
        M_right_up((j - 1) * N/2^i + 1 : j * N/2^i, i) = (1 : N/2^i)' + (j - 1) * N/2^(i - 1);
    end
    M_right_down(:, i) = M_right_up(:, i) + 2^(n - i);
end
M_right_up = M_right_up(:, end : -1 : 1);
M_right_down = M_right_down(:, end : -1 : 1);

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
% Clear and close everything
clear all;
close all;

% Add paths needed
addpath('functions/matGeom/polygons2d/')
addpath('functions/matGeom/geom3d/')
addpath('functions/lips/')

%% CONFIG INFORMATION IS ALL IN HERE

pathplanes = '../input/floorplan_spencer_small.txt';
pathpath = '../input/path_spencer_small_01.txt';

planeheight = 3; %meters

imurate = 200; %hz
lidarrate = 0.5; %hz

totalruntime = 240; %seconds

%% ROBOT TRAJECTORY LOADING AND PROCESSING
% Load our data from file
data = csvread(pathpath,1,0);
%pathtime = linspace(0,totalruntime,size(path_IMUinG,1));
pathtime = totalruntime*data(:,1)/max(data(:,1));
path_IMUinG = data(:,2:end);

% NOTE: Convert from feet poses to meters
path_IMUinG(:,1:3) = 0.3048.*path_IMUinG(:,1:3);

% Convert from the degrees to radians
path_IMUinG(:,4:6) = pi/180.*path_IMUinG(:,4:6);

% Create a spline for each dimension x,y,z and roll,pitch,yaw
% This will be used to evaluate the derivative to give us accelerations
% Can find the value at any point using "ppval"
pp_IMUinG = {};
for k=1:size(path_IMUinG,2)
    pp_IMUinG{end+1} = spline(pathtime,[0; path_IMUinG(:,k); 0]);
    %pp_IMUinG{end+1} = spline(pathtime,path_IMUinG(:,k));
    %pp_IMUinG{end+1} = csapi(pathtime,path_IMUinG(:,k));
    %pp_IMUinG{end+1} = pchip(pathtime,path_IMUinG(:,k));
    %pp_IMUinG{end+1} = csape(pathtime,path_IMUinG(:,k))
end

% Our time for the states is a combo of our different sensors
% Here we will combine the two rates, so everything is in "sync"
timeimu = linspace(0,totalruntime,imurate*totalruntime);
timelidar = linspace(0,totalruntime,lidarrate*totalruntime);

%% PLANES AND ENVIROMENT DATA LOADING AND PROCESSING
% Load our data from file (convert from feet to meters)
planes2d_PinG = 0.3048.*load_2dplanedata(pathplanes);

% Convert to polygons with a height of 8ft
planes3d_PinG = planes2dtopolygons3d(planes2d_PinG, planeheight);

%% VISULIZATION AND PLOTTING OF DATA FOR ANALYSIS
% Plot the polyons (uses geom3d function)
fh = figure(2);

fontsize = 20;
set(gcf,'PaperPositionMode','auto');
set(gcf,'defaultuicontrolfontsize',fontsize);
set(gcf,'defaultuicontrolfontname','Bitstream Charter');
set(gcf,'DefaultAxesFontSize',fontsize);
set(gcf,'DefaultAxesFontName','Bitstream Charter');
set(gcf,'DefaultTextFontSize',fontsize);
set(gcf,'DefaultTextFontname','Bitstream Charter');

clf(fh)
for ii=1:size(planes3d_PinG,2)
	drawPolygon3d(planes3d_PinG{ii}(:,1),planes3d_PinG{ii}(:,2),planes3d_PinG{ii}(:,3),'b');
    hold on;
end

% Plot path points and their IDs
for ii=1:size(path_IMUinG,1)
    drawPoint3d(path_IMUinG(ii,1),path_IMUinG(ii,2),path_IMUinG(ii,3),'or');
    drawCoordinates3d([path_IMUinG(ii,1),path_IMUinG(ii,2),path_IMUinG(ii,3)], [path_IMUinG(ii,4),path_IMUinG(ii,5),path_IMUinG(ii,6)], 1.5)
    hold on;
end
%text(path_IMUinG(:,1),path_IMUinG(:,2),path_IMUinG(:,3),[repmat('  ',size(path_IMUinG,1),1), num2str((1:size(path_IMUinG,1))')])
hold on;

% Plot the spline (use IMU rate so it is smooth)
plot3(ppval(pp_IMUinG{1},timeimu),ppval(pp_IMUinG{2},timeimu),ppval(pp_IMUinG{3},timeimu),'r');
hold on;

% Plot lidar poses over the spline
for time=timelidar
    %drawPoint3d(ppval(pp_IMUinG{1},time),ppval(pp_IMUinG{2},time),ppval(pp_IMUinG{3},time),'og');
    hold on;
    drawCoordinates3d([ppval(pp_IMUinG{1},time),ppval(pp_IMUinG{2},time),ppval(pp_IMUinG{3},time)],...
        [ppval(pp_IMUinG{4},time),ppval(pp_IMUinG{5},time),ppval(pp_IMUinG{6},time)], 0.75)
    hold on;
end

% Draw the global axis
drawCoordinates3d([0,0,0],[0,0,0],2)

% Do the labels and the such...
axis equal
xlabel('x-axis (m)')
ylabel('y-axis (m)')
zlabel('z-axis (m)')
view([0 90])

%% PLOT SAVE TO FILE FOR PAPER FIGURES
save_to_file = 0;
if save_to_file
    set(gcf,'Position',[0 0 1200 400])
    view([-25 25]);
    set(get(gca,'ylabel'),'rotation',-45)
    print(fh,'-dpng','-r500','trajectory_3d.png')
    view([0 90]);
    set(get(gca,'ylabel'),'rotation',90)
    print(fh,'-dpng','-r500','trajectory_2d_top.png')
    view([0 0]);
    print(fh,'-dpng','-r500','trajectory_2d_side.png')
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function c = tapas_hgf_binary_pu_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF)
% for binary inputs in the *presence* of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading exactly to the model introduced in
% Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_hgf_config.
%
% This file refers to UNCERTAIN inputs (Eqs 45-47 in Mathys et al., (2011));
% for inputs without uncertainty, refer to tapas_hgf_binary_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_binary_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.rho        row vector of rhos (representing drift; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%         est.p_prc.al         scalar alpha (perceptual uncertainty)
%         est.p_prc.eta0       scalar eta0 (mean of first input category)
%         est.p_prc.eta1       scalar eta1 (mean of second input category)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels)
%         est.traj.sa          sigma (rows: trials, columns: levels)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_hgf_binary_pu_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - When analyzing a new dataset, take your inputs u and use
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'hgf_binary_pu';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 0, 0];

c.logsa_0mu = [NaN,   log(0.1), log(1)];
c.logsa_0sa = [NaN,          0,      0];

% Rhos
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero to turn off drift.
c.rhomu = [NaN, 0, 0];
c.rhosa = [NaN, 0, 0];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,      0];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,  -3,  -6];
c.omsa = [NaN, 4^2, 4^2];

% Alpha
% Format: scalar.
c.logalmu = log(0.5);
c.logalsa = 1;

% Eta0
% Format: scalar.
c.eta0mu = 0;
c.eta0sa = 0;

% Eta1
% Format: scalar.
c.eta1mu = 1;
c.eta1sa = 0;

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.rhomu,...
    c.logkamu,...
    c.ommu,...
    c.logalmu,...
    c.eta0mu,...
    c.eta1mu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.rhosa,...
    c.logkasa,...
    c.omsa,...
    c.logalsa,...
    c.eta0sa,...
    c.eta1sa,...
         ];

% Check whether we have the right number of priors
expectedLength = 3*c.n_levels+2*(c.n_levels-1)+4;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_hgf_binary_pu;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_binary_pu_transp;

return;


text:
function Conn = hlp_subtractConnectivity(Conn1,Conn2)
    % return a Connectivity object containing the difference in
    % connectivity between two datasets

    Conn = Conn1;

    connmethods = hlp_getConnMethodNames(Conn1);

    for c=1:length(connmethods)
        Conn.(connmethods{c}) = Conn1.(connmethods{c}) - Conn2.(connmethods{c});
    end

    otherfields = setdiff_bc(fieldnames(Conn1),connmethods);
    for f=1:length(otherfields)
        Conn.(otherfields{f}) = Conn1.(otherfields{f});
    end

repo: REST


text:
function [insideNodes,insideNodeStruct] = floodFillFindPerim2(mesh,perimDist,startNode)
% A flood-fill method to find the mesh perimeter
%
% [insideNodes,insideNodeStruct]=floodFillFindPerim2(mesh,perimDist,startNode,busyHandle)
%
% Find expanding rings connected to the start node.
% Stop when all the members of the ring exceed perimDist from the startNode
% But cleverer than this :)
%
% Also stores the average distance of each new set and its offset - May be
% useful for tacking down points later This can break quite easily - for
% example when the floodfill runs around all sides of a bump. But we can
% fix this later  ...
%
% AUTHOR:  Wade
% DATE : 020701 last modified

allWithinDist = 1;
insideNodes   = startNode;
currentNodes  = startNode;
counter       = 0;
insideNodeStruct.offset = 0;
insideNodeStruct.avDist = 0;
nVerts = length(mesh.connectionMatrix);

% Limit to generate no more that 10000 rings....
while ((allWithinDist) && (counter<10000))

   % What nodes are connected to the current ones?
   [newRows connected]=find(mesh.connectionMatrix(currentNodes,:));

   if (~isempty(connected))

      connected=unique(connected(:));
      insideNodes=[insideNodes;connected];

      currentNodes=connected;
	  notCnodes=setdiff(1:nVerts,currentNodes);

	  % This is a cute way of zeroing columns in a sparse matrix. Much
	  % faster than foo(:,currentNodes)=0;
	  diagMat=sparse(notCnodes,notCnodes,ones(length(notCnodes),1),nVerts,nVerts);
	  mesh.connectionMatrix=(mesh.connectionMatrix)*(diagMat);

	  nodeDists=mesh.dist(currentNodes);
      f=find(nodeDists<perimDist);
      allWithinDist=(sum(f(:)));

      % Average distance
      insideNodeStruct.avDist=[insideNodeStruct.avDist;mean(nodeDists(:))];

      % Offset to the list of current rings
      insideNodeStruct.offset=[insideNodeStruct.offset;length(connected)];

   else
         allWithinDist=0;
   end

   counter=counter+1;

end

insideNodes=unique(insideNodes);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%DYADICM Dyadic dataset mapping
%
%   B = [C1,C2]*DYADICM([],P,Q,SIZE)
%   B = [C1,C2]*DYADICM([],FNAME,[],SIZE)
%   W = A*DYADICM([],{U1,U2,FNAME},[],SIZE)
%   B = A*DYADICM([],{V1,V2,FNAME},[],SIZE)
%
% INPUT
%  C1,C2   Datasets / datafiles to be combined
%  A       Input dataset used for training or classification
%  P,Q     Scalar multiplication factors (default 1) to compute P*A1+Q*A2
%          or string (name of a routine)
%  FNAME   String with the name of a function to combine two datasets or
%          mappings, 'plus' for PLUS(U1,U2)
%  U1,U2   Untrained mappings to be combined
%  V1,V2   Trained or fixe mappings to be combined
%  SIZE    Desired images size of output dataset objects
%
% OUTPUT
%  B       Dataset
%  W       Combined mapping
%
% DESCRIPTION
% This special mapping is a low-level routine to facilitate dyadic
% operations on datafiles, datasets and mappings. Datasets to be combined
% should have the same number of objects. Image objects should have the
% same image size. Datafiles should be preprocessed or postprocessed
% versions of the same original datafile. Mappings should have the same
% input and output size.
%
% This routine has been written for use by PRTools programmers only. Users
% are discouraged to call it directly. The routine is called by the
% dyadic operations of the classes 'prmapping' and 'prdatafile'.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, DATAFILES

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = dyadicm(a,p,q,fsize)

  if nargin < 4, fsize = []; end
  if nargin < 3, q = []; end
  if nargin < 2, p = []; end

  if nargin < 1 | isempty(a)
    % note that this routine is sometimes externally set to a combiner
    b = prmapping(mfilename,'fixed',{p,q,fsize});
    b = setsize_out(b,prod(fsize));
    b = setname(b,'dyadicm');
    return
  end

  if isempty(p)
    p = 1;
  end

  if iscell(a) & ismapping(a{1})
    % combining mappings. We are just here because the dyadic mapping
    % operations like mapping/plus.m are programmed like this.
    if isfixed(a{1})
      b = prmapping(mfilename,'fixed',{{a{1},a{2},p}});
    elseif isuntrained(a{1})
      b = prmapping(mfilename,'untrained',{{a{1},a{2},p},[],[]});
    elseif istrained(a{1})
      if isempty(fsize), fsize = getsize(a{1},2); end
      b = prmapping(mfilename,'trained',{a{1},a{2},p},getlabels(a{1}),getsize(a{1},1));
    else
      error('Wrong type')
    end

  elseif isdatafile(a)
    nodatafile; % forces to store this routine as postprocessing

  elseif iscell(p) & isuntrained(p{1}) % train mappings stored in p
    v1 = p{1};     % untrained mapping
    v2 = p{2};     % fixed or trained mapping or constant
    p = p{3};      % operation
    v1 = a*v1;     % train first mapping
    [kin,kout] = size(v1);
    if ismapping(v2) & (isfixed(v2) | istrained(v2))
      ;            % leave it and hope for the best
    elseif ismapping(v2) & isuntrained(v2)
      v2 = a*v2;   % train it
    else           % should be scalar constant
      ;            % leave it
    end
    % store a standard trained mapping
    b = prmapping(mfilename,'trained',{v1,v2,p},getlabels(v1),kin,kout);

  elseif iscell(p) % execute mappings stored in p
    v1 = p{1};     % fixed or trained mapping
    v2 = p{2};     % fixed or trained mapping or constant
    p = p{3};      % operation
    a1 = a*v1;     % prepare datasets
    if ismapping(v2)
      a2 = a*v2;
    else
      a2 = setdata(a1,v2*ones(size(a1)));
    end
    b = feval(p,a1,a2);

  elseif ismapping(p) % execute standard PRTools trained mapping
    v1 = getdata(p,1);     % fixed or trained mapping
    v2 = getdata(p,2);     % fixed or trained mapping or constant
    p = getdata(p,3);      % operation
    a1 = a*v1;     % prepare datasets
    if ismapping(v2)
      a2 = a*v2;
    else
      a2 = setdata(a1,v2*ones(size(a1)));
    end
    b = feval(p,a1,a2);

  else % the basic call by data and parameters

    % The dataset A has horizontally to be split in two datasets A1, A2.
    % If P and Q are scalars or if Q = [], this is done half-half.
    % If P is a string (name of a routine) it is used to combine A1 and A2.

    [a1,a2,fsize] = split_dataset(a,p,q,fsize);

    if ischar(p) % function name in p
      b = feval(p,a1,a2);
    else % addition
      if isempty(q), q = 1; end
      b = a1*p + a2*q;
    end

    if isdataset(b)
      b = setfeatsize(b,fsize);
    end

  end

return

function [a1,a2,fsize] = split_dataset(a,p,q,fsize)

  if iscell(a)

    a1 = a{1};
    a2 = a{2};

  %elseif ischar(p) & ~isempty(q)
%   elseif ~isempty(q)
%     % split as defined by q
%     fsize = q;
%     a1 = a(:,1:q);
%     a2 = a(:,q+1:end);

  else
    % 50-50 split, but preserve possible image structure
    if isempty(fsize) | fsize == 0
      if isdataset(a)
        fsize = getfeatsize(a);
      else % may be a is set of images
        a = double(a);
        fsize = size(a);
      end
      while (fsize(end) == 1) & (length(fsize) > 1)
        fsize = fsize(1:end-1);
      end
      fsize(end) = fsize(end)/2;
    end
    k = size(a,2);

    if k ~= 2*floor(k/2)
      error('Feature size of dataset should be multiple of 2')
    end

    a1 = a(:,1:k/2);
    a2 = a(:,k/2+1:k);
    if isdataset(a)
      a1 = setfeatsize(a1,fsize);
      a2 = setfeatsize(a2,fsize);
    end

  end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit