text:
% ------------------------------------------------------------------------
% Copyright (C) 2016 University of Southern California, SAIL, U.S.
% Author: Maarten Van Segbroeck
% Mail: mvansegbroeck@gmail.com
% Date: 2016-20-12
% ------------------------------------------------------------------------
% vadout = apply_vad(audiodir,p1,p2);
%
%      Apply Voice Activity Detection to all files in a specified audio directory
%
%      --IN--
%      audiodir: directory of audio files (WAV format)
%      p1: speech/non-speech threshold [default:0.1]
%      p2: speech region smoothing [default:20]
%
%      --OUT--
%      vadout: VAD labels at frame level (10 ms)
%
%      When using please cite:
%
%      @article{van2013robust,
%       title={A Robust Frontend for VAD: Exploiting Contextual, Discriminative and Spectral Cues of Human Voice},
%       author={Van Segbroeck, Maarten and Tsiartas, Andreas and Narayanan, Shrikanth},
%       year={2013}
%     }

function vadout=apply_vad(audiodir,p1,p2)

if ~exist('audiodir','var')
  fprintf('Please specify an audio directory\n')
  return
end

% ---vad parameters (to tune)---
% speech noise threshold
if ~exist('p1','var')
  p1=0.1;
end
% speech region smoothing
if ~exist('p2','var')
  p2=20;
end

% set path
addpath mfiles/

filenames=dir(sprintf('%s/*.wav',audiodir));
fs=8000;

% ---feature---
% GT
NbCh=64;
% Gabor
nb_mod_freq=2;
% LTSV
R=50; % context
M=10; % smoothing
ltsvThr=0.5;
ltsvSlope=0.2;
% vprob2 and ltsv2
K=30; order=4;
% -------------

% ---vad model---
load('models/model.mat')

% ---visualize---
visualize=true;

for k = 1 : length(filenames)

 % read in audio
 [sam,fs_orig]=audioread(fullfile(audiodir,filenames(k).name));
 sam_8k=downsample(sam(:,1),fs_orig/fs);

 % [1] extract cochleagram
 gt=FE_GT(sam_8k,fs,NbCh);

 % [2] Gabor filtering applied on mel
 gbf=FE_GBF(sam_8k,fs,nb_mod_freq,false);
 gbf= [gbf gbf(:,ones(1,10)*size(gbf,2))];
 gbf = gbf(:,1:size(gt,2));

 % [3] LTSV
 ltsv=FE_LTSV(sam_8k,fs,R,M,gt,ltsvThr,ltsvSlope);
 ltsv2 = convert_to_context_stream(ltsv, K, order);
 ltsv2= [ltsv2 ltsv2(:,ones(1,10)*size(ltsv,2))];
 ltsv2 = ltsv2(:,1:size(gt,2));

 % [4] vprob prob
 vprob=voicingfeature(sam_8k,fs);
 vprob2 = convert_to_context_stream(vprob, K, order);
 vprob2 = [vprob2 vprob2(:,ones(1,10)*size(vprob,2))];
 vprob2 = vprob2(:,1:size(gt,2));

 % feature for VAD
 test_x = [gt;gbf;ltsv2;vprob2];
 test_x_norm = mvn(test_x);

 % VAD decoding
 [~,~,output] = nntest(dnn, test_x_norm');
 outprob=double(output(:,1));
 vadout=medfilt1(outprob.^2,p2)>p1;

 if visualize
  imagesc(mvn(gt));axis xy;hold on;
  plot(10*vadout,'m','LineWidth',3); zoom xon; hold off
 end
end

author: mvansegbroeck-zz

repo: vad


text:
%% (Internal) Create new QRS detections based on other lead/algorithms detections
%
%   artificial_annotations = combine_anns(time_serie, estimated_labs, header)
%
% Arguments:
%
%      + time_serie: all the QRS detections
%
%      + estimated_labs: estimated labels (TP, FP, FN) of each detection
%
%      + header: ECG header struct.
%
% Output:
%
%      + artificial_annotations: all the QRS detections plus the new ones
%      created.
%
% Example:
%
%         aux_idx = best_detections_idx(1:min(10, length(best_detections_idx)));
%         artificial_annotations = combine_anns(all_annotations(aux_idx), estimated_labs(aux_idx), ECG_struct.header);
%
% See also QRScorrector
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function artificial_annotations = combine_anns(time_serie, estimated_labs, header)

    % cantidad de leads artificiales que generar�
    min_artificial_leads = 3;
    win_size = 20; % seconds
    win_size = round(win_size*header.freq); % milliseconds

%     lreferences = length(time_serie);
    for ii = 1:min_artificial_leads
        artificial_annotations(ii).time = [];
    end

    start_sample = min(cell2mat(cellfun(@(a)(min(a)),time_serie, 'UniformOutput', false)));
    end_sample = max(cell2mat(cellfun(@(a)(max(a)),time_serie, 'UniformOutput', false))) + 1;

%     aux_seq = (start_sample+win_size):round(win_size/2):end_sample;
    aux_seq = min(start_sample+win_size, end_sample):win_size:end_sample;

    if( isempty(aux_seq) )
        str_aux = disp_string_framed(0, 'Recording too short for combining detections'  );
        fprintf(2, '\n%s\n', str_aux);
        return
    end

    if(aux_seq(end) ~= end_sample )
        aux_seq = [aux_seq (aux_seq(end)+win_size) ];
    end

    laux_seq = length(aux_seq);

    aux_idx = arrayfun(@(a)( ...
                                cellfun( @(b)( ...
                                                findStartEnd(  b >= (a - win_size) & b < a ) ...
                                            ), time_serie, 'UniformOutput', false  )...
                            ), aux_seq, 'UniformOutput', false);

    aux_q = cellfun( @(a)( cell2mat(cellfun( @(b, c)( ...
                                    calc_q_val( b, c ) ...
                                ), estimated_labs, a, 'UniformOutput', false  ))...
                ), aux_idx, 'UniformOutput', false  );

    [~, aux_q_idx] = cellfun( @(a)( sort(a, 'descend') ), aux_q, 'UniformOutput', false );

    % TODO: hacer una estrategia que divida los segmentos de  aux_seq hasta
    % que no haya problemas de transiciones de los segmentos solapados.
%     while()
%         aux_idx = arrayfun( @(ii)( cell2mat(arrayfun( @(jj)( find_disagreements(aux_idx, aux_q_idx, ii, jj ) ), 1:laux_seq , 'UniformOutput', false)) ), 1:3, 'UniformOutput', false );
%     end

    aux_val = arrayfun( @(ii)( cell2mat(cellfun( @(a, q_idx)( build_combined_series(time_serie, a, q_idx, ii ) ), aux_idx, aux_q_idx, 'UniformOutput', false)) ), 1:min_artificial_leads, 'UniformOutput', false );

    cant_artificial_leads = min(min_artificial_leads, sum(cellfun(@(a)(~isempty(a)), aux_val )) );

    if( cant_artificial_leads > 0 )

        for ii = 1:cant_artificial_leads
            % avoid annotations very close each other.
            aux_time_serie = aux_val{ii};
            aux_time_serie(find( diff(sort(aux_time_serie)) <= round(0.15 * header.freq) ) +1) = [];
            if( ~isempty(aux_time_serie) )
                artificial_annotations(ii).time = colvec(aux_time_serie);
            end
        end

        if( cant_artificial_leads < min_artificial_leads )
            for ii = cant_artificial_leads+1:min_artificial_leads
                artificial_annotations(ii).time = artificial_annotations(cant_artificial_leads).time;
            end
        end

    end

function start_end_aux = findStartEnd( bAux )

    start_aux = find(  bAux, 1, 'first' );
    end_aux = find(  bAux, 1, 'last' );
    start_end_aux = [start_aux end_aux];

function this_q = calc_q_val(this_labs, strt_end)
%         win_size in milliseconds

    if( isempty(strt_end) || isempty(this_labs) )
        this_q = [];
    else
        this_labs = this_labs(strt_end(1):strt_end(2));

        this_se = sum(this_labs == 3) / sum(this_labs == 3 | this_labs == 1) ;
        this_pp = sum(this_labs == 3) / sum(this_labs == 3 | this_labs == 2) ;

        % paper metric
        %this_q = (2*this_se + this_pp)/3;
        % F1 score
        this_q = 2/(1/this_se + 1/this_pp);
    end

function new_str_end = find_disagreements(str_end, q_idx, ii, jj)

    this_q_idx = q_idx{jj};
    next_q_idx = q_idx{jj+1};

    if( isempty(this_q_idx) || ii > length(q_idx) )
        new_str_end = [];
    else
        if( q_idx(ii) == q_idx(ii+1) )
            new_str_end = [];
        else

        end
    end

function combined_series = build_combined_series(ts, str_end, q_idx, ii)

    if( isempty(q_idx) || ii > length(q_idx) )
        combined_series = [];
    else
        q_idx = q_idx(ii);
        aux_idx = str_end{q_idx};
        aux_val = ts{q_idx};
        if( isempty(aux_idx) || isempty(aux_val)  )
            combined_series = [];
        else
            combined_series = rowvec(aux_val(aux_idx(1):aux_idx(2)));
        end
    end

author: marianux

repo: ecg-kit

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


text:
%% 1 Param, 1 State
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1);
z0 = 1; %ic

%Generate measurement data
p = 1.5;                    %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements
zm(end) = zm(end)*1.1;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'bounds',0,2,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 1 State, End Point Weighted [No Sensitivity]
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1);
z0 = 1; %ic

%Generate measurement data
p = 1.5;                    %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements
zm(end) = zm(end)*1.1;

weights = ones(size(zm)); weights(end) = 100;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','none');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'bounds',0,2,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 1 State, End Point Weighted [With Sensitivity]
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1);
z0 = 1; %ic

%Generate measurement data
p = 1.5;                    %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements
zm(end) = zm(end)*1.1;

weights = ones(size(zm)); weights(end) = 100;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'bounds',0,2,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 1 State, Reversed Timestamps
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1);
z0 = 1; %ic

%Generate measurement data
p = 1.5;                    %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements
zm = flipud(zm); tm = fliplr(tm);

zm(3) = zm(3)*1.5;
weights = ones(size(zm)); weights(3) = 0;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'bounds',0,2,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 1 State, End Point Weighted [No Sensitivity] as NLP
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1);
z0 = 1; %ic

%Generate measurement data
p = 1.5;                    %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements
zm(end) = zm(end)*1.1;

weights = ones(size(zm)); weights(end) = 100;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','none');
opts = optiset('solver','ipopt','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'bounds',0,2,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 1 State, End Point Weighted [With Sensitivity] as NLP
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1);
z0 = 1; %ic

%Generate measurement data
p = 1.5;                    %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements
zm(end) = zm(end)*1.1;

weights = ones(size(zm)); weights(end) = 100;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('solver','ipopt','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'bounds',0,2,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 2 State
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = 2.345;                  %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements
zm(4,1) = zm(4,1)*1.5;

weights = ones(size(zm)); weights(4,1) = 0;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 2 State [only fitting first state]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); z(1)]; %note z1 in 2nd ode allows us to estimate p1
z0 = [1;3]; %ic

%Generate measurement data
p = 2.345;                  %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements

%State to Measure
state = 1;
zm = zm(:,state); zm(4) = zm(4)*1.5;

weights = ones(size(zm)); weights(4) = 0;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('stateIndex',state,'sensitivity','nd');
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 2 State [only fitting second state]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); z(1)]; %note z1 in 2nd ode allows us to estimate p1
z0 = [1;3]; %ic

%Generate measurement data
p = 2.345;                  %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements

%State to Measure
state = 2;
zm = zm(:,state); zm(4) = zm(4)*1.5;

weights = ones(size(zm)); weights(4) = 0;

%Build OPTI Object
p0 = 1; %inital parameter guess
dopts = optidynset('stateIndex',state,'sensitivity','nd');
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 2 State [different measurement times]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = 2.345;                  %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm1 = 0:0.2:1;              %measurement times
tm2 = 0:0.15:1;
[~,zm] = ode45(oi,tm1,z0); zm1 = zm(:,1); %measurements
[~,zm] = ode45(oi,tm2,z0); zm2 = zm(:,2); %measurements

zm1(4) = zm1(4)*1.5;

tm = {tm1;tm2};
zm = {zm1;zm2};

wts1 = ones(size(zm1)); wts1(4) = 0;
wts2 = ones(size(zm2));
weights = {wts1;wts2};

%Build OPTI Object
theta0 = 1; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',theta0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 2 Param, 2 State
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); p(2)*z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = [2.345;1.1];            %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements

tm = fliplr(tm);
zm = flipud(zm);
zm(3,1) = zm(3,1)*1.5;

weights = ones(size(zm));
weights(3,1) = 0;

%Build OPTI Object
p0 = [1;0.1]; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm(:),'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 2 Param, 2 State [Different Measurement Times]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); p(2)*z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = [2.345;1.1];            %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm1 = 0:0.2:1;              %measurement times
tm2 = 0:0.15:1;
[~,zm] = ode45(oi,tm1,z0); zm1 = zm(:,1); %measurements
[~,zm] = ode45(oi,tm2,z0); zm2 = zm(:,2); %measurements

zm2(4) = zm2(4)*0.1;
tm = {tm1;tm2};
zm = {zm1;zm2};

wts1 = ones(6,1);
wts2 = ones(size(zm2)); wts2(4) = 0;
weights = {wts1;wts2};

%Build OPTI Object
p0 = [1;0.1]; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm(:),'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 2 Param, 2 State [Only fit 2nd state]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); p(2)*z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = [2.345;1.1];            %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm] = ode45(oi,tm,z0);   %measurements

%State to Measure
state = 2;
zm = zm(:,state); zm(4) = zm(4)*1.5;

weights = ones(size(zm)); weights(4) = 0;

%Build OPTI Object
p0 = [1;0.1]; %inital parameter guess
dopts = optidynset('stateIndex',state,'sensitivity','nd');
opts = optiset('solver','lmder','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm(:),'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% Lorenz System [Analytical Derivatives + Different Measurement Times + Non-Zero Initial Time + Estimate IC]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*(z(2) - z(1));
                z(1)*(p(2) - z(3)) - z(2);
                z(1)*z(2) - p(3)*z(3)];
z0 = [5.7;10.50;30.58]; %ic
%Analytical Derivatives
dfdz = @(t,z,p) [-p(1), p(1), 0;
                 p(2) - z(3), -1, -z(1);
                 z(2), z(1), -p(3)];
dfdp = @(t,z,p) [z(2) - z(1), 0, 0;
                 0, z(1), 0;
                 0, 0, -z(3)];

%Generate measurement data
p = [10,46,8/3];            %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm1 = 0:0.2:2;
tm2 = 0.2:0.1:4;            %measurement times
tm3 = 0.5:0.5:3;
[~,zm] = ode45(oi,tm1,z0); zm1 = zm(:,1); %measurements
[~,zm] = ode45(oi,[0 tm2],z0); zm2 = zm(2:end,2); %measurements
[~,zm] = ode45(oi,[0 tm3],z0); zm3 = zm(2:end,3); %measurements

zm3(3) = zm3(3)*1.2;

tm = {tm1;tm2;tm3};
zm = {zm1;zm2;zm3};

wts1 = ones(size(zm1));
wts2 = ones(size(zm2));
wts3 = ones(size(zm3)); wts3(3) = 0;
weights = {wts1;wts2;wts3};

%Given states 2 + 3 start from non-zero time, we should estimate
z0(2:3) = NaN;

%Build OPTI Object
p0 = [10+0.1,46+0.2,8/3,10.4,30.55]; %inital parameter guess
dopts = optidynset('dfdz',dfdz,'dfdp',dfdp);
opts = optiset('solver','nl2sol','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm,zm,'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 1 State [Repeated Measurements + z0]
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1);
z0 = 1; %ic

%Generate measurement data
p = 1.5;                    %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm1] = ode45(oi,tm,z0);        %measurements RUN 1
[~,zm2] = ode45(oi,tm,z0*0.95);   %measurements RUN 2

%Concatenate measurements
tm_m = [tm tm]';
zm_m = [zm1;zm2];

wts1 = 2*ones(size(zm1));
wts2 = ones(size(zm2));
weights = [wts1;wts2];

%Estimate initial condition
z0 = NaN;

%Build OPTI Object
theta0 = [1;0.5]; %inital parameter guess
dopts = optidynset('sensitivity','nd');
opts = optiset('solver','nl2sol','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_m,zm_m,'x0',theta0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 2 Param, 1 State [Repeated Measurements + z0]
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1) + p(2);
z0 = 1; %ic

%Generate measurement data
p = [2.5; 5.5];             %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm1] = ode45(oi,tm,z0);        %measurements RUN 1
[~,zm2] = ode45(oi,tm,z0*0.5);   %measurements RUN 2

%Concatenate measurements
tm_m = [tm tm]';
zm_m = [zm1;zm2];

wts1 = 2*ones(size(zm1));
wts2 = ones(size(zm2));
weights = [wts1;wts2];

%Replace z0 with NaN to indicate to estimate it
z0 = NaN;

%Build OPTI Object
theta0 = [1;1;0.5]; %inital parameter guess + initial state guess
dopts = optidynset('integrator','ode45','sensitivity','nd');
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_m,zm_m,'x0',theta0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 1 Param, 2 State, Repeated Measurements + Both z0
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = 2.345;                  %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm1] = ode45(oi,tm,z0);        %measurements RUN 1
[~,zm2] = ode45(oi,tm,z0*0.9);   %measurements RUN 2

%Concatenate measurements
tm_m = [tm tm]';
zm_m = [zm1;zm2];

