text:
%% brute_force_tune
% Code to tune parameters by brute force. This is for heading init.
% Obviously theoretically completely silly.
%
% Works sorta like RANSAC I guess?
%
% Adam Werries 2016, see Apache 2.0 license.

k_max = 60;
% Specify ranges
roll = linspace(-45,45,100);
pitch = linspace(-45,45,100);
yaw = linspace(-180,180,100);
% Repeat arrays
roll = repmat(roll, [1 k_max]);
pitch = repmat(pitch, [1 k_max]);
yaw = repmat(yaw, [1 k_max]);
% Generate random selections of each vector
num_items = length(roll);
roll_i = randperm(num_items);
pitch_i = randperm(num_items);
yaw_i = randperm(num_items);
rms_error_filter = zeros(1,num_items);
max_error_filter = zeros(1,num_items);
parfor i = 1:num_items
    fprintf('Iteration: %d, Roll: %08.5f, Pitch: %08.5f, Yaw: %08.5f\n', i, roll(roll_i(i)), pitch(pitch_i(i)), yaw(yaw_i(i)));
    temp_cond = init_cond;
    temp_cond(7:9) = [deg2rad(roll(roll_i(i))) deg2rad(pitch(pitch_i(i))) deg2rad(yaw(yaw_i(i)))];
    [out_profile,out_IMU_bias_est,out_KF_SD] = Loosely_coupled_INS_GNSS(temp_cond, filter_time, epoch, lla, novatel, imu, LC_KF_config, est_IMU_bias);
    xyz = out_profile(:,2:4);
    llh = ecef2lla(xyz);
    [x,y] = deg2utm(llh(:,1),llh(:,2));
    x = x-min_x;
    y = y-min_y;

    distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
    rms_error_filter(i) = rms(distance);
    max_error_filter(i) = max(distance);
end

[minmax, i] = min(max_error_filter);
fprintf('\nBest max: %08.5f, rms is %08.5f\n', minmax, rms_error_filter(i));
fprintf('Best iteration for max: %d, Roll: %08.5f, Pitch: %08.5f, Yaw: %08.5f\n', i, roll(roll_i(i)), pitch(pitch_i(i)), yaw(yaw_i(i)));
[minrms, i] = min(rms_error_filter);
fprintf('Best rms: %08.5f, max is %08.5f\n', minrms, max_error_filter(i));
fprintf('Best iteration for rms: %d, Roll: %08.5f, Pitch: %08.5f, Yaw: %08.5f\n', i, roll(roll_i(i)), pitch(pitch_i(i)), yaw(yaw_i(i)));
[minrms, i] = min((rms_error_filter+max_error_filter)/2);
fprintf('Best average of RMS and max: %08.4f, rms is  %08.4f, max is %08.4f\n', minrms, rms_error_filter(i), max_error_filter(i));
fprintf('Best iteration for rms: %d, Roll: %08.5f, Pitch: %08.5f, Yaw: %08.5f\n', i, roll(roll_i(i)), pitch(pitch_i(i)), yaw(yaw_i(i)));

author: awerries

repo: kalman-localization


text:
function kern = whiteblockKernExpandParam(kern, params)

% WHITEBLOCKKERNEXPANDPARAM Fill WHITEBLOCK kernel structure with params.
% FORMAT
% DESC returns a white noise block kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : whiteblockKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

kern.variance = params;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_bug1508

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqanalysis

% Stan reported a strange error caused by the following:
% if the cfg.channel in a call to ft_freqanalysis contains channels that
% are not in the data, the function crashes on the second trial

%% Try to reproduce
data = [];
data.trial = {randn(3,100) randn(3,100)};
data.time  = {1:100 1:100};
data.label = {'chan1';'chan2';'chan3'};

cfg = [];
cfg.method  = 'mtmfft';
cfg.channel = {'chan4'};
cfg.output  = 'pow';
cfg.taper   = 'hanning';

try
  freq = ft_freqanalysis(cfg, data);
catch err
  if strcmp(err.message,'no channels were selected')
    % error is expected handling
  else
    error('test doesn''t result in expected error message')
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function result = tetra_unit_sum ( func, norder, xtab, ytab, ztab, weight )

%*****************************************************************************80
%
%% TETRA_UNIT_SUM carries out a quadrature rule in the unit tetrahedron in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that
%
%      0 <= X,
%      0 <= Y,
%      0 <= Z, and
%      X + Y + Z <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function of three variables which is to be integrated,
%    of the form:
%      function value = func ( x, y, z )
%
%    Input, integer NORDER, the order of the rule.
%
%    Input, real XTAB(NORDER), YTAB(NORDER), ZTAB(NORDER), the
%    abscissas.
%
%    Input, real WEIGHT(NORDER), the weights.
%
%    Output, real RESULT, the approximate integral of the function.
%
  quad = 0.0E+00;

  for i = 1 : norder
    quad = quad + weight(i) * feval ( func, xtab(i), ytab(i), ztab(i) );
  end

  volume = tetra_unit_volume ( );
  result = quad * volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% TEST_MAXWELL_CUBE_G_NMNN: data function for Neumann boundary condition.

function g = test_maxwell_cube_g_nmnn (x, y, z, ind)

  g = zeros ([3, size(x)]);
  switch (ind)
    case 1
      g(2,:,:) = -(exp(x) .* cos(y) - z .* cos(y));
      g(3,:,:) = sin(y) + exp(z) .* sin(x);
    case 2
      g(2,:,:) = exp(x) .* cos(y) - z .* cos(y);
      g(3,:,:) = -(sin(y) + exp(z) .* sin(x));
    case 3
      g(1,:,:) = exp(x) .* cos(y) - z .* cos(y);
    case 4
      g(1,:,:) = -(exp(x) .* cos(y) - z .* cos(y));
    case 5
      g(1,:,:) = -(sin(y) + exp(z) .* sin(x));
    case 6
      g(1,:,:) = sin(y) + exp(z) .* sin(x);
    otherwise
      error ('g_nmnn: unknown reference number')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function []=addColormapButton(varargin)
%% Viewing figures: changing colormap
% addColormapButton
% addColormapButton(hf)

% This script creates a new pushtool in a figure's toolbar to manipulate
% the color scheme displayed in the figure images
%%
switch nargin
    case 0
        hf=gcf;
    case 1
        hf=varargin{1};
end

% Create icon
S=zeros(16,16,3);
S(1,:,:)=parula(16);
S(2,:,:)=parula(16);
S(3,:,:)=parula(16);
S(4,:,:)=parula(16);
S(5,:,:)=jet(16);
S(6,:,:)=jet(16);
S(7,:,:)=jet(16);
S(8,:,:)=jet(16);
S(9,:,:)=gray(16);
S(10,:,:)=gray(16);
S(11,:,:)=gray(16);
S(12,:,:)=gray(16);
S(13,:,:)=coldwarm(16);
S(14,:,:)=coldwarm(16);
S(15,:,:)=coldwarm(16);
S(16,:,:)=coldwarm(16);

hb = findall(hf,'Type','uitoolbar');

% Create a uipushtool in the toolbar
uipushtool(hb(1),'TooltipString','colormap','CData',S,'Tag','colormap_button','ClickedCallback',{@colormapFunc,{hf}});

end

%% Colormap function  colormapFunc

function colormapFunc(~,~,inputCell)

hf = inputCell{1};
prompt = {'Colormap'};
title = '';
format = struct('type','list');
format.style = 'popupmenu';
format.items = {' -- Select a color scheme -- ','coldwarm','parula','jet','hsv','hot','cool','gray','bone'};
default = cell(size(prompt,1),1);
default{1,1} = 1; % Default selection will always be instruction item
default = cell2struct(default,prompt,1);
prompt = repmat(prompt,1,2);
options.AlignControls = 'on';
choice = inputsdlg(prompt,title,format,default,options);

% Change color scheme for all existing images accordingly
if choice.Colormap == 2
    colormap(0.8*coldwarm);
elseif choice.Colormap == 3
    colormap(parula);
elseif choice.Colormap == 4
    colormap(jet);
elseif choice.Colormap == 5
    colormap(hsv);
elseif choice.Colormap == 6
    colormap(hot);
elseif choice.Colormap == 7
    colormap(cool);
elseif choice.Colormap == 8
    colormap(gray);
elseif choice.Colormap == 9
    colormap(bone);
end
end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% Modified by Rana Odabas 2018
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function lissajous ( )

%*****************************************************************************80
%
%% LISSAJOUS uses MATLAB to draw a closed planar curve.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LISSAJOUS:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Make a plane curve by connecting a series of points.\n' );

  xy = load ( 'lissajous.txt' );
%
%  Plot the data.
%
  plot ( xy(:,1), xy(:,2), 'LineWidth', 2, 'Color', 'm' );
  hold on
%
%  To avoid clutter, only plot every 10th point.
%
  plot ( xy(1:10:end,1), xy(1:10:end,2), 'k.' );
  grid on
  axis ( [ -1.2, +1.2, -1.2, +1.2 ] )
  axis equal
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( 'Lissajous, x=sin(3t+pi/2), y=sin(4t)' );
  hold off

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LISSAJOUS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function enter_eul_angs(varargin)
% update because euler angles were changed
global ts rhs lr rhs1 lr1 xs xst ys yst zs zst ks phi theta psi rad deg x y z gamma delta alpha  vna vnat N Nt sxyz sav sN arcR arcth arctxtR phia1 phia2 phiat sa thetaa1 thetaa2 thetaat psia1 psia2 psiat gammaa1 gammaa2 gammaat saa  deltaa1 deltaa2 deltaat arcvsh alphaa1 alphaa2 alphaat arcR1 arcth1 arctxtR1

naner=false;

phin=str2num(get(phi,'String'));
if length(phin)==0
    naner=true;
end

thetan=str2num(get(theta,'String'));
if length(thetan)==0
    naner=true;
end

psin=str2num(get(psi,'String'));
if length(psin)==0
    naner=true;
end

if naner
    nan_error;
