text:
% [INPUT]
% ds = A structure representing the dataset.
% sn = A string representing the serial number of the result file.
% temp = A string representing the full path to the Excel spreadsheet used as template for the result file.
% out = A string representing the full path to the Excel spreadsheet to which the results are written, eventually replacing the previous ones.
% bw = An integer [21,252] representing the dimension of each rolling window (optional, default=252).
% sel = A string representing the time series selection method (optional, default='F'):
%   - 'F' for firms;
%   - 'G' for groups.
% rr = A float [0,1] representing the recovery rate in case of default (optional, default=0.4).
% pw = A string representing the probabilities of default weighting method (optional, default='W'):
%   - 'A' for plain average;
%   - 'W' for progressive average.
% md = A string representing the multivariate distribution used by the CIMDO model (optional, default='N'):
%   - 'N' for normal distribution;
%   - 'T' for Student's T distribution.
% analyze = A boolean that indicates whether to analyse the results and display plots (optional, default=false).
%
% [OUTPUT]
% result = A structure representing the original dataset inclusive of intermediate and final calculations.
% stopped = A boolean that indicates whether the process has been stopped through user input.
%
% [NOTES]
% If the number of time series is greater than 10, systemic risk measures are calculated as the average of the results of multiple reduced portfolios.
% Each reduced portfolio contains 6 time series and they are created by ranking firms or groups according to the following criteria:
%   1) top 3 and bottom 3 entities by CDS spreads;
%   2) top 3 and bottom 3 entities by variance of returns;
%   3) top 3 and bottom 3 entities by market capitalization (if data is available);
%   4) top 3 and bottom 3 entities by liabilities (if data is available).

function [result,stopped] = run_cross_entropy(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('ds',@(x)validateattributes(x,{'struct'},{'nonempty'}));
        ip.addRequired('sn',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addRequired('temp',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addRequired('out',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addOptional('bw',252,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 21 '<=' 252 'scalar'}));
        ip.addOptional('sel','F',@(x)any(validatestring(x,{'F' 'G'})));
        ip.addOptional('rr',0.4,@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0 '<=' 1 'scalar'}));
        ip.addOptional('pw','W',@(x)any(validatestring(x,{'A' 'W'})));
        ip.addOptional('md','N',@(x)any(validatestring(x,{'N' 'T'})));
        ip.addOptional('analyze',false,@(x)validateattributes(x,{'logical'},{'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    ds = validate_dataset(ipr.ds,'CrossEntropy');
    sn = ipr.sn;
    temp = validate_template(ipr.temp);
    out = validate_output(ipr.out);
    bw = ipr.bw;
    sel = validate_selection(ipr.sel,ds.Groups);
    rr = ipr.rr;
    pw = ipr.pw;
    md = ipr.md;
    analyze = ipr.analyze;

    nargoutchk(1,2);

    [result,stopped] = run_cross_entropy_internal(ds,sn,temp,out,bw,sel,rr,pw,md,analyze);

end

function [result,stopped] = run_cross_entropy_internal(ds,sn,temp,out,bw,sel,rr,pw,md,analyze)

    result = [];
    stopped = false;
    e = [];

    ds = initialize(ds,sn,bw,sel,rr,pw,md);
    t = ds.T;

    bar = waitbar(0,'Initializing cross-entropy measures...','CreateCancelBtn',@(src,event)setappdata(gcbf(),'Stop',true));
    setappdata(bar,'Stop',false);
    cleanup = onCleanup(@()delete(bar));

    pause(1);
    waitbar(0,bar,'Calculating cross-entropy measures...');
    pause(1);

    try

        k = size(ds.Portfolios,1);
        tk = t * k;

        futures(1:t,1:k) = parallel.FevalFuture;
        futures_max = 0;
        futures_results = cell(t,k);

        for j = 1:k
            windows_r = extract_rolling_windows(ds.Portfolios{j,2},ds.BW);
            windows_pods = extract_rolling_windows(ds.Portfolios{j,3},ds.BW);

            for i = 1:t
                futures(i,j) = parfeval(@main_loop,1,windows_r{i},windows_pods{i},ds.PW,ds.MD);
            end
        end

        for i = 1:tk
            if (getappdata(bar,'Stop'))
                stopped = true;
                break;
            end

            [future_index,value] = fetchNext(futures);
            [future_i,future_j] = ind2sub([t k],future_index);
            futures_results{future_i,future_j} = value;

            futures_max = max([future_index futures_max]);
            waitbar((futures_max - 1) / tk,bar);

            if (getappdata(bar,'Stop'))
                stopped = true;
                break;
            end
        end

    catch e
    end

    try
        cancel(futures);
    catch
    end

    if (~isempty(e))
        delete(bar);
        rethrow(e);
    end

    if (stopped)
        delete(bar);
        return;
    end

    pause(1);
    waitbar(1,bar,'Finalizing cross-entropy measures...');
    pause(1);

    try
        ds = finalize(ds,futures_results);
    catch e
        delete(bar);
        rethrow(e);
    end

    pause(1);
    waitbar(1,bar,'Writing cross-entropy measures...');
    pause(1);

    try
        write_results(ds,temp,out);
        delete(bar);
    catch e
        delete(bar);
        rethrow(e);
    end

    if (analyze)
        analyze_result(ds)
    end

    result = ds;

end

%% PROCESS

function ds = initialize(ds,sn,bw,sel,rr,pw,md)

    n = ds.N;
    t = ds.T;

    if (strcmp(sel,'G'))
        g = ds.Groups;
        gd = ds.GroupDelimiters;
        gs = cell(g,1);

        for i = 1:g
            if (i == 1)
                gs{i} = 1:gd(1);
            elseif (i == g)
                gs{i} = (gd(i - 1) + 1):n;
            else
                gs{i} = (gd(i - 1) + 1):gd(i);
            end
        end

        cds_ref = ds.CDS;
        cds = zeros(t,g);

        r_ref = ds.Returns;
        r = zeros(t,g);

        for i = 1:g
            gs_i = gs{i};
            gs_i_len = numel(gs_i);

            r_i = r_ref(:,gs_i);
            w_i = 1 ./ (repmat(gs_i_len,t,1) - sum(isnan(r_i),2));
            r(:,i) = sum(r_i .* repmat(w_i,1,gs_i_len),2,'omitnan');

            cds_i = cds_ref(:,gs_i);
            w_i = 1 ./ (repmat(gs_i_len,t,1) - sum(isnan(cds_i),2));
            cds(:,i) = sum(cds_i .* repmat(w_i,1,gs_i_len),2,'omitnan');
        end

        cp_ref = ds.Capitalizations;

        if (isempty(cp_ref))
            cp = [];
        else
            cp = zeros(t,g);

            for i = 1:g
                cp(:,i) = sum(cp_ref(:,gs{i}),2,'omitnan');
            end
        end

        lb_ref = ds.Liabilities;

        if (isempty(lb_ref))
            lb = [];
        else
            lb = zeros(t,g);

            for i = 1:g
                lb(:,i) = sum(lb_ref(:,gs{i}),2,'omitnan');
            end
        end

        n = g;
    else
        cds = ds.CDS;
        cp = ds.Capitalizations;
        lb = ds.Liabilities;
        r = ds.Returns;
    end

    pods = cds ./ rr;

    if (n <= 10)
        nc = n;

        if (strcmp(sel,'G'))
            pfc = ds.GroupShortNames.';
        else
            pfc = ds.FirmNames;
        end

        pf = {'Unique' r pods []};
    else
        nc = 6;
        nch = nc / 2;

        pfc = arrayfun(@(x)sprintf('C%d',x),1:nc,'UniformOutput',false);
        pf = [repmat({''},4,1) cell(4,3)];

        rw = extract_rolling_windows(ds.Returns,bw);

        pf_indices = zeros(t,nc);
        pf_pods = zeros(t,nc);
        pf_r = zeros(t,nc);

        for i = 1:t
            [cds_i,indices] = sort(cds(i,:),'ascend');

            indices(isnan(cds_i)) = [];
            indices = [fliplr(indices(end-nch+1:end)) fliplr(indices(1:nch))];

            pf_r(i,:) = r(i,indices);
            pf_pods(i,:) = pods(i,indices);
            pf_indices(i,:) = indices;
        end

        pf(1,1:4) = {'CDS Spreads' pf_r pf_pods pf_indices};

        pf_indices = zeros(t,nc);
        pf_pods = zeros(t,nc);
        pf_r = zeros(t,nc);

        for i = 1:t
            [v_i,indices] = sort(var(rw{i}),'ascend');

            indices(isnan(v_i)) = [];
            indices = [fliplr(indices(end-nch+1:end)) fliplr(indices(1:nch))];

            pf_r(i,:) = r(i,indices);
            pf_pods(i,:) = pods(i,indices);
            pf_indices(i,:) = indices;
        end

        pf(2,1:4) = {'Returns Variance' pf_r pf_pods pf_indices};

        if (isempty(cp))
            offset = 3;
        else
            pf_indices = zeros(t,nc);
            pf_pods = zeros(t,nc);
            pf_r = zeros(t,nc);

            for i = 1:t
                [cp_i,indices] = sort(cp(i,:),'ascend');

                indices(isnan(cp_i)) = [];
                indices = [fliplr(indices(end-nch+1:end)) fliplr(indices(1:nch))];

                pf_r(i,:) = r(i,indices);
                pf_pods(i,:) = pods(i,indices);
                pf_indices(i,:) = indices;
            end

            pf(3,:) = {'Capitalization' pf_r pf_pods pf_indices};
            offset = 4;
        end

        if (~isempty(lb))
            pf_indices = zeros(t,nc);
            pf_pods = zeros(t,nc);
            pf_r = zeros(t,nc);

            for i = 1:t
                [lb_i,indices] = sort(lb(i,:),'ascend');

                indices(isnan(lb_i)) = [];
                indices = [fliplr(indices(end-nch+1:end)) fliplr(indices(1:nch))];

                pf_r(i,:) = r(i,indices);
                pf_pods(i,:) = pods(i,indices);
                pf_indices(i,:) = indices;
            end

            pf(offset,:) = {'Liabilities' pf_r pf_pods pf_indices};
            offset = offset + 1;
        end

        pf(:,offset:end) = [];
    end

    ds.Result = 'CrossEntropy';
    ds.ResultDate = now(); %#ok<TNOW1>
    ds.ResultAnalysis = @(ds)analyze_result(ds);
    ds.ResultSerial = sn;

    ds.BW = bw;
    ds.LGD = 1 - rr;
    ds.MD = md;
    ds.PW = pw;
    ds.RR = rr;
    ds.SEL = sel;

    ds.PortfolioComponents = pfc;
    ds.Portfolios = pf;

    all_label = [' (RR=' num2str(ds.RR * 100) '%, ' ds.PW ', ' ds.MD ')'];

    ds.LabelsMeasuresSimple = {'SI' 'SV' 'CoJPoDs'};
    ds.LabelsMeasures = {['SI' all_label] ['SV' all_label] ['CoJPoDs' all_label]};

    ds.LabelsIndicatorsSimple = {'JPoD' 'FSI' 'PCE'};
    ds.LabelsIndicators = {['JPoD' all_label] ['FSI' all_label] ['PCE' all_label]};

    ds.LabelsSheetsSimple = [{'Indicators' 'Average DiDe'} ds.LabelsMeasuresSimple];
    ds.LabelsSheets = [{'Indicators' 'Average DiDe'} ds.LabelsMeasures];

    ds.Indicators = NaN(t,numel(ds.LabelsIndicators));

    ds.AverageDiDe = NaN(nc);
    ds.DiDe = cell(t,1);
    ds.SI = NaN(t,nc);
    ds.SV = NaN(t,nc);

    ds.CoJPoDs = NaN(t,nc);

    ds.ComparisonReferences = {'Indicators' [] strcat({'CE-'},ds.LabelsIndicatorsSimple)};

end

function window_results = main_loop(r,pods,pw,md)

    window_results = struct();

    nan_indices = any(isnan(r),1);
    n = numel(nan_indices);

    r(:,nan_indices) = [];
    pods(:,nan_indices) = [];

    if (strcmp(pw,'A'))
        pods = mean(pods,1).';
    else
        [t,n] = size(pods);
        w = repmat(fliplr(((1 - 0.98) / (1 - 0.98^t)) .* (0.98 .^ (0:1:t-1))).',1,n);
        pods = sum(pods .* w,1).';
    end

    [g,p] = cimdo(r,pods,md);

    if (any(isnan(p)))
        window_results.JPoD = NaN;
        window_results.FSI = NaN;
        window_results.PCE = NaN;
        window_results.DiDe = NaN(n);
        window_results.SI = NaN(1,n);
        window_results.SV = NaN(1,n);
        window_results.CoJPoDs = NaN(1,n);
    else
        opods = pods;
        pods = NaN(n,1);
        pods(~nan_indices) = opods;

        [jpod,fsi,pce,dide,si,sv,cojpods] = cross_entropy_metrics(pods,g,p);
        window_results.JPoD = jpod;
        window_results.FSI = fsi;
        window_results.PCE = pce;
        window_results.DiDe = dide;
        window_results.SI = si;
        window_results.SV = sv;
        window_results.CoJPoDs = cojpods;
    end

end

function ds = finalize(ds,results)

    t = ds.T;
    k = size(results,2);

    if (k == 1)
        for i = 1:t
            result = results{i};

            ds.Indicators(i,:) = [result.JPoD result.FSI result.PCE];

            ds.DiDe{i} = result.DiDe;
            ds.SI(i,:) = result.SI;
            ds.SV(i,:) = result.SV;

            ds.CoJPoDs(i,:) = result.CoJPoDs;
        end
    else
        jpod_avg = mean(cellfun(@(x)x.JPoD,results,'UniformOutput',true),2,'omitnan');
        fsi_avg = mean(cellfun(@(x)x.FSI,results,'UniformOutput',true),2,'omitnan');
        pce_avg = mean(cellfun(@(x)x.PCE,results,'UniformOutput',true),2,'omitnan');
        ds.Indicators = [jpod_avg fsi_avg pce_avg];

        dide = cellfun(@(x)x.DiDe,results,'UniformOutput',false);
        si = cellfun(@(x)x.SI,results,'UniformOutput',false);
        sv = cellfun(@(x)x.SI,results,'UniformOutput',false);

        cojpods = cellfun(@(x)x.CoJPoDs,results,'UniformOutput',false);

        for i = 1:t
            ds.DiDe{i} = mean(cat(3,dide{i,:}),3);
            ds.SI(i,:) = mean(cat(3,si{i,:}),3);
            ds.SV(i,:) = mean(cat(3,sv{i,:}),3);

            ds.CoJPoDs(i,:) = mean(cat(3,cojpods{i,:}),3);
        end
    end

    ds.AverageDiDe = sum(cat(3,ds.DiDe{:}),3) ./ numel(ds.DiDe);

    si_vec = ds.SI(:);
    si_max = max(si_vec,[],'omitnan');
    si_min = min(si_vec,[],'omitnan');
    ds.SI = (ds.SI - si_min) ./ (si_max - si_min);

    sv_vec = ds.SV(:);
    sv_max = max(sv_vec,[],'omitnan');
    sv_min = min(sv_vec,[],'omitnan');
    ds.SV = (ds.SV - sv_min) ./ (sv_max - sv_min);

end

function write_results(ds,temp,out)

    [out_path,~,~] = fileparts(out);

    try
        if (exist(out_path,'dir') ~= 7)
            mkdir(out_path);
        end

        if (exist(out,'file') == 2)
            delete(out);
        end
    catch
        error('A system I/O error occurred while writing the results.');
    end

    copy_result = copyfile(temp,out,'f');

    if (copy_result == 0)
        error('The results file could not be created from the template file.');
    end

    dates_str = cell2table(ds.DatesStr,'VariableNames',{'Date'});

    labels_all = ds.PortfolioComponents;

    if (size(ds.Portfolios,1) == 1)
        if (strcmp(ds.SEL,'G'))
            header = {'Groups'};
        else
            header = {'Firms'};
        end
    else
        header = {'Components'};
    end

    labels = ds.LabelsIndicatorsSimple;
    tab = [dates_str array2table(ds.Indicators,'VariableNames',labels)];
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{1},'WriteRowNames',true);

    vars = [labels_all.' num2cell(ds.AverageDiDe)];
    tab = cell2table(vars,'VariableNames',[header labels_all]);
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{2},'WriteRowNames',true);

    tab = [dates_str array2table(ds.SI,'VariableNames',labels_all)];
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{3},'WriteRowNames',true);

    tab = [dates_str array2table(ds.SV,'VariableNames',labels_all)];
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{4},'WriteRowNames',true);

    tab = [dates_str array2table(ds.CoJPoDs,'VariableNames',labels_all)];
    writetable(tab,out,'FileType','spreadsheet','Sheet',ds.LabelsSheetsSimple{5},'WriteRowNames',true);

    worksheets_batch(out,ds.LabelsSheetsSimple,ds.LabelsSheets);

end

%% PLOTTING

function analyze_result(ds)

    if (size(ds.Portfolios,1) > 1)
        safe_plot(@(id)plot_portfolios_coverage(ds,id));
    end

    safe_plot(@(id)plot_indicators(ds,id));
    safe_plot(@(id)plot_dide(ds,id));
    safe_plot(@(id)plot_sequence_dide(ds,id));
    safe_plot(@(id)plot_sequence_cojpods(ds,id));

end

function plot_portfolios_coverage(ds,id)

    if (strcmp(ds.SEL,'G'))
        n = ds.Groups;
        labels = ds.GroupShortNames.';
    else
        n = ds.N;
        labels = ds.FirmNames;
    end

    nc = numel(ds.PortfolioComponents);

    k = size(ds.Portfolios,1);

    seq = 1:n;

    counts = cell(k + 1,2);
    count_total = zeros(n,1);

    for i = 1:k
        pf_indices = ds.Portfolios{i,4};
        pf_indices = pf_indices(:);

        [values_unique,~,indices_unique] = unique(pf_indices);

        count = zeros(n,1);
        count(values_unique) = accumarray(indices_unique,1);

        count_total = count_total + count;

        [count,order] = sort(count);
        counts(i + 1,:) = {(count ./ (sum(count) / nc)) labels(order)};
    end

    [count_total,order] = sort(count_total);
    counts(1,:) = {(count_total ./ (sum(count_total) / nc)) labels(order)};

    if (k == 4)
        spp = [2 3];
        spo = {[1 4] 2 3 5 6};
    elseif (k == 3)
        spp = [3 2];
        spo = {[1 5] 2 4 6};
    else
        spp = [2 2];
        spo = {[1 3] 2 4};
    end

    f = figure('Name','Cross-Entropy Measures > Reduced Portfolios Coverage','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    subs = gobjects(k + 1,1);

    sub_1 = subplot(spp(1),spp(2),spo{1});
    bar(sub_1,seq,counts{1,1},'FaceColor',[0.749 0.862 0.933]);
    set(sub_1,'XLim',[0 (ds.N + 1)],'XTick',seq,'XTickLabel',counts{1,2},'XTickLabelRotation',90);
    set(sub_1,'YLim',[0 1],'YTick',0:0.1:1,'YTickLabels',arrayfun(@(x)sprintf('%.f%%',x),(0:0.1:1) .* 100,'UniformOutput',false));
    title('Overall');

    subs(1) = sub_1;

    for i = 1:k
        sub = subplot(spp(1),spp(2),spo{i + 1});
        bar(sub,seq,counts{i + 1,1},'FaceColor',[0.749 0.862 0.933]);
        set(sub,'XLim',[0 (ds.N + 1)],'XTick',seq,'XTickLabel',counts{i + 1,2},'XTickLabelRotation',90);
        set(sub,'YLim',[0 1],'YTick',0:0.1:1,'YTickLabels',arrayfun(@(x)sprintf('%.f%%',x),(0:0.1:1) .* 100,'UniformOutput',false));
        title(ds.Portfolios{i,1});

        subs(i + 1) = sub;
    end

    set(subs,'YGrid','on');

    figure_title(f,'Reduced Portfolios Coverage');

    maximize_figure(f);

end

function plot_indicators(ds,id)

    nc = numel(ds.PortfolioComponents);

    jpod = smooth_data(ds.Indicators(:,1));
    fsi = smooth_data(ds.Indicators(:,2));
    pce = smooth_data(ds.Indicators(:,3));

    f = figure('Name','Cross-Entropy Measures > Indicators','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(2,2,1:2);
    plot(sub_1,ds.DatesNum,jpod);
    set(sub_1,'YLim',plot_limits(jpod,0.1,0));
    t1 = title(sub_1,ds.LabelsIndicatorsSimple{1});
    set(t1,'Units','normalized');
    t1_position = get(t1,'Position');
    set(t1,'Position',[0.4783 t1_position(2) t1_position(3)]);

    sub_2 = subplot(2,2,3);
    plot(sub_2,ds.DatesNum,fsi);
    set(sub_2,'YLim',[1 nc],'YTick',1:nc);
    title(sub_2,ds.LabelsIndicatorsSimple{2});

    sub_3 = subplot(2,2,4);
    plot(sub_3,ds.DatesNum,pce);
    set(sub_3,'YLim',[0 1],'YTick',0:0.1:1,'YTickLabels',arrayfun(@(x)sprintf('%.f%%',x),(0:0.1:1) .* 100,'UniformOutput',false));
    title(sub_3,ds.LabelsIndicatorsSimple{3});

    set([sub_1 sub_2 sub_3],'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set([sub_1 sub_2 sub_3],'XGrid','on','YGrid','on');

    if (ds.MonthlyTicks)
        date_ticks([sub_1 sub_2 sub_3],'x','mm/yyyy','KeepLimits','KeepTicks');
    else
        date_ticks([sub_1 sub_2 sub_3],'x','yyyy','KeepLimits');
    end

    figure_title(f,['Indicators (RR=' num2str(ds.RR * 100) '%, ' ds.PW ', ' ds.MD ')']);

    maximize_figure(f);

end

function plot_dide(ds,id)

    nc = numel(ds.PortfolioComponents);

    dide = ds.AverageDiDe;
    didev = dide(:);

    [dide_x,dide_y] = meshgrid(1:nc,1:nc);
    dide_x = dide_x(:) + 0.5;
    dide_y = dide_y(:) + 0.5;

    dide_txt = cellstr(num2str(didev,'~%.4f'));

    for i = 1:nc^2
        didev_i = didev(i);

        if (didev_i == 0)
            dide_txt{i} = '0';
        elseif (didev_i == 1)
            dide_txt{i} = '';
        end
    end

    lt_indices = (dide < 0.2) & (dide ~= 1);
    ge_indices = (dide >= 0.2) & (dide ~= 1);

    dide(lt_indices) =  0;
    dide(ge_indices) =  1;
    dide = dide - (eye(nc) .* 0.5);
    dide = padarray(dide,[1 1],'post');

    didev = dide(:);
    didev_ones = any(didev == 1);
    didev_zeros = any(didev == 0);

    f = figure('Name','Cross-Entropy Measures > Average Distress Dependency','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    pcolor(dide);

    if (didev_ones && didev_zeros)
        colormap([1 1 1; 0.65 0.65 0.65; 0.749 0.862 0.933]);
    else
        if (didev_ones)
            colormap([0.65 0.65 0.65; 0.749 0.862 0.933]);
        else
            colormap([1 1 1; 0.65 0.65 0.65]);
        end
    end

    text(dide_x,dide_y,dide_txt,'HorizontalAlignment','center');
    axis image;

    ax = gca();
    set(ax,'TickLength',[0 0]);
    set(ax,'XAxisLocation','top','XTick',1.5:(nc + 0.5),'XTickLabels',ds.PortfolioComponents,'XTickLabelRotation',45);
    set(ax,'YDir','reverse','YTick',1.5:(nc + 0.5),'YTickLabels',ds.PortfolioComponents,'YTickLabelRotation',45);

    figure_title(f,'Average Distress Dependency');

    maximize_figure(f);

end

function plot_sequence_dide(ds,id)

    nc = numel(ds.PortfolioComponents);
	nc_ones = ones(1,nc);

    t = ds.T;
    dn = ds.DatesNum;
    mt = ds.MonthlyTicks;

    ts_si = ds.SI;
    ts_sv = ds.SV;

    data = [repmat({dn},1,nc); mat2cell(ts_si,t,nc_ones); mat2cell(ts_sv,t,nc_ones)];

    sequence_titles = ds.PortfolioComponents;

    plots_title = [repmat(ds.LabelsMeasures(1),1,nc); repmat(ds.LabelsMeasures(2),1,nc)];

    x_limits = [dn(1) dn(end)];

    y_limits = [0 1];
    y_ticks = 0:0.1:1;

    core = struct();

    core.N = nc;
    core.Data = data;
    core.Function = @(subs,data)plot_function(subs,data);

    core.OuterTitle = 'Cross-Entropy Measures > Distress Dependency Time Series';
    core.InnerTitle = 'Distress Dependency Time Series';
    core.SequenceTitles = sequence_titles;

    core.PlotsAllocation = [2 1];
    core.PlotsSpan = {1 2};
    core.PlotsTitle = plots_title;

    core.XDates = {mt mt};
    core.XGrid = {true true};
    core.XLabel = {[] []};
    core.XLimits = {x_limits x_limits};
    core.XRotation = {[] []};
    core.XTick = {[] []};
    core.XTickLabels = {[] []};

    core.YGrid = {true true};
    core.YLabel = {[] []};
    core.YLimits = {y_limits y_limits};
    core.YRotation = {[] []};
    core.YTick = {y_ticks y_ticks};
    core.YTickLabels = {[] []};

    sequential_plot(core,id);

    function plot_function(subs,data)

        x = data{1};
        si = data{2};
        sv = data{3};

        d = min(find(isnan(si),1,'first'),find(isnan(sv),1,'first'));

        if (isempty(d))
            xd = [];
        else
            xd = x(d) - 1;
        end

		sub_1 = subs(1);
        plot(sub_1,x,si,'Color',[0.000 0.447 0.741]);

        if (~isempty(xd))
            hold(sub_1,'on');
                plot(sub_1,[xd xd],get(sub_1,'YLim'),'Color',[1.000 0.400 0.400]);
            hold(sub_1,'off');
        end

		sub_2 = subs(2);
        plot(sub_2,x,sv,'Color',[0.000 0.447 0.741]);

        if (~isempty(xd))
            hold(sub_2,'on');
                plot(sub_2,[xd xd],get(sub_2,'YLim'),'Color',[1.000 0.400 0.400]);
            hold(sub_2,'off');
        end

    end

end

function plot_sequence_cojpods(ds,id)

    nc = numel(ds.PortfolioComponents);

    t = ds.T;
    dn = ds.DatesNum;
    mt = ds.MonthlyTicks;

    ts = smooth_data(ds.CoJPoDs);

    data = [repmat({dn},1,nc); mat2cell(ts,t,ones(1,nc))];

    sequence_titles = ds.PortfolioComponents;

    plots_title = repmat(ds.LabelsMeasures(3),1,nc);

    x_limits = [dn(1) dn(end)];
    y_limits = plot_limits(ts,0.1,0);

    y_tick_labels = @(x)sprintf('%.2f%%',x .* 100);

    core = struct();

    core.N = nc;
    core.Data = data;
    core.Function = @(subs,data)plot_function(subs,data);

    core.OuterTitle = 'Cross-Entropy Measures > CoJPoDs Time Series';
    core.InnerTitle = 'CoJPoDs Time Series';
    core.SequenceTitles = sequence_titles;

    core.PlotsAllocation = [1 1];
    core.PlotsSpan = {1};
    core.PlotsTitle = plots_title;

    core.XDates = {mt};
    core.XGrid = {true};
    core.XLabel = {[]};
    core.XLimits = {x_limits};
    core.XRotation = {45};
    core.XTick = {[]};
    core.XTickLabels = {[]};

    core.YGrid = {true};
    core.YLabel = {[]};
    core.YLimits = {y_limits};
    core.YRotation = {[]};
    core.YTick = {[]};
    core.YTickLabels = {y_tick_labels};

    sequential_plot(core,id);

    function plot_function(subs,data)

        x = data{1};
        y = data{2};

        d = find(isnan(y),1,'first');

        if (isempty(d))
            xd = [];
        else
            xd = x(d) - 1;
        end

		sub_1 = subs(1);
        plot(sub_1,x,y,'Color',[0.000 0.447 0.741]);

        if (~isempty(xd))
            hold(sub_1,'on');
                plot(sub_1,[xd xd],get(sub_1,'YLim'),'Color',[1.000 0.400 0.400]);
            hold(sub_1,'off');
        end

    end

end

%% VALIDATION

function out = validate_output(out)

    [path,name,extension] = fileparts(out);

    if (~strcmpi(extension,'.xlsx'))
        out = fullfile(path,[name extension '.xlsx']);
    end

end

function sel = validate_selection(sel,groups)

    if (strcmp(sel,'G') && (groups == 0))
        error('The selection cannot be set to groups because their are not defined in the dataset.');
    end

end

function temp = validate_template(temp)

    sheets = {'Indicators' 'Average DiDe' 'SI' 'SV' 'CoJPoDs'};
    file_sheets = validate_xls(temp,'T');

    if (~all(ismember(sheets,file_sheets)))
        error(['The template must contain the following sheets: ' sheets{1} sprintf(', %s',sheets{2:end}) '.']);
    end

    worksheets_batch(temp,sheets);

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function h = imcredit(s)
%IMCREDIT Display image credit string
%   H = IMCREDIT(S) displays the string S at the bottom right of the image
%   displayed in the current axes.  S is wrapped to 30 characters and
%   displayed in a light gray small font.  H is the handle of the text
%   object created by IMCREDIT.
%
%   This function assumes the typical display created using imshow: a
%   single axes containing a single image, with no axes ticks or tick
%   labels.  The credit string is placed just below the bottom edge of the
%   axes box.
%
%   Note: IMCREDIT requires LINEWRAP, which is available on the MATLAB
%   Central File Exchange.
%
%   Example
%       imshow(rand(256,256))
%       imcredit('Image courtesy of Joe and Frank Hardy, MIT, 1993')
%
%   Possible enhancements:
%   * Adjust credit string position if axes tick labels are visible
%   * Adjust text color based on figure background color
%   * Additional syntaxes, including linewrap and text property options
%
%   See also LINEWRAP.

% Steven L. Eddins
% $Revision: 1.2 $  $Date: 2006/02/07 02:39:05 $

ax_handle = gca;
x = getRight(ax_handle);
y = getBottom(ax_handle);
text_handle = text(x, y, linewrap(s, 30), ...
   'VerticalAlignment', 'top', ...
   'HorizontalAlignment', 'right', ...
   'FontSize', 7, ...
   'Color', [.5 .5 .5]);

if nargout > 0
   h = text_handle;
end

%==========================================================================
function right = getRight(ax_handle)
% right = getRight(ax_handle) returns the x-coordinate, in axes data space,
% of the right edge of the axes box.

x_limits = get(ax_handle, 'XLim');
if strcmp(get(ax_handle, 'XDir'), 'normal')
   right = x_limits(2);
else
   right = x_limits(1);
end
%--------------------------------------------------------------------------

%==========================================================================
function bottom = getBottom(ax_handle)
% bottom = getBottom(ax_handle) returns the y-coordinate, in axes data
% space, of the bottom edge of the axes box.
y_limits = get(ax_handle, 'YLim');
if strcmp(get(ax_handle, 'YDir'), 'normal')
   bottom = y_limits(1);
else
   bottom = y_limits(2);
end
%--------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function [ lines, ori_lines ] = combineEdgesN(  edge )
%COMBINEEDGES Combine some small line segments, should be very conservative
%   lines: combined line segments
%   ori_lines: original line segments
%   line format: [nx ny nz projectPlaneID umin umax LSfov score]
arcList = [];
for i = 1:length(edge)
    panoLst = edge(i).panoLst;
    if size(panoLst,1) == 0
        continue;
    end
    arcList = [arcList; panoLst];
end

%% ori lines
numLine = size(arcList,1);
ori_lines = zeros(numLine,8);
areaXY = abs(sum(arcList(:,1:3).*repmat([0 0 1], [size(arcList,1) 1]),2));
areaYZ = abs(sum(arcList(:,1:3).*repmat([1 0 0], [size(arcList,1) 1]),2));
areaZX = abs(sum(arcList(:,1:3).*repmat([0 1 0], [size(arcList,1) 1]),2));
[~, planeIDs] = max([areaXY areaYZ areaZX], [], 2); % 1:XY 2:YZ 3:ZX

for i = 1:numLine
    ori_lines(i,1:3) = arcList(i,1:3);
    ori_lines(i,4) = planeIDs(i);
    coord1 = arcList(i,4:6);
    coord2 = arcList(i,7:9);
    uv = xyz2uvN([coord1; coord2], planeIDs(i));
    umax = max(uv(:,1))+pi;
    umin = min(uv(:,1))+pi;
    if umax-umin>pi
        ori_lines(i,5:6) = [umax umin]/2/pi;
%         ori_lines(i,7) = umin + 1 - umax;
    else
        ori_lines(i,5:6) = [umin umax]/2/pi;
%         ori_lines(i,7) = umax - umin;
    end
    ori_lines(i,7) = acos(dot(coord1, coord2, 2)./(norm(coord1)*norm(coord2)));
    ori_lines(i,8) = arcList(i,10);
end
% valid = ori_lines(:,3)<0;
% ori_lines(valid,1:3) = -ori_lines(valid,1:3);

%% additive combination
lines = ori_lines;
% panoEdge = paintParameterLine( lines, 1024, 512);
% figure; imshow(panoEdge);
for iter = 1:3
    numLine = size(lines,1);
    valid_line = true(numLine,1);
    for i = 1:numLine
%         fprintf('%d/%d\n', i, numLine);
        if ~valid_line(i)
            continue;
        end
        dotProd = dot(lines(:,1:3), repmat(lines(i,1:3), [numLine 1]), 2);
        valid_curr =  (abs(dotProd) > cos(1*pi/180)) & valid_line;
        valid_curr(i) = false;
        valid_ang = find(valid_curr);
        for j = valid_ang'
            range1 = lines(i,5:6);
            range2 = lines(j,5:6);
            valid_rag = intersection(range1, range2);
            if ~valid_rag
                continue;
            end

            % combine
            [~,I] = max(abs(lines(i,1:3)));
            if lines(i,I)*lines(j,I)>0
                nc = lines(i,1:3)*lines(i, 7) + lines(j,1:3)*lines(j, 7);
            else
                nc = lines(i,1:3)*lines(i, 7) - lines(j,1:3)*lines(j, 7);
            end
            nc = nc / norm(nc);

            if insideRange(range1(1), range2)
                nrmin = range2(1);
            else
                nrmin = range1(1);
            end
            if insideRange(range1(2), range2)
                nrmax = range2(2);
            else
                nrmax = range1(2);
            end

            u = [nrmin;nrmax]*2*pi - pi;
            v = computeUVN( nc, u, lines(i,4));
            xyz = uv2xyzN([u v], lines(i,4));
            len = acos(dot(xyz(1,:), xyz(2,:), 2));
            scr = (lines(i,7)*lines(i,8) + lines(j,7)*lines(j,8))/(lines(i,7)+lines(j,7));

            newLine = [nc lines(i,4) nrmin nrmax len scr];
            lines(i,:) = newLine;
            valid_line(j) = false;
        end
    end
    lines(~valid_line,:) = [];
    fprintf('iter: %d, before: %d, after: %d\n', iter, length(valid_line), sum(valid_line));
%     panoEdge = paintParameterLine( lines, 1024, 512);
%     figure; imshow(panoEdge);
end
%% previous method, bin voting
% %% build up voting space
% numDivision = 15;%round(max(width, height)/4);
%
% normal = arcList(:,1:3);
% valid = normal(:,3)<0;
% normal(valid,:) = -normal(valid,:);
% uv = xyz2uvN(normal, 1);
% thetas = uv(:,1);   phis = uv(:,2);
%
% uBinSize = 2*pi/(4*numDivision);
% vBinSize = pi/2/(numDivision);
% m = min(floor( (thetas-(-pi))/uBinSize) + 1, numDivision*4);
% n = min(floor( phis/vBinSize) + 1, numDivision);
% normalInd = sub2ind([4*numDivision numDivision], m, n);
%
% uniqueNormal = unique(normalInd);
% % decide voting dimension
% [m,n] = ind2sub([4*numDivision numDivision], uniqueNormal);
% u = -pi + (m-1)*uBinSize + uBinSize/2;
% v = 0   + (n-1)*vBinSize + vBinSize/2;
% uniNormal = [cos(v).*sin(u) cos(v).*cos(u) sin(v)];
% areaXY = abs(sum(uniNormal.*repmat([0 0 1], [size(uniNormal,1) 1]),2));
% areaYZ = abs(sum(uniNormal.*repmat([1 0 0], [size(uniNormal,1) 1]),2));
% areaZX = abs(sum(uniNormal.*repmat([0 1 0], [size(uniNormal,1) 1]),2));
% [~, planeIDs] = max([areaXY areaYZ areaZX], [], 2); % 1:XY 2:YZ 3:ZX
%
% subVoteBinNum = 1024;
% uvBin = false(length(uniqueNormal), subVoteBinNum);
% uBinSize = 2*pi/subVoteBinNum;
% for i = 1:length(uniqueNormal)
%     normIDs = find(normalInd==uniqueNormal(i));
%     norms = normal(normIDs,:);
%     rectNorm = sum(norms,1);
%     uniNormal(i,:) = rectNorm./norm(rectNorm);
%     dimIndc = planeIDs(i);
%     for j = normIDs'
%         coord1 = arcList(j,4:6);
%         coord2 = arcList(j,7:9);
%         xx = linspace(coord1(1),coord2(1),1000);
%         yy = linspace(coord1(2),coord2(2),1000);
%         zz = linspace(coord1(3),coord2(3),1000);
%         uv = xyz2uvN([xx' yy' zz'], dimIndc);
%
%         m = min(floor( (uv(:,1)-(-pi))/uBinSize) + 1, subVoteBinNum);
%         uvBin(i,m) = true;
%     end
% end
%
% %% extract line segments
% % numLines = 0;
% lines = [];
% for i = 1:length(uniqueNormal)
% %     fprintf('%d\n',i);
%     bins = int32(uvBin(i,:));
%     changePt(2:length(bins)) = bins(2:end) - bins(1:end-1);
%     changePt(1) = bins(1) - bins(end);
%     startPt = find(changePt==1);
%     endPt   = find(changePt==-1) - 1;
%     endPt = rem(endPt + subVoteBinNum+1, subVoteBinNum+1);
%
%     if endPt(1)>=startPt(1)
%         mtEndPt = endPt;
%     else
%         mtEndPt = [endPt(2:end) endPt(1)];
%     end
%
%     lines(end+1:end+length(startPt),:) = ...
%         [repmat([uniNormal(i,:) planeIDs(i)],[length(startPt) 1]) ...
%         (startPt'-1)/subVoteBinNum mtEndPt'/subVoteBinNum];
% end
end

function b = intersection(range1, range2)
    if range1(2)<range1(1)
        range11 = [range1(1) 1];
        range12 = [0 range1(2)];
    else
        range11 = range1;
        range12 = [0 0];
    end
    if range2(2)<range2(1)
        range21 = [range2(1) 1];
        range22 = [0 range2(2)];
    else
        range21 = range2;
        range22 = [0 0];
    end
    b = max(range11(1),range21(1))<min(range11(2),range21(2));
    if b
        return;
    end
    b2 = max(range12(1),range22(1))<min(range12(2),range22(2));
    b = b||b2;
end

function b = insideRange(pt, range)
if range(2)>range(1)
    b = pt>=range(1) && pt<=range(2);
else
    b1 = pt>=range(1) && pt<=1;
    b2 = pt>=0 && pt<=range(2);
    b = b1 || b2;
end

end

author: zouchuhang

repo: LayoutNet


text:
function pass = test_trigs(pref)
% Test some trigonometric functions

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 100*pref.cheb3Prefs.chebfun3eps;

f = {@cos, @sin, @tan, @cosh, @sinh, @tanh, @tand};

g = chebfun3(@(x,y,z) x.*y.^2.*z.^3);
for jj = 1:numel(f)
    h = f{jj}(g);
    exact = chebfun3(@(x,y,z) f{jj}(x.*y.^2.*z.^3));
    pass(jj) = norm(h - exact) < tol;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [y,D,nz] = ksvddenoise(params,msgdelta)
%KSVDDENOISE K-SVD denoising.
%  [Y,D] = KSVDDENOISE(PARAMS) denoises the specified (possibly
%  multi-dimensional) signal using K-SVD denoising. Y is the denoised
%  signal and D is the trained dictionary produced by K-SVD.
%
%  [Y,D] = KSVDDENOISE(PARAMS,MSGDELTA) specifies the frequency of message
%  printing during the process. MSGDELTA should be a positive number
%  representing the interval in seconds between messages. A zero or
%  negative value cancels all messages. Default is MSGDELTA=5.
%
%  [Y,D,NZ] = KSVDDENOISE(...) also returns the average number of non-zero
%  coefficients in the representations of the denoised blocks.
%
%
%  Required fields in PARAMS:
%  --------------------------
%
%    'x' - Noisy signal.
%      The signal to denoise (can be multi-dimensional). Should be of type
%      double, and (for PSNR computations) with values within [0,1] (to
%      specify a different range, see parameter 'maxval' below).
%
%    'blocksize' - Size of block.
%      Indicates the size of the blocks to operate on. Should be either an
%      array of the form [N1 N2 ... Np], where p is the number of
%      dimensions of x, or simply a scalar N, representing the square block
%      [N N ... N]. See parameter 'stepsize' below to specify the amount of
%      overlap between the blocks.
%
%    'dictsize' - Size of dictionary to train.
%      Specifies the number of dictionary atoms to train by K-SVD.
%
%    'psnr' / 'sigma' - Noise power.
%      Specifies the noise power in dB (psnr) or the noise standard
%      deviation (sigma), used to determine the target error for
%      sparse-coding each block. If both fields are present, sigma is used
%      unless the field 'noisemode' is specified (below). When specifying
%      the noise power in psnr, make sure to set the 'maxval' parameter
%      as well (below) if the signal values are not within [0,1].
%
%    'trainnum' - Number of training blocks.
%      Specifies the number of training blocks to extract from the noisy
%      signal for K-SVD training.
%
%
%  Optional fields in PARAMS:
%  --------------------------
%
%    'initdict' - Initial dictionary.
%      Specifies the initial dictionary for the K-SVD training. Should be
%      either a matrix of size NxL where N=(N1*N2*...*Np), the string
%      'odct' to specify the overcomplete DCT dictionary, or the string
%      'data' to initialize using random signal blocks. When a matrix is
%      specified for 'initdict', L must be >= dictsize, and in this case
%      the dictionary is initialized using the first dictsize columns from
%      initdict. By default, initdict='odct'.
%
%    'stepsize' -  Interval between neighboring blocks.
%      Specifies the interval (in pixels/voxels) between neighboring blocks
%      to denoise in the OMP denoising step. By default, all overlapping
%      blocks are denoised and averaged. This can be changed by specifying
%      an alternate stepsize, as an array of the form [S1 S2 ... Sp] (where
%      p is the number of dimensions of x). This sets the distance between
%      neighboring blocks to be Si in the i-th dimension. Stepsize can also
%      be a scalar S, corresponding to the step size [S S ... S]. Each Si
%      must be >= 1, and, to ensure coverage of the entire noisy signal,
%      size(x,i)-Ni should be a multiple of Si for all i. The default
%      stepsize is 1.
%
%    'iternum' - Number of K-SVD iterations.
%      Specifies the number of K-SVD training iterations to perform. If not
%      specified, the default is 10.
%
%    'maxval' - Maximal intensity value.
%      Specifies the range of the signal values. Used to determine the
%      noise standard deviation when the noise power is specified in psnr.
%      By default, the signal values are assumed to be within [0,1]. When
%      'maxval' is specified, this range changes to [0,maxval].
%
%    'memusage' - Memory usage.
%      This parameter controls memory usage of the function. 'memusage'
%      should be one of the strings 'high', 'normal' (default) or 'low'.
%      When 'memusage' is specified, both KSVD and OMPDENOISE are invoked
%      using their corresponding memusage setting. Additionallly, when set
%      to 'low', the use of the functions OMPDENOISE2 and OMPDENOISE3 to
%      accelerate denoising of 2-D and 3-D signals is disabled, reverting
%      to the function OMPDENOISE. Note that the 'low' setting will
%      significantly increase runtime.
%
%
%  Optional fields in PARAMS - advanced:
%  -------------------------------------
%
%    'noisemode' - Noise power mode.
%      Specifies whether the 'psnr' or 'sigma' fields should be used to
%      determine the noise power. This is useful when both fields are
%      present in PARAMS. 'noisemode' should be one of the string 'psnr' or
%      'sigma'. If it is not present, and both fields are specified,
%      'sigma' is used.
%
%    'gain' - Noise gain.
%      A positive value (usually near 1) controlling the target error for
%      sparse-coding each block. When gain=1, the target error is precisely
%      the value derived from the psnr/sigma fields. When gain is different
%      from 1, the target error is multiplied by this value. The default
%      value is gain = 1.15.
%
%    'lambda' - Weight of the input signal.
%      Specifies the relative weight attributed to the noisy input signal
%      in determining the output. The default value is 0.1*(maxval/sigma),
%      where sigma is the standard deviation of the noise. See function
%      OMPDENOISE for more information.
%
%    'maxatoms' - Maximal number of atoms.
%      This parameter can be used to specify a hard limit on the number of
%      atoms used to sparse-code each block. Default value is
%      prod(blocksize)/2, i.e. half the number of samples in a block. See
%      function OMP2 for more information.
%
%    'exact' - Exact K-SVD update.
%      Specifies whether the exact or approximate dictionary update should
%      be used in the K-SVD training. By default, the approximate
%      computation is used. However, specifying a nonzero value for 'exact'
%      causes the exact computation to be used instead. See function KSVD
%      for more information.
%
%
%   Summary of all fields in PARAMS:
%   --------------------------------
%
%   Required:
%     'x'                      signal to denoise
%     'blocksize'              size of block to process
%     'dictsize'               size of dictionary to train
%     'psnr' / 'sigma'         noise power in dB / standard deviation
%     'trainnum'               number of training signals
%
%   Optional (default values in parenthesis):
%     'initdict'               initial dictionary ('odct')
%     'stepsize'               distance between neighboring blocks (1)
%     'iternum'                number of training iterations (10)
%     'maxval'                 maximal intensity value (1)
%     'memusage'               'low, 'normal' or 'high' ('normal')
%     'noisemode'              'psnr' or 'sigma' ('sigma')
%     'gain'                   noise gain (1.15)
%     'lambda'                 weight of input signal (0.1*maxval/sigma)
%     'maxatoms'               max # of atoms per block (prod(blocksize)/2)
%     'exact'                  exact update instead of approximate (0)
%
%
%  References:
%  [1] M. Elad and M. Aharon, "Image Denoising via Sparse and Redundant
%      representations over Learned Dictionaries", the IEEE Trans. on Image
%      Processing, Vol. 15, no. 12, pp. 3736-3745, December 2006.
%
%  See also KSVD, OMPDENOISE, OMP2.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  August 2009

%%%%% parse input parameters %%%%%
addpath('ompbox');

x = params.x;
blocksize = params.blocksize;
trainnum = params.trainnum;
dictsize = params.dictsize;

p = ndims(x);
if (p==2 && any(size(x)==1) && length(blocksize)==1)
  p = 1;
end

% blocksize %
if (numel(blocksize)==1)
  blocksize = ones(1,p)*blocksize;
end

% maxval %
if (isfield(params,'maxval'))
  maxval = params.maxval;
else
  maxval = 1;
  params.maxval = maxval;
end

% gain %
if (isfield(params,'gain'))
  gain = params.gain;
else
  gain = 1.15;
  params.gain = gain;
end

% msgdelta %
if (nargin<2)
  msgdelta = 5;
end

verbose = 't';
if (msgdelta <= 0)
  verbose='';
  msgdelta = -1;
end

% initial dictionary %

if (~isfield(params,'initdict'))
  params.initdict = 'odct';
end

if (isfield(params,'initdict') && ischar(params.initdict))
  if (strcmpi(params.initdict,'odct'))
    params.initdict = odctndict(blocksize,dictsize,p);
  elseif (strcmpi(params.initdict,'data'))
    params = rmfield(params,'initdict');    % causes initialization using random examples
  else
    error('Invalid initial dictionary specified.');
  end
end

if (isfield(params,'initdict'))
  params.initdict = params.initdict(:,1:dictsize);
end

% noise mode %
if (isfield(params,'noisemode'))
  switch lower(params.noisemode)
    case 'psnr'
      sigma = maxval / 10^(params.psnr/20);
    case 'sigma'
      sigma = params.sigma;
    otherwise
      error('Invalid noise mode specified');
  end
elseif (isfield(params,'sigma'))
  sigma = params.sigma;
elseif (isfield(params,'psnr'))
  sigma = maxval / 10^(params.psnr/20);
else
  error('Noise strength not specified');
end

params.Edata = sqrt(prod(blocksize)) * sigma * gain;   % target error for omp
params.codemode = 'error';

params.sigma = sigma;
params.noisemode = 'sigma';

% make sure test data is not present in params
if (isfield(params,'testdata'))
  params = rmfield(params,'testdata');
end

%%%% create training data %%%

ids = cell(p,1);
if (p==1)
  ids{1} = reggrid(length(x)-blocksize+1, trainnum, 'eqdist');
else
  [ids{:}] = reggrid(size(x)-blocksize+1, trainnum, 'eqdist');
end
params.data = sampgrid(x,blocksize,ids{:});
params.data = remove_dc(params.data,'columns');

%%%%% KSVD training %%%%%

if (msgdelta>0)
  disp('KSVD training...');
end
D = ksvd(params,verbose,msgdelta);

%%%%%  denoise the signal  %%%%%

if (~isfield(params,'lambda'))
  params.lambda = maxval/(10*sigma);
end

params.dict = D;

if (msgdelta>0)
  disp('OMP denoising...');
end

% call the appropriate ompdenoise function
if (p>3  || (isfield(params,'memusage') && strcmpi(params.memusage,'low')))
  [y,nz] = ompdenoise(params,msgdelta);
elseif (p==1)
  [y,nz] = ompdenoise1(params,msgdelta);
elseif (p==2)
  [y,nz] = ompdenoise2(params,msgdelta);
else
  [y,nz] = ompdenoise3(params,msgdelta);
end

return;

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function cuboidhyp = omapobj2cuboidhyp(rpo, reglabel, regori, vp, OMAP_FACTOR)

vps{1} = OMAP_FACTOR * vp{1};
vps{2} = OMAP_FACTOR * vp{2};
vps{3} = OMAP_FACTOR * vp{3};

regoriobj = regori;
statsobj = regionprops(reglabel,'PixelIdxList','PixelList');

%%
rm = regionmask(vps, size(reglabel));
sr = segregion(rm, statsobj);

convexpair = find_convex_pair(rpo, sr, regoriobj);

%% get all bounding quadrilaterals of convex pairs
[r c v] = find(convexpair);

rl = unique([r; c]);
for i = 1:length(rl)
    boundingquad(rl(i)) = fit_bounding_quad(statsobj(rl(i)).PixelList, regori(rl(i)), vps);
end

% % display quadrilaterals
% figure;
% img = zeros(size(reglabel));
% for i = 1:length(rl)
%     img(statsobj(rl(i)).PixelIdxList) = 1;
%     imshow(img,[]);
%     hold on;
%     for j = 1:4
%     plot(boundingquad(rl(i)).pt(j,1), boundingquad(rl(i)).pt(j,2), ...
%         'x', 'Markersize', 10, 'Color', [1 0 0]);
%     text(boundingquad(rl(i)).pt(j,1), boundingquad(rl(i)).pt(j,2), ...
%         num2str(j), 'FontSize', 20, 'Color', [1 0 0]);
%     end
%     pause;
%     hold off
%     img(statsobj(rl(i)).PixelIdxList) = 0;
% end

if ~exist('boundingquad', 'var')
    cuboidhyp = [];
    return;
end

%%

partcubhyp = convexpair2partialcuboidhyp(convexpair, boundingquad);

% upscale
for i = 1:length(partcubhyp)
    for j = 1:8
        if ~isempty(partcubhyp(i).junc3(j).pt)
            partcubhyp(i).junc3(j).pt = partcubhyp(i).junc3(j).pt / OMAP_FACTOR;
        end
    end
end

%% get cuboidhyp from partial cuboids
for i = 1:length(partcubhyp)
    prm = cuboid_pts2prm_incomp(partcubhyp(i), vp);
    pts = cuboid_prm2pts(prm, vp);
    cuboidhyp(i) = cuboid_pts2hyp(pts, vp);
end

%% display...
% global omapmore;
% % figure;
% % img = zeros(size(reglabel));
% [row col] = find(convexpair);
% for i = 1:length(row)
%
%     rimg = disp_region(omapmore, reglabel, [row(i) col(i)], OMAP_FACTOR);
%
%     for j = 1:8
%         figure; imshow(rimg); hold on;
%
%         pt = cat(1, boundingquad(row(i)).pt);
%         for k=1:4, text(pt(k,1),pt(k,2),num2str(k), 'Color',[0 1 1]); end
%         pt = pt([1 2 4 3 1], :) / OMAP_FACTOR;
%         plot(pt(:,1), pt(:,2), 'm', 'LineWidth',2);
%
%         pt = cat(1, boundingquad(col(i)).pt);
%         for k=1:4, text(pt(k,1),pt(k,2),num2str(k), 'Color',[0 1 1]); end
%         pt = pt([1 2 4 3 1], :) / OMAP_FACTOR;
%         plot(pt(:,1), pt(:,2), 'm', 'LineWidth',2);
%
%
%         disp_cubes(cuboidhyp((i-1)*8+j), []);
%         pause;
%         close;
%     end
% end

%% discard invalid cuboids
valid = check_validity(cuboidhyp, vp);
cuboidhyp = cuboidhyp(valid);

%% drop thin cuboids
cuboidhyp = drop_thin_cuboids(cuboidhyp, vp);

%% this does not work yet.
% ch = convexpair2cuboidhyp(convexpair, statsobj, vp);

%%
function convexpair = find_convex_pair(rel, sr, regoriobj)
% rel{1}(i,j) = 1 if i is up, j is down
% rel{2}(i,j) = 1 if i is left, j is right
% rel{3}(i,j) = 1 if i is front, j is back
% sr[n x 4]: sr(i,j)=1 if segment i is in region j
% regoriobj(i): orientation of region
% convexpair = 1,2,3,4, or 5 depending on configuration of convex pair
%    = 1 if top & front & region3or4
%    = 2 if top & left & region4
%    = 3 if top & right & region3
%    = 4 if left & front & region4
%    = 5 if right & front & region3

convexpair = zeros(size(rel{1}));

[r c] = find(rel{1} | rel{2} | rel{3});
row = [r; c];
col = [c; r];

for i = 1:length(row)
    si = row(i);
    sj = col(i);
    orii = regoriobj(si);
    orij = regoriobj(sj);
    regioni = sr(si,:);
    regionj = sr(sj,:);
    if orii==1 && orij==3 && ... % top & front & region3or4
       rel{1}(si,sj)==1 && rel{3}(sj,si)==1 && ...
       (regioni(3)==1 || regioni(4)==1) && ...
       (regionj(3)==1 || regionj(4)==1) && ...
       (regioni(1)==0 && regioni(2)==0)
           convexpair(si,sj) = 1;
    elseif orii==1 && orij==2 && ... % top & left & region4
        rel{1}(si,sj)==1 && rel{2}(sj,si)==1 && ...
        regioni(4)==1 && regionj(4)==1 && ...
        regioni(1)==0 && regioni(2)==0
            convexpair(si,sj) = 2;
    elseif orii==1 && orij==2 && ... % top & right & region3
        rel{1}(si,sj)==1 && rel{2}(si,sj)==1 && ...
        regioni(3)==1 && regionj(3)==1 && ...
        regioni(1)==0 && regioni(2)==0
            convexpair(si,sj) = 3;
    elseif orii==2 && orij==3 && ... % left & front & region4
        rel{2}(si,sj)==1 && rel{3}(sj,si)==1 && ...
        regioni(4)==1 && regionj(4)==1
            convexpair(si,sj) = 4;
    elseif orii==2 && orij==3 && ... & right & front & region3
        rel{2}(sj,si)==1 && rel{3}(sj,si)==1 && ...
        regioni(3)==1 && regionj(3)==1
            convexpair(si,sj) = 5;
    end
end

%%
function sr = segregion(rm, stats)
% sr[n x 4]: sr(i,j)=1 if segment i is in region j

sr = false(length(stats), 4);
for i = 1:length(stats)
    sr(i,1) = any( rm{1}(stats(i).PixelIdxList) );
    sr(i,2) = any( rm{2}(stats(i).PixelIdxList) );
    sr(i,3) = any( rm{3}(stats(i).PixelIdxList) );
    sr(i,4) = any( rm{4}(stats(i).PixelIdxList) );
end

%%
function rel = convert_rpo_to_udlrfb(rpo, vp)
% rpo{1,2,3}(j,k) = 1 if region j is closer to camera than k (k is closer to vp than j)
% rel{1}(i,j) = 1 if i is up, j is down
% rel{2}(i,j) = 1 if i is left, j is right
% rel{3}(i,j) = 1 if i is front, j is back

rel = rpo;
if vp{1}(2) < 0
    rel{1} = rel{1}';
end
if vp{2}(1) < 0
    rel{2} = rel{2}';
end

%%
function rm = regionmask(vp, imgsize)
% %    1 | 2
% %   ---+---
% %    3 | 4
% ***** TODO: fix this. need to check various conditions of vp *****

if vp{3}(2) < 1
    upmask = false(imgsize(1),imgsize(2));
elseif vp{3}(2) > imgsize(1)
    upmask = true(imgsize(1),imgsize(2));
else
    [p1ext p2ext] = extline(vp{2}, vp{3}, imgsize(2), imgsize(1));

    if p1ext(1) < p2ext(1)
        x = [p1ext(1) p2ext(1) imgsize(2)+100 -100];
    else
        x = [p1ext(1) p2ext(1) -100 imgsize(2)+100];
    end
    y = [p1ext(2) p2ext(2) -100 -100];

    upmask = poly2mask(x,y,imgsize(1),imgsize(2));
end
if vp{3}(1) < 1
    leftmask = false(imgsize(1),imgsize(2));
elseif vp{3}(1) > imgsize(2)
    leftmask = true(imgsize(1),imgsize(2));
else
    [p1ext p2ext] = extline(vp{1}, vp{3}, imgsize(2), imgsize(1));
    if p1ext(2) < p2ext(2)
        y = [p1ext(2) p2ext(2) imgsize(1)+100 -100];
    else
        y = [p1ext(2) p2ext(2) -100 imgsize(1)+100];
    end
    x = [p1ext(1) p2ext(1) -100 -100];

    leftmask = poly2mask(x,y,imgsize(1),imgsize(2));
end

rm{1} =  upmask &  leftmask;
rm{2} =  upmask & ~leftmask;
rm{3} = ~upmask &  leftmask;
rm{4} = ~upmask & ~leftmask;


text:
function prior = priorParamInit(prior, varargin)

% PRIORPARAMINIT Prior model's parameter initialisation.
% FORMAT
% DESC wrapper function for initialising prior distributions parameters.
% ARG prior : structure to initialise.
% RETURN prior : initialised prior structure.
%
% SEEALSO : priorCreate
%
% COPYRIGHT : Neil D. Lawrence, 2003
% MODIFICATIONS: Andreas Damianou, Carl Henrik Ek, 2013

% PRIOR

prior = feval([prior.type 'PriorParamInit'], prior, varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%{
Authors:
Jonatas Lopes de Paiva
Claudio Fabiano Motta Toledo
Helio Pedrini

%}

function pop = arrangePop(pop, desc)

pop = nestedSortStruct(pop, 'fitness');

if(nargin > 1)
    if(desc)
        pop = fliplr(pop);
    end
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function dopnoist
%DOPNOIST Unit test for the function DOPNOISE.

%	O. Lemoine - March 1996.

N=128; Fs=100; D=12;

% The only value we can test is iflaw since y is a random variable.

% Pure tone for a fixed target
F0=25; V=0;
[y,iflaw]=dopnoise(N,Fs,F0,D,V);
d=abs(iflaw-iflaw(1));
if sum(any(d>sqrt(eps)))~=0,
 error('dopnoise test 1 failed');
end;

% Null signal for a non-emitting target
F0=0; V=50;
[y,iflaw]=dopnoise(N,Fs,F0,D,V);
d=abs(iflaw-iflaw(1));
if sum(any(d>sqrt(eps)))~=0,
 error('dopnoise test 2 failed');
end;

N=111; Fs=51; D=9;

% The only value we can test is iflaw since y is a random variable.

% Pure tone for a fixed target
F0=12; V=0;
[y,iflaw]=dopnoise(N,Fs,F0,D,V);
d=abs(iflaw-iflaw(1));
if sum(any(d>sqrt(eps)))~=0,
 error('dopnoise test 3 failed');
end;

% Null signal for a non-emitting target
F0=0; V=51;
[y,iflaw]=dopnoise(N,Fs,F0,D,V);
d=abs(iflaw-iflaw(1));
if sum(any(d>sqrt(eps)))~=0,
 error('dopnoise test 4 failed');
end;

author: HeLiangHIT


text:
ind_active = 1;
n_ima = 1;
active_images = 1;

author: strawlab

repo: MultiCamSelfCal


text:
function [N info energy V] = ibr_fuse_depths(D1, D2, vals)
%IBR_FUSE_DEPTHS  Fuse two disparity images using QPBO
%
%   [N info energy V] = ibr_fuse_depths(D0, D1, vals)
%
% Given two disparity maps and input data (image sequence, projection
% matrices and configuration parameters) this function will fuse the two
% maps into a single, lower energy disparity map, using QPBO.
%
%IN:
%   D0 - MxN disparity map to splice, assumed to be current best solution.
%   D1 - MxN disparity map to splice, assumed to be new proposal.
%   vals - structure containing the following input data and parameters:
%      R - MxNxC reference image.
%      I - Lx1 cell array of input images.
%      P - 4x3xL array of transposed projection matrices from reference to
%          input images.
%      SEI - (O+1)xV matrix of indices for smoothness cliques (columns),
%            where O is the order of the smoothness prior (1 or 2).
%      ephoto - handle to function which, when given TxC colour
%               differences, returns Tx1 values for Ephoto.
%      esmooth - handle to function which, when given Sx1 disparity
%                derivatives, returns Sx1 values for Esmooth.
%      d_min - scalar giving the minimum disparity in the scene.
%      d_step - scalar giving the range of disparities in the scene.
%      visibility - boolean indicating whether visibility constraints are
%                   to be used.
%      occl_val - scalar value for occluded pixels in the graph.
%      improve - method to use on unlabelled pixels: 0, fix to D0; 1, use
%                QPBOI; 2, use optimal splice on independent regions; 3,
%                fix to whichever of D0, D1 gives the lowest energy; 4,
%                transform labelling of method 2 using QPBOI.
%      contract - iterations of QPBOP to do.
%      independent - boolean indicating whether to use independent, or
%                    merely strongly-connected, regions for improve methods
%                    2 & 4.
%      show_output - figure handle to display output to, otherwise 0.
%
%OUT:
%   N - MxN logical array indicating fusing of D0 and D1.
%   info - structure containing the following useful info on the
%          optimization:
%      timings - 4x1 vector of cumulative times for [data_term_eval;
%                smoothness_term_eval; qpbo_fuse_time; finish_time];
%      numbers - 3x1 vector of values for [disps_set2D1; num_unlabelled_by_
%                qpbo; independent_unlabelled_regions];
%   energy - scalar value for the energy of the output disparity.
%   V - (M*N)xL logical array of visibilities of the input pixels, given
%       output disparity.

% $Id: ibr_fuse_depths.m,v 1.3 2008/11/17 11:27:35 ojw Exp $
t_start = cputime;

% Initialize values
Kinf = int32(0);
scale_factor = 1e5 / vals.ephoto(1e6);
occl_cost = cast(scale_factor * vals.occl_val, class(Kinf));
Kinf = occl_cost + 1; % Smaller value seems to avoid errors in QPBO
info.timings = zeros(3, 1, numel(vals.improve));
info.numbers = zeros(4, 1, numel(vals.improve), 'uint32');
energy = zeros(1, numel(vals.improve));
num_in = numel(vals.I);
sp = size(D1);
tp = numel(D1);
planar = size(vals.SEI, 1) == 3;
oobv = cast(-1000, class(vals.R));
out_unlabel = vals.improve(end) < 0;

% Calculate the homogenous coordinates of our two labellings
X = repmat(1:sp(2), [sp(1) 1]);
Y = repmat((1:sp(1))', [1 sp(2)]); % Faster than meshgrid
WC = [repmat([X(:) Y(:)], [2 1]) ones(2*tp, 1) [D1(:); D2(:)]];
clear X Y

% Initialise arrays for the data terms
if vals.visibility
    EI = reshape(repmat(uint32(tp+(0:num_in-1)*2*tp), [2*tp 1]), 1, []);
    EI = [repmat(uint32(1:tp), [1 2*num_in]); EI+repmat(uint32(1:2*tp), [1 num_in])];
    E = repmat(cat(3, [occl_cost 0 0 0], [0 0 occl_cost 0]), [tp 1 1 num_in]);
else
    % Data edges are not needed
    EI = zeros(2, 0, 'uint32');
    E = zeros(4, 0, class(Kinf));
end
U = zeros(tp, 2, class(Kinf));
TEI = zeros(2, 0, 'uint32');
TE = zeros(4, 0, class(Kinf));

% For each input image...
for a = 1:num_in
    % Calculate the coordinates in the input image
    T = WC * vals.P(:,:,a);
    N = 1 ./ T(:,3);
    T(:,1) = T(:,1) .* N;
    T(:,2) = T(:,2) .* N;

    % Calculate photoconsistency
    M = vgg_interp2(vals.I{a}, T(:,1), T(:,2), 'linear', oobv);
    M = squeeze(M) - vals.R;
    IA = reshape(cast(scale_factor * vals.ephoto(M), class(Kinf)), tp, 2);
    clear M N

    % Find interactions
    T(:,3) = T(:,3) ./ WC(:,4);
    [T M] = sortrows(T);
    N = find_interactions(T, 0.5); % Optimized version
    N = M(N); % Unsort
    N = uint32(N(:,abs(diff(N))~=tp)); % Remove interactions between the same node

    % Add the pixel interactions to the graph
    M = N(1,:) > tp;
    TEI = [TEI [N(1,:)-uint32(tp*M); (tp*2*a-tp)+N(2,:)]];
    T = zeros(4, numel(M));
    T(2,~M) = Kinf;
    T(4,M) = Kinf;
    TE = [TE T];

    if vals.visibility
        % Set up photoconsistency edges
        E(:,2,1,a) = IA(:,1);
        E(:,4,2,a) = IA(:,2);

        if vals.compress_graph
            % Determine the photoconsistency nodes which have no interactions
            M = ones(tp, 2, class(Kinf));
            M(N(2,:)) = 0;

            % Add those photoconsistency terms to the unaries
            U = U + M .* IA;
        end
    else
        % Add the photoconsistency terms to the unaries
        U = U + IA;
    end
    clear T M N
end
clear WC IA

EI_ = EI;
if vals.visibility
    E = reshape(permute(E, [2 1 3 4]), 4, []);
    if vals.compress_graph
        % The unary and pairwise energies as they stand are entirely
        % correct, i.e. will give the correct labelling. However, it can be
        % compressed into a smaller but equivalent graph, which will be
        % faster to solve, by removing superfluous nodes and edges.
        [U EI EI_ E N T] = compress_graph(U', EI, E, TEI, TE, tp, num_in);
    else
        U = zeros(2, tp+tp*2*num_in, class(Kinf));
        T = TE;
        N = TEI;
    end
    % Concatenate data and visibility edges
    E = [E T];
    EI_ = [EI_ N];
    clear T N
else
    U = U';
end
TE = TE(2,:) ~= 0;
info.timings(1,:) = cputime - t_start; % Time data term evaluation

% Add surface smoothness constraints
SE = (double([D1(:)'; D2(:)']) - vals.d_min) / vals.d_step;
SE = SE(:,vals.SEI);
if planar
    % Planar prior - Finite differences 2nd derivative of disparity
    SE = reshape(SE, 6, []);
    SE = reshape(SE([1 3 5; 1 3 6; 1 4 5; 1 4 6; 2 3 5; 2 3 6; 2 4 5; 2 4 6]',:), 3, 8, []);
    SE = diff(SE, 2);
else
    % Fronto-parallel prior - Finite differences 1st derivative of
    % disparity
    SE = reshape(SE, 4, []);
    SE = diff(reshape(SE([1 3; 1 4; 2 3; 2 4]',:), 2, 4, []));
end
% Apply our smoothness weighting
SE = reshape(cast(scale_factor * vals.esmooth(SE(:)), class(E)), [], size(vals.SEI, 2));
if ~planar
    E = [E SE];
    EI_ = [EI_ vals.SEI];
end
info.timings(2,:) = cputime - t_start; % Time smoothness term evaluation

for a = 1:numel(vals.improve)
    t_start = cputime;
    % Fuse the two labellings, using contract and/or improve if desired
    qpbo_params = int32([tp ((vals.improve(a)==1)+(vals.improve(a)==4)*2) vals.contract(a) vals.contract(a)>0]);
    if vals.improve(a) == 4
        % Add callback function handle
        qpbo_params = {qpbo_params, @(L) (choose_labels(L, U, E, EI, SE, vals.SEI, TE, TEI, num_in, vals.visibility, 2, vals.independent) > 0)};
    end
    try
        if planar
            [M stats] = vgg_qpbo(U, EI_, E, vals.SEI, SE, qpbo_params);
        else
            [M stats] = vgg_qpbo(U, EI_, E, qpbo_params);
        end
    catch
        % Error probably due to probe failure
        stats = [0 0 Inf];
        M = false(tp, 1);
    end
    clear qpbo_params
    info.numbers(2:4,a) = stats;

    if stats(1) && vals.improve(a) >= 2 && vals.improve(a) <= 3
        if nargout > 2 || vals.show_output
            [M info.numbers(3,a) U_ E_ SE_ V] = choose_labels(M, U, E(:,1:size(EI, 2)), EI, SE, vals.SEI, TE, TEI, num_in, vals.visibility, vals.improve(a), vals.independent);
            energy(a) = sum(U_) + sum(E_) + sum(SE_);
        else
            [M info.numbers(3,a)] = choose_labels(M, U, E(:,1:size(EI, 2)), EI, SE, vals.SEI, TE, TEI, num_in, vals.visibility, vals.improve(a), vals.independent);
        end
        N = M > 0;
    elseif nargout > 2 || vals.show_output
        N = M > 0;
        [U_ E_ SE_ V] = calc_vis_energy(N, U, E(:,1:size(EI, 2)), EI, SE, vals.SEI, TE, TEI, num_in);
        energy(a) = sum(U_) + sum(E_) + sum(SE_);
    end
    info.numbers(1,a) = sum(N);
    info.timings(3,a) = cputime - t_start + info.timings(2,a); % Time optimization
end
clear TEI TE U E SE EI_

if nargout > 3 || vals.show_output
    % Generate output visibilities
    T = (tp * N) + (1:tp)';
    for b = 1:num_in
        V(1:tp,b) = V(T);
        T = T + 2*tp;
    end
    V(tp+1:end,:) = [];
end

if vals.show_output
    % Display the output figures
    U_ = double(U_) + accum(EI(1,:)', E_, [tp 1]);
    U_ = reshape(U_, sp(1), sp(2));
    if vals.visibility
        % Take off the occlusion costs and normalize
        EI = reshape(sum(V(1:tp,:), 2), sp(1), sp(2));
        U_ = U_ - (num_in - EI) * double(occl_cost);
        E_ = EI ~= 0;
        U_(E_) = U_(E_) ./ EI(E_);
    end
    set(0, 'CurrentFigure', vals.show_output);
    subplot('Position', [1/3 0.5 1/3 0.5]);
    D1(N) = D2(N);
    sc(D1, 'contrast', vals.d_min+[0 vals.d_step]);
    subplot('Position', [2/3 0.5 1/3 0.5]);
    sc(U_, 'jet');
    subplot('Position', [0 0 1/3 0.5]);
    T = reshape(sc(reshape(M, sp(1), sp(2)), 'prism'), [], 3);
    I = M < 0;
    T(I,:) = 1 - (1 - T(I,:)) * 0.3; % Lighten unlabelled pixels set to 0
    I = M > 1;
    T(I,:) = T(I,:) .* 0.3; % Darken unlabelled pixels set to 1 by optimal splice
    T(M==0,:) = 1;
    T(M==1,:) = 0;
    sc(reshape(T, sp(1), sp(2), 3), [0 1]);
    subplot('Position', [1/3 0 1/3 0.5]);
    sc(reshape(sum(V, 2), sp(1), sp(2)), [0 num_in], 'contrast');
    subplot('Position', [2/3 0 1/3 0.5]);
    U_ = -accum(vals.SEI(2,:)', SE_, [tp 1]);
    sc(reshape(U_, sp(1), sp(2)));
    drawnow;
end
if out_unlabel
    N = M; % Output unlabelled pixels
end
return

function [M num_regions U_ E_ SE_ V] = choose_labels(M, U, E, EI, SE, SEI, TE, TEI, num_in, visibility, improve, independent)
% Calculate visibilities and regions assuming unlabelled pixels are set
% to 0 then 1.
[U_ E_ SE_ V] = calc_vis_energy(M==1, U, E, EI, SE, SEI, TE, TEI, num_in);
[U2 E2 SE2 V2] = calc_vis_energy(M~=0, U, E, EI, SE, SEI, TE, TEI, num_in);
num_regions = double(-min(M(:)));

if improve == 2
    % We want to do optimal splice.
    sz = [num_regions 1];
    % Merge strongly connected regions that are connected by smoothness
    % cliques
    SEI2 = M(SEI);
    SEI2 = SEI2(:,any(SEI2 < 0));
    SEI2 = SEI2(:,~all(SEI2 >= 0 | ojw_bsxfun(@eq, SEI2, min(SEI2))));
    while independent && ~isempty(SEI2)
        num_regions = num_regions - 1;
        T = SEI2(:,1);
        N = min(T);
        T = T(T < 0 & T ~= N);
        M(M==T(1)) = N;
        SEI2(SEI2==T(1)) = N;
        SEI2 = SEI2(:,~all(SEI2 >= 0 | ojw_bsxfun(@eq, SEI2, min(SEI2))));
    end
    if visibility
        % Merge strongly connected regions that are connected by visibility
        % cliques
        TEI2 = TEI(:,M(TEI(1,:))<0);
        tp = numel(M);
        M = repmat(M, [1+2*num_in 1]);
        M(TEI2(2,:)) = M(TEI2(1,:));
        SEI2 = M(TEI);
        SEI2 = SEI2(:,any(SEI2 < 0));
        SEI2 = SEI2(:,~all(SEI2 >= 0 | ojw_bsxfun(@eq, SEI2, min(SEI2))));
        while independent && ~isempty(SEI2)
            num_regions = num_regions - 1;
            T = SEI2(:,1);
            N = min(T);
            T = T(T < 0 & T ~= N);
            M(M==T(1)) = N;
            SEI2(SEI2==T(1)) = N;
            SEI2 = SEI2(:,~all(SEI2 >= 0 | ojw_bsxfun(@eq, SEI2, min(SEI2))));
        end
        % Go through each region and determine whether a labelling of 1 or 0
        % gives a lower energy, starting with the visibility edges
        EI2 = min(M(EI))';
        T = EI2 < 0;
        E2 = E2(T) - E_(T);
        engy = accum(-EI2(T), E2, sz);
        M = M(1:tp);
    else
        engy = zeros(sz);
    end
    % Go through each region and determine whether a labelling of 1 or 0
    % gives a lower energy
    SEI2 = M(SEI);
    T = any(SEI2 < 0);
    SEI2 = SEI2(:,T);
    SE2 = SE2(T) - SE_(T);
    T = -min(SEI2);
    engy = engy + accum(T(:), SE2, sz);
    T = M < 0;
    U2 = U2(T) - U_(T);
    engy = engy + accum(-M(T), U2, sz);
    update = false;
    for b = 1:sz(1)
        if engy(b) < 0
            M(M==-b) = b + 1;
            update = true;
        end
    end
    if update && nargout > 2
        % Recalculate the visibilities and energies
        [U_ E_ SE_ V] = calc_vis_energy(M > 0, U, E, EI, SE, SEI, TE, TEI, num_in);
    end
else
    % Just choose the label which gives the lowest energy
    if (sum(U_) + sum(E_) + sum(SE_)) >= (sum(U2) + sum(E2) + sum(SE2))
        % Update the labelling
        T = M < 0;
        M(T) = 1 - M(T);
        % Update the energy and visibility
        U_ = U2;
        E_ = E2;
        SE_ = SE2;
        V = V2;
    end
end
return

function [U E SE V] = calc_vis_energy(L, U, E, EI, SE, SEI, TE, TEI, num_in)
% Generate visibility maps
tp = numel(L);
V = true(2*tp, num_in);
V(TEI(2,L(TEI(1,:))~=TE')-tp) = false;

% Calculate energies
U = U((0:tp-1)'*2+L+1);
L = [L; V(:)];
E = E((0:size(EI, 2)-1)'*4+L(EI(1,:))*2+L(EI(2,:))+1);
if size(SEI, 1) == 3
    SE = SE((0:size(SEI, 2)-1)'*8+L(SEI(1,:))*4+L(SEI(2,:))*2+L(SEI(3,:))+1);
else
    SE = SE((0:size(SEI, 2)-1)'*4+L(SEI(1,:))*2+L(SEI(2,:))+1);
end
return

function [U EI EI_ E TEI TE] = compress_graph(U, EI, E, TEI, TE, tp, num_in)
% Count the number of interactions per input sample
SE = accum(TEI(2,:)', (1:size(TEI, 2))', [tp+tp*2*num_in 1], @num_first);
SE = SE(tp+1:end);

% Remove single interactions, attaching the photoconsitency edge
% directly to the interacting pixel
M = find(SE > 0);
L = SE(M);
EI(2,M) = TEI(1,L);
M = M(TE(4,L)~=0);
E(:,M) = E([2 1 4 3],M);
TEI(:,L) = [];
TE(:,L) = [];

% Remove the superfluous edges - photoconsistency edges with no
% interactions, that have already been incorporated into the unary term
M = SE ~= 0;
E = E(:,M);
EI = EI(:,M);

% Compress the node indices
M = zeros(tp+2*tp*num_in, 1, 'uint32');
SE = SE < 0;
L = sum(SE);
M([true(tp, 1); SE]) = uint32(1):uint32(L+tp);
EI_ = EI;
EI_(2,:) = M(EI(2,:));
TEI(2,:) = M(TEI(2,:));
U = [U zeros(2, L, class(U))];
return

function B = accum(I, A, sz, varargin)
% Older versions of Matlab can't accumulate integer arrays!
try
    B = accumarray(I, A, sz, varargin{:});
catch
    if isempty(I)
        B = zeros(sz);
    else
        B = accumarray(double(I), double(A), sz, varargin{:});
    end
end
return

function B = num_first(A)
% Return:        A  if numel(A) == 1
%         -numel(A) otherwise
B = -numel(A);
if B == -1
    B = A;
end
return

repo: SelfExSR


text:
function xfat = r8mat_expand_linear ( m, n, x, mfat, nfat )

%*****************************************************************************80
%
%% R8MAT_EXPAND_LINEAR linearly interpolates new data into an R8MAT.
%
%  Discussion:
%
%    In this routine, the expansion is specified by giving the number
%    of intermediate values to generate between each pair of original
%    data rows and columns.
%
%    The interpolation is not actually linear.  It uses the functions
%
%      1, x, y, and xy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of input data.
%
%    Input, real X(M,N), the original data.
%
%    Input, integer MFAT, NFAT, the number of data values to interpolate
%    between each row, and each column, of original data values.
%
%    Output, real XFAT(M2,N2), the fattened data, where
%    M2 = (M-1)*(MFAT+1)+1,
%    N2 = (N-1)*(NFAT+1)+1.
%
  for i = 1 : m

    if ( i < m )
      ihi = mfat;
    else
      ihi = 0;
    end

    for j = 1 : n

      if ( j < n )
        jhi = nfat;
      else
        jhi = 0;
      end

      if ( i < m )
        ip1 = i+1;
      else
        ip1 = i;
      end

      if ( j < n )
        jp1 = j+1;
      else
        jp1 = j;
      end

      x00 = x(i,j);
      x10 = x(ip1,j);
      x01 = x(i,jp1);
      x11 = x(ip1,jp1);

      for ii = 0 : ihi

        s = ii / ( ihi + 1 );

        for jj = 0 : jhi

          t = jj / ( jhi + 1 );

          iii = 1 + ( i - 1 ) * ( mfat + 1 ) + ii;
          jjj = 1 + ( j - 1 ) * ( nfat + 1 ) + jj;

          xfat(iii,jjj) = ...
                                            x00   ...
              + s     * (       x10       - x00 ) ...
              + t     * (             x01 - x00 ) ...
              + s * t * ( x11 - x10 - x01 + x00 );

        end

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ySqrtPred,PInvSqrtPred,RwPred,RwxPred]=ESRIFDiscPred(ySqrtPrev,PInvSqrtPrev,f,FJacob,SQ,u,Gamma)
%%ESRIFDISCPRED Perform the discrete-time prediction step that comes with
%               the extended square root information filter with additive
%               process noise.
%
%INPUTS: ySqrtPrev The xDimX1 square root information state that is to be
%                  propagated. The previous information state is always
%                  PInvSqrtPrev times the previous target state estimate.
%     PInvSqrtPrev The previous inverse square root information matrix.
%                  If P is the covariance matrix of a Gaussian state x,
%                  then P=PSqrt*PSqrt' and PInvSqrtPrev=inv(PSqrt). This
%                  can be either upper triangular or lower triangular.
%                f A function handle for the state transition function
%                  that takes the state as its parameter.
%           FJacob A function handle for calculating the xDim X xDim state
%                  transition matrix. If an empty matrix is passed, then
%                  FJacob will be found using numerical differentiation
%                  via the numDiff function with default parameters.
%               SQ The lower-triangular square root of the process noise
%                  covariance matrix. In the standard linear dynamic
%                  model, one has
%                  x(k)=F(k-1)*x(k-1)+u(k-1)+noise where x is the state
%                  and Q is the covariance matrix of the noise. However,
%                  if the covariance matrix of the noise in such equation
%                  is singular, then one should rewrite it as
%                  x(k)=F(k-1)*x(k-1)+u(k-1)+Gamma*noise
%                  where gamma is some matrix and the covariance matrix of
%                  the untransformed noise is not singular. This SQ is the
%                  lower triangular square root of the noise in the linear
%                  dynamic equation.
%                u An optional xDim X1 vector that is the control input.
%                  If omitted, a zero control input (no control input) is
%                  used.
%            Gamma An optional matrix that transforms the process noise
%                  to the state domain if the process noise covariance
%                  matrix is singular, as disccused for the input SQ. If
%                  this is omitted an identity matrix is used (i.e. there
%                  is no Gamma).
%
%OUTPUTS:ySqrtPred The xDim X 1 predicted square root information state
%                  vector.
%     PInvSqrtPred The predicted xDim X xDim inverse square root state
%                  covariance matrix, which is upper-triangular.
%   RwPred,RwxPred  Noise matrices which can be saved for use in smoothing
%
%The time prediction step is taken from the algorithmic implementation
%of [1] described in chapters V and VI.
%
%Given a Gaussian predicted state with mean x and covariance matrix P, the
%square root information state is
%ySqrt=PInvSqrt*x
%where
%PSqrt=inv(PInvSqrt)
%and
%P=PSqrt*PSqrt';
%The matrix PInvSqrtPred can be upper or lower triangular, when supplied to
%this function. For example, a lower-triangular matrix can be obtained
%using PInvSqrt=inv(chol(P,'lower')). However, the output of this function,
%PInvSqrtUpdate is always upper triangular.
%
%REFERENCES:
%[1] G. J. Bierman, "Factorization Methods for Discrete Sequential
%    Estimation. Academic Press, New York, 1977.
%
%February 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(ySqrtPrev,1);
if(nargin<6||isempty(u))
    u=zeros(xDim,1);
end
if(nargin<7||isempty(Gamma))
    Gamma=eye(xDim);
end
if(isempty(FJacob))
    FJacob=@(x)numDiff(x,f,xDim);
end

xPrev=PInvSqrtPrev\ySqrtPrev;

F=FJacob(xPrev);
PInvSqrtTilde=PInvSqrtPrev/F;

%The mapping matrix found on p.121 of Bierman
A=[inv(SQ),                zeros(xDim,xDim), SQ\u;
    -PInvSqrtTilde*Gamma,   PInvSqrtTilde,   ySqrtPrev];
[~,T] = qr(A);

PInvSqrtPred=T((xDim+1):end,(end-xDim):(end-1));

%Since f may be nonlinear, the information state output described in
%Bierman is not valid, and the predicted state should be calculated
%independently.
ySqrtPred=PInvSqrtPred*f(xPrev);

if(nargout>2)
    RwPred=T(1:xDim,1:xDim);
    RwxPred=T(1:xDim,xDim+1:end-1);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function rs = corrdim(s, bins)

%tstoolbox/@signal/corrdim
%   Syntax:
%     * rs = corrdim(s, bins)
%
%   Input arguments:
%     * s - data points (row vectors)
%     * bins - maximal number of partition per axis (optional)
%
%   Compute the correlation dimension of a time-delay reconstructed
%   timeseries s for dimensions from 1 to D, where D is the dimension of
%   the input vectors using boxcounting approach. The default number of
%   bins is 100.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,2);

if ndim(s) ~= 2
    error('Signal must contain vector data');
end

if nargin<2
    bins = 100;
end

points = data(s);
[N,dim] = size(points);

% scale data to be within 0 and 1
points = points - min(min(points));
points = points / max(max(points));

% give a sortiment of (integer) partitionsizes with almost exponential behaviour
par = [2 3 4 5 6 7 8 10 12 14 16 20 23 27 32 39 46 54 64 77 91 108 128 153 182 216 256 ...
      305 363 431 512 609 725 862 1024 1218 1449 1723 2048 2436 2897 3445 4096 4871 ...
	  5793 6889 8192 9742 11586 13778 16384 19484 23171 27555 32768 38968 46341 55109 65536];

partitions = par(find(par<=bins));	% use no sizes greater than bins

[dummy,dummy2,e] = boxcount(points, partitions);

e = [zeros(1,dim) ; e];             % add zeros for partition size 1
partitions = [1 ; partitions(:)];	% add partition size 1

rs = signal(core(e), s);
a = achse(-log2(partitions));     		% create axis with arbitrary spacing
a = setname(a, 'ld r');
a2 = setname(achse(unit, 1, 1), 'Embedding dimension');
rs = setaxis(rs, 1, a);
rs = setaxis(rs, 2, a2);
rs = setplothint(rs, 'multigraph');
rs = addhistory(rs,  ['Computed correlation dimension']);
rs = addcommandlines(rs, 's = corrdim(s', bins);
rs = setyname(rs, 'ld C(r)');
rs = setlabel(rs, 'Scaling of D2');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function y = gpReversibleDynamicsSamp(model, X);

% GPREVERSIBLEDYNAMICSSAMP Sample from the dynamics for a given input.

% FGPLVM

persistent oldX
if isempty(oldX)
  Xp = [X zeros(size(X))];
else
  Xp = [X X-oldX];
end
oldX = X;
[mu, var] = gpPosteriorMeanVar(model, Xp);
y = gsamp(mu, diag(var), 1);
y = X + y;


text:
clear;
close all
addpath(genpath(pwd));
for i=1:10
    I1 = imread(['../ir/',num2str(i),'.bmp']);
    I2 = imread(['../vis/',num2str(i),'.bmp']);
    J(:,:,1) = I1;
    J(:,:,2) = I2;
    F_echo_dtf = IJF(I1,I2);
    imwrite(F_echo_dtf,[num2str(i),'.bmp']);
    clear J;
end

Q_echo_dtf = Qp_ABF(I1, I2, F_echo_dtf)

author: thfylsty


text:
function [upper,x_min,info_text,numglobals,timing,p_upper] = solve_upper_in_node(p,p_upper,x,upper,x_min,uppersolver,info_text,numglobals,timing,cutiterations);

if nargin < 10
    cutiterations = 1;
end

% Call nonlinear solver (if there are SDP cones and non-SDP solver, or user
% specified use of nonlinear SDP cut generation, remove SDP cone)
[output,timing] = global_solve_upper(pruneSDPCone(p),pruneSDPCone(removeCuts(p_upper)),x,p.options,uppersolver,timing);
output.Primal(p_upper.integer_variables) = round(output.Primal(p_upper.integer_variables));
output.Primal(p_upper.binary_variables) = round(output.Primal(p_upper.binary_variables));

z = apply_recursive_evaluation(p_upper,output.Primal);

[upper_residual,~,feasibility] = constraint_residuals(p_upper,z);
feasible = ~isempty(z) & ~any(isnan(z)) & all(upper_residual(1:p_upper.K.f)>=-p.options.bmibnb.eqtol) & all(upper_residual(1+p_upper.K.f:end)>=-p.options.bmibnb.pdtol);
for i = 1:length(p.evalMap)
    if any(p.evalMap{i}.properties.forbidden)
        if any(z(p.evalMap{i}.variableIndex) >p.evalMap{i}.properties.forbidden(1) & z(p.evalMap{i}.variableIndex) < p.evalMap{i}.properties.forbidden(2))
            feasible = 0;
            break
        end
    end
end
if feasible
    this_upper = p_upper.f+p_upper.c'*z+z'*p_upper.Q*z;
    if isequal(p_upper.K.f,0) && isequal(p_upper.K.l,0) && isequal(p_upper.K.q,0) && isequal(p_upper.K.s,0)
        % This is a simple bound-constrained problem. We can try a
        % cross-over and see if it improves the objective
        xx = output.Primal;
        j = find(abs(xx-p_upper.lb)<=p.options.bmibnb.vartol);
        xx(j) = p_upper.lb(j);
        j = find(abs(xx-p_upper.ub)<=p.options.bmibnb.vartol);
        xx(j) = p_upper.ub(j);
        alt_z = apply_recursive_evaluation(p_upper,xx);
        alt_upper = p_upper.f + p_upper.c'*alt_z + alt_z'*p_upper.Q*alt_z;
        if alt_upper < this_upper
            this_upper = alt_upper;
            z = alt_z;
        end
    end
    if this_upper < upper
        x_min = z;
        if length(info_text) == 0
            if isinf(upper)
                info_text = 'Solution found by upper solver';
            else
                info_text = 'Improved solution found by upper solver';
            end
        else
            if isinf(upper)
                info_text = [info_text ' | ' 'Solution found by upper solver'];
            else
                info_text = [info_text ' | ' 'Improved solution found by upper solver'];
            end
         end
        upper = this_upper;
        numglobals = numglobals + 1;
    end
elseif output.problem == 0 && p_upper.K.s(1)>0 && ~p_upper.solver.uppersolver.constraint.inequalities.semidefinite.linear && cutiterations > 0
    % Piggy-back some old code. Generate elementwise cuts from
    % vioalated SDP constraints (if there are any) and add these to the
    % global upper bound model.
    p_u = p_upper;
    p_u.lpcuts = [];
    p_u.cutState = [];
    p_u=createsdpcut(p_u,z);
    p_upper.F_struc = [p_upper.F_struc(1:p_upper.K.f,:);p_u.lpcuts;p_upper.F_struc(1+p_upper.K.f:end,:)];
    p_upper.K.l  = p_upper.K.l + size(p_u.lpcuts,1);
%    if ~isempty(p_u.socpcuts)
%        p_upper.F_struc = [p_upper.F_struc(1:p_upper.K.f+p_upper.K.l,:);p_u.socpcuts;p_upper.F_struc(1+p_upper.K.f+p_upper.K.l:end,:)];
%         if isequal(p_upper.K.q,0)
%             p_upper.K.q  = repmat(3,1,size(p_u.socpcuts,1)/3);
%         else
%             p_upper.K.q  = [repmat(3,1,size(p_u.socpcuts,1)/3) p_upper.K.q];
%         end
%    end

    if ~isempty(p_u.lpcuts)
        n = length(p_u.cutState);
        if length(info_text) == 0
            info_text = ['Added ' num2str(n) ' cut' plurals(n) ' on SDP cone'];
        else
            info_text = [info_text ' | ' 'Added ' num2str(n) ' cut' plurals(n) ' on SDP cone'];
        end
    end

    % Call recursively if we added cuts, and user wants to run more
    % than one iteration
    if ~isempty(p_u.cutState) > 0 && cutiterations > 1
        [upper,x_min,info_text,numglobals,timing,p_upper] = solve_upper_in_node(p,p_upper,x,upper,x_min,uppersolver,info_text,numglobals,timing,cutiterations-1);
    end
end

function p = pruneSDPCone(p)
if any(p.K.s) && p.options.bmibnb.uppersdprelax
    if ~p.solver.uppersolver.constraint.inequalities.semidefinite.linear
        p.F_struc(end-sum(p.K.s.^2)+1:end,:) = [];
        p.K.s = 0;
    end
end

function s = plurals(n)
if n>1
    s = 's';
else
    s = '';
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ word, done ] = word_next_read ( s, done )

%*****************************************************************************80
%
%% WORD_NEXT_READ "reads" words from a string, one at a time.
%
%  Special cases:
%
%    The following characters are considered to be a single word,
%    whether surrounded by spaces or not:
%
%      " ( ) { } [ ]
%
%    Also, if there is a trailing comma on the word, it is stripped off.
%    This is to facilitate the reading of lists.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string, presumably containing words
%    separated by spaces.
%
%    Input, logical DONE.
%    TRUE, if we are beginning a new string;
%    FALSE, if we are continuing to process the current string.
%
%    Output, string WORD.
%    If DONE is FALSE, then WORD contains the "next" word read.
%    If DONE is TRUE, then WORD is blank, because there was no more to read.
%
%    Output, logical DONE.
%      FALSE if another word was read,
%      TRUE if no more words could be read.
%
  persistent lenc;
  persistent next;

  tab = char ( 9 );
%
%  We "remember" LENC and NEXT from the previous call.
%
%  An input value of DONE = TRUE signals a new line of text to examine.
%
  if ( done )

    next = 1;
    done = 0;
    lenc = s_len_trim ( s );

    if ( lenc <= 0 )
      done = 1;
      word = ' ';
      return
    end

  end
%
%  Beginning at index NEXT, search the string for the next nonblank,
%  which signals the beginning of a word.
%
  ilo = next;
%
%  ...S(NEXT:) is blank.  Return with WORD = ' ' and DONE = TRUE.
%
  while ( 1 )

    if ( lenc < ilo )
      word = ' ';
      done = 1;
      next = lenc + 1;
      return
    end
%
%  If the current character is blank, skip to the next one.
%
    if ( s(ilo) ~= ' ' & s(ilo) ~= tab )
      break
    end

    ilo = ilo + 1;

  end
%
%  ILO is the index of the next nonblank character in the string.
%
%  If this initial nonblank is a special character,
%  then that's the whole word as far as we're concerned,
%  so return immediately.
%
  if ( s(ilo) == '"' | ...
       s(ilo) == '(' | ...
       s(ilo) == ')' | ...
       s(ilo) == '{' | ...
       s(ilo) == '}' | ...
       s(ilo) == '[' | ...
       s(ilo) == ']' )

    word = s(ilo);
    next = ilo + 1;
    return

  end
%
%  Now search for the last contiguous character that is not a
%  blank, TAB, or special character.
%
  next = ilo + 1;

  while ( next <= lenc )

    if ( s(next) == ' ' )
      break;
    elseif ( s(next) == tab )
      break;
    elseif ( s(next) == '"' )
      break;
    elseif ( s(next) == '(' )
      break;
    elseif ( s(next) == ')' )
      break;
    elseif ( s(next) == '{' )
      break;
    elseif ( s(next) == '}' )
      break;
    elseif ( s(next) == '[' )
      break;
    elseif ( s(next) == ']' )
      break;
    end

    next = next + 1;

  end

  if ( s(next-1) == ',' )
    word = s(ilo:next-2);
  else
    word = s(ilo:next-1);
  end

  return
end


text:
classdef GenericCylinder < Shape
	% GenericCylinder is a Shape for a cylinder whose cross section is any 2D
	% shape.

	properties (SetAccess = immutable)
		lsf2d  % level set function in lateral domain
	end

	methods
        function this = GenericCylinder(normal_axis, lsf2d, lprim_cell, dl_max)
			chkarg(istypesizeof(normal_axis, 'Axis'), '"normal_axis" should be instance of Axis.');

			% lsf2d is a level set function defined in 2D.  It takes an argument
			% r = [h, v], where, e.g., (h, v) is (x, y) for narmal_axis: z.
			chkarg(istypeof(lsf2d, 'function_handle'), '"lsf2d" should be function handle.');

			chkarg(istypesizeof(lprim_cell, 'realcell', [1 Axis.count], [1 0]), ...
				'"lprim_cell" should be length-%d row cell array whose each element is row vector with real elements.', Axis.count);

			[h, v, n] = cycle(normal_axis);
			bound_n = [min(lprim_cell{n}), max(lprim_cell{n})];
			sn = diff(bound_n) / 2;  % semiside in normal direction
			cn = mean(bound_n);  % center in normal direction

			% For rhv = r([h, v]), and rn = r(n), the level set
			% function is basically min(lsf2d(rhv), 1 - abs(rn-cn)./sn, [], 2),
			% but it is vectorized, i.e., modified to handle r = [x y z] with
			% column vectors x, y, z.
			function level = lsf(x, y, z)
				chkarg(istypeof(x, 'real'), '"x" should be array with real elements.');
				chkarg(istypeof(y, 'real'), '"y" should be array with real elements.');
				chkarg(istypeof(z, 'real'), '"z" should be array with real elements.');
				chkarg(isequal(size(x), size(y), size(z)), '"x", "y", "z" should have same size.');

				loc = {x, y, z};
				level = min(lsf2d(loc{h}, loc{v}), 1 - abs(loc{n}-cn)./sn);  % intersection of regions defined by two level set functions
			end

			if nargin < 4  % no dl_max
				super_args = {lprim_cell, @lsf};
			else
				super_args = {lprim_cell, @lsf, dl_max};
			end

			this = this@Shape(super_args{:});
			this.lsf2d = lsf2d;
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function element_num = grid_t6_element_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_T6_ELEMENT_NUM counts the elements in a grid of 6 node triangles.
%
%  Example:
%
%    Input:
%
%      NELEMX = 3, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NUM = 2 * NELEMX * NELEMY = 12
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NUM, the number of elements in the grid.
%
  element_num = 2 * nelemx * nelemy;

  return
end


text:
function value = zsign2 ( z1, z2 )

%*****************************************************************************80
%
%% ZSIGN2 is a complex transfer-of-sign function.
%
%  Discussion:
%
%    The L2 norm is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, complex Z1, Z2, the arguments.
%
%    Output, complex VALUE,  a complex value, with the magnitude of
%    Z1, and the argument of Z2.
%
  if ( zabs2 ( z2 ) == 0.0 )
    value = 0.0;
  else
    value = zabs2 ( z1 ) * ( z2 / zabs2 ( z2 ) );
  end

  return
end


text:
function Dy = computeDy(Xmin, Xmax, ny)
    Dy_tmp = cell(2*ny);

    for i=1:2*ny
        if i <= ny
            Dy_tmp{i,1} = Xmax;
        else
            Dy_tmp{i,1} = -Xmin;
        end
    end

    Dy =cell2mat(Dy_tmp);
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


%**************************************************************************
% SARSA Learning applied to Cart-Pole balancing problem.
% The environment of the learning system is a black box, from which it has
% several lines and a reinforcement line. Its task is to learn to give
% responses which maximize the scalar signals on its reinforcement line.
% The Q-value update is the following, if the system takes action a from
%   state s at time t, and arrives at state ss with feedback r at time t+1:
%
%   Q(t+1, s, a) = Q(t, s, a)
%                  + alpha (r + gamma *Q(t,ss, b) - Q(t, s, a))

% get_box:  Given the current state, returns a number from 1 to 162
%           designating the region of the state space encompassing the current
%           state.
%           Returns a value of -1 if a failure state is encountered.

% cart_pole2: The cart and pole dynamics; given the Force and
%             current state, estimates next state

%Code written by: Savinay Nagendra
%email id:        nagsavi17@gmail.com
%**************************************************************************
clc;
clear all;
close all;
% Initialization

NUM_BOXES = 163;
ALPHA = 0.4;             % Learning rate parameter
GAMMA = 0.999;           % Discount factor for future reinf
Q = zeros(NUM_BOXES,2);  % State-Action Values
action = [30 -30];
MAX_FAILURES = 1000;
MAX_STEPS = 150000;
epsilon = 0;
steps = 0;
failures = 0;
thetaPlot = 0;
xPlot = 0;
%Pendulum state initialization
theta = 0;
thetaDot = 0;
x = 0;
xDot = 0;
box = getBox4(theta,thetaDot,x,xDot);
if(rand>epsilon)       % exploit
    [~,actionMax] = max(Q(box,:));
    currentAction = action(actionMax);
else                   % explore
    currentAction = datasample(action,1);
end
actionIndex1 = find(action == currentAction);

while(steps<=MAX_STEPS && failures<+MAX_FAILURES)
    steps = steps + 1;
    [thetaNext,thetaDotNext,thetaacc,xNext,xDotNext] = cart_pole2(currentAction,theta,thetaDot,x,xDot);
    thetaPlot(end + 1) = thetaNext;
    xPlot(end + 1) = xNext;
    newBox = getBox4(thetaNext,thetaDotNext,xNext,xDotNext);
    theta = thetaNext;
    thetaDot = thetaDotNext;
    x = xNext;
    xDot = xDotNext;
    if(newBox == 163)
        r = -1;
        Q(newBox,:) = 0;
        figure(2);
        plot((1:length(thetaPlot)),thetaPlot,'-ob');
        figure(3);
        plot((1:length(xPlot)),xPlot,'-og');

        thetaPlot = 0;
        xPlot = 0;
        %Swing Up. Find the box.
        theta = 0;
        thetaDot = 0;
        x = 0;
        xDot = 0;

        newBox = getBox4(theta,thetaDot,x,xDot);
        failures = failures + 1;
        fprintf('Trial %d was %d steps. \n',failures,steps);
        figure(1);
        plot(failures,steps,'-or');
        hold on;
        steps = 0;
    else
        r = 0;
    end
    if(rand > epsilon)  %exploit
        [~,newActionMax] = max(Q(newBox,:));
        newAction = action(newActionMax);
    else                %explore
        newAction = datasample(action,1);
    end
    actionIndex2 = find(action == newAction);
    Q(box,actionIndex1) = Q(box,actionIndex1) + ALPHA*(r + GAMMA*Q(newBox,actionIndex2) - Q(box,actionIndex1));
    box = newBox;
    currentAction = newAction;
    actionIndex1 = actionIndex2;
end

if(failures == MAX_FAILURES)
    fprintf('Pole not balanced. Stopping after %d failures.',failures);
else
    fprintf('Pole balanced successfully for at least %d steps\n', steps);
    figure(1);
    plot(failures+1,steps,'-or');
    hold on;
    figure(2);
    plot((1:length(thetaPlot)),thetaPlot,'-ob');
    figure(3);
    plot((1:length(xPlot)),xPlot,'-og');
    figure(4);
    plot((1:301),thetaPlot(1:301),'-ob');
    hold on;
    figure(5);
    plot((1:301),xPlot(1:301),'-og');
    hold on;
end


text:
% This function returns the oblique shock wave angle (beta) for a given
% deflection angle (theta) in degrees and ratio of specific heats (gamma).
% and Mach number M. Specify 0 for the weak oblique shock
% or 1 for the strong shock.
%
% Syntax:
% beta(M,theta,gamma,n)     where n specifies weak or strong shock returned
%
% NOTE: Angles supplied and returned from this function are in DEGREES.
%
% Based on an analytical solution to the theta-beta-Mach relation given in
% the following reference:  Rudd, L., and Lewis, M. J., "Comparison of
% Shock Calculation Methods", AIAA Journal of Aircraft, Vol. 35, No. 4,
% July-August, 1998, pp. 647-649.
%
% By Chris Plumley, undergraduate, University of Maryland.

function Beta=beta(M,theta,gamma,n)

theta=theta*pi/180;             % convert to radians
mu=asin(1/M);                   % Mach wave angle
c=tan(mu)^2;
a=((gamma-1)/2+(gamma+1)*c/2)*tan(theta);
b=((gamma+1)/2+(gamma+3)*c/2)*tan(theta);
d=sqrt(4*(1-3*a*b)^3/((27*a^2*c+9*a*b-2)^2)-1);
Beta=atan((b+9*a*c)/(2*(1-3*a*b))-(d*(27*a^2*c+9*a*b-2))/(6*a*(1-3*a*b))*tan(n*pi/3+1/3*atan(1/d)))*180/pi;


text:
function [mn] = calcSpectralMoment(angFreq,S_f,order)
% Parameters
% ------------
%     angFreq : [1 n] float vector
%         Vector of wave angular frequencies
%
%     S_f : [1 n] float vector
%         Vector of wave spectra
%
%     order : float
%         Order of the spectral moment
%
% Returns
% ------------
%     mn : float vector
%         Spectral moment of the nth order
%

% Check the dimensions of angFreq and S_f
[raf,caf]   = size(angFreq);
[rsf,csf]   = size(S_f);
if raf ~= rsf && caf ~= csf
    if raf == csf && caf == rsf
        S_f = S_f.';
    else
        disp('The vectors for wave angular frequency and spectral moment must be the same size')
        mn = NaN;
        return
    end
end

% Check the dimensions of order
[ro,co]     = size(order);
if ro ~= 1 || co ~=1
    disp('The order of the spectal moment calculation must be a scale value')
    mn = NaN;
    return
end

% Spectral moment calculation
mn = trapz(angFreq,angFreq.^(order).*S_f);
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function varargin = conv(varargout)
%CONV   Operation not supported for TRIGTECH objects.
%
%   Use the function CIRCCONV.
%
% See also CIRCCONV.

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

error('CHEBFUN:TRIGTECH:conv:NotSupported', ...
    'CONV is not supported for trigtech objects. Use the function CIRCCONV.');

end