wts1 = 2*ones(size(zm1));
wts2 = ones(size(zm2));
weights = [wts1;wts2];

%Replace z0 with NaN to indicate to estimate it
z0(1:2) = NaN;

%Build OPTI Object
theta0 = [1;0.5;2.5]; %inital parameter guess + initial state guess
dopts = optidynset('integrator','ode45','sensitivity','nd');
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_m,zm_m,'x0',theta0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 2 Param, 2 State + Solve for both z0 + Repeated Measurements [Only fit 1st state]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); p(2)*z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = [2.345;1.1];            %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm1] = ode45(oi,tm,z0);        %measurements RUN 1
[~,zm2] = ode45(oi,tm,z0*0.9);   %measurements RUN 2

%Concatenate measurements
tm_m = [tm tm]';
zm_m = [zm1;zm2];

wts1 = 2*ones(size(zm1));
wts2 = ones(size(zm2));
weights = [wts1;wts2];

%Replace z0 with NaN to indicate to estimate it
z0(1:2) = NaN;

%States to Measure
state = 1;
zm_m = zm_m(:,state);
weights = weights(:,state);

%Build OPTI Object
p0 = [1;0.1;0.1;0.1]; %inital parameter guess + initial state guess
dopts = optidynset('stateIndex',state,'sensitivity','nd');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_m,zm_m,'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% 2 Param, 2 State + Solve for both z0 + Repeated Measurements [Only fit 2nd state]
clc
%ODE to fit
ode = @(t,z,p) [p(1)*z(1) + z(2); p(2)*z(1)];
z0 = [1;3]; %ic

%Generate measurement data
p = [2.345;1.1];            %true parameter value
oi = @(t,z) ode(t,z,p);     %ode integrator function
tm = 0:0.2:1;               %measurement times
[~,zm1] = ode45(oi,tm,z0);        %measurements RUN 1
[~,zm2] = ode45(oi,tm,z0*0.9);   %measurements RUN 2

%Concatenate measurements
tm_m = [tm tm]';
zm_m = [zm1;zm2];

wts1 = 2*ones(size(zm1));
wts2 = ones(size(zm2));
weights = [wts1;wts2];

%Replace z0 with NaN to indicate to estimate it
z0(1:2) = NaN;

%States to Measure
state = 2;
zm_m = zm_m(:,state);
weights = weights(:,state);

%Build OPTI Object
p0 = [1;0.1;0.1;0.1]; %inital parameter guess + initial state guess
dopts = optidynset('stateIndex',state,'sensitivity','nd');
opts = optiset('solver','auto','display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_m,zm_m,'x0',p0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

%% DIW MAPLE ODEs [Analytical Derivatives + IC Estimate + Different Measurement Times + Repeated Points + Non-Zero Initial Time]
clc
ode = @(t,z,p) [-p(1)*z(1) + 4;
                2*z(1) - p(1)*z(2) + 5;
                -4*z(1) - 2*z(3) - p(2)];
z0 = [-1.5;1.25;1]; %ic

%Analytical Derivatives
dfdz = @(t,z,p) [-p(1), 0, 0;
                 2, -p(1), 0;
                 -4, 0, -2];
dfdp = @(t,z,p) [-z(1), 0
                 -z(2), 0
                 0, -1];

% Measurement Times for each State
 tm1  = 0.5:0.1:2;              %state 1
 tm2  = 1:0.1:2;                %state 2
 tm3  = 0.2:0.2:2;              %state 3

 p = [2.345;1.1];            %true parameter value
% Solve ODEs and index Measurement Data
% Note 0 is required below just to match initial condition in this example
odeInt = @(t,z) ode(t,z,p);     %ode integrator function
[~,zm1] = ode45(odeInt,[0 tm1],z0); zm1 = zm1((2:end),1);
[~,zm1b] = ode45(odeInt,[0 tm1],z0*0.5); zm1b = zm1b((2:end),1);
[~,zm2] = ode45(odeInt,[0 tm2],z0); zm2 = zm2((2:end),2);
[~,zm3] = ode45(odeInt,[0 tm3],z0); zm3 = zm3((2:end),3);

% Group measurements and time stamps in cell arrays
 tm_multi = {[tm1 tm1];tm2;tm3};
 zm_multi = {[zm1;zm1b];zm2;zm3};

 wts1 = ones(size(zm1));
 wts1b = 10*ones(size(zm1));
 wts2 = ones(size(zm2));
 wts3 = ones(size(zm3));
 weights = {[wts1;wts1b];wts2;wts3};

% We will need to estimate all initial conditions in this problem
 z0 = [NaN;NaN;NaN];

% New Initial Guess Vector [p;z0];
 theta0 = [1;0.5;0.5;0.5;0.5];

%Build OPTI Object
dopts = optidynset('dfdz',dfdz,'dfdp',dfdp,'initialT',0,'sensitivity','user');
opts = optiset('solver','auto','display','iter','derivCheck','on','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_multi,zm_multi,'x0',theta0,'z0',z0,'weights',weights,'options',opts)

[x,f,e,i] = solve(Opt)
plot(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Draw various graphs obtained from the same point set.
%
%   output = drawSampleGraphs(input)
%
%   Example
%   drawSampleGraphs
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2020-06-04,    using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2020 INRAE.

%% Read input points

% read a set of 16 points within a [10 28]^2 bounding box
pts = load('sedgewick_points.txt');

% display the sample points
figure; axis equal; axis([10 28 10 28]); hold on;
drawPoint(pts, 'bo');
title('Sample points');

%% Delaunay Graph

% compute the graph
[nodes, edges] = delaunayGraph(pts);

% display the sample points
figure; axis equal; axis([10 28 10 28]); hold on;
drawPoint(pts, 'bo');

% display the graph
drawGraphEdges(nodes, edges);
title('Delaunay Graph');

%% Minimal Spanning Tree

% compute the graph
edges = euclideanMST(pts);

% display the sample points
figure; axis equal; axis([10 28 10 28]); hold on;
drawPoint(pts, 'bo');

% display the graph
drawGraphEdges(pts, edges);
title('Minimal Spanning Tree');

%% Relative Neighborhood Graph

% compute the graph
edges = relativeNeighborhoodGraph(pts);

% display the sample points
figure; axis equal; axis([10 28 10 28]); hold on;
drawPoint(pts, 'bo');

% display the graph
drawGraphEdges(pts, edges);
title('Relative Neighborhood Graph');

%% Gabriel Graph

% compute the graph
edges = gabrielGraph(pts);

% display the sample points
figure; axis equal; axis([10 28 10 28]); hold on;
drawPoint(pts, 'bo');

% display the graph
drawGraphEdges(pts, edges);
title('Gabriel Graph');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [X1, Y1] = Direction2CubeMap(D, r, c)
%
%        [X1, Y1] = Direction2CubeMap(D, r, c)
%
%
%        Input:
%           -D: 3D directions of the img format
%        Output:
%           -X1: X coordinates in the CubeMap format
%           -Y1: Y coordinates in the CubeMap format
%
%     Copyright (C) 2011  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

[rD, cD, ~] = size(D);

X1 = zeros(rD, cD);
Y1 = zeros(rD, cD);
totD = rD * cD;

face_cm = round((r / 4 + c / 3) / 2);

%Calculating faces' directions
Mul = [1, 1, 1, -1, -1, -1];

T = [2, 1, 3, 2, 1, 3];
A = [1, 3, 1, 1, 3, 1];
B = [3, 2, 2, 3, 2, 2];

X1_1 = [1.5,2.5, 1.5, 1.5, 0.5, 1.5];
X1_2 = [0.5,0.5, 0.5,-0.5, 0.5,-0.5];
Y1_1 = [2.5,1.5, 3.5, 0.5, 1.5, 1.5];
Y1_2 = [0.5,0.5,-0.5, 0.5,-0.5,-0.5];

for i=1:6
    indx = find(    (Mul(i) * D(:,:,T(i)) > 0)&...
                    (Mul(i) * D(:,:,T(i)) >= abs(D(:,:,A(i))))&...
                    (Mul(i) * D(:,:,T(i)) >= abs(D(:,:,B(i)))));

    indx_shifted = indx + (T(i) - 1) * totD;

    X1(indx) = X1_1(i) + X1_2(i) * D(indx + (A(i) - 1) * totD) ./ D(indx_shifted);
    Y1(indx) = Y1_1(i) + Y1_2(i) * D(indx + (B(i) - 1) * totD) ./ D(indx_shifted);
end

X1 = RemoveSpecials(X1) * face_cm;
Y1 = flipud(RemoveSpecials(Y1) * face_cm);

end

author: banterle


text:
function [edges,weights] = spm_vb_edgeweights(vxyz,img)
% Compute edge set and edge weights of a graph
% FORMAT [edges,weights]= spm_vb_edgeweights(vxyz,img)
%
% vxyz     list of neighbouring voxels (see spm_vb_neighbors)
% img      image defined on the node set, e.g. wk_ols. The edge weights
%          are uniform if this is not given, otherwise they are a function
%          of the distance in physical space and that between the image
%          at neighoring nodes

% edges    [Ne x 2] list of neighboring voxel indices
% weights  [Ne x 1] list of edge weights
% Ne       number of edges (cardinality of edges set)
% N        number of nodes (cardinality of node set)
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging

% Lee Harrison
% $Id: spm_vb_edgeweights.m 6079 2014-06-30 18:25:37Z spm $

N       = size(vxyz,1);
[r,c,v] = find(vxyz');
edges   = [c,v];
% undirected graph, so only need store upper [lower] triangle
i       = find(edges(:,2) > edges(:,1));
edges   = edges(i,:);
if nargin < 2,
    weights = ones(size(edges,1),1);
    return
else
    ka  = 16;
    M   = mean(img,2)*ones(1,N);
    C   = (1/N)*(img-M)*(img-M)';
    Hf  = inv(C);
    A   = spm_vb_incidence(edges,N);
    dB  = img*A'; % spatial gradients of ols estimates
    dg2 = sum((dB'*Hf).*dB',2); % squared norm of spatial gradient of regressors
    ds2 = 1 + dg2; % squared distance in space is 1 as use only nearest neighbors
    weights = exp(-ds2/ka); % edge weights
end

author: spm


text:
classdef MAD < Algorithms.CDAlg
    properties
        nMADUsed = uint8(255);
    end
    methods
        function obj = MAD(k)
            obj@Algorithms.CDAlg('Multivariate Alteration Detection', ...
                ['Nielsen, Allan & Conradsen, Knut & Simpson, James. (1998). ',...
                'Multivariate Alteration Detection (MAD) and MAF Postprocessing',...
                'in Multispectral, Bitemporal Image Data: ',...
                'New Approaches to Change Detection Studies.',...
                'Remote Sensing of Environment.',...
                '64. 1-19. 10.1016/S0034-4257(97)00162-4. '] ...
                );
            if exist('k', 'var'), obj.nMADUsed=uint8(k); end
        end
        [DI, k] = detectChange(obj, t1, t2);
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function res = ctranspose(a)
a.adjoint = xor(a.adjoint,1);
res = a;

author: thomaskuestner


text:
%ANGDIFF Difference of two angles
%
% ANGDIFF(TH1, TH2) is the difference between angles TH1 and TH2, ie. TH1-TH2
% on the circle.  The result is in the interval [-pi pi).  Either or both
% arguments can be a vector:
% - If TH1 is a vector, and TH2 a scalar then return a vector where TH2 is modulo
%   subtracted from the corresponding elements of TH1.
% - If TH1 is a scalar, and TH2 a vector then return a vector where the
%   corresponding elements of TH2 are modulo subtracted from TH1.
% - If TH1 and TH2 are vectors then return a vector whose elements are the modulo
%   difference of the corresponding elements of TH1 and TH2, which must be the
%   same length.
%
% ANGDIFF(TH) as above but TH=[TH1 TH2].
%
% ANGDIFF(TH) is the equivalent angle to the scalar TH in the interval [-pi pi).
%
% Notes::
% - The MathWorks Robotics Systems Toolbox defines a function with the same name
%   which computes TH2-TH1 rather than TH1-TH2.
% - If TH1 and TH2 are both vectors they should have the same
%   orientation, which the output will assume.
%

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function d = angdiff(th1, th2)

    switch nargin
        case 1
            if length(th1) == 2
                d = th1(1) - th1(2);
            else
                d = th1;
            end
        case 2
            if length(th1) > 1 && length(th2) > 1
                % if both arguments are vectors, they must be the same
                assert(all(size(th1) == size(th2)), 'SMTB:angdiff:badarg', 'vectors must be same shape');
            end
            % th1 or th2 could be scalar
            d = th1 - th2;
    end

    % wrap the result into the interval [-pi pi)
    d = mod(d+pi, 2*pi) - pi;
end

% Simplistic version of the code, easy to see what it does, but slow...
%
% for very negative angles keep adding 2pi
%     while true
%         k = find(d < -pi);
%         if isempty(k)
%             break;
%         end
%         d(k) = d(k) + 2*pi;
%     end
%
%     % for very positive angles keep subtracting 2pi
%     while true
%         k = find(d > pi);
%         if isempty(k)
%             break;
%         end
%         d(k) = d(k) - 2*pi;
%     end

author: petercorke

repo: spatialmath-matlab


text:
function [cut, map, short] = vl_aibcut(parents, n)
% VL_AIBCUT  Cut VL_AIB tree
%  CUT = VL_AIBCUT(PARENTS, N) cuts the binary merge tree PARENTS and
%  returns a cut CUT of N nodes. The format of PARENTS is the same
%  used by the VL_AIB() function.
%
%  A cut is a set of N nodes such that no node is a descendant of any
%  other node in the CUT and such that any leaf descend from a node in
%  the cut. CUT lists the nodes of the binary merge tree PARENT that
%  form the cut.
%
%  Nodes with null parent (as defined by PARENTS) are comprised in the
%  cut if the other nodes are not enough to fill a cut of N elements.
%
%  [CUT, MAP] = VL_AIBCUT(...) returns a vector MAP with the same size as
%  PARENTS. MAP assign each node below or in the cut to the
%  corresponding element in the CUT vector (each element above the cut
%  or with null parent is mapped to 0). To get the index of the
%  corresponding cut nodes use CUT(MAP). MAP by itself is useful to
%  quantize the leaves in a sequences of N contiguous indexes,
%  starting from one (see also VL_AIBCUTPUSH()).
%
%  [CUT, MAP, SHORT] = VL_AIBCUT(...) returns also a vector SHORT that
%  represents a merge tree compressed to the cut. This is obtained by
%  mapping all nodes below to the cut element above them. Nodes in or
%  above the cut are mapped to themselves. Null parents are left
%  unchanged, except if the corresponding node is in the cut (in which
%  case the map-to-itself rule has the precedence).
%
%  See also: VL_HELP(), VL_AIB().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% --------------------------------------------------------------------
%                                           Determine nodes in the cut
% --------------------------------------------------------------------

if n > 1
  root    = max(parents) ;

  % count number of null nodes
  z = sum(parents(1:root) == 0) ;

  % determine number of leves
  nleaves = (root - z + 1) / 2 ;

  % find first node of the cut
  mu   = root - min(n, nleaves) + 1 ;

  % correction for presence of null nodes
  nz   = find(parents(1:mu) > 0) ;
  mu   = nz(end) ;

  % find node belnoging to the cut
  cut  = find(parents(1:mu) > mu) ;

  % In the presence of null nodes, the cut size might exceed
  % nleaves, which is the maximum cut size we can obtain with the
  % specified tree. The additional nodes have to be picked up from
  % the null nodes.

  if length(cut) < n
    sel_z = find(parents == 0) ;
    cut = [sel_z(1:n-length(cut)) cut] ;
  end

  % aesthetic reasons only
  cut = sort(cut) ;

else
  mu   = max(parents) ;
  cut  = mu ;
end

% --------------------------------------------------------------------
%                                       Short-circuit nodes to the cut
% --------------------------------------------------------------------

stop = [cut find(parents == 0)] ;
short = 1:length(parents) ;

while 1
  [drop,sel] = setdiff(short(1:mu), stop)  ;
  sel = setdiff(sel, stop) ;
  if isempty(sel), break ; end
  short(sel) = parents(short(sel))  ;
end

short(setdiff(find(parents == 0), cut)) = 0 ;

% --------------------------------------------------------------------
%                                                  Build quantizer map
% --------------------------------------------------------------------

map             = 1:numel(parents) ;
map(cut)        = 1:n ;
map(short >  0) = map(short(short > 0)) ;
map(short == 0) = 0 ;
map(mu+1:end)   = 0 ;

author: jianxiongxiao

repo: ProfXkit


text:
classdef(Abstract) AbstractLaunchVehicleElectricalPowerSrcSnk < matlab.mixin.SetGet & matlab.mixin.Heterogeneous
    %AbstractLaunchVehicleElectricalPowerSrc Summary of this class goes here
    %   Detailed explanation goes here

    properties(Abstract)
        id
    end

    properties(Dependent)
        lvdData
    end

    methods
        name = getName(obj);

        stage = getAttachedStage(obj);

        pwrRate = getElectricalPwrRate(obj, elemSet, steeringModel, hasSunLoS, body2InertDcm, elemSetSun);

        newState = createDefaultInitialState(obj, stageState)

        useTF = openEditDialog(obj);

        tf = isInUse(obj);

        newObj = copy(obj);

        summStr = getSummaryStr(obj);

        function lvdData = get.lvdData(obj)
            lvdData = obj.getAttachedStage().launchVehicle.lvdData;
        end
    end

    methods(Sealed)
        function tf = eq(a,b)
            tf = [a.id] == [b.id];
        end
    end
end