else

    if get(deg,'Value')
        phin=pi*phin/180;
        thetan=pi*thetan/180;
        psin=pi*psin/180;
    end

    an=ea_bounding(phin,thetan,psin); % correct angles
    if an(1)
        phin=an(2);
        thetan=an(3);
        psin=an(4);
        if get(deg,'Value')
            set(phi,'String',num2str(180*phin/pi));
            set(theta,'String',num2str(180*thetan/pi));
            set(psi,'String',num2str(180*psin/pi));
        else
            set(phi,'String',num2str(phin));
            set(theta,'String',num2str(thetan));
            set(psi,'String',num2str(psin));
        end
    end

    Ms=matrices(phin,thetan,psin);
    M=Ms{1}*Ms{2}*Ms{3};

    als=0.5; % transparensy of xyz
    lcs=[0.4 0.4 0.4]; %labels color

    if get(sxyz,'value')
        % x
        xsv=M*[1;0;0];
        xsv1=ks*xsv;
        %xs=arrowa(0,0,0,xsv1(1),xsv1(2),xsv1(3),rhs1*0.8,lr1*0.8,[1 0 0],als,hpar);
        arrow_update(xs,0,0,0,xsv1(1),xsv1(2),xsv1(3),rhs1*0.8,lr1*0.8);
        %xst=text('parent',hpar,'position',xsv1+ts*xsv,'string','x','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(xst,'position',xsv1+ts*xsv);

        % y
        ysv=M*[0;1;0];
        ysv1=ks*ysv;
        %ys=arrowa(0,0,0,ysv1(1),ysv1(2),ysv1(3),rhs1*0.8,lr1*0.8,[0 1 0],als,hpar);
        arrow_update(ys,0,0,0,ysv1(1),ysv1(2),ysv1(3),rhs1*0.8,lr1*0.8);
        %yst=text('parent',hpar,'position',ysv1+ts*ysv,'string','y','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(yst,'position',ysv1+ts*ysv);

        % z
        zsv=M*[0;0;1];
        zsv1=ks*zsv;
        %zs=arrowa(0,0,0,zsv1(1),zsv1(2),zsv1(3),rhs1*0.8,lr1*0.8,[0 0 1],als,hpar);
        arrow_update(zs,0,0,0,zsv1(1),zsv1(2),zsv1(3),rhs1*0.8,lr1*0.8);
        %zst=text('parent',hpar,'position',zsv1+ts*zsv,'string','z','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(zst,'position',zsv1+ts*zsv);

        arrow_visible_off_on(xs,true);
        set(xst,'visible','on');

        arrow_visible_off_on(ys,true);
        set(yst,'visible','on');

        arrow_visible_off_on(zs,true);
        set(zst,'visible','on');

    else

        arrow_visible_off_on(xs,false);
        set(xst,'visible','off');

        arrow_visible_off_on(ys,false);
        set(yst,'visible','off');

        arrow_visible_off_on(zs,false);
        set(zst,'visible','off');

    end

    if (nargin>=1)&&varargin{1}
        vn(1,1)=str2num(get(x,'string'));
        vn(2,1)=str2num(get(y,'string'));
        vn(3,1)=str2num(get(z,'string'));
        gamman=varargin{2};
        deltan=varargin{3};
        alphan=varargin{4};
    else
        % axis v
        axan=euler2axan(phin,thetan,psin,M);
        % {{gamma,delta,alpha},v}
        gda=axan{1};
        gamman=gda{1};
        deltan=gda{2};
        alphan=gda{3};
        vn=axan{2};
        if get(deg,'Value')
            set(gamma,'string',num2str(180*gamman/pi));
            set(delta,'string',num2str(180*deltan/pi));
            set(alpha,'string',num2str(180*alphan/pi));
        else
            set(gamma,'string',num2str(gamman));
            set(delta,'string',num2str(deltan));
            set(alpha,'string',num2str(alphan));
        end

        set(x,'string',num2str(vn(1)));
        set(y,'string',num2str(vn(2)));
        set(z,'string',num2str(vn(3)));
    end

    if get(sav,'value')
        vn1=ks*vn;
        %vna=arrow(0,0,0,vn(1),vn(2),vn(3),rhs1*0.8,lr1*0.8,[0.9 0.2 1],hpar);
        arrow_update(vna,0,0,0,vn(1),vn(2),vn(3),rhs1*0.8,lr1*0.8);
        %vnat=text('parent',hpar,'position',vn+ts*vn1,'string','v','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(vnat,'position',vn+ts*vn1);

        arrow_visible_off_on(vna,true);
        set(vnat,'visible','on');
    else
        arrow_visible_off_on(vna,false);
        set(vnat,'visible','off');
    end

    % line of nodes
    if get(sN,'value')
        Nv=Ms{1}*[1;0;0];
        Nv1=ks*Nv;
        %N=arrow(0,0,0,Nv(1),Nv(2),Nv(3),rhs1*0.8,lr1*0.8,[0.4 0.4 4],hpar);
        arrow_update(N,-Nv(1),-Nv(2),-Nv(3),2*Nv(1),2*Nv(2),2*Nv(3),rhs1*0.8,lr1*0.8);
        %Nt=text('parent',hpar,'position',Nv+ts*Nv1,'string','N','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(Nt,'position',Nv+ts*Nv1);

        arrow_visible_off_on(N,true);
        set(Nt,'visible','on');

    else

        arrow_visible_off_on(N,false);
        set(Nt,'visible','off');

    end

    % angles:
    if get(sa,'value')
        % phi:
        an=arc_data(phin,arcR,arcth,arctxtR);
        Ml=an{1};
        Mt=an{2};
        txv=an{3};
        %phia1=plot3(Ml(1,:),Ml(2,:),zeros(1,length(Ml(1,:))),'-k','parent',hpar);
        set(phia1,'XData',Ml(1,:),'YData',Ml(2,:),'ZData',zeros(1,length(Ml(1,:))));
        %phia2=plot3(Mt(1,:),Mt(2,:),zeros(1,length(Mt(1,:))),'-k','parent',hpar);
        set(phia2,'XData',Mt(1,:),'YData',Mt(2,:),'ZData',zeros(1,length(Mt(1,:))));
        %phiat=text('parent',hpar,'position',[txv 0],'string','\phi','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(phiat,'position',[txv 0]);

        set(phia1,'visible','on');
        set(phia2,'visible','on');
        set(phiat,'visible','on');

        % theta:
        an=arc_data(thetan,arcR,arcth,arctxtR);
        Ml1=an{1};
        Ml=Ms{1}*[zeros(1,length(Ml1(1,:))); -Ml1(2,:); Ml1(1,:)];
        Mt1=an{2};
        Mt=Ms{1}*[zeros(1,length(Mt1(1,:))); -Mt1(2,:); Mt1(1,:)];
        txv1=an{3};
        txv=Ms{1}*[0; -txv1(2); txv1(1)];
        %thetaa1=plot3(Ml(1,:),Ml(2,:),Ml(3,:),'-k','parent',hpar);
        set(thetaa1,'XData',Ml(1,:),'YData',Ml(2,:),'ZData',Ml(3,:));
        %thetaa2=plot3(Mt(1,:),Mt(2,:),Mt(3,:),'-k','parent',hpar);
        set(thetaa2,'XData',Mt(1,:),'YData',Mt(2,:),'ZData',Mt(3,:));
        %thetaat=text('parent',hpar,'position',txv,'string','\theta','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(thetaat,'position',txv);

        set(thetaa1,'visible','on');
        set(thetaa2,'visible','on');
        set(thetaat,'visible','on');

        % psi:
        an=arc_data(psin,arcR,arcth,arctxtR);
        Ml1=an{1};
        Ml=Ms{1}*Ms{2}*[Ml1(1,:); Ml1(2,:); zeros(1,length(Ml1(1,:)))];
        Mt1=an{2};
        Mt=Ms{1}*Ms{2}*[Mt1(1,:); Mt1(2,:); zeros(1,length(Mt1(1,:)))];
        txv1=an{3};
        txv=Ms{1}*Ms{2}*[txv1(1); txv1(2); 0];
        %psia1=plot3(Ml(1,:),Ml(2,:),Ml(3,:),'-k','parent',hpar);
        set(psia1,'XData',Ml(1,:),'YData',Ml(2,:),'ZData',Ml(3,:));
        %psia2=plot3(Mt(1,:),Mt(2,:),Mt(3,:),'-k','parent',hpar);
        set(psia2,'XData',Mt(1,:),'YData',Mt(2,:),'ZData',Mt(3,:));
        %psiat=text('parent',hpar,'position',txv,'string','\psi','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(psiat,'position',txv);

        set(psia1,'visible','on');
        set(psia2,'visible','on');
        set(psiat,'visible','on');

    else
        set(phia1,'visible','off');
        set(phia2,'visible','off');
        set(phiat,'visible','off');

        set(thetaa1,'visible','off');
        set(thetaa2,'visible','off');
        set(thetaat,'visible','off');

        set(psia1,'visible','off');
        set(psia2,'visible','off');
        set(psiat,'visible','off');
    end

    % axis-angle angles
    if get(saa,'value')
        % gamma:
        an=arc_data(gamman,arcR,arcth,arctxtR);
        Ml=an{1};
        Mt=an{2};
        txv=an{3};
        %gammaa1=plot3(Ml(1,:),Ml(2,:),zeros(1,length(Ml(1,:))),'-k','parent',hpar);
        set(gammaa1,'XData',Ml(1,:),'YData',Ml(2,:),'ZData',zeros(1,length(Ml(1,:))));
        %gammaa2=plot3(Mt(1,:),Mt(2,:),zeros(1,length(Mt(1,:))),'-k','parent',hpar);
        set(gammaa2,'XData',Mt(1,:),'YData',Mt(2,:),'ZData',zeros(1,length(Mt(1,:))));
        %gammaat=text('parent',hpar,'position',[txv 0],'string','\gamma','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(gammaat,'position',[txv 0]);

        set(gammaa1,'visible','on');
        set(gammaa2,'visible','on');
        set(gammaat,'visible','on');

        % delta:
        an=arc_data(deltan,arcR,arcth,arctxtR);
        csgm=cos(gamman);
        sngm=sin(gamman);
        Mrot=[csgm, -sngm, 0;
              sngm, csgm,  0;
              0,    0,     1];
        Ml1=an{1};
        Ml=Mrot*[Ml1(1,:); zeros(1,length(Ml1(1,:))); Ml1(2,:)];
        Mt1=an{2};
        Mt=Mrot*[Mt1(1,:); zeros(1,length(Mt1(1,:))); Mt1(2,:)];
        txv1=an{3};
        txv=Mrot*[txv1(1); 0; txv1(2)];
        %deltaa1=plot3(Ml(1,:),Ml(2,:),Ml(3,:),'-k','parent',hpar);
        set(deltaa1,'XData',Ml(1,:),'YData',Ml(2,:),'ZData',Ml(3,:));
        %deltaa2=plot3(Mt(1,:),Mt(2,:),Mt(3,:),'-k','parent',hpar);
        set(deltaa2,'XData',Mt(1,:),'YData',Mt(2,:),'ZData',Mt(3,:));
        %deltaat=text('parent',hpar,'position',txv,'string','\delta','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(deltaat,'position',txv);

        set(deltaa1,'visible','on');
        set(deltaa2,'visible','on');
        set(deltaat,'visible','on');

        an=arc_data(alphan,arcR1,arcth1,arctxtR1);
        csdl=cos(deltan);
        sndl=sin(deltan);
        Mrot1=[csdl, 0,  -sndl;
               0,    1,  0    ;
               sndl, 0,  csdl ];
        Ml1=an{1};
        Ml=Mrot*Mrot1*[zeros(1,length(Ml1(1,:))); -Ml1(2,:); Ml1(1,:)];
        Mt1=an{2};
        Mt=Mrot*Mrot1*[zeros(1,length(Mt1(1,:))); -Mt1(2,:); Mt1(1,:)];
        txv1=an{3};
        txv=Mrot*Mrot1*[0; -txv1(2); txv1(1)];
        %alphaa1=plot3(Ml(1,:)+vn(1)*arcvsh,Ml(2,:)+vn(2)*arcvsh,Ml(3,:)+vn(3)*arcvsh,'-k','parent',hpar);
        set(alphaa1,'XData',Ml(1,:)+vn(1)*arcvsh,'YData',Ml(2,:)+vn(2)*arcvsh,'ZData',Ml(3,:)+vn(3)*arcvsh);
        %alphaa2=plot3(Mt(1,:)+vn(1)*arcvsh,Mt(2,:)+vn(2)*arcvsh,Mt(3,:)+vn(3)*arcvsh,'-k','parent',hpar);
        set(alphaa2,'XData',Mt(1,:)+vn(1)*arcvsh,'YData',Mt(2,:)+vn(2)*arcvsh,'ZData',Mt(3,:)+vn(3)*arcvsh);
        %alphaat=text('parent',hpar,'position',txv+vn*arcvsh,'string','\alpha','HorizontalAlignment','center','VerticalAlignment','middle','color',lcs);
        set(alphaat,'position',txv+vn*arcvsh);

        set(alphaa1,'visible','on');
        set(alphaa2,'visible','on');
        set(alphaat,'visible','on');
    else
        set(gammaa1,'visible','off');
        set(gammaa2,'visible','off');
        set(gammaat,'visible','off');

        set(deltaa1,'visible','off');
        set(deltaa2,'visible','off');
        set(deltaat,'visible','off');

        set(alphaa1,'visible','off');
        set(alphaa2,'visible','off');
        set(alphaat,'visible','off');
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = pwq_yalmip(varargin)
%PWQ_YALMIP Defines a piecewise quadratic function using data from MPT
%
%Only intended for internal use in YALMIP
%
% Currently only a container for PWQ functions. Can not be
% used in actual optmization problem.

switch class(varargin{1})

    case {'struct','cell'} % Should only be called internally

        if isa(varargin{2},'double')
            % Called from YALMIP to get double
            pwastruct = varargin{1};
            x = varargin{2};
            index = varargin{5};

            val = inf;
            for i = 1:length(pwastruct)
                [ii,jj] = isinside(pwastruct{i}.Pn,x);
                if ii
                    for k = 1:length(jj)
                        Q = pwastruct{i}.Ai{jj(k)};
                        if index>1 | min(size(pwastruct{i}.Bi{jj(k)}))>1
                            % FIX: Why?? Where is this feature used
                            val = min(val,x'*Q*x + reshape(pwastruct{i}.Bi{jj(k)}(index,:),1,[])*x+pwastruct{i}.Ci{jj(k)}(index));
                        else
                            val = min(val,x'*Q*x + reshape(pwastruct{i}.Bi{jj(k)},1,[])*x+pwastruct{i}.Ci{jj(k)});
                        end
                    end
                end
            end
            if isinf(val)
                val = nan;
            end
            varargout{1} = min(val);
            return
        end

        if nargin<3
            pwaclass = 'general'
        end

        if isa(varargin{1},'struct')
            varargin{1} = {varargin{1}};
        end

        % Put in standard format
        if ~isfield(varargin{1}{1},'Bi')
            if ~isfield(varargin{1}{1},'Fi')
                error('Wrong format on input to PWQ (requires Bi or Fi)');
            else
                for i = 1:length(varargin{1})
                    varargin{1}{1}.Ai = cell(1, varargin{1}{i}.Fi);
                    varargin{1}{i}.Bi = varargin{1}{i}.Fi
                    varargin{1}{i}.Ci = varargin{1}{i}.Gi
                end
            end
        end

        if ~isfield(varargin{1}{1},'Pfinal')
            error('Wrong format on input to PWQ (requires field Pn)');
        end

        % This will be a container for binary variables in the furture
        varargin{end+1} = binvar(length(varargin{1}{1}.Pn),1);

        % Create one variable for each row
        % Inefficient but the only way currently in YALMIP
        varargout{1} = [];
        varargin{end+1} = 1;

        for i = 1:length(varargin{1}{1}.Ci{1})
            varargin{end} = i;
            varargout{1} = [varargout{1};yalmip('define',mfilename,varargin{:})];
        end

    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph

        switch varargin{1}

            case 'graph'

                varargout{1} = [];
                varargout{2} = struct('convexity','none','monotoncity','none','definiteness','none');
                varargout{3} = [];

            case {'integer','exact'}

                % FIX : Should create case for overlapping convex PWAs,
                % used in a nonconvex fashion...

                % Can only generate the first class of PWA functions
                t     = varargin{2};     % The YALMIP variables modelling this pwa
                pwq_struct = varargin{3};% MPT structure
                x     = varargin{4};     % Argument
                flag  = varargin{5};     % Type of PWA function
                d     = varargin{6};     % Binary for nonconvex cases
                index = varargin{7};     % Which row in Bix+Ci

                switch flag

                    case {'general','convex'}

                        if length(d)==1
                            % Don't introduce any binary variables when
                            % there is only one quadratic cost
                            F = ([]);
                            cost = x'*pwq_struct{1}.Ai{1}*x+pwq_struct{1}.Bi{1}*x+pwq_struct{1}.Ci{1};
                        else
                            n = length(x);
                            m = length(d);
                            z = sdpvar(n,m,'full');
                            F = (sum(d) == 1);
                            cost = 0;
                            [Mm,mm] = derivebounds(x);
                            try
                                [aux,mx,Mx] = bounding_box(pwq_struct{1}.Pfinal);
                            catch
                                Mx = 10000*ones(length(x),1);
                                mx = -10000*ones(length(x),1);
                            end
                            Mx = min([Mx Mm],[],2);
                            mx = max([mx mm],[],2);
                            for i = 1:m
                                %    bounds(z(:,i),mx,Mx);
                                F = F + (-(Mx-mx)*(1-d(i)) <= z(:,i)-x <= (Mx-mx)*(1-d(i)));
                                F = F + (mx*d(i) <= z(:,i) <= Mx*d(i));
                                [H,K] = double(pwq_struct{1}.Pn(i));
                                [M,m]  = derivebounds(H*x-K);
                                F = F + (H*x-K <= M*(1-d(i)));
                                cost = cost + z(:,i)'*pwq_struct{1}.Ai{i}*z(:,i)+pwq_struct{1}.Bi{i}(:)'*z(:,i)+d(i)*pwq_struct{1}.Ci{i};
                            end
                        end

                        varargout{1} = F;
                        varargout{2} = struct('convexity','none','monotoncity','none','definiteness','none','model','integer');
                        varargout{2}.replacer = cost;
                        varargout{3} = x;

                    case {'convexoverlapping'}

                        n = length(x);
                        cost = 0;
                        r = binvar(length(pwq_struct),1);
                        d = {};

                        % Derive bounds from polytopes
                        [aux,mx,Mx] = bounding_box(pwq_struct{1}.Pfinal);
                        for i = 2:length(pwq_struct)
                            [aux,L,U] = bounding_box(pwq_struct{i}.Pfinal);
                            Mx = max([Mx U],[],2);
                            mx = min([mx L],[],2);
                        end
                        bounds(x,mx,Mx);

                        % Some overlapping function is active
                        % (it will automatically be the smallest one)
                        F = (sum(r) == 1);
                        for j = 1:length(pwq_struct)

                            if length(pwq_struct{j}.Pn) == 1
                                d{j} = r(j);
                            else
                                d{j} = binvar(length(pwq_struct{j}.Pn),1);
                                F = F + (sum(d{j}) == r(j));
                            end
                            m = length(d{j});
                            z = sdpvar(n,m,'full');

                            for i = 1:m
                                bounds(z(:,i),mx,Mx);
                                F = F + (-(Mx-mx)*(1-d{j}(i)) <= z(:,i)-x <= (Mx-mx)*(1-d{j}(i)));
                                F = F + (mx*d{j}(i) <= z(:,i) <= Mx*d{j}(i));
                                [H,K] = double(pwq_struct{j}.Pn(i));
                                [M,m]  = derivebounds(H*x-K);
                                F = F + (H*x-K <= M*(1-d{j}(i)));
                                cost = cost + z(:,i)'*pwq_struct{j}.Ai{i}*z(:,i)+pwq_struct{j}.Bi{i}*z(:,i)+d{j}(i)*pwq_struct{j}.Ci{i};
                            end
                        end

                        varargout{1} = F;
                        varargout{2} = struct('convexity','none','monotoncity','none','definiteness','none','model','integer');
                        varargout{2}.replacer = cost;
                        varargout{3} = x;

                    otherwise

                        varargout{1} = [];
                        varargout{2} = struct('convexity','convex','monotoncity','none','definiteness','none');
                        varargout{3} = [];
                        return
                end

            otherwise
                error('PWA_YALMIP called with CHAR argument?');
        end

    otherwise
        error('Strange type on first argument in PWQ_YALMIP');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function procNoiseParam=processNoiseSuggest(algorithm,maxVal,T,sigmaw2,manDur)
%%PROCESSNOISESUGGEST Use one of a number methods for choosing the scaling
%             parameter for the process noise covariance in a number of
%             continuous-time and discrete-time dynamic models. This only
%             provides the scaling parameter for one dimension of motion.
%             In a multidimensional system, if the measurement noise
%             variance varies between dimensions (such as tracking in range
%             and angle with a range-angle state), then one might use
%             different values in different dimensions. The method only
%             considers scalar measurements in the respective dimensions
%             and is thus only good for a rough approximation.
%
%INPUTS: algorithm A string specifying the algorithm and model used to
%           estimate a good process noise parameter. Possible values are:
%                  'PolyKal-ROT' Use a rule-of-thumb method for determining
%                                the q0 term for the QPolyKal function, for
%                                discrete-time models, and the DPoly
%                                function, for continuous-time models. This
%                                covers the continuous white noise
%                                acceleration (CWNA) model (order=1), the
%                                discretized continuous white noise
%                                acceleration (DCWNA) model (order=1), the
%                                continuous Wiener process acceleration
%                                (CWPA) model  (order=2), and the
%                                discretized continuous Wiener process
%                                acceleration (DCWPA) model (order=2),
%                                among others.
%           'PolyKalDirectDisc-ROT' Use a rule-of-thumb method for
%                                determining the sigmaV2 parameter for the
%                                QPolyKalDirectDisc function. This covers
%                                the discrete white noise acceleration
%                                (DWNA) model (order=1), among others.
%           'PolyKalDirectAlt-ROT' Use a rule-of-thumb method for
%                                determining the sigmaV2 parameter for the
%                                QPolyKalDirectAlt function. This covers
%                                the discrete Wiener process acceleration
%                                (DWPA) model (order=2), among others.
%                 'CWNA-OptMMSE' Use the method of Blair for determining
%                                the asymptotically optimal value of q0 in
%                                terms of MSE for the CWNA and DCWNA
%                                dynamic models for a maneuver having a
%                                fixed maximum acceleration.
%                 'DWNA-OptMMSE' Use the method of Blair for determining
%                                the asymptotically optimal value of
%                                sigmaV2 for the DWNA model for a maneuver
%                                having a fixed maximum acceleration.
%               'CWNA-ConstMeas' Use the method of Blair for determining
%                                the asymptotically optimal value of q0
%                                for the CWNA and DCWNA models such that
%                                the MSE of a Kalman filter under a maximum
%                                acceleration maneuver is not worse than
%                                the measurement accuracy.
%               'DWNA-ConstMeas' Use the method of Blair for determining
%                                the asymptotically optimal value of
%                                sigmaV2 for the DWNA model such that
%                                the MSE of a Kalman filter under a maximum
%                                acceleration maneuver is not worse than
%                                the measurement accuracy.
%    maxVal All of the methods require a maximum bound related to how the
%           target can maneuver. For algorithms CWNA-OptMMSE,
%           DWNA-OptMMSE, CWNA-ConstMeas, and DWNA-ConstMeas, maxVal is
%           the maximum acceleration of the target. For PolyKal-ROT and
%           PolyKalDirectDisc-ROT, and PolyKalDirectAlt-ROT maxVal is the
%           maximum value of a moment one order higher than the maximum
%           order of the dynamic model. For example, if order=1, then
%           maxVal is a maximum acceleration.  If order=2, then maxVal is
%           a maximum jerk.
%         T For all of the algorithms except PolyKalDirectDisc-ROT, this
%           parameter is required and is the typical time between
%           measurements of the target.
%   sigmaw2 For algorithms CWNA-OptMMSE, DWNA-OptMMSE, CWNA-ConstMeas, and
%           DWNA-ConstMeas, this parameter is required and is the variance
%           of the 1D measurement, which should be the lowest order moment
%           of the state. In other words, position.
%    manDur For algorithms CWNA-OptMMSE, DWNA-OptMMSE, CWNA-ConstMeas, and
%           DWNA-ConstMeas, this parameter is required and is the
%           maximum number of samples of duration T that a maximum
%           acceleration maneuver is expected to take. This can take
%           values, 3, 4, 5, 6 and infinity. The default if omitted or an
%           empty matrix is passed is Inf.
%
%OUTPUTS: procNoiseParam The value of q0 or sigmaV2 for the specified
%                        dynamic model, chosen according to the selected
%                        ad-hoc parameter.
%
%The output of this function goes into functions like QPolyKal. It is meant
%that one computes the process noise parameter once (for QPolyKal, this
%would correspond to the q0 parameter), and then one calls the other
%functions with different samples prediction intervals. This function
%should not be called to change the process noise parameter as the sampling
%period changes.
%
%All of the rule-of-thumb methods for choosing the process noise values are
%from Chapter 6.2 and 6.3 of [1] and modified slightly, as described in the
%comments for the implementation below.
%
%The solutions for the ConstMeas and MMSE methods are from [2], which is an
%extension of the work in [3] and [4].
%
%The rule-of thumb parameters might be suitable for other dynamic models.
%For example, when using the Singer model, given with aGaussMarkov and
%DPoly in continuous-time and with FGaussMarkov and QGaussMarkov in
%discrete-time, there does not appear to be a clear way to set the scaling
%value for the process noise in the literature. However, if tau=infinity in
%the model, then in continuous time, it reduces to the CWPA model. Thus,
%methods for choosing the process noise parameter in the CWPA and the DCWPA
%models might be good starting points for setting the process noise
%parameter in the continuous and discrete Singer models.
%
%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.
%[2] W. D. Blair, "Design of nearly constant velocity filters for radar
%    tracking of maneuvering targets," in Proceedings of the IEEE Radar
%    Conference, Atlanta, GA, 7-11 May 2012, pp. 1008-1013.
%[3] W. D. Blair, "Design of nearly constant velocity track filter for
%    brief maneuvers," in Proceedings of the 14th International Conference
%    on Information Fusion, Chicago, IL, 5-8 Jul. 2011.
%[4] W. D. Blair, "Design of nearly constant velocity track filters for
%    tracking maneuvering targets," in Proceedings of the 11th
%    International Conference on Information Fusion, Cologne, Germany, 30
%    Jun. - 3 Jul. 2008.
%
%April 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(manDur))
    manDur=Inf;
end

switch(algorithm)
    case 'PolyKal-ROT'
        %The rule of thumb is generalized from the rules of thumb in
        %Chapters 6.2.2 and 6.2.3 of Bar-Shalom's book, when considering
        %the CWNA, DCWNA, CWPA and DCWPA. For the CWNA and DCWNA models, it
        %is suggested that sqrt(q0*T) be on the order of the changes in
        %velocity over a sampling interval. In other words, it is the
        %average acceleration. It might be more useful to express things in
        %terms of a maximum acceleration, though. In Chapter 6.3.2 for the
        %DWNA model, it is suggested that  0.5*aMax<=sigmaV<=aMax. Using a
        %similar logic, one might try  0.5*aMax<=sqrt(q0*T)<=aMax. Using
        %the same logic for the CWPA and DCWPA models of Chapter 6.2.3, one
        %might choose q0 such that  0.5*jerkMax<=sqrt(q0*T)<=jerkMax, where
        %jerkMax is the maximum possible jerk (derivative of acceleration).
        %This function generalizes the rule  to any order and chooses the
        %point midway in that specified range.

        procNoiseParam=(0.75*maxVal)^2/T;
        return;
    case 'PolyKalDirectDisc-ROT'
        %The rule-of-thumb is generalized from that given in Chapter 6.3.2
        %of Bar-Shalom's book. In the book, it is suggested that
        %0.5*aMax<=sigmaV<=aMax for the DWNA model (order=1). This function
        %generalizes it to any order and chooses the point midway in that
        %specified range.

        procNoiseParam=(0.75*maxVal)^2;
        return;
    case 'PolyKalDirectAlt-ROT'
        %The rule-of-thumb is generalized from that given in Chapter 6.3.3
        %of Bar-Shalom's book. In the book, it is suggested that
        %0.5*deltaAMax<=sigmaV<=deltaAMax for the DWPA model (order=2),
        %where deltaAMax is the maximum change in acceleration over an
        %interval of duration T. This can be related to a maximum jerk as
        %0.5*jerkMax*T<=sigmaV<=jerkMax*T. This function generalizes it to
        %any order and chooses the point midway in that range.

        procNoiseParam=(0.5*maxVal*T)^2;
        return;
    case 'CWNA-OptMMSE'
        optType=0;
        isDiscrete=false;
    case 'DWNA-OptMMSE'
        optType=0;
        isDiscrete=true;
    case 'CWNA-ConstMeas'
        optType=1;
        isDiscrete=false;
    case 'DWNA-ConstMeas'
        optType=1;
        isDiscrete=true;
    otherwise
        error('Unknown algorithm chosen')
end

%Equation 3, the deterministic maneuvering index.
GammaD=maxVal*T^2/sqrt(sigmaw2);

if(0.01>GammaD||GammaD>10)
   warning('The deterministic maneuvering index is outside of the tabulated range of values. The process noise parameter might be inaccurate.')
end
%Blair's papers do not say the base, but if one uses base 10, one gets the
%results in the plots. If one uses base e, the results do not match the
%plots in the paper.
logVal=log10(GammaD);

%optType=0 means MMSE (so use max)
%optType=1 means limit the maximum error to the measurement noise (so use
%          min)
switch(manDur)
    case Inf
        if(optType)
            curRow=2;
        else
            curRow=1;
        end
    case 3
        if(optType)
            curRow=4;
        else
            curRow=3;
        end
    case 4
        if(optType)
            curRow=6;
        else
            curRow=5;
        end
    case 5
    %If one wants a five-sample maneuver, then we will average the values
    %from the length four and the length-six maneuvers.
        if(optType)
            curRow=[6;8];
        else
            curRow=[5;7];
        end
    case 6
        if(optType)
            curRow=8;
        else
            curRow=7;
        end
    otherwise
        error('An untabulated value for the maximum maneuver duration has been given')
end

%Table I in Blair's paper
k1CoeffTable=[1.677, -0.726,  0.230, -0.012, 0.005,  0.000;%k^{max}_{\inf}
              0.872, -0.102, -0.019,  0.010, 0.006,  0.001;%k^{min}_{\inf}
              1.539, -0.189, -0.651,  0.187, 0.284,  0.059;%k^{max}_{3}
              0.707,  0.346, -0.318, -0.097, 0.088,  0.029;%k^{min}_{3}
              1.630, -0.473, -0.403,  0.272, 0.173,  0.023;%k^{max}_{4}
              0.803,  0.197, -0.385, -0.009, 0.138,  0.035;%k^{min}_{4}
              1.671, -0.697,  0.025,  0.250,-0.0324,-0.028;%k^{max}_{6}
              0.869,  0.008, -0.320,  0.088, 0.123,  0.023];%k^{min}_{6}

a0=k1CoeffTable(curRow,1);
a1=k1CoeffTable(curRow,2);
a2=k1CoeffTable(curRow,3);
a3=k1CoeffTable(curRow,4);
a4=k1CoeffTable(curRow,5);
a5=k1CoeffTable(curRow,6);

%Equations 6 and 7. The mean command is for averaging if one wanted a
%length 5 maneuver. Note that equation 6 has a repeated cubed term. That
%appears to be a mistake.
k1=mean(a0+logVal*(a1+logVal*(a2+logVal*(a3+logVal*(a4+a5*logVal)))));

if(isDiscrete)%Discrete-time
    procNoiseParam=(k1*maxVal)^2;
else%Continuous-time
    procNoiseParam=T*(k1*maxVal)^2;
end

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 display(rot,varargin)
% standart output

displayClass(rot,inputname(1),varargin{:});
if length(rot)~=1, disp(['  size: ' size2str(rot)]); end

if length(rot) <= 19 && ~isempty(rot)

  Euler(rot);

elseif ~getMTEXpref('generatingHelpMode')  && ~isempty(rot)

  disp(' ')
  s = setappdata(0,'data2beDisplayed',rot);
  disp(['  show Euler angles'])
  disp(' ')

else

  disp(' ')

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [dists, wh] = canlab_fast_euclidean_distance(X, C, varargin)
% dists = canlab_fast_euclidean_distance(X, C)
%
% X = obs x variables (dimensions, n) for Set 1
% C = obs x variables (dimensions, n) for Set 2
%
% dists
%
% Adapted (largely just borrowed) from:
%%=========================================================
%  Fast Euclidean Distance Calculation
%
%   This script demonstrates the use of matrix multiplication to quickly
%   calculate the Euclidean distance between a large number of vectors.
%
% $Author: ChrisMcCormick $    $Date: 2014/08/22 22:00:00 $    $Revision: 1.0 $
%
% NOTE: (Tor):
% Dot product is SLOWER if n < 6, but gets dramatically faster
% with n = 6 or more. But in some tests it is about equally fast...
% Matrix with fewer observations should be C when using loop (sum sq diffs)
% approach.
%
%%=========================================================
%
% Optional arguments:
%
% 'tracktime'   Track and report time using both matrix and SSD method
% 'target_distance' Distance that counts as "close enough"
% 'squared_distance' Distance returned is squared - saves computation time
%                    Also assumes that target_distance is squared

dotracktime = true;
dosqrt = true;          % save more time by comparing squares to squared dist target
% don't do final sqrt if false, return squared distance

[m, n] = size(X);
k = size(C, 1);

%%  Sum-of-squared-differences approach
%%=========================================================

if dotracktime
    % Measure the time.
    tic();
end

% Create a matrix to hold the distances between each data point and
% each model.
dists = zeros(m, k);

% For each model...
for (i = 1 : k)

    % Subtract model i from all data points.
    diffs = bsxfun(@minus, X, C(i, :));

    % Take the square root of the sum of squared differences.
    dists(:, i) = sqrt(sum(diffs.^2, 2));
end

if dotracktime
    elapsed = toc();
    fprintf('Sum-of-squared-differences took %.3f seconds.\n', elapsed);
    if (exist('OCTAVE_VERSION')) fflush(stdout); end

    dists1 = dists;

end

%%=========================================================
%%  Matrix-multiply approach
%%=========================================================

if dotracktime
    tic();
end

% Calculate the sum of squares for all input vectors and
% for all cluster centers / models.
%
% Matrix dimensions:
%   X  [m  x  n]
%  XX  [m  x  1]
%   C  [k  x  n]
%  CC  [1  x  k]
XX = sum(X.^2, 2);
CC = sum(C.^2, 2)';

% Calculate the dot product between each input vector and
% each cluster center / model.
%
% Matrix dimensions:
%   X  [m  x  n]
%   C  [k  x  n]
%   C' [n  x  k]
%  XC  [m  x  k]
XC = X * C';

% Calculate the Euclidean distance between all input vectors in X
% and all clusters / models in C using the following equation:
%
%   z = sqrt(||x||^2 - 2xc' + ||c||^2)
%
%  Step 1: Subtract the column vector XX from every column of XC.
%  Step 2: Add the row vector CC to every row of XC.
%
% Matrix dimensions:
%     XX  [m  x  1]
%     XC  [m  x  k]
%     CC  [1  x  k]
%  dists  [m  x  k]
%
dists = sqrt(bsxfun(@plus, CC, bsxfun(@minus, XX, 2*XC)));

if dotracktime

    elapsed = toc();
    fprintf('Dot-product took %.3f seconds.\n', elapsed);

    dists2 = dists;

    % Make sure the resulting distances are nearly identical.
    fprintf('Difference in result: %f\n', sum(abs(dists1(:) - dists2(:))));
    if (exist('OCTAVE_VERSION')) fflush(stdout); end

end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function newimg = move3dimage(img,Vy,Vx,Vz,method,offsets,jacobian_modulation)
%
% Calculate the moved image:
%    newimg = move3dimage3(img,Vy,Vx,Vz,method,offsets)
%
% Input:
%	Vy, Vx, Vz	- the motion field
%	method		- interpolation method, default value is 'linear'
%   zoffset		- Offset of Z for the motion field dimension respective to
%				  the img dimension
%
% In version 3, the input method is allowed to be larger than the dimension
% of motion fields. This actually allows better recostruction of the moved
% images because the motion fields could extend larger than the original
% dimension. An additional parameter has been added, the 'offsets'
% parameter.
%
%
if ~exist('method','var') || isempty(method)
	method = 'linear';
end

if ~exist('offsets','var') || isempty(offsets)
	offsets = [0 0 0];
elseif length(offsets) == 1
	offsets = [0 0 offsets];
end

if ~exist('jacobian_modulation','var') || isempty(jacobian_modulation)
	jacobian_modulation = 0;
end

% Computer Jacobian
if jacobian_modulation ~= 0
	jac = compute_jacobian(Vy,Vx,Vz);
end

dimimg = mysize(img);
dimmotion = mysize(Vy);
x0 = single([1:dimmotion(2)])+offsets(2);
y0 = single([1:dimmotion(1)])+offsets(1);
z0 = single([1:dimmotion(3)])+offsets(3);
[xx,yy,zz] = meshgrid(x0,y0,z0);	% xx, yy and zz are the original coordinates of image pixels

Vy = max((yy-Vy),1); clear yy; Vy = min(Vy,dimimg(1));
Vx = max((xx-Vx),1); clear xx; Vx = min(Vx,dimimg(2));
Vz = max((zz-Vz),1); clear zz; Vz = min(Vz,dimimg(3));

if dimimg(3) > 1
	newimg = zeros(dimmotion,'single');
	if offsets(3) == 0 && size(img,3) == size(Vy,3)
		spacing = 20;
		if mod(dimimg(3),spacing) == 1
			spacing = 19;
		end

		N = ceil(dimimg(3)/spacing);
		fprintf('Moving image');
		for k = 1:N
			fprintf('.');
			zmin = (k-1)*spacing+1;
			zmax = min(dimimg(3),k*spacing);
			zmin2 = floor(min(min(min(Vz(:,:,zmin:zmax)))));
			zmax2 = ceil(max(max(max(Vz(:,:,zmin:zmax)))));
			%newimg(:,:,zmin:zmax) = interp3(xx(:,:,zmin2:zmax2),yy(:,:,zmin2:zmax2),zz(:,:,zmin2:zmax2),img(:,:,zmin2:zmax2),Vx(:,:,zmin:zmax),Vy(:,:,zmin:zmax),Vz(:,:,zmin:zmax),method);
			newimg(:,:,zmin:zmax) = interp3(x0,y0,z0(zmin2:zmax2),img(:,:,zmin2:zmax2),Vx(:,:,zmin:zmax),Vy(:,:,zmin:zmax),Vz(:,:,zmin:zmax),method);
		end
		fprintf('\n');
	else
		newimg = interp3(img,Vx,Vy,Vz,method,0);
	end
	%newimg = interp3(img,Vx,Vy,Vz,method,0);
else
	newimg = interp2(img,Vx,Vy,method,0);
end

%newimg = interpn(img,yy-Vy,xx-Vx,zz-Vz,method,0);

if jacobian_modulation ~= 0
	% Apply Jacobin intensity modulation
	jac(isinf(jac)) = 1;
	jac(jac<0) = 1;
	newimg = newimg .* jac;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function moveWaitForDTWhenInterrupted( t_Kuka ,...
     Pos, VEL,joints_indices,max_torque,min_torque, w)
 %% Perform line motion the motion is interruptile through a collision with the environment
 % After a collision the robot waits for a double touch on the end-effector
 % along the Z axes, upon which the robot returns to its motion again.

 % Copyright: Mohammad SAFEEA, 11th-April-2018

 while true
    [ res ] = movePTP_ConditionalTorque_LineEEF( t_Kuka ,...
     Pos, VEL,joints_indices,max_torque,min_torque);
    if(res==1)
        disp('Motion completed successfully');
        break;
    elseif(res==0)
        disp('Motion interrupted');
        detectDoubleTouch(t_Kuka,w);
        pause(0.5);
    elseif(res==-1)
        disp('Error while performing the motion');
        break;
    end
end
end

function x=detectDoubleTouch(t_Kuka,w)
%% Detect double touch at the EEF along the Z axes
buffer=zeros(50,1);
daSum=0;
while true
    [ f ] = getEEF_Force( t_Kuka );
    normForce=(f{3}*f{3})^0.5;
    buffer(1)=[];
    buffer=[buffer;normForce];
    pause(0.04);
    temp=buffer>(5+w);
    for i=2:50
        if (temp(i)-temp(i-1))==1
            daSum=daSum+1;
        end
    end
    if daSum==2
        x=1;
        return;
    end
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% Author: Guosheng Lin (guosheng.lin@gmail.com)
% example code for multi-scale predictions, fusion and accuracy evaluation

function demo_predict_mscale_sunrgbd()

rng('shuffle');

addpath('./my_utils');
dir_matConvNet='../libs/matconvnet/matlab';
run(fullfile(dir_matConvNet, 'vl_setupnn.m'));

run_config=[];
ds_config=[];

run_config.use_gpu=true;
% run_config.use_gpu=false;
run_config.gpu_idx=1;

%-------------------------------------------------------------------------------------------------------------------------
% settings for using trained model:

ds_name_subfix='sunrgbd';

% result dir:
result_name=['result_' datestr(now, 'YYYYmmDDHHMMSS') '_predict_custom_data'];
result_dir=fullfile('../cache_data', ['test_examples_' ds_name_subfix], result_name);

% using a trained model:
run_config.trained_model_path='../model_trained/refinenet_res101_sunrgbd.mat'; % resnet101 based refinenet
% run_config.trained_model_path='../model_trained/refinenet_res152_sunrgbd.mat'; % resnet152 based refinenet

% provide class_info for the trained model:
ds_config.class_info=gen_class_info_sunrgbd();

% control the size of input images to avoid excessively small or large images
run_config.input_img_short_edge_min=450;
run_config.input_img_short_edge_max=800;
%-------------------------------------------------------------------------------------------------------------------------

% specify the folder that contains testing images:
img_data_dir='../datasets/example_imgs_sunrgbd';

% providing ground-truth mask folder for evaluation, or set to empty if not available
% if the folder of ground-truth masks are provided (an example is shown below),
% this code will perform evaluation after multi-scale prediction and fusion.
gt_mask_dir=[];

%-------------------------------------------------------------------------------------------------------------------------
% settings for accuracy evaluation after multi-scale prediction and fusion:
% the following example specifies images for testing and the ground-truth masks for evaluation:

% Notes: for the evaluation purpose, the image folder should only contain testing images, excluding any training images:
% Notes: the ground-truth masks should follow the class label definition in the variable specified above: 'ds_config.class_info'

% uncomment the following lines to enable accuracy evaluation, should replace by your folders:

% img_data_dir='../datasets/sunrgbd/images_testonly';
% gt_mask_dir='../datasets/sunrgbd/gt_testonly';
%-------------------------------------------------------------------------------------------------------------------------

%-------------------------------------------------------------------------------------------------------------------------
% settings for multi-scale prediction, you can consider 3 scales or 5 scales:

prediction_scales=[0.4 0.6 0.8 1 1.2]; % 5 scales
% prediction_scales=[0.6 0.8 1]; % 3 scales
% prediction_scales=[0.8]; % or only use 1 scale
%-------------------------------------------------------------------------------------------------------------------------

ds_config.img_data_dir=img_data_dir;
result_evaluate_param=[];
result_evaluate_param.gt_mask_dir=gt_mask_dir;
ds_config.result_evaluate_param=result_evaluate_param;

run_config.prediction_scales=prediction_scales;
run_config.model_name=result_name;
run_config.root_cache_dir=result_dir;

do_predict_and_evaluate(run_config, ds_config);

end

function do_predict_and_evaluate(run_config, ds_config)

run_config.gen_net_opts_fn=@gen_net_opts_model_type1;

run_config.run_evaonly=true;
ds_config.use_custom_data=true;
ds_config.use_dummy_gt=true;
run_config.use_dummy_gt=ds_config.use_dummy_gt;

ds_config.ds_name='tmp_data';
mkdir_notexist(run_config.root_cache_dir);

diary_dir=run_config.root_cache_dir;
mkdir_notexist(diary_dir);
diary(fullfile(diary_dir, 'output.txt'));
diary on

run_dir_name=fileparts(mfilename('fullpath'));
[~, run_dir_name]=fileparts(run_dir_name);
run_config.run_dir_name=run_dir_name;
run_config.run_file_name=mfilename();

ds_info=gen_dataset_info(ds_config);
my_diary_flush();

if run_config.use_gpu
	gpu_num=gpuDeviceCount;
	if gpu_num>=1
		gpuDevice(run_config.gpu_idx);
    else
        error('no gpu found!');
	end
end

train_opts=run_config.gen_net_opts_fn(run_config, ds_info.class_info);

imdb=my_gen_imdb(train_opts, ds_info);
data_norm_info=[];
data_norm_info.image_mean=128;
imdb.ref.data_norm_info=data_norm_info;

[net_config, net_exp_info]=prepare_running_model(train_opts);

prediction_scales=run_config.prediction_scales;
scale_num=length(prediction_scales);
predict_result_dirs=cell(scale_num, 1);

for s_idx=1:scale_num

	one_scale=prediction_scales(s_idx);
	one_result_dir=fullfile(run_config.root_cache_dir, sprintf('predict_result_%d', s_idx));
	predict_result_dirs{s_idx}=one_result_dir;

	fprintf('\n\n--------------------------------------------------\n\n');
	fprintf('conduct prediction using image scale: %1.2f  (current scale / total scales: %d/%d) \n', one_scale, s_idx, scale_num);

	train_opts.root_cache_dir=one_result_dir;
	train_opts.input_img_scale=one_scale;
	train_opts.eva_param=update_eva_param_mscale(train_opts.eva_param, train_opts);

	my_net_tool(train_opts, imdb, net_config, net_exp_info);

	fprintf('\n\n--------------------------------------------------\n\n');
	disp('results are saved in:');
	disp(train_opts.root_cache_dir);

	my_diary_flush();

end

if length(predict_result_dirs)>1

    fprintf('\n\n--------------------------------------------------\n\n');
    disp('fusing multiscale predictions');

    fuse_param=[];
    fuse_param.predict_result_dirs=predict_result_dirs;
    fuse_param.fuse_result_dir=fullfile(run_config.root_cache_dir, 'predict_result_final_fused');
    fuse_param.cache_fused_score_map=true;
    fuse_multiscale_results(fuse_param, ds_config.class_info);

    fprintf('\n\n--------------------------------------------------\n\n');
    disp('final fused results are saved in:');
    disp(fuse_param.fuse_result_dir);

    my_diary_flush();

    final_prediction_dir=fuse_param.fuse_result_dir;

else

    final_prediction_dir=predict_result_dirs{1};
end

% evaluate the final fused result, if ground-truth is provided
result_evaluate_param=ds_config.result_evaluate_param;
if ~isempty(result_evaluate_param) && ~isempty(result_evaluate_param.gt_mask_dir)
	result_evaluate_param.predict_result_dir=fullfile(final_prediction_dir, 'predict_result_mask');
	fprintf('\n\n--------------------------------------------------\n\n');
	disp('perform evaluation of the predicted masks');
	result_info=evaluate_predict_results(result_evaluate_param, ds_config.class_info);
    disp_evaluate_result(result_info, ds_config.class_info);

    eva_result_cached_file=fullfile(final_prediction_dir, 'eva_result_info.mat');
    fprintf('saving evaluation result to: %s\n', eva_result_cached_file);
    save(eva_result_cached_file, 'result_info');
end

my_diary_flush();
diary off

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
%% clear the workspace and select data
clear; clc; close all;

%% choose multiple datasets or just one
neuron = Sources2D();
nams = {'./data_2p.tif'};          % you can put all file names into a cell array; when it's empty, manually select files
nams = neuron.select_multiple_files(nams);  %if nam is [], then select data interactively

%% parameters
% -------------------------    COMPUTATION    -------------------------  %
pars_envs = struct('memory_size_to_use', 8, ...   % GB, memory space you allow to use in MATLAB
    'memory_size_per_patch', 0.3, ...   % GB, space for loading data within one patch
    'patch_dims', [30, 40],...  %GB, patch size
    'batch_frames', 1000);           % number of frames per batch

% -------------------------      SPATIAL      -------------------------  %
gSig = 0.5;           % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
gSiz = 10;          % pixel, neuron diameter
ssub = 1;           % spatial downsampling factor
with_dendrites = false;   % with dendrites or not
if with_dendrites
    % determine the search locations by dilating the current neuron shapes
    updateA_search_method = 'dilate';  %#ok<UNRCH>
    updateA_bSiz = 20;
    updateA_dist = neuron.options.dist;
else
    % determine the search locations by selecting a round area
    updateA_search_method = 'ellipse'; %#ok<UNRCH>
    updateA_dist = 5;
    updateA_bSiz = neuron.options.dist;
end
spatial_constraints = struct('connected', true, 'circular', false);  % you can include following constraints: 'circular'

% -------------------------      TEMPORAL     -------------------------  %
Fs = 5;             % frame rate
tsub = 1;           % temporal downsampling factor
deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
    'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
    'smin', -3, ...         % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
    'optimize_pars', true, ...  % optimize AR coefficients
    'optimize_b', true);% optimize the baseline);
nk = 1;             % detrending the slow fluctuation. usually 1 is fine (no detrending)
                    % when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'local_min';  % compute the local minimum as an estimation of trend.

% -------------------------     BACKGROUND    -------------------------  %
bg_model = 'svd';  % model of the background {'ring', 'svd'(default), 'nmf'}
nb = 1;             % number of background sources for each patch (only be used in SVD and NMF model)
bg_neuron_factor = 0.5;
ring_radius = round(bg_neuron_factor * gSiz);  % when the ring model used, it is the radius of the ring used in the background model.
                    %otherwise, it's just the width of the overlapping area

% -------------------------      MERGING      -------------------------  %
show_merge = false;  % if true, manually verify the merging step
merge_thr = 0.7;     % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
method_dist = 'max';   % method for computing neuron distances {'mean', 'max'}
dmin = 10;       % minimum distances between two neurons. it is used together with merge_thr
dmin_only = 2;  % merge neurons if their distances are smaller than dmin_only.

% -------------------------  INITIALIZATION   -------------------------  %
K = [];             % maximum number of neurons per patch. when K=[], take as many as possible
min_corr = 0.7;     % minimum local correlation for a seeding pixel
min_pnr =10;       % minimum peak-to-noise ratio for a seeding pixel
min_pixel = 2^2;      % minimum number of nonzero pixels for each neuron
bd = 0;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
frame_range = [];   % when [], uses all frames
save_initialization = false;    % save the initialization procedure as a video.
use_parallel = true;    % use parallel computation for parallel computing
show_init = true;   % show initialization results
choose_params = false; % manually choose parameters
center_psf = false;  % set the value as true when the background fluctuation is large (usually 1p data)
                    % set the value as false when the background fluctuation is small (2p)

% ----------------------   MANUAL INTERVENTION  --------------------  %
with_manual_intervention = false;

% -------------------------  FINAL RESULTS   -------------------------  %
save_demixed = true;    % save the demixed file or not
kt = 1;                 % frame intervals

% -------------------------    UPDATE ALL    -------------------------  %
neuron.updateParams('gSig', gSig, ...       % -------- spatial --------
    'gSiz', gSiz, ...
    'ring_radius', ring_radius, ...
    'ssub', ssub, ...
    'search_method', updateA_search_method, ...
    'bSiz', updateA_bSiz, ...
    'dist', updateA_bSiz, ...
    'spatial_constraints', spatial_constraints, ...
    'tsub', tsub, ...                       % -------- temporal --------
    'deconv_options', deconv_options, ...    '
    'nk', nk, ...
    'detrend_method', detrend_method, ...
    'background_model', bg_model, ...       % -------- background --------
    'nb', nb, ...
    'ring_radius', ring_radius, ...
    'merge_thr', merge_thr, ...             % -------- merging ---------
    'dmin', dmin, ...
    'method_dist', method_dist, ...
    'min_corr', min_corr, ...               % ----- initialization -----
    'min_pnr', min_pnr, ...
    'min_pixel', min_pixel, ...
    'bd', bd, ...
    'center_psf', center_psf);
neuron.Fs = Fs;

%% distribute data and be ready to run source extraction
neuron.getReady_batch(pars_envs);

%% initialize neurons in batch mode
neuron.initComponents_batch(K, save_initialization, use_parallel);

%% udpate spatial components for all batches
neuron.update_spatial_batch(use_parallel);

%% udpate temporal components for all bataches
neuron.update_temporal_batch(use_parallel);

%% update background
neuron.update_background_batch(use_parallel);

%% delete neurons

%% merge neurons

%% save workspace
neuron.save_workspace_batch();

%% get the correlation image and PNR image for all neurons
neuron.correlation_pnr_batch();

%% concatenate temporal components
neuron.concatenate_temporal_batch();
neuron.viewNeurons([], neuron.C_raw);

author: zhoupc


text:
%qrscorrect() -  Correct for false positive and negative QRS peaks
%   and align peaks for maximum correlation.
%
% Usage:
%   >>Peaks=qrscorrect(Peaks,Ecg,fs)
%
% Inputs:
%   Peaks: Indeces of located QRS peaks
%   Ecg: vector of ecg data
%   sampling frequency.
%
% Ouput:
%   Peaks: Corrected peaks
%
% Note: progress bar designed to work as part of fmrib_qrsdetect.
%
% Author: Rami Niazy, FMRIB Centre, University of Oxford.
%
%
%
% Copyright (c) 2004 Rami Niazy, FMRIB Centre, University of Oxford.
%

% Copyright (C) 2004 Rami K. Niazy, FMRIB Centre, University of Oxford
%                    rami@fmrib.ox.ac.uk
%
% 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 2 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
%

%   Dec 23, 2004
%   use binary prcorr2
%   instead of corrcoef

%   Dec 17, 2004
%   use corrcoef instead of
%   corr2

%   Dec 15, 2004
%   Median of all peaks used
%   for correction instead of
%   section median.

%   Nov 4, 2004
%   Changed order of correction
%   Added additional alignment

%   Nov 3, 2004
%   Added check and
%   adjustment of edge peaks clearance

%   Oct 26, 2004
%   Fixed bug for non-int fs case

function Peaks=qrscorrect(Peaks,Ecg,fs)

% init
%------
P=zeros(1,length(Ecg));
P(Peaks)=1;
dP_std=std(diff(Peaks));
dP_med=median(diff(Peaks));
pad=round(5*fs);
secL=round(10*fs);
sections=floor(length(P)/secL)+1;

if length(P((sections-1)*secL+1:end))<pad
    sections=sections-1;
end

barth=2;
barth_step=barth;
if isempty(javachk('desktop'))
    clearcom='clc';
else
    clearcom='home';
end

% Correct for FP
%---------------
for s=1:sections

    %wait bar
    if s==1
        barth=5;
        barth_step=barth;
        Flag25=0;
        Flag50=0;
        Flag75=0;
        fprintf('\nStage 2 of 5: Correcting for False Positive detection.\n');

    end
    %-----------------calc------------------------------------------

    f_flag=0;
    if s==1
        sec=P(1:secL+pad);
    elseif s==sections
        sec=P((s-1)*secL+1-pad:end);
    else
        sec=P((s-1)*secL+1-pad:s*secL+pad);
    end

    sec_P=find(sec==1);
    sec_dP=diff(sec_P);
%     dP_med=median(sec_dP);

    false_p=1;
    while ~isempty(false_p)
        false_p=find(sec_dP<(dP_med-3*dP_std));%-3*dP_std
        if ~isempty(false_p)
            f_flag=1;
            sec(sec_P(false_p(1)+1))=0;
            sec_P=find(sec==1);
            sec_dP=diff(sec_P);
            false_p=find(sec_dP<(dP_med-3*dP_std));
        end
    end

    if f_flag==1
        if s==1
            P(1:secL+pad)=sec;
        elseif s==sections
            P((s-1)*secL+1-pad:end)=sec;
        else
            P((s-1)*secL+1-pad:s*secL+pad)=sec;
        end
    end

    %---------------------------------------------------------------
    %update wait bar
    percentdone=floor(s*100/sections);
    if floor(percentdone)>=barth
        if percentdone>=25 & Flag25==0
            fprintf('25%% ')
            Flag25=1;
        elseif percentdone>=50 & Flag50==0
            fprintf('50%% ')
            Flag50=1;
        elseif percentdone>=75 & Flag75==0
            fprintf('75%% ')
            Flag75=1;
        elseif percentdone==100
            fprintf('100%%\n')
        else
            fprintf('.')
        end

        while barth<=percentdone
            barth=barth+barth_step;
        end
        if barth>100
            barth=100;
        end
    end
end

%Check Edge Peaks
%----------------
Peaks=find(P==1);
pc1=2;
pc2=1;

while (pc1>1 | pc2>0)
	dP=diff(Peaks);
	hwinl=round(mean(dP)/2);
	searchw=round(0.33*mean(dP));

	pc1=1;
	if length(Peaks)>5
        for p=1:5
            try
                tmp=Ecg(Peaks(p)-hwinl-searchw:Peaks(p)+hwinl);
                break;
            catch
                pc1=pc1+1;
            end
        end
	end

	if pc1>1
        Peaks(1:pc1-1)=[];
	end

	pc2=0;
	if length(Peaks)>5
        for p=length(Peaks):-1:length(Peaks-4)
            try
                tmp=Ecg(Peaks(p)-hwinl:Peaks(p)+hwinl+searchw);
                break;
            catch
                pc2=pc2+1;
            end
        end
	end

	if pc2>0
        Peaks(end-pc2:end)=[];
	end
end

%Finding corrected peaks
%-----------------------
searchwL=length(searchw);
peaksL=length(Peaks);
qrs=zeros(hwinl*2+1,peaksL);
qrs2=zeros(hwinl*2+1,peaksL);
barth=2;

% Find average Heart Beat for use in alignment
%---------------------------------------------
for p=1:peaksL
    qrs(:,p)=Ecg(Peaks(p)-hwinl:Peaks(p)+hwinl);
end

m_qrs=mean(qrs,2);
% figure;plot(qrs)

% Align Peaks (1)
%----------------
for p=1:peaksL

    %wait bar
    if p==1
        barth=5;
        barth_step=barth;
        Flag25=0;
        Flag50=0;
        Flag75=0;
        fprintf('\nStage 3 of 5: Aligning QRS Peaks  (1)\n');
    end

    %-----------------calc------------------------------------------

    ppn=1;
    for B=Peaks(p)-searchw:Peaks(p)+searchw
        C(ppn)=prcorr2(Ecg(B-hwinl:B+hwinl),m_qrs);;
        ppn=ppn+1;
	end
	[CV,CP]=max(C);
	Beta=CP-(searchw+1);
	Peaks(p)=Peaks(p)+Beta;

    %---------------------------------------------------------------
    %update wait bar
    percentdone=floor(p*100/(peaksL));
    if floor(percentdone)>=barth
        if percentdone>=25 & Flag25==0
            fprintf('25%% ')
            Flag25=1;
        elseif percentdone>=50 & Flag50==0
            fprintf('50%% ')
            Flag50=1;
        elseif percentdone>=75 & Flag75==0
            fprintf('75%% ')
            Flag75=1;
        elseif percentdone==100
            fprintf('100%%\n')
        else
            fprintf('.')
        end

        while barth<=percentdone
            barth=barth+barth_step;
        end
        if barth>100
            barth=100;
        end
    end
end
barth=2;
P=zeros(1,length(Ecg));
P(Peaks)=1;

% Correct for FN
%---------------
for s=1:sections

    %wait bar
    if s==1
        barth=5;
        barth_step=barth;
        Flag25=0;
        Flag50=0;
        Flag75=0;
        fprintf('\nStage 4 of 5: Correcting for False Negative detection.\n');
    end

    %-----------------calc------------------------------------------

    f_flag=0;
    if s==1
        sec=P(1:secL+pad);
    elseif s==sections
        sec=P((s-1)*secL+1-pad:end);
    else
        sec=P((s-1)*secL+1-pad:s*secL+pad);
    end

    sec_P=find(sec==1);
    sec_dP=diff(sec_P);
    dP_med=median(sec_dP);

    false_n=1;
    while ~isempty(false_n)
        false_n=find(sec_dP>(1.5*dP_med));
        if ~isempty(false_n)
            f_flag=1;
            sec(sec_P(false_n(1))+round(dP_med))=1;%sec_dP(false_n(1)-1)
            sec_P=find(sec==1);
            sec_dP=diff(sec_P);
            false_n=find(sec_dP>(1.5*dP_med));
        end
    end

    if f_flag==1
        if s==1
            P(1:secL+pad)=sec;
        elseif s==sections
            P((s-1)*secL+1-pad:end)=sec;
        else
            P((s-1)*secL+1-pad:s*secL+pad)=sec;
        end
    end

    %---------------------------------------------------------------
    %update wait bar
    percentdone=floor(s*100/sections);
    if floor(percentdone)>=barth
        if percentdone>=25 & Flag25==0
            fprintf('25%% ')
            Flag25=1;
        elseif percentdone>=50 & Flag50==0
            fprintf('50%% ')
            Flag50=1;
        elseif percentdone>=75 & Flag75==0
            fprintf('75%% ')
            Flag75=1;
        elseif percentdone==100
            fprintf('100%%\n')
        else
            fprintf('.')
        end

        while barth<=percentdone
            barth=barth+barth_step;
        end
        if barth>100
            barth=100;
        end
    end
end

barth=2;
Peaks=find(P==1);
peaksL=length(Peaks);
barth=2;
% Align Peaks (2)
%----------------
for p=2:peaksL-1

    %wait bar
    if p==2
        barth=5;
        barth_step=barth;
        Flag25=0;
        Flag50=0;
        Flag75=0;
         fprintf('\nStage 5 of 5: Aligning QRS Peaks  (2)\n');
    end

    %-----------------calc------------------------------------------

    ppn=1;
    for B=Peaks(p)-searchw:Peaks(p)+searchw
        C(ppn)=prcorr2(Ecg(B-hwinl:B+hwinl),m_qrs);
        ppn=ppn+1;
	end
	[CV,CP]=max(C);
	Beta=CP-(searchw+1);
	Peaks(p)=Peaks(p)+Beta;

    %---------------------------------------------------------------
    %update wait bar
    percentdone=floor((p-1)*100/(peaksL-2));
    if floor(percentdone)>=barth
        if percentdone>=25 & Flag25==0
            fprintf('25%% ')
            Flag25=1;
        elseif percentdone>=50 & Flag50==0
            fprintf('50%% ')
            Flag50=1;
        elseif percentdone>=75 & Flag75==0
            fprintf('75%% ')
            Flag75=1;
        elseif percentdone==100
            fprintf('100%%\n')
        else
            fprintf('.')
        end

        while barth<=percentdone
            barth=barth+barth_step;
        end
        if barth>100
            barth=100;
        end
    end
end

% for p=2:peaksL-1
%     qrs2(:,p-1)=Ecg(Peaks(p)-hwinl:Peaks(p)+hwinl);
% end
% figure;plot(qrs2);
return;

author: PatternRecognition

repo: OpenBMI


text:
function button = rxMakeButton(name,pos,direc);
%  button = rxMakeButton(name,pos,[direc]);
%
% Make a radio button for mrRx, with
% the specified name and position.
%
% If a cell of names is specified,
% will make a vertical set of radio
% buttons, one for each entry, filling
% up the box specified by pos.
%
% direc: direction flag--only of relevance
% for multiple buttons. If 0 [default],
% stack buttons vertically; otherwise, place
% horizontally.
%
% Returns a vector of handles to the buttons.
%
% ras 02/05.
if ieNotDefined('direc')
    direc = 0;
end

if iscell(name)
    % recursively make buttons
    nButtons = length(name);
    for i = 1:nButtons
        subpos = pos;
        if direc==0
            % stack vertically
            subpos(2) = pos(2) + (i-1)*pos(4)/nButtons;
            subpos(4) = pos(4)/nButtons;
        else
            % place horizontally
            subpos(1) = pos(1) + (i-1)*pos(3)/nButtons;
            subpos(3) = pos(3)/nButtons;
        end
        button(i) = rxMakeButton(name{i},subpos,direc);

        % new callback:
        cbstr = sprintf('selectButton(get(gcbo,''UserData''),%i);',i);
        cbstr = sprintf('%s \n rxRefresh;',cbstr);
        set(button(i),'Callback',cbstr);
    end
    set(button,'UserData',button,'Style','radiobutton');
    return
end

button = uicontrol('Style','checkbox',...
                   'Units','Normalized',...
                   'Position',pos,...
                   'String',name,...
                   'Min',0,'Max',1,...
                   'Value',0,...
                   'BackgroundColor',get(gcf,'Color'),...
                   'ForegroundColor',[.3 .3 .3],...
                   'Callback','rxRefresh;');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function calcAllAertsFeatures_batchHN_temp(pathData,pathText,namePT,nameCT,nameROI,outcomes,featType,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% function calcAllSeparateTextures_batchHN(pathData,pathText,namePT,nameCT,nameROI,outcomes,featType,scale_mat,Ng_mat,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes SEPARATE texture features for all patients, for
% all different combinations of the following texture extraction parameters:
% - Scale: Resolution at which the ROI is isotropically resampled.
% - Ng: Number of gray-levels in the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to test all possible combinations. This function is
% used for SEPARATE scans specifically. See Ref. [1,2] and 'prepareVolume.m'
% for more details.
%
% Texture features are computed for all head and neck (HN) DICOM imaging
% data downloaded from The Cancer Imaging Archive (TCIA) website at:
% <http://dx.doi.org/xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx>, and first organized
% in a 'DATA' directory using the function readAllDICOM_HN.m.  Results are
% then saved in a folder 'TEXTURES' in the HN WORKSPACE.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathData: Full path to the HN sData files directory.
%              --> Ex: '/myProject/WORKSPACE/DATA'
% 2. pathText: Full path to the HN non texture features directory.
%              --> Ex: '/myProject/WORKSPACE/FEATURES/TEXTURES'
% 3. namePT: Cell of strings of all PET sData files to read
%            --> Ex: {'HGJ_001_PT.PTscan.mat';'HGJ_022_PT.PTscan.mat'}
% 4. namePT: Cell of strings of all CT sData files to read
%            --> Ex: {'HGJ_001_CT.CTscan.mat';'HGJ_022_CT.CTscan.mat'}
% 5. nameROI: Cell of strings specifying the ROI names to analyze for the
%             patients defined by "namePT" and "nameCT"
%             --> Ex: {'GTV';'GTV-P'}
% 6. outcomes: Structure specifying the status (1 or 0) for different
%              outcomes in HN cancer. Contains: outcomes.Failure,
%              outcomes.Locoregional, outcomes.Distant. See ref.[1] for
%              more details.
% 7. featType: Either 'GTVp' for primary GTV, or 'GTVtot' for primaty GTV +
%              nodal GTVs
%              --> Ex: 'GTVp'
% 8. scale_mat: Array vector specifying the different 'Scale' values to test.
%               --> Ex: [1,2,3,4,5]
% 9. Ng_mat: Array vector specifying the different 'Ng' values to test.
%            --> Ex: [8,16,32,64]
% 10. nBatch: Number of parallel batch.
%             --> Ex: 8
% 11.  matlabPATH: Full path to the MATLAB excutable on the system.
%      --> Ex: 'matlab' (symbolic link)
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
ind = strfind(namePT{1},'_'); cohortID = namePT{1}(1:ind-1);
cd(pathText), mkdir(['batchLog_',cohortID,'_',featType,'_AertsFeatures']), cd(['batchLog_',cohortID,'_',featType,'_AertsFeatures']), pathBatch = pwd;
time = 60; % Number of seconds to wait before checking if parallel computations are done
nameOutcomes = fieldnames(outcomes); nOutcomes = numel(nameOutcomes);
scans = {'PT','CT'}; nScans = numel(scans);

% PRODUCE BATCH COMPUTATIONS
nPatient = numel(namePT); valid = ones(nPatient,1);
for i = 1:nPatient
    if isempty(nameROI{i})
        valid(i) = 0;
    end
end
indValid = find(valid); nPatient = numel(indValid);
if nPatient < nBatch
    nBatch = nPatient;
end
[patients] = batchPatients(nPatient,nBatch);
save('workspace','pathData','pathText','namePT','nameCT','nameROI','patients','indValid','featType'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    fprintf(fid,['calcAllAertsFeatures_HN(pathData,pathText,namePT(indValid(patients{',num2str(i),'})),nameCT(indValid(patients{',num2str(i),'})),nameROI(indValid(patients{',num2str(i),'})),featType)\n']);
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

% % GROUPING RESULTS FROM ALL BATCH
% nPatient = numel(namePT);
% names = {namePT,nameCT};
% for scan = 1:nScans
%     cd(pathText)
%     sign = struct;
%     signTemp = zeros(nPatient,4); % 4-feature signature
%     tempText = cell(1,nPatient); % Cell used to load patient textures only once
%     for p = 1:nPatient
%         ind = strfind(names{scan}{p},'.'); namePatientScan = names{scan}{p}(1:ind(1)-1);
%         if exist([namePatientScan,'_',featType,'_sign.mat'],'file')
%             load([namePatientScan,'_',featType,'_sign']) % Variable 'signature' is now in MATLAB workspace
%         else
%             load([namePatientScan,'_','GTVp','_sign']) % Variable 'signature' is now in MATLAB workspace
%         end
%         tempText{p} = signature;
%     end
%     for p = 1:nPatient
%         signTemp(p,:) = tempText{p};
%     end
%     sign.Data = signTemp;
%     for o = 1:nOutcomes
%         [sign.Spearman.(nameOutcomes{o}).rs,sign.Spearman.(nameOutcomes{o}).p] = corr(sign.Data,outcomes.(nameOutcomes{o}),'type','Spearman','rows','pairwise');
%     end
%     cd .., save(['AertsSign_',cohortID,'_',scans{scan},'_',featType],'sign')
% end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%% Copyright (C) 2014, 2016, 2018-2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym children (@var{f})
%% Return "children" (terms, lhs/rhs, etc) of symbolic expression.
%%
%% For a scalar expression, return a row vector of sym expressions:
%% @example
%% @group
%% syms x y
%% f = 2*x*y + sin(x);
%% C = children(f)
%%   @result{} C = (sym) [2⋅x⋅y  sin(x)]  (1×2 matrix)
%%
%% children(C(1))
%%   @result{} ans = (sym) [2  x  y]  (1×3 matrix)
%% children(C(2))
%%   @result{} ans = (sym) x
%% @end group
%% @end example
%%
%% A symbol/number/boolean has itself as children:
%% @example
%% @group
%% children(x)
%%   @result{} ans = (sym) x
%% @end group
%% @end example
%%
%% For matrices/vectors, return a cell array where each entry is
%% a row vector.  The cell array is the same shape as the input.
%% @example
%% @group
%% A = [x*y 2; 3 x]
%%   @result{} A = (sym 2×2 matrix)
%%       ⎡x⋅y  2⎤
%%       ⎢      ⎥
%%       ⎣ 3   x⎦
%%
%% C = children (A)
%%   @result{} C = @{ ... @}
%% @end group
%%
%% @group
%% class (C), size (C)
%%   @result{} ans = cell
%%   @result{} ans =
%%        2   2
%%
%% C@{:@}
%%   @result{} ans = (sym) [x  y]  (1×2 matrix)
%%   @result{} ans = (sym) 3
%%   @result{} ans = (sym) 2
%%   @result{} ans = (sym) x
%% @end group
%% @end example
%%
%%
%% For sets, @code{children} can be used to extract
%% an matrix (array) containing the set elements, @pxref{finiteset}.
%% This is useful for accessing the elements of a set.
%%
%% @seealso{@@sym/lhs, @@sym/rhs, @@sym/eq, @@sym/lt, finiteset}
%% @end defmethod

function r = children(f)

  cmd = {
    'f, = _ins'
    'f = sympify(f)'  % mutable -> immutable
    'def scalarfcn(a):'
    '    if not hasattr(a, "args") or len(a.args) == 0:'
    '        return sympy.Matrix([a])'  % children(x) is [x]
    '    return sympy.Matrix([a.args])'
    '# note, not for MatrixExpr'
    'if isinstance(f, sp.MatrixBase):'
    '    r = [scalarfcn(a) for a in f.T]'  % note transpose
    'else:'
    '    r = scalarfcn(f)'
    'return r,' };

  r = pycall_sympy__ (cmd, f);

  if (~isscalar(f))
    r = reshape(r, size(f));
  end

end

%!test
%! % basics, sum
%! syms x y
%! f = 2*x + x*x + sin(y);
%! assert (isempty (setxor (children(f), [2*x x*x sin(y)])))

%!test
%! % basics, product
%! syms x y
%! f = 2*x*sin(y);
%! assert (isempty (setxor (children(f), [2 x sin(y)])))

%!test
%! % basics, product and powers
%! syms x y
%! f = 2*x^2*y^3;
%! assert (isempty (setxor (children(f), [2 x^2 y^3])))

%!test
%! % eqn, ineq
%! syms x y
%! lhs = 2*x^2; rhs = y^3 + 7;
%! assert (isequal (children(lhs == rhs), [lhs rhs]))
%! assert (isequal (children(lhs < rhs),  [lhs rhs]))
%! assert (isequal (children(lhs >= rhs), [lhs rhs]))

%!test
%! % matrix
%! syms x y
%! f = [4 + y  1 + x;  2 + x  3 + x];
%! c = children(f);
%! ec = {[4 y], [1 x]; [2 x], [3 x]};
%! assert (isequal (size(c), size(ec)))
%! for i=1:length(c)
%!   assert (isempty (setxor (c{i}, ec{i})))
%! end

%!test
%! % matrix, sum/prod
%! syms x y
%! f = [x + y; x*sin(y); sin(x)];
%! ec = {[x y]; [x sin(y)]; [x]};
%! c = children(f);
%! assert (isequal (size(c), size(ec)))
%! for i=1:length(c)
%!   assert (isempty (setxor (c{i}, ec{i})))
%! end

%!test
%! % scalar symbol
%! syms x
%! assert (isequal (children(x), x))

%!test
%! % scalar number
%! x = sym(6);
%! assert (isequal (children(x), x))

%!test
%! % symbolic size matrix
%! syms n m integer
%! A = sym('a', [n m]);
%! C = children (A);
%! assert (isequal (C(2), n))
%! assert (isequal (C(3), m))

%!xtest
%! % symbolic size matrix, fails on newer SymPy Issue #1089
%! syms n m integer
%! A = sym('a', [n m]);
%! assert (isequal (children (A), [sym('a') n m]))

repo: octsympy


text:
function varargout=msentropy(varargin)
%
% [y,scale,info]=msentropy(x,dn,dm,dr,N,N0,minM,maxM,maxScale,minR,maxR)
%
%    Wrapper to the Multiscale Entropy C code written by Madalena Costa (mcosta@fas.harvard.edu):
%         http://physionet.org/physiotools/mse/mse-1.htm
%
% Calculates the multi scale entropy of a signal 'x'. A tutorial on Mulsticale
% entropy is available at:
% http://www.physionet.org/physiotools/mse/tutorial/
%
%
% Please cite these publications when referencing this material:
%     Costa M., Goldberger A.L., Peng C.-K. Multiscale entropy analysis of biological signals. Phys Rev E 2005;71:021906.
%     Costa M., Goldberger A.L., Peng C.-K. Multiscale entropy analysis of physiologic time series. Phys Rev Lett 2002; 89:062102.
%
% Also include the standard citation for PhysioNet:
%     Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG,
%     Mietus JE, Moody GB, Peng C-K, Stanley HE. PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. Circulation 101(23):e215-e220 [Circulation Electronic Pages; http://circ.ahajournals.org/cgi/content/full/101/23/e215]; 2000 (June 13)
%
% Readers of may also wish to read:
%     Costa M, Peng C-K, Goldberger AL, Hausdorff JM. Multiscale entropy analysis of human gait dynamics. Physica A 2003;330:53-60.
%
% Required Parameters:
%
% x
%       Nx1 vector of doubles in which to caculate the multiscale entropy.
%
% Optional Parameters are:
% dn
%       1x1 double. Sets the scale increment to dn (1-40; default: 1).
% dm
%       1x1 double. Sets the m increment to dm (1-10; default: 1).
% dr
%       1x1 double. Sets the scale increment to dr (>0; default: 0.05).
% N
%       1x1 integer. Stop the analysis with row N.
%       By default, analysis ends at row 39999, or at the end of the data set if there are fewer rows.
% N0
%       1x1 integer. Begin the analysis with row N0.
%       By default, analysis begins with row 1.
% minM
%       1x1 integer betwee 1-10. Set the minimum pattern length for SampEn to minN (1-10; default: 2).
% maxM
%        1x1 integer betwee 1-10. Set the maximum m to maxM (1-10; default: 2).
% maxScale
%        1x1 integer betwee 1-40. Set the maximum scale for coarse-graining to maxScale (1-40; default: 20).
% minR
%        1x1 double >0. Set the minimum similarity criterion for SampEn to minR (>0; default: 0.15).
% maxR
%        1x1 double > 0. Set the maximum m to maxR (>0; default: 0.15).
%
%
% Outputs:
% y
%       A 1xM vector of doubles corresponding to estimated sample entropies at each scale.
% scale
%       A 1xM vector of integers specifying the scales in which 'y' was
%       estimated.
%
% info
%       An optional 3x1 cell array of strings providing loggin and verbose information from
%       the calculation.
%
% Wrapper written by Ikaro Silva, 2013
% Last Modified: March 20, 2014
% Version 0.0.1
%
% Since 0.9.5
%
% %Example
% N=30000;
% noise=randn(N,1);
% maxScale=10;
%[entropyNoise,scale1]=msentropy(noise,[],[],[],[],[],[],[],maxScale);
% %Simulate determistic system with noise-like 2nd order statistics
% nlinear=zeros(N,1);nlinear(1)=0.2;u=4;
% for n=2:N;nlinear(n)=u*nlinear(n-1)*(1-nlinear(n-1));end
%[entropyDeterm,scale2]=msentropy(nlinear,[],[],[],[],[],[],[],maxScale);
%subplot(2,1,1);
%plot(noise(1:1000));hold on;grid on;plot(nlinear(1:1000),'r');legend('Stochastic','Deterministic')
%subplot(2,1,2);
%plot(scale1,entropyNoise);hold on;grid on;plot(scale2,entropyDeterm,'r');legend('Stochastic','Deterministic')
%
%
% See also SURROGATE, DFA, WFDBDESC, PHYSIONETDB, RDANN, ANN2RR, MAPRECORD

%endOfHelp
persistent javaWfdbExec config
if(isempty(javaWfdbExec))
    [javaWfdbExec,config]=getWfdbClass('mse');
end

%Set default pararamter values
inputs={'x','dn','dm','dr','N','N0','minM','maxM','maxScale','minR','maxR'};
outputs={'y','scale','info'};
dn=[];
dm=[];
dr=[];
N=[];
N0=[];
minM=[];
maxM=[];
maxScale=[];
minR=[];
maxR=[];
wfdb_argument={};
info=[];
scale=[];
y=[];
x=[];
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};'])
    end
end
if(~isempty(dn))
    wfdb_argument{end+1}='-a';
    wfdb_argument{end+1}=[num2str(dn)];
end
if(~isempty(dm))
    wfdb_argument{end+1}='-b';
    wfdb_argument{end+1}=[num2str(dm)];
end
if(~isempty(dr))
    wfdb_argument{end+1}='-c';
    wfdb_argument{end+1}=[num2str(dr)];
end
if(~isempty(N0))
    wfdb_argument{end+1}='-i';
    wfdb_argument{end+1}=[num2str(N0-1)];
end
if(~isempty(N))
    wfdb_argument{end+1}='-I';
    wfdb_argument{end+1}=[num2str(N-1)];
end
if(~isempty(minM))
    wfdb_argument{end+1}='-m';
    wfdb_argument{end+1}=[num2str(minM)];
end
if(~isempty(maxM))
    wfdb_argument{end+1}='-M';
    wfdb_argument{end+1}=[num2str(maxM)];
end
if(~isempty(maxScale))
    wfdb_argument{end+1}='-n';
    wfdb_argument{end+1}=[num2str(maxScale)];
end
if(~isempty(minR))
    wfdb_argument{end+1}='-r';
    wfdb_argument{end+1}=[num2str(minR)];
end
if(~isempty(maxR))
    wfdb_argument{end+1}='-R';
    wfdb_argument{end+1}=[num2str(maxR)];
end
javaWfdbExec.setArguments(wfdb_argument);

if(config.inOctave)
    x=cellstr(num2str(x));
    x=java2mat(javaWfdbExec.execWithStandardInput(x));
    Nx=x.size;
    out=cell(Nx,1);
    for n=1:Nx
        out{n}=x.get(n-1);
    end
else
    out=cell(javaWfdbExec.execWithStandardInput(x).toArray);
end
M=length(out);
if(M<4)
    error(['Error calculating MSE:' out{:}])
end
info=out(1:3);
out(1:4)=[];
M=M-4;
scale=zeros(M,1)+NaN;
y=zeros(M,1)+NaN;
for m=1:M
    str=out{m};
    sep=regexp(str,'\s');
    scale(m)=str2num(str(1:sep));
    y(m)=str2num(str(sep(1):sep(2)));
end

for n=1:nargout
    eval(['varargout{n}=' outputs{n} ';'])
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox