text:
function [c, r, t, cpulse, acq_codes] = tapas_physio_read_physlogfiles_brainproducts(log_files, ...
    cardiac_modality, verbose)
% reads out physiological time series (ECG, PMU, resp belt) and timing vector for BrainProducts .eeg file
%
%   [cpulse, rpulse, t, c, acq_codes] = tapas_physio_read_physlogfiles_brainproducts(logfiles, ...
%                               verbose)
%
%
% IN    log_files
%       .cardiac        contains ECG or pulse oximeter time course
%                           for BrainProducts: usually the same as respiration
%       .respiration    contains breathing belt amplitude time course
%                           for BrainProducts: usually the same as cardiac
%       .sampling_interval  is ignored here, read from logfile
%
%   cardiac_modality
%                       'ecg1_filtered'     filtered 1st ECG channel signal
%                                           (Default)
%                       'ecg2_filtered'     filteered 2nd ECG channel
%                                           (sometimes less gradient artifacts)
%                       'ecg1_raw'          raw 1st ECG channel
%
%       verbose
%       .level              debugging plots are created if level >=3
%       .fig_handles        appended by handle to output figure
%
% OUT
%   r                   respiratory time series
%   c                   cardiac time series (ECG or pulse oximetry)
%   t                   vector of time points (in seconds)
%   cpulse              time events of R-wave peak in cardiac time series (seconds)
%                       for Biopac: usually empty, kept for compatibility
%   acq_codes           slice/volume start events marked by number <> 0
%                       for time points in t
% EXAMPLE
%   [ons_secs.cpulse, ons_secs.rpulse, ons_secs.t, ons_secs.c] =
%       tapas_physio_read_physlogfiles_GE(logfiles);
%
%   See also tapas_physio_main_create_regressors

% Author: Lars Kasper
% Created: 2017-01-27
% Copyright (C) 2017 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% user input

% volume trigger event
event_type = 'Response';
event_value = 'R128';

switch lower(cardiac_modality)
    case {'ecg','ecg1_filtered'}
        % chose one ecg channel (only for visualisation)
        ecg_ch = 1;
    case 'ecg2_filtered'
       % chose one ecg channel (only for visualisation)
        ecg_ch = 2;
end

% is ecg data flipped?
ecg_is_flipped = 1;

%% data

% read data and header info using field trip (included in SPM)
if ~exist('ft_read_header', 'file')
   pathSpm = fileparts(which('spm'));
   addpath(genpath(fullfile(pathSpm, 'external', 'fieldtrip')));
end

hdr = ft_read_header(log_files.cardiac);
data = ft_read_data(log_files.cardiac);

fs = hdr.Fs; % sampling frequency in Hz
N = hdr.nSamples; % number of samples
dt = 1/fs; % sampling interval in seconds
t = linspace(0,dt*(N-1), N); % time vector in seconds

fh = [];

% plot first 10 seconds of the raw data
fh(end+1,1) = tapas_physio_get_default_fig_params();
plot_end = 10 * fs;
plot(t(1:plot_end), data(:,(1:plot_end)));
xlabel('time in seconds');
ylabel(hdr.chanunit);
legend(hdr.label);

% extract ECG data
if ecg_is_flipped
    s = -data(ecg_ch,:);
else
    s = data(ecg_ch,:);
end

% plot first 10 seconds again for check
fh(end+1,1) = tapas_physio_get_default_fig_params();
plot_end = 10 * fs;
plot(t(1:plot_end), s(:,(1:plot_end)));
xlabel('time in seconds');
ylabel(hdr.chanunit{ecg_ch});
legend(hdr.label{ecg_ch});
%% events
% display all events in the data on the command line
cfg = [];
cfg.dataset = log_files.cardiac;
cfg.trialdef.eventtype  = '?';
ft_definetrial(cfg);

% define volume trigger
cfg.trialdef.eventtype = event_type;
cfg.trialdef.eventvalue = event_value;
cfg = ft_definetrial(cfg);
trigger_pos = cfg.trl(:,1);
n_trigger = length(trigger_pos);

% plot raw data and events
fh(end+1,1) = tapas_physio_get_default_fig_params();
max_value = max(s);
plot(t, s);
hold all;
stem(t(trigger_pos), max_value * ones(1,n_trigger));

% find segments
% compute the number of samples between each trigger (pos-1)
diff_trigger_pos = diff(trigger_pos);
% find the positons where a change is happening
diff_diff_trigger_pos = find(diff(diff_trigger_pos));
% add first and last event
start_segment = [1; diff_diff_trigger_pos(2:2:end)+1];
end_segment = [diff_diff_trigger_pos(1:2:end)+1; n_trigger];
% number of samples per trials
n_trial_samples = diff_trigger_pos(start_segment);
% position of start segment in samples
sample_start_segment = trigger_pos(start_segment);
% position of end segment in samples (not including the whole trial)
sample_end_segment = trigger_pos(end_segment);
% number of segments
n_segments = length(start_segment);
% number of trials per segment
n_trials = end_segment - start_segment + 1;
% TR per trial
TR_trials = n_trial_samples*dt;

% plot segments
stem(t(sample_start_segment), max_value * ones(1,n_segments), '--', 'LineWidth', 5);
stem(t(sample_end_segment), max_value * ones(1,n_segments), '--', 'LineWidth', 5);
legend('signal', 'all trigger', 'start segment', 'end segment');

%% load only first segment for sanity check
cfg = [];
cfg.dataset = log_files.cardiac;
cfg.trialdef.eventtype = event_type;
cfg.trialdef.eventvalue = event_value;
cfg.trialfun = 'ft_trialfun_segment';
cfg.trialdef.prestim = 0;
cfg.trialdef.poststim = TR_trials(1);
cfg.trialdef.segment_start = start_segment(1);
cfg.trialdef.segment_end = end_segment(1);
cfg = ft_definetrial(cfg);

% plot segments
fh(end+1,1) = tapas_physio_get_default_fig_params();
max_value = max(s);
plot(t, s);
hold all;
stem(t(cfg.trl(:,1)), max_value * ones(1,n_trials(1)));
stem(t(cfg.trl(:,2)), max_value * ones(1,n_trials(1)));
legend('signal', 'start trial', 'end trial');

%% save results
[~,name] = fileparts(file_name);
save_name = fullfile(file_path, [name, '_segments.mat']);
save(save_name, 'sample_start_segment', 'sample_end_segment', ...
    'n_segments', 'n_trials', 'TR_trials', 'start_segment', 'end_segment', ...
    'event_type', 'event_value', 'n_trial_samples', 'fs', ...
    'ecg_is_flipped', 'ecg_ch', 'file_path', 'file_name');

author: translationalneuromodeling

repo: tapas


text:
% function [Results, Gt, F1, Fs2, Fs3, Ft2, Ft3] = GenerativeTriTL(Train_Data,Test_Data,Parameter_Setting)
function [Results, pzd_t] = GenerativeTriTL(TrainData, TestData, TrainLabel, TestLabel, numIdentical, numAlike, numDistinct, numIter, numSource, numTrain, numTarget, numTest)
% function [Results, pz_d] = TriTL(Train_Data,Test_Data,Parameter_Setting)

% The common program for CD_PLSA, which can deal with multiple classes,
% multiple source domains and multiple target domains

%%%% Input:
% The parameter Train_data stores the file pathes of training data and the
% corresponding labels
% The parameter Test_data stores the file pathes of test data and the
% corresponding labels
% The parameter Parameterfile stores the parameter setting information

%%%% Output
% The variable Results is a matrix with size numIteration x numTarget, where
% numIteration is the number of iterations, numTarget is the number of
% target domains. Results record the detailed results of each iteration.

% The variable pz_d is a matrix with size n x c, where n is the number of
% instances in all target domains, specifically, n = n_1 + ... + nt (n_t is
% the number of instances in t-th target domain), c is the number of
% classes
%
% Note that if you want to deal with large data set, you should set larget
% memory for Matlab. You can set it in the file C:\boot.ini (This may not
% be true in your system), change '/fastdetect' to '/fastdetect /3GB'.
%
% Be good luck for your research, if you have any questions, you can
% contact the email: zhuangfz@ics.ict.ac.cn

%read the parameters
numK_1 = double(numIdentical);
numK_2 = double(numAlike);
numK_3 = double(numDistinct);
numIteration = double(numIter);

TrainX = TrainData;
TrainY = TrainLabel;
TestX = TestData;
TestY = TestLabel;
labelset = union(TestY,[]);

numC = length(labelset);
numFeature = size(TestX,1);

start = 1;
% if numK_3 == 0
%     start = 0;
% end
if start == 1
    DataSetX = [TrainX TestX];
    Learn.Verbosity = 1;
    Learn.Max_Iterations = 20;
    Learn.heldout = .1; % for tempered EM only, percentage of held out data
    Learn.Min_Likelihood_Change = 1;
    Learn.Folding_Iterations = 20; % for TEM only: number of fiolding in iterations
    Learn.TEM = 0; %tempered or not tempered
    [Pw_z,Pz_d,Pd,Li,perp,eta] = pLSA(DataSetX,[],numK_1+numK_2,Learn); %start PLSA
    %xlswrite(strcat('pwz_','common_selected','.xls'),Pw_z);
    csvwrite(strcat('pzw_','common_selected','.plsa'),Pw_z);
end
%pwy = xlsread(strcat('pwz_','common_selected','.xls'));

%% Following are Initializaitons
pzw = csvread(strcat('pzw_','common_selected','.plsa'));
pwy_a = pzw(:,1:numK_1); % the common topics using the same words
pwy_b_s = []; % the common topics using different words
pwy_b_t = []; % the common topics using different words
for i = 1:numSource
    pwy_b_s = [pwy_b_s, pzw(:,1+numK_1:numK_1+numK_2)];
end
for i = 1:numTarget
    pwy_b_t = [pwy_b_t, pzw(:,1+numK_1:numK_1+numK_2)];
end
% pwy_c_s = []; % different topics using different words
% pwy_c_t = []; % different topics using different words
pwy_c_s = ones(numFeature,numK_3*numSource)/numFeature;
pwy_c_t = ones(numFeature,numK_3*numTarget)/numFeature;
clear pzw;
pdz_s = zeros(sum(numTrain),numC);
for i = 1:size(pdz_s,1)
    pdz_s(i,TrainY(i)) = 1;
end
for i = 1:numSource
    pos = 0;
    if i > 1
        for t = 1:i-1
            pos = pos + numTrain(t);
        end
    end
    if i == 1
        pos = 0;
    end
    for j = 1:numC
        pdz_s(pos+1:pos+numTrain(i),j) = pdz_s(pos+1:pos+numTrain(i),j)/sum(pdz_s(pos+1:pos+numTrain(i),j));
    end
end

% In our paper, pdz_t is assigned as the predicted results by supervised
% classifiers
% The initialization of the target-domain label
pdz_t = zeros(sum(numTest),numC);
flag = 1;
if flag == 1
    w_models = [];
    for i = 1:numSource
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTrain(t);
            end
        end
        if i == 1
            pos = 0;
        end
        TempTrainX = TrainX(:,pos+1:pos+numTrain(i));
        TempTrainY = TrainY(:,pos+1:pos+numTrain(i));
        for v = 1:length(TempTrainY)
            if TempTrainY(v) > 1
                TempTrainY(v) = -1;
            end
        end

        TempTrainXY = scale_cols(TempTrainX,TempTrainY);
        fprintf('.....................................\n');
        w00 = zeros(size(TempTrainXY,1),1);
        lambda = exp(linspace(-0.5,6,20));
        wbest = [];
        f1max = -inf;
        for j = 1:length(lambda)
            w_0 = train_cg(TempTrainXY,w00,lambda(j));
            f1 = logProb(TempTrainXY,w_0);
            if f1 > f1max
                f1max = f1;
                wbest = w_0;
                %se_lambda = lambda(j);
            end
        end
        w_models = [w_models wbest];
        clear TempTrainX;
        clear TempTrainY;
        clear TempTrainXY;
    end
%     csvwrite(strcat('lg_models/','model_lg.model'),w_models);
end

TempGt = zeros(size(pdz_t));
% w_models = csvread(strcat('lg_models/','model_lg.model'));
for i = 1:numTarget
    pos = 0;
    if i > 1
        for t = 1:i-1
            pos = pos + numTest(t);
        end
    end
    if i == 1
        pos = 0;
    end
    TempTestX = TestX(:,pos+1:pos+numTest(i));
    for j = 1:numSource
        wbest = w_models(:,j);
        ptemp = 1./(1 + exp(-wbest'*TempTestX));
    end
    TempGt(pos+1:pos+numTest(i),:) = TempGt(pos+1:pos+numTest(i),:) + [(ptemp'+0.5)/2 ((1-ptemp)'+0.5)/2];
    clear TempTestX;
end
TempGt = TempGt/numSource;
pdz_t = TempGt; % not yet normalize
%% The initialization pyz
pyz_a = ones(numK_1,numC)/numK_1;
pyz_b = ones(numK_2,numC)/numK_2;
pyz_c_s = ones(numK_3,numC*numSource)/numK_3;
pyz_c_t = ones(numK_3,numC*numTarget)/numK_3;

%% the initial accuracy
iter_results = [];
for i = 1:numTarget
    pos = 0;
    if i > 1
        for t = 1:i-1
            pos = pos + numTest(t);
        end
    end
    if i == 1
        pos = 0;
    end
    pzd = pdz_t(pos+1:pos+numTest(i),:);
    nCorrect = 0;
    for j = 1:size(pzd,1)
        [va vi] = max(pzd(j,:));
        if labelset(vi) == TestY(pos+j)
            nCorrect = nCorrect + 1;
        end
    end
    iter_results(1,i+1) = nCorrect/(numTest(i));
    iter_results(1,1) = 0;
end

% the normalization of pdz_t
for i = 1:numTarget
    pos = 0;
    if i > 1
        for t = 1:i-1
            pos = pos + numTest(t);
        end
    end
    if i == 1
        pos = 0;
    end
    for j = 1:numC
        pdz_t(pos+1:pos+numTest(i),j) = pdz_t(pos+1:pos+numTest(i),j)/sum(pdz_t(pos+1:pos+numTest(i),j));
    end
end

pzr_s = ones(1,numC*numSource)/numC;
pzr_t = ones(1,numC*numTarget)/numC;
pr = ones(1,numSource+numTarget)/(numSource+numTarget);

for i = 1:numSource
    pos = 0;
    if i > 1
        for t = 1:i-1
            pos = pos + numTrain(t);
        end
    end
    if i == 1
        pos = 0;
    end
    pr(i) = sum(sum(TrainX(:,pos+1:pos+numTrain(i))));
end
for i = 1:numTarget
    pos = 0;
    if i > 1
        for t = 1:i-1
            pos = pos + numTest(t);
        end
    end
    if i == 1
        pos = 0;
    end
    pr(numSource+i) = sum(sum(TestX(:,pos+1:pos+numTest(i))));
end
pr = pr/sum(pr);
%% stepLen
stepLen = 1000;
% fprintf('the 0 iteration,the value of objective is %g\n',fvalue);
%% Start to interate
% update all variables
% pwy_a; pwy_b_s; pwy_b_t; pwy_c_s; pwy_c_t; pdz_s; pdz_t;
% pyz_a; pyz_b; pyz_c_s; pyz_c_t; pzr_s; pzr_t; pr;
for iterID = 1:numIteration

    % update pwy_a; pwy_b_s; pwy_b_t; pwy_c_s; pwy_c_t;
    temp_pwy_a = zeros(size(pwy_a));
    temp_pwy_b_s = [];
    temp_pwy_b_t = [];
    temp_pwy_c_s = [];
    temp_pwy_c_t = [];
    % update pdz_t;
    temp_pdz_t = [];
    % update pzr_s; pzr_t; pr
    temp_pzr_s = [];
    temp_pzr_t = [];
    temp_pr = [];
    for i = 1:numSource
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTrain(t);
            end
        end
        if i == 1
            pos = 0;
        end
        A = pyz_a;
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_s(1,j+(i-1)*numC);
        end

        tempsum2 = pwy_a*A*pdz_s(pos+1:pos+numTrain(i),:)';
        tempsum2 = tempsum2*pr(i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        temp_pwy_a = temp_pwy_a + pwy_a.*(MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,pdz_s(pos+1:pos+numTrain(i),:),stepLen)*A'*pr(i));
        I = sum(MatrixProduce(pwy_a',MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,pdz_s(pos+1:pos+numTrain(i),:),stepLen),stepLen).*pyz_a)*pr(i);

        %%-----------------%%
        A = pyz_b;
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_s(1,j+(i-1)*numC);
        end

        tempsum2 = pwy_b_s(:,(i-1)*numK_2+1:i*numK_2)*A*pdz_s(pos+1:pos+numTrain(i),:)';
        tempsum2 = tempsum2*pr(i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        B = pwy_b_s(:,(i-1)*numK_2+1:i*numK_2).*(MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,pdz_s(pos+1:pos+numTrain(i),:),stepLen)*A'*pr(i));
        temp_pwy_b_s = [temp_pwy_b_s B];
        I = I + sum(MatrixProduce(pwy_b_s(:,(i-1)*numK_2+1:i*numK_2)',MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,pdz_s(pos+1:pos+numTrain(i),:),stepLen),stepLen).*pyz_b)*pr(i);

        A = pyz_c_s(:,(i-1)*numC+1:i*numC);
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_s(1,j+(i-1)*numC);
        end

        tempsum2 = pwy_c_s(:,(i-1)*numK_3+1:i*numK_3)*A*pdz_s(pos+1:pos+numTrain(i),:)';
        tempsum2 = tempsum2*pr(i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        B = pwy_c_s(:,(i-1)*numK_3+1:i*numK_3).*(MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,pdz_s(pos+1:pos+numTrain(i),:),stepLen)*A'*pr(i));
        temp_pwy_c_s = [temp_pwy_c_s B];
        I = I + sum(MatrixProduce(pwy_c_s(:,(i-1)*numK_3+1:i*numK_3)',MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,pdz_s(pos+1:pos+numTrain(i),:),stepLen),stepLen).*pyz_c_s(:,(i-1)*numC+1:i*numC))*pr(i);
        temp_pzr_s = [temp_pzr_s I];
        temp_pr = [temp_pr sum(I)];
    end
    for i = 1:numTarget
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTest(t);
            end
        end
        if i == 1
            pos = 0;
        end
        A = pyz_a;
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_t(1,j+(i-1)*numC);
        end

        tempsum2 = pwy_a*A*pdz_t(pos+1:pos+numTest(i),:)';
        tempsum2 = tempsum2*pr(numSource+i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        temp_pwy_a = temp_pwy_a + pwy_a.*(MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,pdz_t(pos+1:pos+numTest(i),:),stepLen)*A'*pr(numSource+i));
        H = (MatrixProduce([TestX(:,pos+1:pos+numTest(i))./tempsum2]',pwy_a,stepLen)*A*pr(numSource+i));
        I = sum(MatrixProduce(pwy_a',MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,pdz_t(pos+1:pos+numTest(i),:),stepLen),stepLen).*pyz_a)*pr(numSource+i);
        %%-----------------%%
        A = pyz_b;
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_t(1,j+(i-1)*numC);
        end

        tempsum2 = pwy_b_t(:,(i-1)*numK_2+1:i*numK_2)*A*pdz_t(pos+1:pos+numTest(i),:)';
        tempsum2 = tempsum2*pr(numSource+i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        B = pwy_b_t(:,(i-1)*numK_2+1:i*numK_2).*(MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,pdz_t(pos+1:pos+numTest(i),:),stepLen)*A'*pr(numSource+i));
        temp_pwy_b_t = [temp_pwy_b_t B];
        H = H + (MatrixProduce([TestX(:,pos+1:pos+numTest(i))./tempsum2]',pwy_b_t(:,(i-1)*numK_2+1:i*numK_2),stepLen)*A*pr(numSource+i));
        I = I + sum(MatrixProduce(pwy_b_t(:,(i-1)*numK_2+1:i*numK_2)',MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,pdz_t(pos+1:pos+numTest(i),:),stepLen),stepLen).*pyz_b)*pr(numSource+i);

        A = pyz_c_t(:,(i-1)*numC+1:i*numC);
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_t(1,j+(i-1)*numC);
        end

        tempsum2 = pwy_c_t(:,(i-1)*numK_3+1:i*numK_3)*A*pdz_t(pos+1:pos+numTest(i),:)';
        tempsum2 = tempsum2*pr(numSource+i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        B = pwy_c_t(:,(i-1)*numK_3+1:i*numK_3).*(MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,pdz_t(pos+1:pos+numTest(i),:),stepLen)*A'*pr(numSource+i));
        temp_pwy_c_t = [temp_pwy_c_t B];
        H = H + (MatrixProduce([TestX(:,pos+1:pos+numTest(i))./tempsum2]',pwy_c_t(:,(i-1)*numK_3+1:i*numK_3),stepLen)*A*pr(numSource+i));
        H = pdz_t(pos+1:pos+numTest(i),:).*H;
        temp_pdz_t = [temp_pdz_t; H];
        I = I + sum(MatrixProduce(pwy_c_t(:,(i-1)*numK_3+1:i*numK_3)',MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,pdz_t(pos+1:pos+numTest(i),:),stepLen),stepLen).*pyz_c_t(:,(i-1)*numC+1:i*numC))*pr(numSource+i);
        temp_pzr_t = [temp_pzr_t I];
        temp_pr = [temp_pr sum(I)];
    end

    % update pyz_a; pyz_b; pyz_c_s; pyz_c_t;
    temp_pyz_a = zeros(size(pyz_a));
    temp_pyz_b = zeros(size(pyz_b));
    temp_pyz_c_s = [];
    temp_pyz_c_t = [];
    for i = 1:numSource
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTrain(t);
            end
        end
        if i == 1
            pos = 0;
        end
        A = pdz_s(pos+1:pos+numTrain(i),:);
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_s(1,j+(i-1)*numC);
        end

        tempsum2 = pwy_a*pyz_a*A';
        tempsum2 = tempsum2*pr(i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        temp_pyz_a = temp_pyz_a + pyz_a.*(MatrixProduce(pwy_a',MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,A,stepLen),stepLen)*pr(i));

        %%-----------------%%
        tempsum2 = pwy_b_s(:,(i-1)*numK_2+1:i*numK_2)*pyz_b*A';
        tempsum2 = tempsum2*pr(i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        temp_pyz_b = temp_pyz_b + pyz_b.*(MatrixProduce(pwy_b_s(:,(i-1)*numK_2+1:i*numK_2)',MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,A,stepLen),stepLen)*pr(i));

        tempsum2 = pwy_c_s(:,(i-1)*numK_3+1:i*numK_3)*pyz_c_s(:,(i-1)*numC+1:i*numC)*A';
        tempsum2 = tempsum2*pr(i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        B = pyz_c_s(:,(i-1)*numC+1:i*numC).*(MatrixProduce(pwy_c_s(:,(i-1)*numK_3+1:i*numK_3)',MatrixProduce(TrainX(:,pos+1:pos+numTrain(i))./tempsum2,A,stepLen),stepLen)*pr(i));
        temp_pyz_c_s = [temp_pyz_c_s B];
    end
    for i = 1:numTarget
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTest(t);
            end
        end
        if i == 1
            pos = 0;
        end
        A = pdz_t(pos+1:pos+numTest(i),:);
        for j = 1:numC
            A(:,j) = A(:,j)*pzr_t(1,j+(i-1)*numC);
        end
        tempsum2 = pwy_a*pyz_a*A';
        tempsum2 = tempsum2*pr(numSource+i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        temp_pyz_a = temp_pyz_a + pyz_a.*(MatrixProduce(pwy_a',MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,A,stepLen),stepLen)*pr(numSource+i));

        %%-----------------%%
        tempsum2 = pwy_b_t(:,(i-1)*numK_2+1:i*numK_2)*pyz_b*A';
        tempsum2 = tempsum2*pr(numSource+i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        temp_pyz_b = temp_pyz_b + pyz_b.*(MatrixProduce(pwy_b_t(:,(i-1)*numK_2+1:i*numK_2)',MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,A,stepLen),stepLen)*pr(numSource+i));

        tempsum2 = pwy_c_t(:,(i-1)*numK_3+1:i*numK_3)*pyz_c_t(:,(i-1)*numC+1:i*numC)*A';
        tempsum2 = tempsum2*pr(numSource+i);
        [xs ys] = find(tempsum2 < 10^(-20));
        for q = 1:size(xs,1)
            tempsum2(xs(q,1),ys(q,1)) = 1;
        end
        B = pyz_c_t(:,(i-1)*numC+1:i*numC).*(MatrixProduce(pwy_c_t(:,(i-1)*numK_3+1:i*numK_3)',MatrixProduce(TestX(:,pos+1:pos+numTest(i))./tempsum2,A,stepLen),stepLen)*pr(numSource+i));
        temp_pyz_c_t = [temp_pyz_c_t B];
    end

    % normalize all the variables
    pwy_a = temp_pwy_a;
    pwy_b_s = temp_pwy_b_s;
    pwy_b_t = temp_pwy_b_t;
    pwy_c_s = temp_pwy_c_s;
    pwy_c_t = temp_pwy_c_t;
    % pdz_s = temp_pdz_s;
    pdz_t = temp_pdz_t;
    pyz_a = temp_pyz_a;
    pyz_b = temp_pyz_b;
    pyz_c_s = temp_pyz_c_s;
    pyz_c_t = temp_pyz_c_t;
    pzr_s = temp_pzr_s;
    pzr_t = temp_pzr_t;
    pr = temp_pr;

    for t = 1:numK_1
        pwy_a(:,t) = pwy_a(:,t)/sum(pwy_a(:,t));
    end
    for t = 1:numC
        pyz_a(:,t) = pyz_a(:,t)/sum(pyz_a(:,t));
        pyz_b(:,t) = pyz_b(:,t)/sum(pyz_b(:,t));
    end
    pr = pr/sum(pr);
    for t = 1:numK_2*numSource
        pwy_b_s(:,t) = pwy_b_s(:,t)/sum(pwy_b_s(:,t));
    end
    for t = 1:numK_3*numSource
        pwy_c_s(:,t) = pwy_c_s(:,t)/sum(pwy_c_s(:,t));
    end
    for t = 1:numC*numSource
        pyz_c_s(:,t) = pyz_c_s(:,t)/sum(pyz_c_s(:,t));
    end
    for t = 1:numK_2*numTarget
        pwy_b_t(:,t) = pwy_b_t(:,t)/sum(pwy_b_t(:,t));
    end
    for t = 1:numK_3*numTarget
        pwy_c_t(:,t) = pwy_c_t(:,t)/sum(pwy_c_t(:,t));
    end
    for t = 1:numC*numTarget
        pyz_c_t(:,t) = pyz_c_t(:,t)/sum(pyz_c_t(:,t));
    end
    for i = 1:numTarget
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTest(t);
            end
        end
        if i == 1
            pos = 0;
        end
        for t = 1:numC
            pdz_t(pos+1:pos+numTest(i),t) = pdz_t(pos+1:pos+numTest(i),t)/sum(pdz_t(pos+1:pos+numTest(i),t));
        end
    end
    for i = 1:numSource
        pzr_s(1,(i-1)*numC+1:i*numC) = pzr_s(1,(i-1)*numC+1:i*numC)/sum(pzr_s(1,(i-1)*numC+1:i*numC));
    end
    for i = 1:numTarget
        pzr_t(1,(i-1)*numC+1:i*numC) = pzr_t(1,(i-1)*numC+1:i*numC)/sum(pzr_t(1,(i-1)*numC+1:i*numC));
    end

    %%%%%% The output results
    pzd_t = [];
    for i = 1:numTarget
        A = pdz_t(pos+1:pos+numTest(i),:);
        for t = 1:numC
            A(:,t) = A(:,t)*pzr_t(1,(i-1)*numC+t);
        end
        A = A*pr(numSource+i);
        pzd_t = [pzd_t; A];
    end
    iter_results(iterID+1,1) = iterID;
    for i = 1:numTarget
        pos = 0;
        if i > 1
            for t = 1:i-1
                pos = pos + numTest(t);
            end
        end
        if i == 1
            pos = 0;
        end
        pzd = pzd_t(pos+1:pos+numTest(i),:);
        nCorrect = 0;
        for j = 1:size(pzd,1)
            [va vi] = max(pzd(j,:));
            if labelset(vi) == TestY(pos+j)
                nCorrect = nCorrect + 1;
            end
        end
        iter_results(iterID+1,i+1) = nCorrect/(numTest(i));
    end
end
%% output
Results = iter_results
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [Tc,dT] = nnInterMex(T,omega,x,varargin)
%
% This wrapper runs a CPP file which codes the same scheme as nnInter
% but is hopefully faster, see nnInter for details
%==============================================================================

function [Tc,dT] = nnInterMex(T,omega,x,varargin)

Tc = mfilename('fullpath'); dT = [];
if nargin == 0
  testOneImgModel(mfilename);
    return
elseif nargin == 1 && isempty(T)
    return
end

% flag for computing the derivative
doDerivative = (nargout>1);
matrixFree   = 0;
boundary     = 0; % boundary condition (0: zero padding; 1: replicate)

for k=1:2:length(varargin) % overwrite default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end

dim = length(omega)/2;
m   = size(T);         if dim == 1, m = numel(T); end;
n   = length(x)/dim;

% call CPP subroutine
try
    Tc = nnInterMexC(double(T(:)),omega,m,x(:),boundary==1);
catch err
    FAIRerror(err);
end
if doDerivative
    if matrixFree, dT = zeros(n,dim); else dT = sparse(n,dim*n); end
end
%==============================================================================


text:
function len = length(f)
%LENGTH   Length of a SINGFUN.
%   LENGTH(F) is the length of the smooth part contained in F.
%
% See also SIZE.

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

% The length of a SINGFUN object is the length of its smooth part:
len = length(f.smoothPart);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function scanOut3M = interpScanToScan(baseScan,movScan,addToPlanCFlag,planC)
% function scanOut3M = interpScanToScan(baseScan,movScan,addToPlanCFlag,planC)
%
% This function interpolates movScan to baseScan grid by taking into
% account DICOM imageOrientation. baseScan and movScan must have sam eframe
% of reference for interpolation to be valid.
%
% baseScan = 1;
% movScan = 2;
% addToPlanCFlag = true;
% global planC % or load from file
% scanOut3M = interpScanToScan(baseScan,movScan,addToPlanCFlag,planC);
%
% APA, 3/7/2023

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

% Check FOR
baseFOR = planC{indexS.scan}(baseScan).scanInfo(1).frameOfReferenceUID;
movFOR = planC{indexS.scan}(movScan).scanInfo(1).frameOfReferenceUID;

if ~isequal(baseFOR,movFOR)
    scanOut3M = [];
    return;
end

% Get physical (x,y,z) coordinates of baseScan in DICOM LPS
basePositionMatrix = planC{indexS.scan}(baseScan).Image2PhysicalTransM;
baseSizV = size(planC{indexS.scan}(baseScan).scanArray);
[rowM,colM] = meshgrid(1:baseSizV(1), 1:baseSizV(2));
rowColM = [rowM(:),colM(:)] - 1;
onesV = ones(size(rowColM,1),1);
numSlcVoxels = length(onesV);
xyzM = zeros(numSlcVoxels*baseSizV(3),4);
for slc = 1:baseSizV(3)
    rowColSlcM = [rowColM,onesV*(slc-1),onesV];
    xyzM(1+(slc-1)*numSlcVoxels:slc*numSlcVoxels,1:4) = (basePositionMatrix * rowColSlcM')';
end

% Get moving scan voxel indices (i,j,k) corresponding to (x,y,z) coordinates of the base scan
movSizV = size(planC{indexS.scan}(movScan).scanArray);
movPositionMatrix = planC{indexS.scan}(movScan).Image2PhysicalTransM;
ijkM = movPositionMatrix \ xyzM'; % col,row,slc
ijkM = round(ijkM);
nanV = sum(ijkM(1:3,:) < 0) > 0;
nanV = nanV | ijkM(2,:) > movSizV(1)-1 | ijkM(1,:) > movSizV(2)-1 | ijkM(3,:) > movSizV(3)-1;
indV = sub2ind(movSizV,ijkM(2,~nanV)+1,ijkM(1,~nanV)+1,ijkM(3,~nanV)+1);

% Get moving scan values at (i,j,k)
scanOut3M = zeros(baseSizV,'single');
movScan3M = single(planC{indexS.scan}(movScan).scanArray) - ...
    planC{indexS.scan}(movScan).scanInfo(1).CTOffset;
% slice order in cerr is opposite to dicom. (dicom 1st slice is inf, cerr 1st slice is sup)
% row, col order is same in cerr and dicom. i.e. (1st row, 1st col) is at the top left
% of the slice.
% Hence, flip moving scan along slice dim.
movScan3M = flip(movScan3M,3);
scanOut3M(~nanV) = movScan3M(indV);
scanOut3M = flip(scanOut3M,3);

% Add scan to planC
if addToPlanCFlag
    [xVals, yVals, zVals] = getScanXYZVals(planC{indexS.scan}(baseScan));
    deltaXYZv = getScanXYZSpacing(baseScan,planC);
    zV = zVals(:);
    regParamsS.horizontalGridInterval = deltaXYZv(1);
    regParamsS.verticalGridInterval = deltaXYZv(2);
    regParamsS.coord1OFFirstPoint = xVals(1);
    regParamsS.coord2OFFirstPoint   = yVals(end);

    regParamsS.zValues  = zV;
    regParamsS.sliceThickness = [planC{indexS.scan}(baseScan).scanInfo(:).sliceThickness];

    assocTextureUID = '';
    scanType = [planC{indexS.scan}(movScan).scanType,'_interp_',...
        planC{indexS.scan}(baseScan).scanType];

    planC = scan2CERR(scanOut3M,scanType,'Passed',regParamsS,assocTextureUID,planC);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [MAT,VARNAME,tag_contents]=xml2mat(XML)

%XML2MAT converts XML string into matlab structure variable
%Syntax: [MAT,VARNAME,tag_contents]=xml2mat(XML)
%Description
% XML is an XML formated string using rules compliant with
%     the proceedure implemented in MAT2XML.
%     It can also be a file name with MbML text.
% VARNAME is the variable name.
%     If recovering the MAT variable with the original name is desired than
%     the followinf line will do the trick:
%
%          [MAT,VARNAME]=xml2mat(XML);eval([VARNAME,'=MAT'])
%
% See Also: mat2xml
% Jonas Almeida, 20 Aug 2002, XML4MAT Tbox

if strncmp(XML,'%60;',4)  % XML provided as encoded XML string
    XML=spcharout(XML);
end

% is this a xml string or xml filename ?
if XML(1)~='<'
    XML=strrep(file2str(XML),'''','''''');
    % Remove non-content lines if they exist
    XML=regexprep(XML,'<[?!].*?>','');
end

%Analise XML line
tag_ini=find(XML=='<');
tag_end=find(XML=='>');
n=length(tag_ini); % number of enclosed tag structures

% extract tag_names properties and contents
if n>0
for i=1:n
    tag_close(i)=(XML(tag_ini(i)+1)=='/'); % 1 for closing contents and 0 for opening
    tag_contents{i,1}=XML(tag_ini(i)+1+tag_close(i):tag_end(i)-1);  % first column contains names
end
tag_path=[0]; % first name is root name
to_do=zeros(n,1); % 1 needs doing
to_do(1)=1;
tag_open=~tag_close;
i=1;tag_contents{i,2}=xml2whos(tag_contents{i,1});tag_contents{i,2}.fields=[];tag_contents{i,3}=tag_path;tag_path=[tag_path,i];
for i=2:n
    tag_contents{i,2}=xml2whos(tag_contents{i,1});tag_contents{i,2}.fields=[]; % second column contains WHO properties
    if tag_open(i)
        tag_contents{i,3}=tag_path;  % third column contains structre path
        tag_path=[tag_path,i];
        to_do(i)=1; % do this one
        to_do(tag_contents{i,3}+(tag_contents{i,3}==0))=0; % do host later
        tag_contents{tag_path(end-1),2}.fields=[tag_contents{tag_path(end-1),2}.fields,i];
    else
        tag_path(end)=[]; % move back one level
    end
end

% RECOVER DATA
todo_list=find(to_do==1)';do_i=1;
while do_i<=length(todo_list);
    i=todo_list(do_i);
    % for each case extract value
    tag_contents{i,4}=XML(tag_end(i)+1:tag_ini(i+1)-1);  % 4th column contains value as string
    w=tag_contents{i,2};
    % recover number format if appropriated
    if ~isfield(w,'class')
        w.class='char';
        tag_contents{i,4}=spcharout(tag_contents{i,4});
        w.size=size(tag_contents{i,4}); %correct size tag as well
        %tag_contents{i,2}=w;
    elseif strcmp(w.class,'char')
        tag_contents{i,4}=spcharout(tag_contents{i,4});
        w.size=size(tag_contents{i,4}); %correct size tag as well
        %tag_contents{i,2}=w;
    elseif strcmp(w.class,'struct')&(~isfield(w,'size'))
        n_fields=length(w.fields);
        for i=1:n_fields
            field_names{i}=tag_contents{w.fields(i),2}.name;
        end
        unique_names=unique(field_names);
        n_unique=length(unique_names);
        n_certo=(n_fields/n_unique);
        w.size=[1 n_certo];
    else % it is a numeric type, say "double" or "single"
        tag_contents{i,4}=str2num(tag_contents{i,4});
        if ~isfield(w,'size');
            w.size=size(tag_contents{i,4});
        end
    end
    tag_contents{i,2}=w;
    if (length(w.size>2)|(w.size(1)>1))
        iis='i1';for ii=2:length(w.size);iis=[iis,',i',num2str(ii)];end
        nn=prod(w.size); %number of elements
        eval(['[',iis,']=ind2sub(w.size,[1:nn]);']); % generation of indexes
        iis='i1(ind)';for ii=2:length(w.size);iis=[iis,',i',num2str(ii),'(ind)'];end % indexes of indexes
        for ind=1:nn
            eval(['valor(',iis,')=tag_contents{i,4}(ind);'])
        end
        if exist('valor')==1;tag_contents{i,4}=valor;clear valor;end
    end
    do_i=do_i+1;
end

% RECOVER STRUCTURE
j=1;tags=find(tag_open);
eval_i={'MAT'};to_eval={};
while j<=length(tags)
    i=tags(j);
    [to_eval,eval_i,j]=tag2eval(tag_contents,to_eval,eval_i,j,tags);
    %w=tag_contents{i,2}
    j=j+1;
end

for i=1:length(to_eval)
    %disp(to_eval{i})
    eval(to_eval{i})
end

%MAT.tag_contents=tag_contents;
VARNAME=tag_contents{1,2}.name;
else
    MAT=XML;
    VARNAME=[];
    tag_contents=[];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ics_competing_talker(filename, clfile, t_outfile, m_outfile,thrd)

% filename - mixture filename
% clfile - clean target filename
% t_outfile - output file: Target
% m_outfile - output file: Competing talker
% thrd - SNR threshold in dB
%
%   Copyright (c) 2011 by Philipos C. Loizou

[x,Srate,nb] = wavread(filename); % mixture signal
cl = wavread(clfile); % clean speech signal

% =============== Initialize variables ===============
%
len = floor(20*Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC = 50; % window overlap in percent of frame size
len1 = floor(len*PERC/100);
len2 = len-len1;
win = hanning(len);
nFFT = len;

%--- allocate memory and initialize various variables

x_old = zeros(len1,1);
Nframes = floor(length(x)/len2)-1;
xfinal = zeros(Nframes*len2,1); % for target signal
x_old_m = zeros(len1,1);
xfinal_m = zeros(Nframes*len2,1); % for masker signal

%% ===========  Start Processing    ============
k = 1;
thrd = 10^(thrd/10);
m=1;
nFFT2=floor(nFFT/2);

for n = 1:Nframes
    insign = win.*x(k:k+len-1);  % mixture
    cl_sign = win.*cl(k:k+len-1);  % clean target
    tn = insign - cl_sign;   % masker signal

    %--- Take fourier transform of  frames

    spec = fft(insign,nFFT);
    sig = abs(spec); % compute the magnitude
    sig2 = sig.^2;

    cl_spec = fft(cl_sign,nFFT);
    cl_sig = abs(cl_spec); % compute the magnitude
    cl_sig2 = cl_sig.^2;

    tn_spec = fft(tn,nFFT);
    tn_sig = abs(tn_spec); % compute the magnitude
    tn_sig2 = tn_sig.^2;

    %% Compute IDEAL BINARY MASK
    %
    ksi_IEC = cl_sig2./tn_sig2;  % Ideal ksi

    hw_true = zeros(len,1);
    indp = find(ksi_IEC>= thrd); % ideal IBM
    hw_true(indp)=1;  % true ksi
    % ideal IBM to recover the target signal

    hw_true_m = zeros(len,1);
    indp = find(ksi_IEC< thrd); % ideal IBM
    hw_true_m(indp)= 1;  % true ksi
    % complementary ideal IBM to recover the mask signal

    x_hw_m = hw_true_m;
    x_hw = hw_true;

    % Synthesize target signal
    xi_w = ifft( x_hw .* spec);    %
    xi_w = real( xi_w);
    xfinal(k:k+ len1-1) = x_old+ xi_w(1:len1);
    x_old = xi_w( len2+ 1: len);

    % Synthesize masker signal
    xi_w_m = ifft( x_hw_m .* spec);
    xi_w_m = real( xi_w_m);
    xfinal_m( k: k+ len1-1) = x_old_m+ xi_w_m( 1:len1);
    x_old_m = xi_w_m( len2+ 1: len);

    k = k+len2;
end

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

if max(abs(xfinal))>1
    xfinal = xfinal*0.6/max(abs(xfinal));
    fprintf('Max amplitude exceeded 32768 for file %s\n',filename);
end

wavwrite(xfinal, Srate, 16, t_outfile);

if max(abs(xfinal_m))>1
    xfinal_m = xfinal_m* 0.6/max(abs(xfinal_m));
    fprintf('Max amplitude exceeded 32768 for file %s\n',filename);
end

wavwrite(xfinal_m, Srate, 16, m_outfile);

author: jtkim-kaist

repo: Speech-enhancement


text:
function flood_fill(xc,yc,x,y,new,old,width,height)
% flood fills area of matrix
% (xc,yc) = start point of flood fill
% (x,y) = current test location
% new = new value to be filled
% old = old value to replace
%
% example use:
% global fill
% width = 20;
% height = 20;
% fill = zeros(width,height);
% new = 1;
% old = 0;
% xc = 10;
% yc = 10;
% flood_fill(xc,yc,xc,yc,new,old,width,height);
%
% Results:
% The flood fill will alter all values of 0 to one flood filling from a
% start point (xc,yc)
% If the flood fill starts in an enclosed space it will fill up to the
% boundary
%
% Author: James Goodwin, June 2003
% email: JamesRichardGoodwin@lyocs.com

global fill;
global call;

call = call + 1;

if (x<2 || x>=width)
    x = xc;
end

if (y<2 || y>=height)
    y = yc;
end

if (fill(y,x) == old)
    fill(y,x) = new;
    flood_fill(xc,yc,x+1,y,new,old,width,height);
    flood_fill(xc,yc,x,y+1,new,old,width,height);
    flood_fill(xc,yc,x-1,y,new,old,width,height);
    flood_fill(xc,yc,x,y-1,new,old,width,height);
end

author: Sable

repo: mcbench-benchmarks


text:
%MDL_IR140 Create model of the ABB IRB 140 manipulator
%
% MDL_IRB140_MOD is a script that creates the workspace variable irb140 which
% describes the kinematic characteristics of an ABB IRB 140 manipulator
% using modified DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%
% Reference::
% - ABB IRB 140 data sheet
% - "The modeling of a six degree-of-freedom industrial robot for
%   the purpose of efficient path planning",
%   Master of Science Thesis, Penn State U, May 2009,
%   Tyler Carter
%
% See also mdl_irb140, mdl_puma560, mdl_stanford, mdl_twolink, SerialLink.
%
% Notes::
% - SI units of metres are used.
% - The tool frame is in the centre of the tool flange.
% - Zero angle configuration has the upper arm vertical and lower arm
%   horizontal.

% MODEL: ABB, IRB140, 6DOF, modified_DH

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

clear L

% joint angle limits from
% A combined optimization method for solving the inverse kinematics problem...
% Wang & Chen
% IEEE Trans. RA 7(4) 1991 pp 489-
L(1) = Revolute('d', 0.352, 'a', 0, 'alpha', 0, 'offset', 0, 'modified');
L(2) = Revolute('d', 0, 'a', 0.070, 'alpha', pi/2, 'offset', 0, 'modified');
L(3) = Revolute('d', 0, 'a', 0.360, 'alpha', 0, 'offset', 0, 'modified');
L(4) = Revolute('d', 0.380, 'a', 0, 'alpha', pi/2, 'offset', 0, 'modified');
L(5) = Revolute('d', 0, 'a', 0, 'alpha', -pi/2, 'offset', 0, 'modified');
L(6) = Revolute('d', 0, 'a', 0, 'alpha', pi/2, 'offset', 0, 'modified');

L(1).m = 34655.36e-3;
L(1).r = [27.87 43.12 -89.03]*1e-3;
L(1).I = [
    512052539.74 1361335.88 51305020.72
    1361335.88 464074688.59 70335556.04
    51305020.72 70335556.04 462745526.12]*1e-9;

L(2).m = 15994.59e-3;
L(2).r = [ 198.29 9.73 92.43]*1e03;
L(2).I = [
    94817914.40 -3859712.77 37932017.01
    -3859712.77 328604163.24 -1088970.86
    37932017.01 -1088970.86 277463004.88]*1e-9;

L(3).m = 20862.05e-3;
L(3).r = [ -4.56 -79.96 -5.86];
L(3).I = [
    500060915.95 -1863252.17 934875.78
    -1863252.17 75152670.69 -15204130.09
    934875.78 -15204130.09 515424754.34]*1e-9;

irb140 = SerialLink(L, 'name', 'IRB 140', ...
    'manufacturer', 'ABB', 'comment', 'modified DH');

%
% some useful poses
%
qz = [0 0 0 0 0 0]; % zero angles, L shaped pose

clear L

author: petercorke

repo: robotics-toolbox-matlab


text:
function [m, A, Eigenfaces] = EigenfaceCore(Training_Data)

%----------------------------Calculate the mean image ------------------------
% ---------------------compute the covariance matrix --------------------------
m = mean(Training_Data,2);
Train_Number = size(Training_Data,2);
temp_m = [];
for i = 1 : Train_Number
    temp_m = [temp_m m];
end
A = double(Training_Data) - temp_m;

%-------------------use the "svd" function to compute the eigenvectors
%--------------------and eigenvalues of the covariance matrix.
disp('Computing...Wait a second please')
L = A'*A;
%size(L)
[V D] = eig(L);

%-----------------------------Sort and eliminate eigenvalues---------------
Eig_vec = [];
for i = 1 : size(V,2)
    if( D(i,i)>1000 ) % Set Threshold value whatever you like
        Eig_vec = [Eig_vec V(:,i)];
    end
end

Eigenfaces = A * Eig_vec;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function varargout = circumCenter(a, b, c)
%CIRCUMCENTER  Circumcenter of three points
%
%   CC = circumCenter(P1, P2, P3)
%
%   Example
%     A = [10 10]; B = [30 10]; C = [10 20];
%     circumCenter(A, B, C)
%     ans =
%         20    15
%
%     % works also for multiple input points
%     circumCenter([A;A;A], [B;B;B], [C;C;C])
%     ans =
%         20    15
%         20    15
%         20    15
%
%
%   See also
%     points2d, circles2d, circumCircle, centroid
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-12-09,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% pre-compute some terms
ah = sum(a .^ 2, 2);
bh = sum(b .^ 2, 2);
ch = sum(c .^ 2, 2);

dab = a - b;
dbc = b - c;
dca = c - a;

% common denominator
D  = .5 ./ (a(:,1) .* dbc(:,2) + b(:,1) .* dca(:,2) + c(:,1) .* dab(:,2));

% center coordinates
xc =  (ah .* dbc(:,2) + bh .* dca(:,2) + ch .* dab(:,2) ) .* D;
yc = -(ah .* dbc(:,1) + bh .* dca(:,1) + ch .* dab(:,1) ) .* D;

if nargout <= 1
    varargout = {[xc yc]};
else
    varargout = {xc, yc};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
classdef ParadigmMKLFBCSP < ParadigmBase
    % Multiple Kernel Learning Filter-Bank Common Spatial Patterns (mklFBCSP)
    %
    % This paradigm implements a generalization of mklCSP [1] to multiple frequency bands as in
    % FBCSP [2] and to multiple time windows. See also ParadigmFBCSP and ParadigmMKLCSP.
    %
    % References:
    % [1] Samek, W., Binder, A., & Muller, K. R.
    %     "Multiple kernel learning for brain-computer interfacing."
    %     In Engineering in Medicine and Biology Society (EMBC) pp. 7048-7051 (2013)
    % [2] Kai K. Ang, Zhang Y. Chin, Haihong Zhang, Cuntai Guan,
    %     "Filter Bank Common Spatial Pattern (FBCSP) in Brain-Computer Interface"
    %     In 2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence) (June 2008), pp. 2390-2397.
    %
    % Name:
    %   Multiple Kernel Learning Filter-Bank Common Spatial Patterns
    %
    %                            Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                            2014-02-05

    methods

        function defaults = preprocessing_defaults(self)
            defaults = {'EpochExtraction',[0.5 3.5],'Resampling',200};
        end

        function defaults = machine_learning_defaults(self)
            % set up the default parameters for machine learning
            defaults = {'logreg', 'Variant','lars'};
        end

        function model = calibrate(self,varargin)
            % calibrate an mklCSP model from a corpus of training sets
            args = arg_define(varargin, ...
                arg_norep({'collection','Collection'}), ...
                arg_norep({'goal_identifier','GoalIdentifier'}), ...
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 10000]),'Number of CSP patterns (times two).','cat','Feature Extraction','type','expression','shape','row'),...
                arg({'shrinkage','ShrinkageLevel'},0,[0 1],'Shrinkage level. The amount of shrinkage (regularization) to apply during covariance estimation.'), ...
                arg({'freqwnds','FreqWindows'},[0.5 3; 4 7; 8 12; 13 30; 31 42],[0 0.5 200 1000],'Frequency bands of interest. Matrix containing one row for the start and end of each frequency band from which CSP patterns shall be computed. Values in Hz.','cat','Feature Extraction'), ...
                arg({'timewnds','TimeWindows'},[],[],'Time windows of interest. Matrix containing one row for the start and end of each time window from which CSP patterns shall be computed. Values in seconds. If both this and the freqwnds parameter are non-empty, they should have the same number of rows.','cat','Feature Extraction'), ...
                arg({'winfunc','WindowFunction'},'rect',{'barthann','bartlett','blackman','blackmanharris','bohman','cheb','flattop','gauss','hamming','hann','kaiser','nuttall','parzen','rect','taylor','triang','tukey'},'Type of window function. Typical choices are rect (rectangular), hann, gauss, blackman and kaiser.'),...
                arg({'winparam','WindowParameter','param'},[],[],'Parameter of the window function. This is mandatory for cheb, kaiser and tukey and optional for some others.','shape','scalar'), ...
                arg({'verbose','Verbose'},true,[],'Verbose output.'), ...
                arg({'hotpatching','HotPatching'},false,[],'Hot-patch the data. This can be enabled to ensure that a long-running computation survives bad data.'), ...
                arg_sub({'flt','SignalProcessing'}, self.preprocessing_defaults(), @flt_pipeline, 'Signal processing stages. These parameters control filter stages that run on the signal level; they can be enabled, disabled and configured for the given paradigm. The prediction operates on the outputs of this stage.'), ...
                arg_sub({'ml','MachineLearning'},{'Learner',self.machine_learning_defaults()},@ml_train,'Machine learning stage of the paradigm. Operates on the feature vectors that are produced by the feature-extraction stage.'),...
                arg({'arg_dialogsel','ConfigLayout'},self.dialog_layout_defaults(),[],'Parameters displayed in the config dialog. Cell array of parameter names to display (dot-notation allowed); blanks are translated into empty rows in the dialog. Referring to a structure argument lists all parameters of that struture, except if it is a switchable structure - in this case, a pulldown menu with switch options is displayed.','type','cellstr','shape','row'));

            if ~isempty(args.freqwnds) && ~isempty(args.timewnds) && size(args.freqwnds,1) ~= size(args.timewnds,1)
                error('If both time and frequency windows are specified, both arrays must have the same number of rows (together they define the windows in time and frequency).'); end
            if isempty(args.timewnds)
                args.timewnds = zeros(size(args.freqwnds,1),0); end
            if isempty(args.freqwnds)
                args.freqwnds = zeros(size(args.timewnds,1),0); end

            % pre-parse arguments for flt_window and flt_spectrum (for fast subsequent online use)
            for w = 1:max(size(args.freqwnds,1),size(args.timewnds,1))
                time_args{w} = arg_report('vals',@flt_window,{'time',{args.timewnds(w,:),args.winfunc,args.winparam}});
                freq_args{w} = arg_report('vals',@flt_spectrum,{'freq',args.freqwnds(w,:)});
            end

            if args.verbose
                fprintf('Now training model for: %s...\n',hlp_tostring(args.goal_identifier)); end

            % first solve CSP for each subject in the corpus individually and aggregate CSP filters
            filters = [];
            patterns = [];

            % find the unique subjects in the collection
            try
                corpus = [args.collection{:}];
            catch e
                error('The dataset collection must have the same field names for each recording.');
            end
            if ~isfield(corpus,'subject')
                error('The datasets in the collection must each have a .subject field.'); end
            subjects = {corpus.subject};
            if all(cellfun('isclass',subjects,'char'))
                subjects = unique(subjects);
            elseif all(cellfun('isclass',subjects,'double'))
                subjects = unique([subjects{:}]);
            else
                error('The subject identifiers must either be all strings or all doubles');
            end

            if args.verbose
                fprintf('Pre-processing each of %i recordings (%i subjects) in the corpus and solving CSP...\n',length(args.collection),length(subjects)); end

            % remove actual data from corpus so we can micro-cache it
            for s=1:length(corpus)
                if isfield(corpus(s).streams{1},'tracking')
                    corpus(s).streams{1} = corpus(s).streams{1}.tracking.expression; end
            end

            % for each subject...
            for subj=subjects
                % find all recordings that match that subject
                recordings = corpus(cellfun(@(s)isequal(s,subj),{corpus.subject}));
                % calculate FBCSP filters
                [newfilters,newpatterns,chanlocs] = hlp_microcache('filters',@ParadigmMKLFBCSP.filters_for_subject,recordings, args.flt, time_args, freq_args, args.shrinkage, args.patterns);
                % if you get an error here then your data sets had varying number of channels
                filters = [filters newfilters];
                patterns = [patterns newpatterns];
            end
            model.featuremodel = struct('filters',{filters},'patterns',{patterns}, ...
                'n_subjects',length(subjects),'time_args',{time_args},'freq_args',{freq_args}, ...
                'chanlocs',{chanlocs}, 'hotpatching', {args.hotpatching});
            if args.verbose
                fprintf('Preprocessing and extracting features for reference data...\n'); end
            % get the data of the reference subject
            [reference,remaining] = utl_collection_closest(args.collection,args.goal_identifier); %#ok<ASGLU,NASGU>
            % preprocess each recording in the reference collection and concatenate them across epochs into a single set
            for r=1:length(reference)
                refsets{r} = exp_eval_optimized(flt_pipeline('signal',reference{r}.streams{1}, args.flt)); end
            refdata = exp_eval(set_joinepos(refsets{:}));
            % extract features and get target labels
            features = self.feature_extract(refdata,model.featuremodel);
            targets = set_gettarget(refdata);

            if args.hotpatching && size(features,1) < 5
                fprintf('You have too few trials in the data; hot-patching.\n');
                features = features(1+mod(0:4,size(features,1)),:);
                features = features+0.1*randn(size(features));
                targets = 1+mod(0:size(features,1)-1,2);
                targets = targets(:);
            end

            if args.hotpatching && length(targets) ~= size(features,1)
                fprintf('Your # of target markers does not match the # of extracted features; hot-patching.\n');
                if isempty(targets)
                    targets = 1+mod(0:size(features,1)-1,2);
                else
                    targets = targets(1+mod(0:size(features,1)-1,length(targets)));
                end
                targets = targets(:);
            end

            if args.hotpatching && length(unique(targets))==1
                fprintf('Your reference data has only one class; hot-patching the data.\n');
                for ii=1:min(length(targets),max(2,round(length(targets)/10)))
                    targets(ii) = 3-targets(ii); end
            end

            if args.hotpatching && any(~isfinite(features(:)))
                fprintf('Some of your features are non-finite; hot-patching the data.\n');
                tofix = find(~isfinite(features(:)));
                features(tofix) = randn(1,length(tofix));
            end

            if args.verbose
                fprintf('Training predictive model (this may take a while)...\n'); end
            % train classifier, overriding with the correct feature shape (based on the group size)
            if isfield(args.ml.learner,'shape')
                args.ml.learner.shape = [2*args.patterns,length(subjects)]; end
            try
                model.predictivemodel = ml_train('data',{features,targets}, args.ml);
            catch e
                if args.hotpatching && ~isempty(strfind(e.message,'Null probability for class'))
                    fprintf('One of the classes has a probability of 0; hot-patching the data.\n');
                    targets = 1+mod(0:length(targets)-1,2); targets = targets(:);
                    model.predictivemodel = ml_train('data',{features,targets}, args.ml);
                else
                    rethrow(e);
                end
            end
            % set the filter graph based on the reference data
            model.tracking.filter_graph = refsets{end};
            % also store channel locations for model visualization
            model.chanlocs = refdata.chanlocs;
        end

        function predictions = predict(self,bundle,model)
            % extract features
            features = self.feature_extract(bundle.streams{1},model.featuremodel);
            % apply classifier
            predictions = ml_predict(features, model.predictivemodel);
        end

        function features = feature_extract(self,signal,featuremodel)
            try
                N = featuremodel.n_subjects;
                W = length(featuremodel.freq_args);
                F = size(featuremodel.filters,2)/N;
                T = size(signal.data,3);
                features = zeros(T,F*N);
                for w = 1:W
                    % filter data in time & frequency
                    data = exp_eval(flt_spectrum('signal',flt_window('signal',signal,featuremodel.time_args{w}),featuremodel.freq_args{w}));
                    mask = false(1,F); mask((w-1)*(F/W)+(1:(F/W))) = true; mask = repmat(mask,1,N);
                    for t=1:size(signal.data,3)
                        features(t,mask) = sum((data.data(:,:,t)'*featuremodel.filters(:,mask)).^2,1); end
                end
                features = log(features/size(signal.data,2));
            catch e
                if featuremodel.hotpatching
                    fprintf('Trying to prevent error during feature extraction: %s\n',hlp_handleerror(e));
                    fprintf('size(featuremodel.filters): %s\n',hlp_tostring(size(featuremodel.filters)));
                    fprintf('size(signal.data): %s\n',hlp_tostring(size(signal.data)));
                    fprintf('trying to hot-fix the issue...');
                    featuremodel.filters = featuremodel.filters(1+mod(0:size(signal.data,1)-1,size(featuremodel.filters,1)),:);
                    features = self.feature_extract(signal,featuremodel);
                    fprintf('succeeded.\n');
                else
                    rethrow(e);
                end
            end
        end

        function visualize(self,varargin) %#ok<*INUSD>
            % visualize an mklCSP model
            args = arg_define(varargin, ...
                arg_norep({'model','Model'},[],[],'BCI Model to visualize.'), ...
                arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
                arg({'paper','PaperFigure'},false,[],'Use paper-style font sizes. Whether to generate a plot with font sizes etc. adjusted for paper.'));

            f = figure;

            % mask out unused filters
            mask = args.model.predictivemodel.model.w(:)' ~= 0;
            args.model.featuremodel.patterns = args.model.featuremodel.patterns(:,mask);
            args.model.featuremodel.filters = args.model.featuremodel.filters(:,mask);

            % number of plots, and index of pattern per subplot
            np = nnz(mask);
            horz = ceil(sqrt(np));
            vert = ceil(np/horz);

            % get number of pairs, and index of pattern per subplot
            % for each CSP pattern...
            for p=1:np
                subplot(horz,vert,p,'Parent',f);
                if args.patterns
                    topoplot(args.model.featuremodel.patterns(:,p),args.model.chanlocs);
                else
                    topoplot(args.model.featuremodel.filters(:,p),args.model.chanlocs);
                end
                t = title(['CSP Pattern ' num2str(p)]);
                if args.paper
                    set(t,'FontUnits','normalized');
                    set(t,'FontSize',0.1);
                end
            end
        end

        function layout = dialog_layout_defaults(self)
            % define the default configuration dialog layout
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.FIRFilter.Frequencies', ...
                'SignalProcessing.FIRFilter.Type', 'SignalProcessing.EpochExtraction', '', ...
                'PatternPairs', '', 'MachineLearning.Learner'};
        end

    end

    methods (Static)
        function [filters, patterns, chanlocs] = filters_for_subject(recordings, flt, time_args, freq_args, shrinkage, n_patterns)
            % get the CSP filters for a given subject

            % preprocess and concatenate across trials
            preproc = {};
            for r=1:length(recordings)
                preproc{r} = exp_eval_optimized(flt_pipeline('signal',recordings(r).streams{1}, flt)); end
            preproc = exp_eval(set_joinepos(preproc{:}));

            % for each window...
            filters = [];
            patterns = [];
            for w = 1:length(time_args)
                for k=1:2
                    % filter trial subrange in time and frequency
                    classdata = exp_eval(flt_spectrum('signal',flt_window('signal',set_picktrials(preproc,'rank',k),time_args{w}),freq_args{w}));
                    covar{k} = reshape(classdata.data,size(classdata.data,1),[])*reshape(classdata.data,size(classdata.data,1),[])'/(size(classdata.data,2)*size(classdata.data,3)); % cov(reshape(classdata.data,size(classdata.data,1),[])');
                    covar{k}(~isfinite(covar{k})) = 0; %#ok<*AGROW>
                    covar{k} = (1-shrinkage)*covar{k} + shrinkage*eye(size(covar{k}))*trace(covar{k})/length(covar{k});
                end
                try
                    [V,D] = eig(covar{1},covar{1}+covar{2}); %#ok<ASGLU,NASGU>
                    P = inv(V);
                    % if you get an error here then your data sets had varying number of channels
                    filters = [filters V(:,[1:n_patterns end-n_patterns+1:end])];
                    patterns = [patterns P([1:n_patterns end-n_patterns+1:end],:)'];
                catch e
                    fprintf('Got a degenerate FBCSP solution, replacing by identity matrix:%s\n',e.message);
                    n_chans = preproc.nbchan;
                    if ~n_chans
                        % no epochs, need to determine the number of channels in the filter stage prior to epoching
                        raw = exp_eval_optimized(utl_get_argument(utl_find_filter(preproc,'set_makepos'),'signal'));
                        n_chans = raw.nbchan;
                    end
                    filters = [filters eye(n_chans,2*n_patterns)];
                    patterns = [patterns eye(n_chans,2*n_patterns)];
                end
            end
            chanlocs = preproc.chanlocs;
        end
    end
end

% (turn off a few editor warnings because some actual implementations are missing in this file)
%#ok<*INUSD,*STOUT,*MANU>

repo: REST


text:
function prob_test027 ( )

%*****************************************************************************80
%
%% TEST027 tests BURR_MEAN, BURR_VARIANCE, BURR_SAMPLE;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST027\n' );
  fprintf ( 1, '  For the Burr PDF:\n' );
  fprintf ( 1, '  BURR_MEAN computes the mean;\n' );
  fprintf ( 1, '  BURR_VARIANCE computes the variance;\n' );
  fprintf ( 1, '  BURR_SAMPLE samples;\n' );

  a = 1.0;
  b = 2.0;
  c = 3.0;
  d = 2.0;

  check = burr_check ( a, b, c, d );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST027 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = burr_mean ( a, b, c, d );
  variance = burr_variance ( a, b, c, d );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );
  fprintf ( 1, '  PDF parameter C = %14f\n', c );
  fprintf ( 1, '  PDF parameter D = %14f\n', d );
  fprintf ( 1, '  PDF mean =        %14f\n', mean );
  fprintf ( 1, '  PDF variance =    %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = burr_sample ( a, b, c, d, seed );
  end

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

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function res = bbIntersection(rect1, rect2)
res = zeros(1, 4);
res(1) = max(rect1(1), rect2(1));
res(2) = max(rect1(2), rect2(2));

res(3) = min(rect1(1) + rect1(3), rect2(1) + rect2(3)) - res(1);
res(4) = min(rect1(2) + rect1(4), rect2(2) + rect2(4)) - res(2);

end

author: aosokin


text:
clc
clear all
for jj=2:2
for j=73:75
filename=['.\data\3lab' num2str(jj) num2str(j)];
% filename=['.\data\x' num2str(j)];
if(exist(filename)==2)
M=1000;
N=1;
csi_trace = read_bf_file(filename);
csi=zeros(3,30,M-N);
for i=(N+1):M
    csi_entry = csi_trace{i};
    csientry = get_scaled_csi(csi_entry);
    perm = csi_entry.perm;
    for k=1:3
        if perm(k)==1
            csi(1,:,i-N)=csientry(1,perm(k),:);
%             time_csi(1,:,i-N)=ifft(csi(1,:,i-N),point);
%             time_csi_abs(1,:,i-N)=abs(time_csi(1,:,i-N));
        elseif perm(k)==2
            csi(2,:,i-N)=csientry(1,perm(k),:);
%             time_csi(2,:,i-N)=ifft(csi(2,:,i-N),point);
%             time_csi_abs(2,:,i-N)=abs(time_csi(2,:,i-N));
        elseif perm(k)==3
            csi(3,:,i-N)=csientry(1,perm(k),:);
%             time_csi(3,:,i-N)=ifft(csi(3,:,i-N),point);
%             time_csi_abs(3,:,i-N)=abs(time_csi(3,:,i-N));
        end
    end
end
csi=squeeze(csi);
save(filename,'csi');
end
end
end

repo: DeepFi


text:
function colorStr = getStringFromLineSpecColor(lineSpecColor)
    if(strcmpi(lineSpecColor,'r'))
        colorStr = 'Red';

    elseif((isnumeric(lineSpecColor) && all(lineSpecColor == [178,0,255]/255)) || (ischar(lineSpecColor) && strcmpi(lineSpecColor,'m')))
        colorStr = 'Magenta';

	elseif(lineSpecColor == [255,131,0]/255)
        colorStr = 'Orange';

    elseif((isnumeric(lineSpecColor) && all(lineSpecColor == [255,216,0]/255)) || (ischar(lineSpecColor) && strcmpi(lineSpecColor,'y')))
        colorStr = 'Yellow';

    elseif((isnumeric(lineSpecColor) && all(lineSpecColor == [76,220,0]/255)) || (ischar(lineSpecColor) && strcmpi(lineSpecColor,'g')))
        colorStr = 'Green';

    elseif((isnumeric(lineSpecColor) && all(lineSpecColor == [0,210,255]/255)) || (ischar(lineSpecColor) && strcmpi(lineSpecColor,'c')))
        colorStr = 'Cyan';

	elseif(strcmpi(lineSpecColor,'b'))
        colorStr = 'Blue';

	elseif(strcmpi(lineSpecColor,'k'))
        colorStr = 'Black';

	elseif(strcmpi(lineSpecColor,'w'))
        colorStr = 'White';

	elseif(lineSpecColor == [255,0,220]/255)
        colorStr = 'Pink';

	elseif(lineSpecColor == [139,69,19]/255)
        colorStr = 'Brown';

    else
        colorStr = 'Black';
        warning(['Did not recongize line spec color string: ', num2str(lineSpecColor)]);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
load cameraRTxj.mat

addpath /n/fs/vision/www/pvt/
addpath /n/fs/vision/www/pvt/drawCamera/

scale = 0.0001;

drawCameraPly('~/Downloads/cam1.ply',extrinsicsC2WframeTestDevice{1},scale);
drawCameraPly('~/Downloads/cam2.ply',extrinsicsC2WframeTestDevice{2},scale);
drawCameraPly('~/Downloads/cam3.ply',extrinsicsC2WframeTestDevice{3},scale);
drawCameraPly('~/Downloads/cam4.ply',extrinsicsC2WframeTestDevice{4},scale);

author: jianxiongxiao

repo: ProfXkit


text:
function C = cell2ColorOrder(colors);
%
% C = cell2ColorOrder(colors);
%
% Convert a cell array containing color signifiers
% into an N x 3 matrix, useable for color orders for
% axes. (e.g., in setLineColors or mybar).
%
% color signifiers can be 1x3 matrices of [R G B]
% values (from 0-1 for each channel), or else one
% of the following characters:
%
%   'r': red [1 0 0]
%   'g': green [0 1 0]
%   'b': blue [0 0 1]
%   'k': black [0 0 0]
%   'y': yellow [1 1 0]
%   'm': magenta [1 0 1]
%   'c': cyan [0 1 1]
%   'w': white [1 1 1]
%   'e' or 'a': gray [.6 .6 .6]
%
%
% ras, 06/05.
if ~iscell(colors), error('Requires a cell input'); end

colors = unNestCell(colors);
C = zeros(length(colors),3);

for i = 1:length(colors)
    if isnumeric(colors{i}) & length(colors{i})==3
        C(i,:) = colors{i};
    elseif ischar(colors{i})
        ch = colors{i}(1);
        switch ch
            case 'r', C(i,:) = [1 0 0];
            case 'g', C(i,:) = [0 1 0];
            case 'b', C(i,:) = [0 0 1];
            case 'k', C(i,:) = [0 0 0];
            case 'y', C(i,:) = [1 1 0];
            case 'm', C(i,:) = [1 0 1];
            case 'c', C(i,:) = [0 1 1];
            case 'w', C(i,:) = [1 1 1];
            case {'e','a'}, C(i,:) = [.6 .6 .6];
            otherwise,
                warning(sprintf('Unrecognized character %s',ch));
        end
    else
        error(sprintf('Entry %i is not a color signifier.',i));
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% function [gamma,x,w,sigu,like]=champagne(y,f,sigu,nem,nd,vcs,nupd,gamma0,retx);
% Output:
% gamma(nd,nd,nv) = voxel covariance matrices
% x(nv*nd,nt) = voxel current density
% w(nv*nd,nk) = reconstruction filter
% sigu(nk,nk) = noise covariance matrix
% like(nem,1) = likelihood
%
% Input:
% #s:
% nk = # of sensors
% nv = # of voxels
% nt = # of data time points

% y(nk,nt) = sensor data % preferable to be scaled to pT units
% f(nk,nv*nd) = lead field matrix % preferable to be normalized.
% sigu(nk,nk) = noise covariance matrix,  [] for automatic initialization
% nem = maximum # of iterations
% nd = # of orientations
% vcs = voxel covariance structure: 0 = scalar, 1 = diagonal, 2 = general
% nupd = 2: learn diagonal noise covariance, = 1: learn scalar, = 0: use input noise cov
% gamma0(nd,nd,nv) = initial voxel covariances, [] for automatic initialization
% retx = 1: return voxel current density, = 0: don't

% Two parameters in code that can be adjusted are:
% eps1 = Default is 1e-8 For numerical stability while inverting model covariance matrix.
% eps1z = Default is 1e-8 For numerical stability while calculating
%         voxel variance and auxilliary variable z

function [gamma,x,w,sigu,like]=champagne(y,f,sigu,nem,nd,vcs,nupd,gamma0,retx, fig);

if vcs==2 && nd>1
    [gamma x w sigu like]=champ_mat(y,f,sigu,nem,nd,nupd,gamma0,retx, fig);
else
    [gamma x w sigu like]=champ_vec(y,f,sigu,nem,nd,vcs,nupd,gamma0,retx, fig);
end

return

function [gamma,x,w,sigu,like]=champ_vec(y,f,sigu,nem,nd,vcs,nupd,gamma0,retx, fig);

eps1=1e-8;
eps1z=1e-8;
[nk nvd]=size(f);
nv=nvd/nd;
nt=size(y,2);
nb=1000;

cyy=y*y'/nt;

% Initialize voxel variances

if isempty(gamma0)
    f2=sum(f.^2,1);
    invf2=zeros(1,nvd);
    ff=find(f2>0);
    invf2(ff)=1./f2(ff);
    w=spdiags(invf2',0,nvd,nvd)*f';

    if nt>nb
        inu0=zeros(nvd,1);
        for it=1:nb:nt-nb+1
            inu0=inu0+sum((w*y(:,it:it+nb-1)).^2,2);
        end
        inu0=(inu0+sum((w*y(:,it+nb:nt)).^2,2))/nt;
    else
        inu0=mean((w*y).^2,2);
    end
%    disp(max(inu0-mean((w*y).^2,2)));
    inu0v=mean(reshape(inu0,nd,nv),1);
    inu1=reshape(ones(nd,1)*inu0v,nvd,1);
    vvec=inu1;
else
    vvec=zeros(nvd,1);
    for iv=1:nv
        jv=(iv-1)*nd+1:iv*nd;
        if vcs==0
            vvec(jv)=mean(diag(gamma0(:,:,iv)));
        else
            vvec(jv)=diag(gamma0(:,:,iv));
        end
    end
end

% initialize noise covariance

if nupd>0
    if isempty(sigu)
        ndr=.1;
        sigu=ndr*mean(mean(y.^2))*eye(nk);
    else
        sigu=mean(diag(sigu))*eye(nk);
    end
end

% Learn voxel variances

v=zeros(nv,1);
figure(fig);

like=zeros(nem,1);
for iem=1:nem
    vmat=spdiags(vvec,0,nvd,nvd);
    c=f*vmat*f'+sigu;
    [p d]=svd(c);
    d=max(real(diag(d)),0);
    invd=zeros(nk,1);
    ff=find(d>=eps1);
    invd(ff)=1./d(ff);
%    invd=1./d;
    invc=p*spdiags(invd,0,nk,nk)*p';

%    like(iem)=-.5*(sum(log(d))+nk*log(2*pi))-.5*sum(sum(y.*(invc*y)))/nt;
    like(iem)=-.5*(sum(log(max(d,eps1)))+nk*log(2*pi))-.5*sum(sum(invc.*cyy));
    subplot(2,2,1);plot((1:iem),like(1:iem));
    title(['Likelihood: ' int2str(iem) ' / ' int2str(nem)]);
    xlabel('iteration');
    set(gca(),'XLim',[0 iem]);

    fc=f'*invc;
    w=vmat*fc;
    x2=sum((w*cyy).*w,2);
    z=sum(fc.*f',2);

    if vcs==0
        x20=sum(reshape(x2,nd,nv),1);
        z0=sum(reshape(z,nd,nv),1);
        v0=(sqrt(z0)./max(z0,eps1z)).*sqrt(x20);
        vvec=reshape(ones(nd,1)*v0,nvd,1);
    else
        vvec=(sqrt(z)./max(z,eps1z)).*sqrt(x2);
    end

    if nupd>0
        fw=eye(nk)-f*w;
        sigy1=sum((fw*cyy).*fw,2);
        fgf=sum((fw*f*vmat).*f,2);
        ilam=sigy1+fgf;
        if nupd==1
            sigu=mean(ilam)*eye(nk);
        else
            sigu=diag(ilam);
        end
    end

    v=sum(reshape(vvec,nd,nv),1);
    subplot(2,2,2);plot((1:nv),v);
    title(['Voxel power: ' num2str(nv) ' / ' num2str(nv)]);
    xlabel('voxel index');
    set(gca(),'XLim',[1 nv]);
    drawnow
end

if retx==1
    x=w*y;
else
    x=[];
end

if nd==1
    gamma=reshape(vvec,1,1,nv);
else
    gamma=zeros(nd,nd,nv);
    for iv=1:nv
        gamma(:,:,iv)=diag(vvec((iv-1)*nd+1:iv*nd));
    end
end

return

function [gamma,x,w,sigu,like]=champ_mat(y,f,sigu,nem,nd,nupd,gamma0,retx, fig);

eps1=1e-8;
eps1z=1e-8;
[nk nvd]=size(f);
nv=nvd/nd;
nt=size(y,2);
nb=1000;

cyy=y*y'/nt;

% Initialize voxel covariances

if isempty(gamma0)
    f2=sum(f.^2,1);
    invf2=zeros(1,nvd);
    ff=find(f2>0);
    invf2(ff)=1./f2(ff);
    w=spdiags(invf2',0,nvd,nvd)*f';

    if nt>nb
        inu0=zeros(nvd,1);
        for it=1:nb:nt-nb+1
            inu0=inu0+sum((w*y(:,it:it+nb-1)).^2,2);
        end
        inu0=(inu0+sum((w*y(:,it+nb:nt)).^2,2))/nt;
    else
        inu0=mean((w*y).^2,2);
    end
%    disp(max(inu0-mean((w*y).^2,2)));
    inu0v=mean(reshape(inu0,nd,nv),1);
    inu1=reshape(ones(nd,1)*inu0v,nvd,1);
    vmat=spdiags(inu1,0,nvd,nvd);
else
    vmat=sparse(nvd,nvd);
    for iv=1:nv
        jv=(iv-1)*nd+1:iv*nd;
        vmat(jv,jv)=gamma0(:,:,iv);
    end
end

% initialize noise covariance

if nupd>0
    if isempty(sigu)
        ndr=.1;
        sigu=ndr*mean(mean(y.^2))*eye(nk);
    else
        sigu=mean(diag(sigu))*eye(nk);
    end
end

% Learn voxel covariances

v=zeros(nv,1);
figure;

like=zeros(nem,1);
for iem=1:nem
    c=f*vmat*f'+sigu;
    [p d]=svd(c);
    d=max(real(diag(d)),0);
    invd=zeros(nk,1);
    ff=find(d>=eps1);
    invd(ff)=1./d(ff);
%    invd=1./d;
    invc=p*spdiags(invd,0,nk,nk)*p';

%    like(iem)=-.5*(sum(log(d))+nk*log(2*pi))-.5*sum(sum(y.*(invc*y)))/nt;
    like(iem)=-.5*(sum(log(max(d,eps1)))+nk*log(2*pi))-.5*sum(sum(invc.*cyy));
     figure(fig)
    subplot(2,2,1);plot((1:iem),like(1:iem));
    title(['Likelihood: ' int2str(iem) ' / ' int2str(nem)]);
    xlabel('iteration');
    set(gca(),'XLim',[0 iem]);

    fc=f'*invc;
    w=vmat*fc;
%    x=w*y;
%    x2=mean(x.^2,2);
%    z=sum(fc.*f',2);

    for iv=1:nv
        jv=((iv-1)*nd+1:iv*nd);
        x2=w(jv,:)*cyy*w(jv,:)';
        z=fc(jv,:)*f(:,jv);

        [pz dz]=eig(z);
        dzd=max(real(diag(dz)),0);
        dz5=sqrt(dzd);
        z5=pz*diag(dz5)*pz';
        invdz5=1./sqrt(max(dzd,eps1z));
        invz5=pz*diag(invdz5)*pz';

        [px dx]=eig(z5*x2*z5);
        dx5=sqrt(max(real(diag(dx)),0));
        cx5=px*diag(dx5)*px';
        vmat(jv,jv)=invz5*cx5*invz5;
        v(iv)=sum(diag(vmat(jv,jv)));
    end

    if nupd>0
        fw=eye(nk)-f*w;
        sigy1=sum((fw*cyy).*fw,2);
        fgf=sum((fw*f*vmat).*f,2);
        ilam=sigy1+fgf;
        if nupd==1
            sigu=mean(ilam)*eye(nk);
        else
            sigu=diag(ilam);
        end
    end

    subplot(2,2,2);plot((1:nv),v);
    title(['Voxel power: ' num2str(nv) ' / ' num2str(nv)]);
    xlabel('voxel index');
    set(gca(),'XLim',[1 nv]);
    drawnow
end

if retx==1
    x=w*y;
else
    x=[];
end

gamma=zeros(nd,nd,nv);
for iv=1:nv
    jv=(iv-1)*nd+1:iv*nd;
    gamma(:,:,iv)=vmat(jv,jv);
end

return

author: spm


text:
%  Figure 7.43      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% script to generate Fig. 7.43
clf;
f=[-10 1 0;-16 0 1;0 0 0];
g=[0 0 10]';
h=[1 0 0];
j=0;
[np,dp]=ss2tf(f,g,h,0);
np=[0 0 0 10];
nc=conv([-1 0.718],[1 1.87]);
r=[-0.95+6.17*i -0.95-6.17*i];
dc=poly(r);
num=conv(np,nc);
den=conv(dp,dc);
rlocus(num,den);
axis([-9 4 -9 9])
title('Fig.7.43 Root locus for reduced-order DC Servo design')
grid;


text:
function [AtA,A] = corrMatrix2D(obj,i)
% calucate 2D correlation matrix
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

nCha = size(obj.kCalib{i},3);

% A = [];
if(isreal(obj.kCalib{i}))
    A = zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision);
else
    A = complex(zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision),zeros(prod(obj.calibSize - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision));
end

counter = 1;
for n=1:nCha
	A(:,counter:counter+prod(obj.kernelSize)-1)  = im2col(obj.kCalib{i}(:,:,n),obj.kernelSize,'sliding').';
    counter = counter + prod(obj.kernelSize);
% 	A = [A, tmp];
end

AtA = A'*A;

end

author: thomaskuestner


text:
function [res,tau1,tau2,tau,p] = AccSumHugeN(p)
%ACCSUMHUGEN  Faithful rounding of sum(p) for huge dimension
%
%   res = AccSumHugeN(p)
%
%On return, res is a faithful rounding of sum(p), also in the presence
%  of underflow. Input vector p may be single or double precision.
%
%Implements Algorithm 8.1 from
%  S.M. Rump, T. Ogita, S. Oishi: Accurate Floating-point Summation II:
%    Sign, K-fold Faithful and Rounding to Nearest, Siam J. Sci. Comput.,
%    31(2):1269-1302, 2008.
%Requires (4m+4K+3)n flops for m executions of repeat-until loop in
%  Transform and K executions of the while-loop.
%
%Reference implementation! Slow due to interpretation!
%

% written  03/03/07     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest, complex input
%

  if ~isreal(p)
    res = complex(AccSumHugeN(real(p)),AccSumHugeN(imag(p)));
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if isa(p,'double')
    nmax = 2^50;            % nmax = 1,125,899,906,842,624
  else
    nmax = 2^21;            % nmax = 2,097,152
  end
  if length(p)>nmax
    error(['maximum length of input vector for AccSumHugeN ' int2str(nmax) '.'])
  end

  kPhi = 2;
  [tau1,tau2,p,sigma,Ms] = Transform(p,0,kPhi);   % Ms = 2^M
  if sigma<=realmin             % p_i identical zero
    res = tau1;
    tau = 0*res;                % make sure tau has same precision
    if rndold, setround(rndold); end
    return
  end
  if isa(p,'double'), prec='double'; else prec='single'; end
  u = 0.5*eps(prec);
  tau = zeros(1,16);
  K = 0;
  phi = Ms*u;
  factor = 2*Ms*Ms*u;
  sigmas = phi*sigma;
  while 1
    K = K+1;
    sigma = sigmas;
    [tau(K),p] = ExtractVector(p,sigma);
    sigmas = phi*sigma;
    if ( factor*sigma<=abs(tau1) ) | ( sigma<=realmin )
      taus = tau2 + sum(p);
      for k=K:-1:1
        taus = taus + tau(k);
      end
      res = tau1 + taus;
      if rndold, setround(rndold); end
      return
    end
  end

repo: OpenMAS


text:
classdef CosineSim < dagnn.Filter

  methods
    function outputs = forward(self, inputs, params)
      outputs{1} = vl_nncosinesim(inputs{1}, inputs{2}) ;
    end

    function [derInputs, derParams] = backward(self, inputs, params, derOutputs)
      derInputs{1} = vl_nncosinesim(inputs{1}, inputs{2}, derOutputs{1}) ;
      derParams = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes{1} = ones(1,4) ;
      outputSizes{1}(4) = inputSizes{1}(4) ;
    end

    function obj = CosineSim(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MFT


text:
function  [h, compUp] =  lfmComputeH4(gamma1_p, gamma1_m, sigma2, t1, preFactor, preExp,...
    mode, term )
% LFMCOMPUTEH4 Helper function for computing part of the LFM kernel.
% FORMAT
% DESC computes a portion of the LFM kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : David Luengo, 2007
%
% COPYRIGHT : Mauricio Alvarez, 2008
%
% MODIFICATIONS : Neil D. Lawrence, 2007
%
%
%
% SEEALSO : lfmKernParamInit, lfmXlfmKernCompute

% KERN

% Evaluation of h

if nargin<8
    term =[];
end

if ~mode
    if ~term
        if nargout>1
            compUp = lfmComputeUpsilonVector(gamma1_p,sigma2, t1);
            h = compUp*( preExp/preFactor(1) - conj(preExp)/preFactor(2)).';
        else
            h = lfmComputeUpsilonVector(gamma1_p,sigma2, t1)*( preExp/preFactor(1) - conj(preExp)/preFactor(2)).';
        end
    else
        if nargout>1
            compUp = lfmComputeUpsilonVector(gamma1_p,sigma2, t1);
            h = compUp*(preExp/preFactor(1)).' - conj(compUp)*(preExp/preFactor(2)).';
        else
            upsilon = lfmComputeUpsilonVector(gamma1_p,sigma2, t1);
            h = upsilon*(preExp/preFactor(1)).' - conj(upsilon)*(preExp/preFactor(2)).';
        end
    end
else
    if nargout > 1
        compUp{1} = lfmComputeUpsilonVector(gamma1_p,sigma2, t1);
        compUp{2} = lfmComputeUpsilonVector(gamma1_m,sigma2, t1);
        h =  compUp{1}*( preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
            + compUp{2}*( preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
    else
        h =  lfmComputeUpsilonVector(gamma1_p,sigma2, t1)*( preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
            + lfmComputeUpsilonVector(gamma1_m,sigma2, t1)*( preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ x, know ] = p03_sol ( m, know )

%*****************************************************************************80
%
%% P03_SOL returns known solutions for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution.
%
  if ( know == 0 )
    know = 1;
    x(1:m,1) = [ ...
      0.999980569087140, ...
      0.500000721280566, ...
      0.333341891834645, ...
      0.249997266604697 ]';
  else
    know = 0;
    x = zeros ( m, 1 );
  end

  return
end


text:
function [Ef, Varf, xtnn, xt1, xt2] = gp_cpred(gp,x,y,xt,ind,varargin)
%GP_CPRED Conditional predictions using specific covariates
%
%  Description
%    GP_CPRED(GP,X,Y,XT,IND,OPTIONS) does predictions using only
%    covariates specified in vector IND. Other covariates are fixed to
%    either mean, median or values chosen by user. Returns predictions for
%    latent values, variance and corresponding inputs. If IND=0, only time
%    is used as a covariate for Cox-PH model.
%
%   OPTIONS is optional parameter-value pair
%      method - which value to fix the not used covariates, 'median'
%               (default), 'mean' or 'mode'
%      var    - vector specifying optional values for not used covariates,
%               elements corresponding to mean/median values should
%               be set to NaN.
%      plot   - option for plotting, 'off' (default) or 'on'
%      normdata - a structure with fields xmean, xstd, ymean, and ystd
%               to allow plotting in the original data scale (see
%               functions normdata and denormdata)
%      target - option for choosing what is computed 'mu' (default),
%               'f' or 'cdf'
%      tr     - Euclidean distance treshold for not using grid points when
%               doing predictions with 2 covariates, default 0.25
%      predcf - an index vector telling which covariance functions are
%               used for prediction. Default is all (1:gpcfn).
%               See additional information below.
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, the expected
%               value for the ith case.

ip=inputParser;
ip.FunctionName = 'GP_CPRED';
ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('xt',  @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addRequired('ind', @(x) ~isempty(x) && isvector(x))
ip.addParamValue('var',  [], @(x) isreal(x))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                 isvector(x) && isreal(x) && all(isfinite(x)&x>0))
ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                 (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
ip.addParamValue('method', 'median', @(x)  ismember(x, {'median', 'mean' 'mode'}))
ip.addParamValue('plot', 'off', @(x)  ismember(x, {'on', 'off'}))
ip.addParamValue('tr', 0.25, @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('target', 'mu', @(x) ismember(x,{'f','mu','cdf'}))
ip.addParamValue('prct', [5 50 95], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('normdata', struct(), @(x) isempty(x) || isstruct(x))
ip.addParamValue('xlabels', [], @(x) isempty(x) || iscell(x));
ip.parse(gp, x, y, xt, ind, varargin{:});
zt=ip.Results.zt;
options=struct();
options.predcf=ip.Results.predcf;
options.prct=ip.Results.prct;
options.tstind=ip.Results.tstind;
method = ip.Results.method;
vars = ip.Results.var;
plot_results = ip.Results.plot;
tr = ip.Results.tr;
target = ip.Results.target;
if strcmp(target,'f')
    options = rmfield(options,'prct');
end
yt=ip.Results.yt;
if ~isempty(yt)
  options.yt=yt;
end
z=ip.Results.z;
if ~isempty(z)
  options.z=z;
end
if ~isempty(zt)
  options.zt=zt;
end
if isempty(zt)
  options.zt=z;
end
% normdata
nd=ip.Results.normdata;
ipnd=inputParser;
ipnd.FunctionName = 'normdata';
ipnd.addParamValue('xmean',zeros(1,size(x,2)),@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('xstd',ones(1,size(x,2)),@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('xlog',zeros(1,size(x,2)),@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('ymean',0,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('ystd',1,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('ylog',0,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.parse(nd);
nd=ipnd.Results;

[tmp, nin] = size(x);

if iscell(gp)
  liktype=gp{1}.lik.type;
else
  liktype=gp.lik.type;
end

if isequal(liktype, 'Coxph') && isequal(target,'mu')
    target='f';
    warning('GP_CPRED: Target ''mu'' not applicable for a Cox-PH model. Switching to target ''f''')
end

if ~isempty(vars) && (~isvector(vars) || length(vars) ~= nin)
  error('Vector defining fixed variable values must be same length as number of covariates')
end

xto=xt; [n,tmp]=size(xto);

if length(ind)==1

  if ind~=0
    [xtnn, iu] = unique(xt(:,ind));
  else
    xtnn = xt(1,:);
    iu = 1;
  end
  if ~isempty(z)
    options.zt = options.zt(iu);
  end
  if ~isempty(yt)
    options.yt = options.yt(iu);
  end
  xt = repmat(feval(method,xt), size(xtnn,1), 1);
  if ~isempty(vars)
    xt(:,~isnan(vars)) = repmat(vars(~isnan(vars)), length(xtnn), 1);
  end

  if ind>0
    xt(:,ind) = xtnn;
  end
  if ~strcmp(liktype, 'Coxph')
    switch target
      case 'f'
        [Ef, Varf] = gp_pred(gp, x, y, xt, options);
      case 'mu'
        prctmu = denormdata(gp_predprctmu(gp, x, y, xt, options),nd.ymean,nd.ystd);
        Ef = prctmu; Varf = [];
      case 'cdf'
        cdf = gp_predcdf(gp, x, y, xt, options);
        Ef = cdf; Varf = [];
    end
  else
    [Ef1,Ef2,Covf] = pred_coxph(gp,x,y,xt,'z',options.z,'zt',options.zt);
    nt=size(Ef1,1);
    if ind>0
      % conditional posterior given Ef1=E[Ef1]
      Ef = Ef2;
      Varf = diag(Covf(nt+1:end,nt+1:end)-Covf(nt+1:end,1:nt)*(Covf(1:nt,1:nt)\Covf(1:nt,nt+1:end)));
    else
      % conditional posterior given Ef2=E[Ef2]
      Ef = Ef1;
      Varf = diag(Covf(1:nt, 1:nt)-Covf(1:nt,nt+1:end)*(Covf(nt+1:end,nt+1:end)\Covf(nt+1:end,1:nt)));
      xtnn = gp.lik.xtime;
    end
  end
  if isequal(plot_results, 'on')
    if ind>0
      if ind>=1&numel(nd.xmean)>=ind
        xtnn=denormdata(xtnn,nd.xmean(ind),nd.xstd(ind));
      end
      deltadist=gp_finddeltadist(gp);
      if ~ismember(ind,deltadist)
        switch target
          case 'f'
            plot(xtnn, Ef, 'ob', xtnn, Ef, '-k', xtnn, Ef-1.64*sqrt(Varf), '--b', xtnn, Ef+1.64*sqrt(Varf), '--b')
          case 'mu'
            plot(xtnn, prctmu(:,2), 'ob', xtnn, prctmu(:,2), '-k', xtnn, prctmu(:,1), '--b', xtnn, prctmu(:,3), '--b')
          case 'cdf'
            plot(xtnn, Ef, 'o-b')
        end
      else
        switch target
          case 'f'
            plot(xtnn, Ef, 'ob', [xtnn xtnn]',[Ef-1.64*sqrt(Varf) Ef+1.64*sqrt(Varf)]', '-b')
            xlim([1.5*xtnn(1)-0.5*xtnn(2) 1.5*xtnn(end)-.5*xtnn(end-1)])
            set(gca,'xtick',xtnn)
          case 'mu'
            plot(xtnn, prctmu(:,2), 'ob', [xtnn xtnn]',[prctmu(:,1) prctmu(:,3)]', '-b')
            xlim([1.5*xtnn(1)-0.5*xtnn(2) 1.5*xtnn(end)-.5*xtnn(end-1)])
            set(gca,'xtick',xtnn)
          case 'cdf'
            plot(xtnn, Ef, 'ob')
        end
      end
    else
      % use stairs for piecewise constant baseline hazard
      xtnn = gp.lik.stime;
      [xx,yy]=stairs(xtnn, [Ef;Ef(end)]);
      [xx,yyl]=stairs(xtnn, [Ef-1.64*sqrt(Varf);Ef(end)-1.64*sqrt(Varf(end))]);
      [xx,yyu]=stairs(xtnn, [Ef+1.64*sqrt(Varf);Ef(end)+1.64*sqrt(Varf(end))]);
      plot(xx, yy, '-k', xx, yyl, '--b', xx, yyu, '--b')
    end
  end

elseif length(ind)==2

  uu1=unique(xt(:,ind(1)));
  uu2=unique(xt(:,ind(2)));
  nu1=numel(uu1);
  nu2=numel(uu2);
  if nu1==2 || nu2==2
    % First or second covariate binary

    if nu1>2 && nu2==2
      % switch indeces, so that binary covariate is first
      tmp=ind(1);ind(1)=ind(2);ind(2)=tmp;
      tmp=uu1;uu1=uu2;uu2=tmp;
    end

    xt1=xt(xt(:,ind(1))==uu1(1),:);
    xt2=xt(xt(:,ind(1))==uu1(2),:);
    [xtnn1, iu1] = unique(xt1(:,ind(2)));
    [xtnn2, iu2] = unique(xt2(:,ind(2)));

    options1=options;
    options2=options;
    if ~isempty(z)
      options1.zt = options.zt(iu1);
      options2.zt = options.zt(iu2);
    end

    xt1 = repmat(feval(method,xt1), length(xtnn1), 1);
    xt2 = repmat(feval(method,xt2), length(xtnn2), 1);
    if ~isempty(vars)
      xt1(:,~isnan(vars)) = repmat(vars(~isnan(vars)), length(xtnn1), 1);
      xt2(:,~isnan(vars)) = repmat(vars(~isnan(vars)), length(xtnn2), 1);
    end
    xt1(:,ind(1)) = uu1(1); xt1(:,ind(2)) = xtnn1;
    xt2(:,ind(1)) = uu1(2); xt2(:,ind(2)) = xtnn2;

    if ~strcmp(liktype, 'Coxph')
      switch target
        case 'f'
          [Ef1, Varf1] = gp_pred(gp, x, y, xt1);
          [Ef2, Varf2] = gp_pred(gp, x, y, xt2);
        case 'mu'
          prctmu1 = denormdata(gp_predprctmu(gp, x, y, xt1, options1),nd.ymean,nd.ystd);
          prctmu2 = denormdata(gp_predprctmu(gp, x, y, xt2, options2),nd.ymean,nd.ystd);
      end
    else
      [Ef11,Ef12,Covf] = pred_coxph(gp,x,y,xt1, 'z', z);
      Ef1 = Ef12; Varf1 = diag(Covf(size(Ef11,1)+1:end,size(Ef11,1)+1:end));
      [Ef21,Ef22,Covf] = pred_coxph(gp,x,y,xt2, 'z', z);
      Ef2 = Ef22; Varf2 = diag(Covf(size(Ef21,1)+1:end,size(Ef21,1)+1:end));
    end

    if isequal(plot_results, 'on')
      xtnn1=denormdata(xtnn1,nd.xmean(ind(2)),nd.xstd(ind(2)));
      xtnn2=denormdata(xtnn2,nd.xmean(ind(2)),nd.xstd(ind(2)));
      if nu1>2 && nu2==2
        lstyle10='or';lstyle11='-r';lstyle12='--r';
        lstyle20='ob';lstyle21='-b';lstyle22='--b';
      else
        lstyle10='ob';lstyle11='-b';lstyle12='--b';
        lstyle20='or';lstyle21='-r';lstyle22='--r';
      end
      deltadist=gp_finddeltadist(gp);
      if ~ismember(ind(2),deltadist)
        switch target
          case 'f'
            plot(xtnn1, Ef1, lstyle10, xtnn1, Ef1, lstyle11, xtnn1, Ef1-1.64*sqrt(Varf1), lstyle12, xtnn1, Ef1+1.64*sqrt(Varf1), lstyle12); hold on;
            plot(xtnn2, Ef2, lstyle20, xtnn2, Ef2, lstyle21, xtnn2, Ef2-1.64*sqrt(Varf2), lstyle22, xtnn2, Ef2+1.64*sqrt(Varf2), lstyle22); hold off;
          case 'mu'
            plot(xtnn1, prctmu1(:,2), lstyle10, xtnn1, prctmu1(:,2), lstyle11, xtnn1, prctmu1(:,1), lstyle12, xtnn1, prctmu1(:,3), lstyle12); hold on;
            plot(xtnn2, prctmu2(:,2), lstyle20, xtnn2, prctmu2(:,2), lstyle21, xtnn2, prctmu2(:,1), lstyle22, xtnn2, prctmu2(:,3), lstyle22); hold off;
        end
      else
        delta=(diff(xtnn1(1:2))/10);
        switch target
          case 'f'
            plot(xtnn1-delta, Ef1, lstyle10, [xtnn1 xtnn1]'-delta, [Ef1-1.64*sqrt(Varf1) Ef1+1.64*sqrt(Varf1)]', lstyle11); hold on;
            plot(xtnn2+delta, Ef2, lstyle20, [xtnn2 xtnn2]'+delta, [Ef2-1.64*sqrt(Varf2) Ef2+1.64*sqrt(Varf2)]', lstyle21); hold off;
            xlim([1.5*xtnn1(1)-0.5*xtnn1(2) 1.5*xtnn1(end)-.5*xtnn1(end-1)])
          case 'mu'
            plot(xtnn1-delta, prctmu1(:,2), lstyle10, [xtnn1 xtnn1]'-delta, [prctmu1(:,1) prctmu1(:,3)]', lstyle11); hold on;
            plot(xtnn2+delta, prctmu2(:,2), lstyle20, [xtnn2 xtnn2]'+delta, [prctmu2(:,1) prctmu2(:,3)]', lstyle21); hold off;
            xlim([1.5*xtnn1(1)-0.5*xtnn1(2) 1.5*xtnn1(end)-.5*xtnn1(end-1)])
        end
      end
    end
    switch target
      case 'f'
        Ef = {Ef1  Ef2}; Varf = {Varf1 Varf2}; xtnn={xtnn1 xtnn2};
      case 'mu'
        Ef = {prctmu1 prctmu2}; Varf = {[] []}; xtnn={xtnn1 xtnn2};
    end

  else
    % both the first and the second covariate are non-binary
    deltadist=gp_finddeltadist(gp);
    if ~ismember(ind(1),deltadist)
      xtnn1 = linspace(min(xt(:,ind(1))), max(xt(:,ind(1))), 20);
    else
      xtnn1 = unique(xt(:,ind(1)));
    end
    if ~ismember(ind(2),deltadist)
      xtnn2 = linspace(min(xt(:,ind(2))), max(xt(:,ind(2))), 20);
    else
      xtnn2 = unique(xt(:,ind(2)));
    end
    [XT1, XT2] = meshgrid(xtnn1, xtnn2); XT1=XT1(:); XT2=XT2(:);
    xtnn1=denormdata(xtnn1,nd.xmean(ind(1)),nd.xstd(ind(1)));
    xtnn2=denormdata(xtnn2,nd.xmean(ind(2)),nd.xstd(ind(2)));
    if ~isempty(z)
      options.zt = repmat(options.zt(1), size(XT1));
    end
    xt = repmat(feval(method,xt), length(XT1), 1);
    if ~isempty(vars)
      xt(:,~isnan(vars)) = repmat(vars(~isnan(vars)), length(XT1), 1);
    end
    xt(:,ind) = [XT1 XT2];
    if ~strcmp(liktype, 'Coxph')
      switch target
        case 'f'
          [Ef, Varf] = gp_pred(gp, x, y, xt);
        case 'mu'
          prctmu = gp_predprctmu(gp, x, y, xt, options, 'prct', 50);
          Ef = prctmu; Varf = [];
        case 'cdf'
          cdf = gp_predcdf(gp, x, y, xt, options);
          Ef = cdf; Varf = [];
      end
    else
      [Ef1,Ef2,Covf] = pred_coxph(gp,x,y,xt, options);
      Ef = Ef2; Varf = diag(Covf(size(Ef1,1)+1:end,size(Ef1,1)+1:end));
    end

    indd = zeros(size(Ef));

    for i2=1:n
      for i3=1:numel(XT1)
        if sqrt(sum((xto(i2,ind)-xt(i3,ind)).^2)) < tr
          indd(i3) = 1;
        end
      end
    end

    XT1(indd==0) = NaN; XT2(indd==0) = NaN; Ef(indd==0) = NaN; Varf(indd==0) = NaN;

    if isequal(plot_results, 'on')
      xtnn1=denormdata(xtnn1,nd.xmean(ind(1)),nd.xstd(ind(1)));
      xtnn2=denormdata(xtnn2,nd.xmean(ind(2)),nd.xstd(ind(2)));
      surf(reshape(XT1,numel(xtnn2),numel(xtnn1)), reshape(XT2,numel(xtnn2),numel(xtnn1)), reshape(Ef,numel(xtnn2),numel(xtnn1)))
      view(2)
      axis tight
      shading flat
      colormap(mapcolor(Ef,repmat(nanmedian(Ef(:)),[1 2])))
      colorbar('EastOutside')
    end

    %xtnn = [XT1(indd==1), XT2(indd==1)]; Ef = Ef(indd==1); Varf = Varf(indd==1);
    xtnn = [XT1, XT2];
  end

else
  error('Only 1 or 2 covariates can be defined for predicting')
end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
%Example of a script that can be used to set up different thermodynamic models depending on the pH temp etc

% [modelT,directions]=setupThermoModel(model,metAbbrAlbertyAbbr,metGroupCont,Alberty2006,temp,cpHa,ppHa,epHa,isCyto,isExt,biomassRxnAbbr,Ecoli_symphID_rxnAbbr,Legendre,useKeqData,nStdDevGroupCont,figures)

mkdir one
cd one
metBoundsFile=[];
rxnBoundsFile=[];
Legendre=1;
useKeqData=1;
nStdDevGroupCont=0;
cumNormProbCutoff=0.2;
figures=0;
printToFile=1;
[modelT,directions]=setupThermoModel(iAF1,metAbbrAlbertyAbbr,metGroupCont,Alberty2006,310.15,7.7,7.7,7.7,0.25,0.25,'Ec_biomass_iAF1260_core_59p81M',Ecoli_symphID_rxnAbbr,metBoundsFile,rxnBoundsFile,Legendre,useKeqData,nStdDevGroupCont,cumNormProbCutoff,figures,printToFile);

cd ../
mkdir two
cd two
Legendre=1;
useKeqData=1;
nStdDevGroupCont=1;
figures=1;
[modelT,directions]=setupThermoModel(iAF1,metAbbrAlbertyAbbr,metGroupCont,Alberty2006,310.15,7.7,7.7,7.7,0.25,0.25,'Ec_biomass_iAF1260_core_59p81M',Ecoli_symphID_rxnAbbr,metBoundsFile,rxnBoundsFile,Legendre,useKeqData,nStdDevGroupCont,cumNormProbCutoff,figures,printToFile);

cd ../
mkdir three
cd three
Legendre=1;
useKeqData=1;
nStdDevGroupCont=0;
figures=0;
[modelT,directions]=setupThermoModel(iAF1,metAbbrAlbertyAbbr,metGroupCont,Alberty2006,298.15,7,7,7,0,0,'Ec_biomass_iAF1260_core_59p81M',Ecoli_symphID_rxnAbbr,metBoundsFile,rxnBoundsFile,Legendre,useKeqData,nStdDevGroupCont,cumNormProbCutoff,figures,printToFile);

cd ../
mkdir four
cd four
Legendre=1;
useKeqData=0;
nStdDevGroupCont=0;
figures=0;
[modelT,directions]=setupThermoModel(iAF1,metAbbrAlbertyAbbr,metGroupCont,Alberty2006,310.15,7.7,7.7,7.7,0.25,0.25,'Ec_biomass_iAF1260_core_59p81M',Ecoli_symphID_rxnAbbr,metBoundsFile,rxnBoundsFile,Legendre,useKeqData,nStdDevGroupCont,cumNormProbCutoff,figures,printToFile);

cd ../
mkdir five
cd five
Legendre=1;
useKeqData=0;
nStdDevGroupCont=1;
figures=0;
[modelT,directions]=setupThermoModel(iAF1,metAbbrAlbertyAbbr,metGroupCont,Alberty2006,310.15,7.7,7.7,7.7,0.25,0.25,'Ec_biomass_iAF1260_core_59p81M',Ecoli_symphID_rxnAbbr,metBoundsFile,rxnBoundsFile,Legendre,useKeqData,nStdDevGroupCont,cumNormProbCutoff,figures,printToFile);

cd ../
mkdir six
cd six
Legendre=1;
useKeqData=0;
nStdDevGroupCont=0;
figures=0;
[modelT,directions]=setupThermoModel(iAF1,metAbbrAlbertyAbbr,metGroupCont,Alberty2006,298.15,7,7,7,0,0,'Ec_biomass_iAF1260_core_59p81M',Ecoli_symphID_rxnAbbr,metBoundsFile,rxnBoundsFile,Legendre,useKeqData,nStdDevGroupCont,cumNormProbCutoff,figures,printToFile);

author: opencobra

repo: cobratoolbox