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).
% bws = An integer [1,10] representing the number of steps between each rolling window (optional, default=10).
% fevd = A string representing the FEVD type used by the variance decomposition (optional, default='G'):
%   - 'G' for generalized FEVD;
%   - 'O' for orthogonal FEVD.
% lags = An integer [1,3] representing the number of lags of the VAR model used by the variance decomposition (optional, default=2).
% h = An integer [1,10] representing the prediction horizon used by the variance decomposition (optional, default=4).
% 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.

function [result,stopped] = run_spillover(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('bws',10,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 1 '<=' 10 'scalar'}));
        ip.addOptional('fevd','G',@(x)any(validatestring(x,{'G' 'O'})));
        ip.addOptional('lags',2,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 1 '<=' 3 'scalar'}));
        ip.addOptional('h',4,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 1 '<=' 10 'scalar'}));
        ip.addOptional('analyze',false,@(x)validateattributes(x,{'logical'},{'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    ds = validate_dataset(ipr.ds,'Spillover');
    sn = ipr.sn;
    temp = validate_template(ipr.temp);
    out = validate_output(ipr.out);
    bw = ipr.bw;
    bws = ipr.bws;
    fevd = ipr.fevd;
    lags = ipr.lags;
    h = ipr.h;
    analyze = ipr.analyze;

    nargoutchk(1,2);

    [result,stopped] = run_spillover_internal(ds,sn,temp,out,bw,bws,fevd,lags,h,analyze);

end

function [result,stopped] = run_spillover_internal(ds,sn,temp,out,bw,bws,fevd,lags,h,analyze)

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

    indices = unique([1:bws:ds.T ds.T]);
    ds = initialize(ds,sn,bw,bws,indices,fevd,lags,h);

    rng(double(bitxor(uint16('T'),uint16('B'))));
    cleanup_1 = onCleanup(@()rng('default'));

    bar = waitbar(0,'Initializing spillover measures...','CreateCancelBtn',@(src,event)setappdata(gcbf(),'Stop',true));
    setappdata(bar,'Stop',false);
    cleanup_2 = onCleanup(@()delete(bar));

    pause(1);
    waitbar(0,bar,'Calculating spillover measures...');
    pause(1);

    try

        windows = extract_rolling_windows(ds.Returns,ds.BW);
        windows = windows(indices);
        windows_len = numel(windows);

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

        for i = 1:windows_len
           futures(i) = parfeval(@main_loop,1,windows{i},ds.FEVD,ds.Lags,ds.H);
        end

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

            [future_index,value] = fetchNext(futures);
            futures_results{future_index} = value;

            futures_max = max([future_index futures_max]);
            waitbar((futures_max - 1) / windows_len,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 spillover measures...');
    pause(1);

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

    pause(1);
    waitbar(1,bar,'Writing spillover 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,bws,indices,fevd,lags,h)

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

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

    ds.BW = bw;
    ds.BWS = bws;
    ds.FEVD = fevd;
    ds.H = h;
    ds.Lags = lags;
    ds.Windows = numel(indices);
    ds.WindowsIndices = indices;

    all_label = [' (' fevd ', H=' num2str(ds.H) ', LAGS=' num2str(ds.Lags) ')'];

    ds.LabelsIndicatorsSimple = {'SI'};
    ds.LabelsIndicators = {['SI' all_label]};

    ds.LabelsSheetsSimple = {'From' 'To' 'Net' 'Indicators'};
    ds.LabelsSheets = {['From' all_label] ['To' all_label] ['Net' all_label] 'Indicators'};

    ds.VarianceDecompositions = cell(t,1);

    ds.SpilloversFrom = NaN(t,n);
    ds.SpilloversTo = NaN(t,n);
    ds.SpilloversNet = NaN(t,n);

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

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

end

function window_results = main_loop(r,fevd,lags,h)

    window_results = struct();

    vd = variance_decomposition(r,fevd,lags,h);
    window_results.VarianceDecomposition = vd;

    [sf,st,sn,si] = spillover_metrics(vd);
    window_results.SpilloversFrom = sf;
    window_results.SpilloversTo = st;
    window_results.SpilloversNet = sn;
    window_results.SI = si;

end

function ds = finalize(ds,results)

    window_indices = ds.WindowsIndices;

    for i = 1:numel(results)
        result = results{i};
        index = window_indices(i);

        ds.VarianceDecompositions{index} = result.VarianceDecomposition;

        ds.SpilloversFrom(index,:) = result.SpilloversFrom;
        ds.SpilloversTo(index,:) = result.SpilloversTo;
        ds.SpilloversNet(index,:) = result.SpilloversNet;

        ds.Indicators(index) = result.SI;
    end

    if (ds.BWS > 1)
        n = ds.N;
        t = ds.T;
        x = ds.DatesNum;

        nan_indices = ~ismember((1:t).',window_indices);
        step_indices = find(nan_indices);
        step_indices_len = numel(step_indices);

        vd = ds.VarianceDecompositions;
        vd(cellfun(@isempty,vd)) = {NaN(n,n)};

        for i = 1:n
            for j = 1:n
                vd_ij = cellfun(@(vdf)vdf(i,j),vd);
                vd_ij(nan_indices) = spline(x(~nan_indices),vd_ij(~nan_indices),x(nan_indices));

                for k = 1:step_indices_len
                    step_index = step_indices(k);

                    vd_k = vd{step_index};
                    vd_k(i,j) = vd_ij(step_index);
                    vd{step_index} = vd_k;
                end
            end
        end

        for i = 1:step_indices_len
            step_index = step_indices(i);

            vd_i = vd{step_index};
            vd_i = bsxfun(@rdivide,vd_i,sum(vd_i,2));
            ds.VarianceDecompositions{step_index} = vd_i;

            [sf,st,sn,si] = spillover_metrics(vd_i);
            ds.SpilloversFrom(step_index,:) = sf;
            ds.SpilloversTo(step_index,:) = st;
            ds.SpilloversNet(step_index,:) = sn;
            ds.Indicators(step_index,1) = si;
        end
    end

    ds.SpilloversFrom = distress_data(ds.SpilloversFrom,ds.Defaults);
    ds.SpilloversTo = distress_data(ds.SpilloversTo,ds.Defaults);
    ds.SpilloversNet = distress_data(ds.SpilloversNet,ds.Defaults);

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 output file could not be created from the template file.');
    end

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

    for i = 1:(numel(ds.LabelsSheetsSimple) - 1)
        sheet = ds.LabelsSheetsSimple{i};
        measure = ['Spillovers' strrep(sheet,' ','')];

        tab = [dates_str array2table(ds.(measure),'VariableNames',ds.FirmNames)];
        writetable(tab,out,'FileType','spreadsheet','Sheet',sheet,'WriteRowNames',true);
    end

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

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

end

%% PLOTTING

function analyze_result(ds)

    safe_plot(@(id)plot_index(ds,id));
    safe_plot(@(id)plot_spillovers(ds,id));
    safe_plot(@(id)plot_sequence(ds,id));

end

function plot_index(ds,id)

    si = smooth_data(ds.Indicators(:,1));

    f = figure('Name','Spillover Measures > Spillover Index','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(1,6,1:5);
    plot(sub_1,ds.DatesNum,si,'Color',[0.000 0.447 0.741]);
    set(sub_1,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_1,'XGrid','on','YGrid','on');
    title(sub_1,ds.LabelsIndicators{1});

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

    sub_2 = subplot(1,6,6);
    boxplot(sub_2,si,'Notch','on','Symbol','k.');
    set(findobj(f,'type','line','Tag','Median'),'Color','g');
    set(findobj(f,'-regexp','Tag','\w*Whisker'),'LineStyle','-');
    set(sub_2,'TickLength',[0 0],'XTick',[],'XTickLabels',[]);

    figure_title(f,'Spillover Index');

    maximize_figure(f);

end

function plot_spillovers(ds,id)

    from = smooth_data(ds.SpilloversFrom);
    from = bsxfun(@rdivide,from,sum(from,2,'omitnan'));
    from(isnan(from)) = 0;

    to = smooth_data(ds.SpilloversTo);
    to = bsxfun(@rdivide,to,sum(to,2,'omitnan'));
    to(isnan(to)) = 0;

    net = smooth_data(ds.SpilloversNet);
    net = [min(net,[],2,'omitnan') max(net,[],2,'omitnan')];

    f = figure('Name','Spillover Measures > Spillovers','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(2,2,1);
    area(sub_1,ds.DatesNum,from,'EdgeColor',[0.000 0.447 0.741],'FaceColor','none');
    set(sub_1,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_1,'YLim',[0 1],'YTick',0:0.2:1,'YTickLabels',arrayfun(@(x)sprintf('%.f%%',x),(0:0.2:1) .* 100,'UniformOutput',false));
    title(sub_1,'Spillovers From Others');

    sub_2 = subplot(2,2,3);
    area(sub_2,ds.DatesNum,to,'EdgeColor',[0.000 0.447 0.741],'FaceColor','none');
    set(sub_2,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_2,'YLim',[0 1],'YTick',0:0.2:1,'YTickLabels',arrayfun(@(x)sprintf('%.f%%',x),(0:0.2:1) .* 100,'UniformOutput',false));
    title(sub_2,'Spillovers To Others');

    sub_3 = subplot(2,2,[2 4]);
    fill(sub_3,[ds.DatesNum; flipud(ds.DatesNum)],[net(:,1); fliplr(net(:,2))],[0.65 0.65 0.65],'EdgeColor','none','FaceAlpha',0.35);
    hold on;
        plot(sub_3,ds.DatesNum,mean(net,2),'Color',[0.000 0.447 0.741]);
        plot(sub_3,ds.DatesNum,zeros(ds.T,1),'Color',[1.000 0.400 0.400]);
    hold off;
    set(sub_3,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_3,'YLim',[-1 1],'YTick',-1:0.2:1,'YTickLabels',arrayfun(@(x)sprintf('%.f%%',x),(-1:0.2:1) .* 100,'UniformOutput',false));
    set(sub_3,'XGrid','on','YGrid','on');
    title(sub_3,'Net Spillovers');

    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,['Spillovers (' ds.FEVD ', H=' num2str(ds.H) ', LAGS=' num2str(ds.Lags) ')']);

    maximize_figure(f);

end

function plot_sequence(ds,id)

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

    from_all = smooth_data(ds.SpilloversFrom);
    to_all = smooth_data(ds.SpilloversTo);
    net_all = smooth_data(ds.SpilloversNet);

    data = [repmat({dn},1,n); mat2cell(from_all,t,ones(1,n)); mat2cell(to_all,t,ones(1,n)); mat2cell(net_all,t,ones(1,n))];

    plots_title = [repmat({'From'},1,n); repmat({'To'},1,n); repmat({'Net'},1,n)];

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

    ft = [from_all to_all];
    y_limits_from = [0 1];
    y_limits_to = plot_limits(max(max(ft)),0.1,0);
    y_limits_to(2) = ceil(y_limits_to(2) * 10) / 10;
    y_limits_net = [-1 1];

    y_ticks_from = 0:0.2:1;
    y_ticks_to = 0:0.2:y_limits_to(2);
    y_ticks_net = -1:0.2:1;
    y_tick_labels = @(x)sprintf('%.f%%',x * 100);

    core = struct();

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

    core.OuterTitle = 'Spillover Measures > Spillovers Time Series';
    core.InnerTitle = 'Spillovers Time Series';
    core.SequenceTitles = ds.FirmNames;

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

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

    core.YGrid = {true true true};
    core.YLabel = {[] [] []};
    core.YLimits = {y_limits_from y_limits_to y_limits_net};
    core.YRotation = {[] [] []};
    core.YTick = {y_ticks_from y_ticks_to y_ticks_net};
    core.YTickLabels = {y_tick_labels y_tick_labels y_tick_labels};

    sequential_plot(core,id);

    function plot_function(subs,data)

        x = data{1};
        from = data{2};
        to = data{3};
        net = data{4};

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

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

		sub_1 = subs(1);
        plot(sub_1,x,from,'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,to,'Color',[0.000 0.447 0.741]);
        hold(sub_2,'on');
            plot(sub_2,x,ones(numel(x),1),'Color',[1.000 0.400 0.400]);
            if (~isempty(xd))
                plot(sub_2,[xd xd],get(sub_2,'YLim'),'Color',[1.000 0.400 0.400]);
            end
        hold(sub_2,'off');

		sub_3 = subs(3);
        plot(sub_3,x,net,'Color',[0.000 0.447 0.741]);
        hold(sub_3,'on');
            plot(sub_3,x,zeros(numel(x),1),'Color',[1.000 0.400 0.400]);
            if (~isempty(xd))
                plot(sub_3,[xd xd],get(sub_3,'YLim'),'Color',[1.000 0.400 0.400]);
            end
        hold(sub_3,'off');

    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 temp = validate_template(temp)

    sheets = {'From' 'To' 'Net' 'Indicators'};
    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 rworlds = gen_random_obs(rows,cols,density,nworlds)
%
%Generated set of random obstacles to use to create random worlds
%
nobs = round((density/100)*rows*cols); %find number of obstacles
for i = 1:nworlds
    w(1:rows,1:cols) = -1;
    for j = 1:nobs
        tf = 1;
        while tf==1 %select empty location
            tr = round(1 + ( rows - 1) .* rand);
            tc = round(1 + ( cols - 1) .* rand);
            tf = 0;
            if w(tr,tc) == 1
                tf = 1;
            end
        end
        w(tr,tc) = 1;
        tdir = round(1 + ( 4 - 1) .* rand);%set a random direction
        switch tdir
            case 1
                td = 'N';
            case 2
                td = 'E';
            case 3
                td = 'W';
            case 4
                td = 'S';
        end
        tv = round((min(rows,cols)-1) .* rand);%set a random velocity
        rw(i).obs(j).r = tr;
        rw(i).obs(j).c = tc;
        rw(i).obs(j).vel = tv;
        rw(i).obs(j).dir = td;
    end
end
rworlds = rw;
end

author: Sable

repo: mcbench-benchmarks


text:
clear all; close all;
I=imread('cameraman.tif');
I=im2double(I);
I=imnoise(I, 'salt & pepper', 0.01);
J=ordfilt2(I, 1, ones(4,4));
K=ordfilt2(I, 9, ones(3));
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
function bool = checkFormulae(formulaA, formulaB, exceptions)
% Compares two formulae for the number of elements for all optionally except the
% elements listed in exceptions
%
% USAGE:
%
%    bool = checkFormulae(formulaA, formulaB, exceptions)
%
% INPUT:
%    formulaA:
%    formulaB:
%
% OPTIONAL INPUT:
%    exceptions:    String array of metabolite abbreviation exceptions e.g. {'H'}
%
% OUTPUT:
%    bool:          Boolean where 1 means the formulae match and 0 otherwise
%
% .. Author: - Ronan M. T. Fleming 20 July 2009

if ~exist('exceptions','var')
    exceptions={};
end

if ischar(exceptions)
    ex=exceptions;
    clear exceptions
    exceptions{1}=ex;
end

elements = {'C','H','O','P','S','K','Na','N'};

bool=1;
for n=1:length(elements)
    if ~any(strcmp(elements{n},exceptions))
        numA=numAtomsOfElementInFormula(formulaA,elements{n});
        numB=numAtomsOfElementInFormula(formulaB,elements{n});
        %check to see if number of elements is not the same for each
        %metabolite
        if numA~=numB
            bool=0;
            break;
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function j = jacobi_symbol ( q, p )

%*****************************************************************************80
%
%% JACOBI_SYMBOL evaluates the Jacobi symbol (Q/P).
%
%  Discussion:
%
%    If P is prime, then
%
%      Jacobi Symbol (Q/P) = Legendre Symbol (Q/P)
%
%    Else
%
%      let P have the prime factorization
%
%        P = Product ( 1 <= I <= N ) P(I)**E(I)
%
%      Jacobi Symbol (Q/P) =
%
%        Product ( 1 <= I <= N ) Legendre Symbol (Q/P(I))**E(I)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Daniel Zwillinger,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996, pages 86-87.
%
%  Parameters:
%
%    Input, integer Q, an integer whose Jacobi symbol with
%    respect to P is desired.
%
%    Input, integer P, the number with respect to which the Jacobi
%    symbol of Q is desired.  P should be 2 or greater.
%
%    Output, integer J, the Jacobi symbol (Q/P).
%    Ordinarily, J will be -1, 0 or 1.
%    -2, not enough factorization space.
%    -3, an error during Legendre symbol calculation.
%

%
%  P must be greater than 1.
%
  if ( p <= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_SYMBOL - Fatal error!\n' );
    fprintf ( 1, '  P must be greater than 1.\n' );
    l = -2;
    return
  end
%
%  Decompose P into factors of prime powers.
%
  [ nfactor, factor, power, nleft ] = i4_factor ( p );

  if ( nleft ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_SYMBOL - Fatal error!\n' );
    fprintf ( 1, '  Not enough factorization space.\n' );
    j = -2;
    error ( 'JACOBI_SYMBOL - Fatal error!' );
  end
%
%  Force Q to be nonnegative.
%
  qq = q;

  while ( qq < 0 )
    qq = qq + p;
  end
%
%  For each prime factor, compute the Legendre symbol, and
%  multiply the Jacobi symbol by the appropriate factor.
%
  j = 1;
  for i = 1 : nfactor
    pp = factor(i);
    l = legendre_symbol ( qq, pp );
    if ( l < -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'JACOBI_SYMBOL - Fatal error!\n' );
      fprintf ( 1, '  Error during Legendre symbol calculation.\n' );
      j = -3;
      error ( 'JACOBI_SYMBOL - Fatal error!' );
    end
    j = j * l^power(i);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [c, alld] = getVertexColors(xyz, v, actcolor, varargin)
% :Usage:
% ::
%
%    [c, alld] = getVertexColors(xyz, v, actcolor, [basecolor], [mind], 'vert', [xyz2], [actcolor2], 'vert', [xyz3], [actcolor3])
%
% given a point list of XYZ mm coordinates (3 columns)
% and a list of vertices in an isosurface,
% returns FaceVertexCData color values for brain near points and brain not near points.
% c is vertex color specification, 3 columns indicating RGB values
%
% :Inputs:
%
%   **xyz:**
%        a 3-vol list of vertices to color
%
%   **v:**
%        can be a matrix of vertices
%        or a handle to a patch object containing vertices
%        if it's a handle, this function sets the color to interp
%        and the FaceVertexCData to the color matrix c
%
%   **actcolor:**
%        [r g b] activation color
%
%   **basecolor:**
%        [r g b] baseline color - optional.
%
%   **mind:**
%        optional - min distance to color vertex
%        Vertices within mind of an xyz coordinate will be colored
%
%   **colorscale
%       optional.  followed by vector of values by which to multiply input
%       color
%       these are scaled to be between .3 and one.
%
%       if entered, this will make the colors vary by, for example, Z score
%       so Z-scores are an acceptable input.
%
%       cscale should be in the same coordinate order as xyz
%
%       for ADDITIONAL clusters, repeat the 'colorscale', Z argument pair in the function call
%
%       YOU CAN ALSO pass true RGB values for each xyz coordinate in: 'colorscale', rgblist,
%       IF cscale is a 3-vector, it specifies the ACTUAL colors, and is not scaled to .3 - 1
%
%   Following basecolor and mind:
%
%   additional xyz coordinate lists, with syntax:
%   vert', xyz2 [your xyz input], [r g b] color for xyz plot
%
%   also, you can enter 'ovlcolor' followed by [r g b] for overlaps between xyz sets
%   colors will ONLY appear in the overlap color if they share actual coordinates in common,
%   not necessarily if surface vertices are within the specified distance from both sets of coords.
%
% :Examples:
% ::
%
%    % to get a good brain surface, try this:
%    figure
%    p = patch('Faces', faces, 'Vertices', vertices, 'FaceColor', [.5 .5 .5], ...
%            å  'EdgeColor', 'none', 'SpecularStrength', .2, 'FaceAlpha', 1, 'SpecularExponent', 200);
%    lighting gouraud;
%    camlight right
%    axis image;
%    myLight = camlight(0, 0);
%    set(myLight, 'Tag', 'myLight');
%    set(gcf, 'WindowButtonUpFcn', 'lightFollowView');
%    lightfollowview
%    drawnow
%
% ..
% by Tor Wager  August 25, 2002
% ..

    global do_fsavg_right do_fsavg_left ras

    mind = 3;
    basecolor = [.5 .5 .5];
    alld = [];
    xyza = xyz;
    cscale = [];
    allda = [];
    vv = [];
    doalph = 0;
    cscale = [];
    alphascale = {};
    do_fsavg_left = false;
    do_fsavg_right = false;
    doverbose = true;

    % -----------------------------------------------------------------------
    % * set up input arguments
    % -----------------------------------------------------------------------

    if ~isempty(varargin), basecolor = varargin{1}; end

    if length(varargin) > 1, mind = varargin{2}; end

    ind = 1;
    for i = 3:length(varargin)
        if strcmp(varargin{i}, 'vert')
            vv{ind} = varargin{i+1};

            % intersections
            xyzb{ind, 1} = intersect(xyz, vv{ind}, 'rows');
            for j = ind-1:-1:1
                xyzb{ind, j} = intersect(vv{j}, vv{ind}, 'rows');
                xyza = intersect(xyza, xyzb{ind, j}, 'rows');
            end

            cc{ind} = varargin{i+2};
            ind = ind+1;
        elseif strcmp(varargin{i}, 'ovlcolor')
            ocol = varargin{i+1};

        elseif strcmp(varargin{i}, 'alphascale')
            doalph = 1;
            alphascale{1} = varargin{i + 1};

        elseif strcmp(varargin{i}, 'allcolor')
            acol = varargin{i+1};

        elseif strcmp(varargin{i}, 'colorscale')

            cscale{end+1} = varargin{i+1};

            % now do this ahead of time
            %             if min(size(cscale{end})) == 1
            %                 % scale colors (may be necessary) - only if single vector, not RGB values
            %                 cscale{end} = cscale{end} ./ max(cscale{end});
            %             end
            if any(cscale{end} < 0), error('Some color scale values are less than zero.'), end
        elseif strcmp(varargin{i}, 'fsavg_left')
            ras = load(which('lh.avgMapping_allSub_RF_ANTs_MNI152_orig_to_fsaverage.mat'));
            do_fsavg_left = true;
        elseif strcmp(varargin{i}, 'fsavg_right')
            % uses freesurfer inflated brain with Thomas Yeo group's RF_ANTs mapping
            % from MNI to Freesurfer. (https://doi.org/10.1002/hbm.24213)
            ras = load(which('rh.avgMapping_allSub_RF_ANTs_MNI152_orig_to_fsaverage.mat'));
            do_fsavg_right = true;

        elseif strcmp(varargin{i}, 'noverbose')
            doverbose = false;

        end
    end

    if isempty(cscale)
        cscale{1} = ones(size(xyz, 1), 1);
        for i = 1:length(vv)    % additional vertices
            cscale{end+1} = ones(size(vv{i}, 1), 1);
        end
    end

    %if ~isempty(cscale) don't need this.
    %    if length(cscale) < length(vv)
    %        cscale{length(vv)} = [];
    %    end
    %end

    if ishandle(v)
        p = v;
        v = get(p, 'Vertices');
        c = get(p, 'FaceVertexCData');   % get existing colors from surface
    else
        % uh-oh, not a handle!
        warning('Figure handle missing: Figure was closed?')
        p = findobj('Type', 'patch');
        v = get(p(1), 'Vertices');
        c = get(p, 'FaceVertexCData');   % get existing colors from surface
    end

    % if FaceVertexCData is v x 1 matrix, it is mapped to the figure colormap
    % if it is v x 3, it is true color.
    % We can replace the colormapped FaceVertexCData with its true-color
    % equivalent, making it colormap independent.
    % This will allow us to use the existing values and average in new colors
    % for blobs
    if  ishandle(p) && size(v,1) == size(c, 1) && size(c, 2) == 1
        cm = get(gcf, 'Colormap');
        c2 = round((c./max(c)) .* length(cm));
        c2(c2==0)=1;
        c2 = cm(c2, :);
        c = c2;  % this is the true-color equivalent
        clear cm c2
        set(p, 'FaceVertexCData', c);
    end

    bad = any(size(c) - size(v)); % this could occur for various reasons?

    if bad
        c = repmat(basecolor, size(v, 1), 1);
    end

% c = true-color vertices x 3
% v = xyz coords of vertices x 3
% cscale: mapped true colors or index for  xyz, cell with cscale{1} = xyz coords to color x 3 or x 1

    % -----------------------------------------------------------------------
    % * main xyz color change
    % -----------------------------------------------------------------------

    t1 = clock;
    if doverbose, fprintf('Main color vertices: '), end

    c = change_colors(c, xyz, v, mind, cscale, actcolor, p, alphascale, doverbose);
    drawnow();

    % -----------------------------------------------------------------------
    % * additional optional vertices
    % -----------------------------------------------------------------------
    if exist('vv', 'var')
        for j = 1:length(vv)
            if doverbose, fprintf('\nAdditional vertices: '), end

            % cscale:
            % pass in vector of ones length coords for solid color
            % or scalar vals for color mapping
            % pass in cell array

            % cc{j} is color, cscale{j+1} is scaling vals for coords
            c = change_colors(c, vv{j}, v, mind, cscale(j+1), cc{j}, p, alphascale, doverbose);
        end
        drawnow();
    end

    % -----------------------------------------------------------------------
    % * figure out which vertices should be colored with ocol (overlap color)
    % -----------------------------------------------------------------------

    if exist('ocol', 'var') && exist('xyzb', 'var') && ~isempty(cat(1, xyzb{:}))

        alld = zeros(size(v, 1), 1);  % keeps track of overlap vertices
        xyzb = cat(1, xyzb{:});

        t1 = clock;
        if doverbose, fprintf('\nOverlap vertices: '), end

        cscaletmp = {ones(size(xyzb, 1), 1)};
        c = change_colors(c, xyzb, v, mind, cscaletmp, ocol, p, alphascale, doverbose);
        drawnow();

        if doverbose,  fprintf('%3.0f.done in %3.0f s\n', i, etime(clock, t1)), end
    end

    % -----------------------------------------------------------------------
    % * figure out which vertices should be colored with acol (all color)
    % -----------------------------------------------------------------------
    if exist('acol', 'var') && ~isempty(xyza) && length(vv)>1

        xyzall = xyza;
        for i = 2:length(vv)
            xyzall = intersect(xyzall, vv{i}, 'rows');
        end

        t1 = clock;
        if doverbose, fprintf('\nAll overlap vertices: '), end

        cscaletmp = {ones(size(xyzall, 1), 1)};
        c = change_colors(c, xyzall, v, mind, cscaletmp, acol, p, alphascale, doverbose);
        drawnow();

        if doverbose, fprintf('%3.0f.done in %3.0f s\n', i, etime(clock, t1)), end

    end

    % -----------------------------------------------------------------------
    % * final color change
    % -----------------------------------------------------------------------

    if exist('p', 'var')
        set(p, 'FaceColor', 'interp')
        set(p, 'FaceVertexCData', c)
        drawnow()
    end

    %lightFollowView
end

% -----------------------------------------------------------------------
% * SUB-FUNCTIONS
% -----------------------------------------------------------------------

function c = change_colors(c, coords, v, mind, cscale, actcolor, p, alphascale, doverbose)

% c is list of colors. by default, usually [.5 .5 .5] for all colors
% (basecolor)
    global do_fsavg_left do_fsavg_right ras

    if isempty(coords), if doverbose, disp('Coords is empty. Nothing to plot.'), end, return, end

    % select surface vertices that are close to coords
    % vertices that are far will be omitted
    % smallv is a reduced set of vertices
    % --------------------------------------------------------------------
    cmax = max(coords, [], 1);
    cmin = min(coords, [], 1);
    if doverbose, fprintf('%3.0f vertices.  selecting: ', size(v, 1)); end

    wh = any(v - repmat(cmax, size(v, 1), 1) > mind, 2);
    wh2 = any(repmat(cmin, size(v, 1), 1) - v > mind, 2);

    % list vertices to test and possibly change color
    whverts = (1:size(v, 1))';
    whverts(wh | wh2) = [];     % indices of smallv in big (original) list
    if do_fsavg_left || do_fsavg_right
        smallv = ras.ras';
    else
        smallv = v;
    end
    smallv(wh | wh2,:) = [];     % vertices--restricted list

    if doverbose, fprintf('%3.0f\n', size(whverts, 1)); end

    if isempty(smallv), return, end

    % select coords that are close to surface vertices5
    % coords, cscale{1}, and alphascale{1} all have same indices
    % ---------------------------------------------------------------------
    cmax = max(smallv, [], 1);
    cmin = min(smallv, [], 1);
    if doverbose, fprintf('%3.0f coords.  selecting: ', size(coords, 1)); end

    % omit wh and wh2 - outside scope of this coord set
    wh = any(coords - repmat(cmax, size(coords, 1), 1) > mind, 2);
    wh2 = any(repmat(cmin, size(coords, 1), 1) - coords > mind, 2);

    % if cscale is matrix, must select these values of cscale as well!
    if length(cscale) > 0 && size(cscale{1}, 1) == size(coords, 1)
        cscale{1}(wh | wh2,:) = [];
    end

    if length(alphascale) > 0 && size(alphascale{1}, 1) == size(coords, 1)
        alphascale{1}(wh | wh2,:) = [];
    end

    coords(wh | wh2,:) = [];

    if isempty(coords), return, end

    nc = size(coords, 1);
    if doverbose, fprintf('%3.0f\n', nc); end

    % break up coords into list and run
    % ----------------------------------

    % break up coords into list
    xyz2 = {}; indx = 1;
    for kk = 1:1000:nc
        setwh{indx} = (kk:min(nc, kk + 1000 - 1))';
        xyz2{indx} = coords(setwh{indx},:);

        indx = indx + 1;
    end

    if doverbose, fprintf('Running %3.0f sets of coordinates: 000', length(xyz2)); end

    indxval = 1;
    wh_coords_near_surface = false(size(coords, 1), 1);

    for setno = 1:length(xyz2)
        if doverbose, fprintf('\b\b\b%03d', setno); end

        for i = 1:size(xyz2{setno}, 1)
            % i indexes coordinate; e.g., 1000 iterations for a typical
            % coord set
            % find vertices that are within range of point i in set  setno
            % vertex_indices are indices into the ORIGINAL list
            vertex_indices = find_in_radius(xyz2, setno, i, smallv, mind, whverts);

            % two modes:
            % - if cscale{1} is a matrix, treats as rgb values, and put in
            % color stored in cscale.  actcolor is ignored.
            % - if cscale{1} is a vector, treat it as a scaling value for actcolor
            % In either case, indxval should index location of coordinate in FULL
            % list (corresponding to full list in cscale)

            % cscale{1}(i) and alphascale{1}(i) contain unique values for each coord i
            % these are used to map colors for single point i to multiple
            % nearby vertices vertex_indices

            if length(cscale) > 0 && ~isempty(cscale{1})
                mycscale = cscale{1}(setwh{setno}(i), :);  % color for this point, index into full list of coords/scalevals
            else
                mycscale = [];
            end

            if length(alphascale) > 0 && ~isempty(alphascale{1})
                myalphascale = alphascale{1}(setwh{setno}(i), :);  % alpha for this point, index into full list of coords/scalevals

                % adjust by cube of mind, to adjust for overlap in vertices
                % affected by adjacent coords - done now in cluster_surf
                %myalphascale = myalphascale ./ mind^3;

            else
                myalphascale = [];
            end

            c = color_change_vertices(c, mycscale, actcolor, vertex_indices, myalphascale);

            if ~isempty(vertex_indices), wh_coords_near_surface(indxval) = 1; end

            indxval = indxval + 1;
        end
    end % setno

    if exist('p', 'var')
        set(p, 'FaceColor', 'interp')
        set(p, 'FaceVertexCData', c)
        drawnow
    end

end % change_colors

function z = dist_tmp(w, p)

    %
    % if isstr(w)
    %   switch (w)
    %     case 'deriv',
    %       z = '';
    %     otherwise
    %       error('Unrecognized code.')
    %   end
    %   return
    % end

    % CALCULATION
    if nargin == 1
        p = w;
        w = w';
    end

    [S, R] = size(w);
    [R2, Q] = size(p);
    if (R ~= R2), error('Inner matrix dimensions do not match.'), end

    z = zeros(S, Q);
    if (Q<S)
        p = p';
        copies = zeros(1, S);
        for q=1:Q
            z(:,q) = sum((w-p(q+copies,:)).^2, 2);
        end
    else
        w = w';
        copies = zeros(1, Q);
        for i=1:S
            z(i,:) = sum((w(:,i+copies)-p).^2, 1);
        end
    end
    z = sqrt(z);

end

function [vertex_indices, d] = find_in_radius(xyz2, setno, i, smallv, mind, whverts)
    % output: indices of vertices in BIG list, and distances

    % get vertices v within box -- fast method
    wh = find(all(abs(bsxfun(@minus, xyz2{setno}(i,:), smallv)) <= mind, 2));
    d = dist_tmp(smallv(wh,:), xyz2{setno}(i,:)');

    % convert back to big list
    vertex_indices = whverts(wh(d < mind));
end

function c = color_change_vertices(c, mycscale, actcolor, vertex_indices, myalphascale)

% pass in/out larger list c, but update only a few coordinates
% indexed by vertex_indices with a particular color value.

% cscale is ordered from least color
% change to most for graphical/aesthetic purposes, so this function is
% called repetitively to index lower i (less change) to higher i (more
% change) successively.

% two modes: if cscale{1} is a matrix, treats as rgb values, and put in
% color stored in cscale
% if cscale{1} is a vector, treat it as a scaling value for actcolor
% In either case, the ith index point should correspond to the
% coordinate being worked on.

n = length(vertex_indices);

if n
    if length(mycscale) > 1
        % if we have rgb values rather than scaling values
        % this occurs if heatmap = yes and colorscale = no, we pass in rgb values
        % actcolor is ignored

        if isempty(myalphascale)
            % solid colors
            c(vertex_indices,:) = repmat(mycscale, n, 1);
        else
            w = myalphascale;  % the weight for new color vs. old.
            c(vertex_indices, :) = (1-w) .* c(vertex_indices, :) + w .* repmat(mycscale, n, 1);
        end

    else
        % this will scale the activation color in proportion to cscale
        % for each voxel

        w = mycscale;  % the weight for new color vs. old.  1 = all new, 0 = all old
        c(vertex_indices, :) = (1-w) .* c(vertex_indices, :) + w .* repmat(actcolor, n, 1);
    end
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function complete_run(ssh_path, eddies_save_path, tracks_save_path, fake_eddies, viewer_data_save_path, viewer_path, varargin)
% NOTE: Please provide full path names for all paths.
%
% The purpose of this function is to automate the entire gauntlet of
% running all of the code to go from having only MATLAB SSH data to getting
% the data set up to be viewed on the tracks viewer. This function should
% automate the following steps:
% 1. Scan SSH data for eddies
% 2. Save detected eddies
% 3. Build LNN tracks for these eddies (Possibly with fake eddies, if the
%    user of the script wants fake eddies)
% 4. Format and process saved data so that it can be viewed in Hung's
%    viewer that we include
%
% NOTE: The process of turning NetCDF file data into MATLAB data is not
% accomplished in this function.
% The function set_up_ssh_data is provided that will accomplish this for
% the user, but the user will have to provide the variable names that match
% the function inputs. Refer to the comments at the top of the function
% set_up_ssh_data for a description of what each of the required variables are.
%
% This script will not open the viewer for the user.
% To open the viewer, run the start_track_viewer script inside the
% tracks_viewer directory.

p = inputParser;
defaultLat = NaN;
defaultLon = NaN;
defaultDates = NaN;
defaultAreaMap = NaN;
defaultScanType = 'v2';
defaultMinPixelSize = 9;
defaultThresholdStep = 0.05;
defaultSSHUnits = 'centimeters';
defaultPaddingFlag = true;

addRequired(p, 'ssh_dir');
addRequired(p, 'eddies_save_path');
addRequired(p, 'fake_eddies');

addParameter(p, 'lat', defaultLat);
addParameter(p, 'lon', defaultLon);
addParameter(p, 'dates', defaultDates);
addParameter(p, 'area_map', defaultAreaMap);
addParameter(p, 'scan_type', defaultScanType);
addParameter(p, 'minimumArea', defaultMinPixelSize, @isnumeric);
addParameter(p, 'thresholdStep', defaultThresholdStep, @isnumeric);
addParameter(p, 'isPadding', defaultPaddingFlag);
addParameter(p, 'sshUnits', defaultSSHUnits);

parse(p, ssh_path, eddies_save_path, fake_eddies, varargin{:});

lat = p.Results.lat;
lon = p.Results.lon;
dates = p.Results.dates;
area_map = p.Results.area_map;
scan_type = p.Results.scan_type;
minimumArea = p.Results.minimumArea;
thresholdStep = p.Results.thresholdStep;
isPadding = p.Results.isPadding;
SSH_Units = p.Results.sshUnits;
if ~strcmp(ssh_path(end), '/')
    ssh_path = strcat(ssh_path, '/');
end
if ~strcmp(eddies_save_path(end), '/')
    eddies_save_path = strcat(eddies_save_path, '/');
end
if ~strcmp(tracks_save_path(end), '/')
    tracks_save_path = strcat(tracks_save_path, '/');
end
if ~strcmp(viewer_data_save_path(end), '/')
    viewer_data_save_path = strcat(viewer_data_save_path, '/');
end
ssh_names = get_ssh_names(ssh_path);
if isnan(lat)
    vars = load([ssh_path, 'lat.mat']);
    names = fieldnames(vars);
    lat = vars.(names{1});
end
if isnan(lon)
    vars = load([ssh_path, 'lon.mat']);
    names = fieldnames(vars);
    lon = vars.(names{1});
end
if isnan(dates)
    vars = load([ssh_path, 'dates.mat']);
    names = fieldnames(vars);
    dates = vars.(names{1});
end
if isnan(area_map)
    vars = load([ssh_path, 'area_map.mat']);
    names = fieldnames(vars);
    area_map = vars.(names{1});
end
if ~exist(eddies_save_path, 'dir')
    mkdir(eddies_save_path);
end
old_path = cd('eddyscan/');
disp('Scanning SSH data for eddies');
for i = 1:length(ssh_names)
    disp(['Iteration ', num2str(i)]);
    vars = load(ssh_names{i});
    names = fieldnames(vars);
    ssh_data = vars.(names{1});
    ant_eddies = scan_single(ssh_data, lat, lon, dates(i), 'anticyc', scan_type, area_map, 'sshUnits', SSH_Units,...
                             'thresholdStep', thresholdStep, 'minimumArea', minimumArea, 'isPadding', isPadding);
    par_save([eddies_save_path, 'anticyc_', num2str(dates(i)), '.mat'], ant_eddies);
    cyc_eddies = scan_single(ssh_data, lat, lon, dates(i), 'cyclonic', scan_type, area_map, 'sshUnits', SSH_Units,...
                             'thresholdStep', thresholdStep, 'minimumArea', minimumArea, 'isPadding', isPadding);
    par_save([eddies_save_path, 'cyclonic_', num2str(dates(i)), '.mat'], cyc_eddies);
end
cd(old_path);
% Track eddies and save tracks
if ~exist(tracks_save_path, 'dir')
    mkdir(tracks_save_path);
end
if ~exist(viewer_data_save_path, 'dir')
    mkdir(viewer_data_save_path);
end
old_path = cd('track_lnn/');
disp('Tracking eddies');
if strcmp(fake_eddies, 'yes')
    anticyclonic_tracks = tolerance_track_lnn(eddies_save_path, 'anticyc', 1, 1, minimumArea);
    cyclonic_tracks = tolerance_track_lnn(eddies_save_path, 'cyclonic', 1, 1, minimumArea);
    modified_eddies_path = [viewer_data_save_path, 'modified_eddies/'];
    modified_tracks_path = [viewer_data_save_path, 'modified_tracks/'];
    if ~exist(modified_eddies_path, 'dir')
        mkdir(modified_eddies_path);
    end
    if ~exist(modified_tracks_path, 'dir')
        mkdir(modified_tracks_path);
    end
    disp('Modifying eddy data to include fake eddies (for the sake of the tracks)');
    process_eddies_and_tracks_tolerance('anticyc', dates, eddies_save_path, anticyclonic_tracks, modified_eddies_path, modified_tracks_path);
    process_eddies_and_tracks_tolerance('cyclonic', dates, eddies_save_path, cyclonic_tracks, modified_eddies_path, modified_tracks_path);
    vars = load([modified_tracks_path, 'anticyc_tracks_processed.mat']);
    s = fieldnames(vars);
    anticyclonic_tracks = vars.(s{1});
    vars = load([modified_tracks_path, 'cyclonic_tracks_processed.mat']);
    s = fieldnames(vars);
    cyclonic_tracks = vars.(s{1});
    eddies_save_path = modified_eddies_path;
else
    anticyclonic_tracks = track_lnn(eddies_save_path, 'anticyc', 1); % Possibly change this so it can be an optional argument in input parser
    cyclonic_tracks = track_lnn(eddies_save_path, 'cyclonic', 1); % Same as above
end
save([tracks_save_path, 'anticyclonic_tracks.mat'], 'anticyclonic_tracks');
save([tracks_save_path, 'cyclonic_tracks.mat'], 'cyclonic_tracks');
eddies_t = reformat_track_data_to_chelton(eddies_save_path, dates, anticyclonic_tracks, cyclonic_tracks);%#ok
save([tracks_save_path, 'chelton_structured_tracks.mat'], 'eddies_t');
% Set up data for viewer
disp('Preparing data for viewer');
prepare_eddy_data_for_viewer(viewer_data_save_path, dates, eddies_save_path, cyclonic_tracks, anticyclonic_tracks, ssh_path);
cd(old_path);
% Deleting old viewer files
files = dir([viewer_path, 'SSH/']);
old_path = cd([viewer_path, 'SSH/']);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        if files(i).isdir
            rmdir(file, 's');
        else
            delete(file);
        end
    end
end
files = dir([viewer_path, 'track_data/']);
cd([viewer_path, 'track_data/']);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        if files(i).isdir
            rmdir(file, 's');
        else
            delete(file);
        end
    end
end
files = dir([viewer_path, 'attributes/']);
cd([viewer_path, 'attributes/']);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        if files(i).isdir
            rmdir(file, 's');
        else
            delete(file);
        end
    end
end
files = dir([viewer_path, 'pixels_data/']);
cd([viewer_path, 'pixels_data/']);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        if files(i).isdir
            rmdir(file, 's');
        else
            delete(file);
        end
    end
end
cd(old_path);
% Copying new viewer files into place
files = dir(ssh_path);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        copyfile([ssh_path, file], [viewer_path, 'SSH/', file]);
    end
end
files = dir(tracks_save_path);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        copyfile([tracks_save_path, file], [viewer_path, 'track_data/', file]);
    end
end
files = dir(viewer_data_save_path);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        if ~files(i).isdir
            if strcmp(file, 'anticyclonic_tags.mat') || strcmp(file, 'cyclonic_tags.mat')
                copyfile([viewer_data_save_path, file], [viewer_path, 'track_data/', file]);
            else
                copyfile([viewer_data_save_path, file], [viewer_path, 'attributes/', file]);
            end
        end
    end
end
files = dir([viewer_data_save_path, 'pixel_data/']);
for i = 1:length(files)
    file = files(i).name;
    if ~strcmp(file, '.') && ~strcmp(file, '..')
        copyfile([viewer_data_save_path, 'pixel_data/', file], [viewer_path, 'pixels_data/', file]);
    end
end

end

function par_save(path, eddies)%#ok
save(path, 'eddies');
end

function [ssh_names] = get_ssh_names(path)
% path is the path to the eddies directory
% type is anticyclonic or cyclonic
if ~strcmp(path(end), '/')
    path = strcat(path, '/');
end
files = dir(path);
x = 0;
for i = 1:length(files)
    if ~isempty(strfind(files(i).name, 'ssh_'))
        x = x + 1;
    end
end
ssh_names = cell(x, 1);
x = 1;
for i = 1:length(files)
    if files(i).isdir && ~isequal(files(i).name, '.') && ~isequal(files(i).name, '..')
        rec_names = get_ssh_names([path, files(i).name, '/']);
        for j = 1:length(rec_names)
            ssh_names{x} = rec_names{j};
            x = x + 1;
        end
        continue;
    end
    file = files(i).name;
    [~, name, ext] = fileparts([path, file]);
    if ~isempty(strfind(name, 'ssh_')) && strcmp(ext, '.mat')
        ssh_names{x} = [path, file];
        x = x + 1;
    end
end
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function test_tutorial_timefrequencyanalysis(datadir)

% MEM 2gb
% WALLTIME 00:20:00
% DEPENDENCY ft_freqanalysis ft_preprocessing ft_multiplotTFR ft_singleplotTFR

if nargin==0
  datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/timefrequencyanalysis');
end

load(fullfile(datadir, 'dataFIC.mat'));

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'MEG';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 2:2:30;                         % analysis 2 to 30 Hz in steps of 2 Hz
cfg.t_ftimwin    = ones(length(cfg.foi),1).*0.5;   % length of time window = 0.5 sec
cfg.toi          = -0.5:0.05:1.5;                  % time window "slides" from -0.5 to 1.5 sec in steps of 0.05 sec (50 ms)
TFRhann = ft_freqanalysis(cfg, dataFIC);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim         = [-3e-27 3e-27];
cfg.showlabels   = 'yes';
cfg.layout       = 'CTF151.lay';
figure
ft_multiplotTFR(cfg, TFRhann);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.maskstyle    = 'saturation';
cfg.zlim         = [-3e-27 3e-27];
cfg.channel      = 'MRC15';
figure
ft_singleplotTFR(cfg, TFRhann);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.maskstyle    = 'saturation';
cfg.xlim         = [0.9 1.3];
cfg.zlim         = [-1.5e-27 1.5e-27];
cfg.ylim         = [15 20];
cfg.showlabels   = 'markers';
figure
ft_topoplotTFR(cfg, TFRhann);

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'MRC15';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 2:1:30;
cfg.t_ftimwin    = 7./cfg.foi;  % 7 cycles per time window
cfg.toi          = -0.5:0.05:1.5;
TFRhann7 = ft_freqanalysis(cfg, dataFIC);

cfg              = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.maskstyle    = 'saturation';
cfg.zlim         = [-3e-27 3e-27];
cfg.channel      = 'MRC15';
cfg.interactive  = 'no';
figure
ft_singleplotTFR(cfg, TFRhann7);

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'MRC15';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 2:1:30;
cfg.t_ftimwin    = 4./cfg.foi;
cfg.toi          = -0.5:0.05:1.5;
TFRhann4 = ft_freqanalysis(cfg, dataFIC);

cfg.t_ftimwin    = 5./cfg.foi;
TFRhann5 = ft_freqanalysis(cfg, dataFIC);

cfg.t_ftimwin    = 10./cfg.foi;
TFRhann10 = ft_freqanalysis(cfg, dataFIC);

cfg = [];
cfg.output     = 'pow';
cfg.channel    = 'MEG';
cfg.method     = 'mtmconvol';
cfg.foi        = 1:2:30;
cfg.t_ftimwin  = 5./cfg.foi;
cfg.tapsmofrq  = 0.4 *cfg.foi;
cfg.toi        = -0.5:0.05:1.5;
cfg.pad        = 'maxperlen';
TFRmult = ft_freqanalysis(cfg, dataFIC);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim         = [-3e-27 3e-27];
cfg.showlabels   = 'yes';
cfg.layout       = 'CTF151.lay';
figure
ft_multiplotTFR(cfg, TFRmult)

cfg = [];
cfg.channel    = 'MEG';
cfg.method     = 'wavelet';
cfg.width      = 7;
cfg.output     = 'pow';
cfg.foi        = 1:2:30;
cfg.toi        = -0.5:0.05:1.5;
TFRwave = ft_freqanalysis(cfg, dataFIC);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim         = [-3e-25 3e-25];
cfg.showlabels   = 'yes';
cfg.layout       = 'CTF151.lay';
figure
ft_multiplotTFR(cfg, TFRwave)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%TRAIN_MODEL Summary of this function goes here
%   Function: train face alignment model
%   Detailed explanation goes here
%   Input:
%       dbnames: the names of database
% Configure the parameters for training model
global params;
global Tr_Data;
config_tr;

if size(dbnames) > 1 & sum(strcmp(dbnames, 'COFW')) > 0
    disp('Sorry, COFW cannnot be combined with others')
    return;
end

if sum(strcmp(dbnames, 'COFW')) > 0
    load('../initial_shape/InitialShape_29.mat');
    params.meanshape        = S0;
else
    load('../initial_shape/InitialShape_68.mat');
    params.meanshape        = S0;
end

if params.isparallel
    disp('Attention, if error occurs, plese ensure you used the correct version of parallel initialization.');

    if isempty(gcp('nocreate'))
        parpool(4);
    else
        disp('Already initialized');
    end

    %{
    if matlabpool('size') <= 0
        matlabpool('open','local',4);
    else
        disp('Already initialized');
    end
    %}
end

% load trainning data from hardware
Tr_Data    = [];
% Tr_Bboxes  = [];
for i = 1:length(dbnames)
    % load training samples (including training images, and groundtruth shapes)
    imgpathlistfile = strcat('..\datasets\', dbnames{i}, '\Path_Images.txt');
    tr_data = loadsamples(imgpathlistfile, 2);
    Tr_Data = [Tr_Data; tr_data];
end

% Augmentate data for traing: assign multiple initial shapes to each image
Data = Tr_Data; % (1:10:end);
Param = params;

if Param.flipflag % if conduct flipping
    Data_flip = cell(size(Data, 1), 1);
    for i = 1:length(Data_flip)
        Data_flip{i}.img_gray    = fliplr(Data{i}.img_gray);
        Data_flip{i}.width_orig  = Data{i}.width_orig;
        Data_flip{i}.height_orig = Data{i}.height_orig;
        Data_flip{i}.width       = Data{i}.width;
        Data_flip{i}.height      = Data{i}.height;

        Data_flip{i}.shape_gt    = flipshape(Data{i}.shape_gt);
        Data_flip{i}.shape_gt(:, 1)  =  Data{i}.width - Data_flip{i}.shape_gt(:, 1);

        Data_flip{i}.bbox_gt        = Data{i}.bbox_gt;
        Data_flip{i}.bbox_gt(1)     = Data_flip{i}.width - Data_flip{i}.bbox_gt(1) - Data_flip{i}.bbox_gt(3);

        Data_flip{i}.bbox_facedet        = Data{i}.bbox_facedet;
        Data_flip{i}.bbox_facedet(1)     = Data_flip{i}.width - Data_flip{i}.bbox_facedet(1) - Data_flip{i}.bbox_facedet(3);
    end
    Data = [Data; Data_flip];
end

% choose corresponding points for training
for i = 1:length(Data)
    Data{i}.shape_gt = Data{i}.shape_gt(Param.ind_usedpts, :);
    Data{i}.bbox_gt = getbbox(Data{i}.shape_gt);

    % modify detection boxes
    shape_facedet = resetshape(Data{i}.bbox_facedet, Param.meanshape);
    shape_facedet = shape_facedet(Param.ind_usedpts, :);
    Data{i}.bbox_facedet = getbbox(shape_facedet);
end

Param.meanshape        = S0(Param.ind_usedpts, :);

dbsize = length(Data);

% load('Ts_bbox.mat');

augnumber = Param.augnumber;

for i = 1:dbsize
    % initializ the shape of current face image by randomly selecting multiple shapes from other face images
    % indice = ceil(dbsize*rand(1, augnumber));

    indice_rotate = ceil(dbsize*rand(1, augnumber));
    indice_shift  = ceil(dbsize*rand(1, augnumber));
    scales        = 1 + 0.2*(rand([1 augnumber]) - 0.5);

    Data{i}.intermediate_shapes = cell(1, Param.max_numstage);
    Data{i}.intermediate_bboxes = cell(1, Param.max_numstage);

    Data{i}.intermediate_shapes{1} = zeros([size(Param.meanshape), augnumber]);
    Data{i}.intermediate_bboxes{1} = zeros([augnumber, size(Data{i}.bbox_gt, 2)]);

    Data{i}.shapes_residual = zeros([size(Param.meanshape), augnumber]);
    Data{i}.tf2meanshape = cell(augnumber, 1);
    Data{i}.meanshape2tf = cell(augnumber, 1);

    % if Data{i}.isdet == 1
    %    Data{i}.bbox_facedet = Data{i}.bbox_facedet*ts_bbox;
    % end
    for sr = 1:params.augnumber
        if sr == 1
            % estimate the similarity transformation from initial shape to mean shape
            % Data{i}.intermediate_shapes{1}(:,:, sr) = resetshape(Data{i}.bbox_gt, Param.meanshape);
            % Data{i}.intermediate_bboxes{1}(sr, :) = Data{i}.bbox_gt;
            Data{i}.intermediate_shapes{1}(:,:, sr) = resetshape(Data{i}.bbox_facedet, Param.meanshape);
            Data{i}.intermediate_bboxes{1}(sr, :) = Data{i}.bbox_facedet;

            meanshape_resize = resetshape(Data{i}.intermediate_bboxes{1}(sr, :), Param.meanshape);

            Data{i}.tf2meanshape{1} = fitgeotrans(bsxfun(@minus, Data{i}.intermediate_shapes{1}(1:end,:, 1), mean(Data{i}.intermediate_shapes{1}(1:end,:, 1))), ...
                (bsxfun(@minus, meanshape_resize(1:end, :), mean(meanshape_resize(1:end, :)))), 'NonreflectiveSimilarity');
            Data{i}.meanshape2tf{1} = fitgeotrans((bsxfun(@minus, meanshape_resize(1:end, :), mean(meanshape_resize(1:end, :)))), ...
                bsxfun(@minus, Data{i}.intermediate_shapes{1}(1:end,:, 1), mean(Data{i}.intermediate_shapes{1}(1:end,:, 1))), 'NonreflectiveSimilarity');

            % calculate the residual shape from initial shape to groundtruth shape under normalization scale
            shape_residual = bsxfun(@rdivide, Data{i}.shape_gt - Data{i}.intermediate_shapes{1}(:,:, 1), [Data{i}.intermediate_bboxes{1}(1, 3) Data{i}.intermediate_bboxes{1}(1, 4)]);
            % transform the shape residual in the image coordinate to the mean shape coordinate
            [u, v] = transformPointsForward(Data{i}.tf2meanshape{1}, shape_residual(:, 1)', shape_residual(:, 2)');
            Data{i}.shapes_residual(:, 1, 1) = u';
            Data{i}.shapes_residual(:, 2, 1) = v';
        else
            % randomly rotate the shape
            % shape = resetshape(Data{i}.bbox_gt, Param.meanshape);       % Data{indice_rotate(sr)}.shape_gt
            shape = resetshape(Data{i}.bbox_facedet, Param.meanshape);       % Data{indice_rotate(sr)}.shape_gt

            if params.augnumber_scale ~= 0
                shape = scaleshape(shape, scales(sr));
            end

            if params.augnumber_rotate ~= 0
                shape = rotateshape(shape);
            end

            if params.augnumber_shift ~= 0
                shape = translateshape(shape, Data{indice_shift(sr)}.shape_gt);
            end

            Data{i}.intermediate_shapes{1}(:, :, sr) = shape;
            Data{i}.intermediate_bboxes{1}(sr, :) = getbbox(shape);

            meanshape_resize = resetshape(Data{i}.intermediate_bboxes{1}(sr, :), Param.meanshape);

            Data{i}.tf2meanshape{sr} = fitgeotrans(bsxfun(@minus, Data{i}.intermediate_shapes{1}(1:end,:, sr), mean(Data{i}.intermediate_shapes{1}(1:end,:, sr))), ...
                bsxfun(@minus, meanshape_resize(1:end, :), mean(meanshape_resize(1:end, :))), 'NonreflectiveSimilarity');
            Data{i}.meanshape2tf{sr} = fitgeotrans(bsxfun(@minus, meanshape_resize(1:end, :), mean(meanshape_resize(1:end, :))), ...
                bsxfun(@minus, Data{i}.intermediate_shapes{1}(1:end,:, sr), mean(Data{i}.intermediate_shapes{1}(1:end,:, sr))), 'NonreflectiveSimilarity');

            shape_residual = bsxfun(@rdivide, Data{i}.shape_gt - Data{i}.intermediate_shapes{1}(:,:, sr), [Data{i}.intermediate_bboxes{1}(sr, 3) Data{i}.intermediate_bboxes{1}(sr, 4)]);
            [u, v] = transformPointsForward(Data{i}.tf2meanshape{1}, shape_residual(:, 1)', shape_residual(:, 2)');
            Data{i}.shapes_residual(:, 1, sr) = u';
            Data{i}.shapes_residual(:, 2, sr) = v';
            % Data{i}.shapes_residual(:, :, sr) = tformfwd(Data{i}.tf2meanshape{sr}, shape_residual(:, 1), shape_residual(:, 2));
        end
    end
end

% train random forests for each landmark
randf = cell(Param.max_numstage, 1);
Ws    = cell(Param.max_numstage, 1);

%{
if nargin > 2
    n = size(LBFRegModel_initial.ranf, 1);
    for i = 1:n
      randf(1:n, :) = LBFRegModel_initial.ranf;
    end
end
%}

for s = 1:Param.max_numstage
    % learn random forest for s-th stage
    disp('train random forests for landmarks...');

    %{
    if isempty(randf{s})
        if exist(strcat('randfs\randf', num2str(s), '.mat'))
            load(strcat('randfs\randf', num2str(s), '.mat'));
        else
            tic;
            randf{s} = train_randomfs(Data, Param, s);
            toc;
            save(strcat('randfs\randf', num2str(s), '.mat'), 'randf', '-v7.3');
        end
    end
    %}
    if isempty(randf{s})
        tic;
        randf{s} = train_randomfs(Data, Param, s);
        toc;
    end

    % derive binary codes given learned random forest in current stage
    disp('extract local binary features...');

    tic;
    binfeatures = derivebinaryfeat(randf{s}, Data, Param, s);
    % save(strcat('LBFeats\LBFeats', num2str(s), '.mat'), 'binfeatures', '-v7.3');
    toc;

    % learn global linear regrassion given binary feature
    disp('learn global regressors...');
    tic;
    [W, Data] = globalregression(binfeatures, Data, Param, s);
    Ws{s} = W;
    % save(strcat('Ws\W', num2str(s), '.mat'), 'W', '-v7.3');
    toc;

end

LBFRegModel.ranf = randf;
LBFRegModel.Ws   = Ws;

author: jwyang

repo: face-alignment


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% compile cpp files
% change the include and lib path if necessary
function compile

include = ' -I/usr/local/include/opencv/ -I/usr/local/include/ -I/usr/include/opencv/';
lib = ' -lopencv_core -lopencv_highgui -lopencv_imgproc -lopencv_video';
eval(['mex lk.cpp -O' include lib]);

mex distance.cpp
mex imResampleMex.cpp
mex warp.cpp

disp('Compilation finished.');

author: yuxng


text:
function slide30

	a = [0 0 0;
		1 0 0;
		1 1 0;
		0 1 0;
		0 0 1;
		1 0 1;
		1 1 1;
		0 1 1];

	b = [1 2 6 5;
		2 3 7 6;
		3 4 8 7;
		4 1 5 8;
		1 2 3 4;
		5 6 7 8];

	p1 = patch('faces',b,...
		'vertices',a,...
		'facecolor',[.5 .5 .5],...
		'edgecolor',[1,1,1],...
		'facealpha',0.5);

	view(3)
	axis square

end


text:
function out = Rstar(d, n, e, fun)
%finds rate that achieves a given excess distortion at a given blocklength for
%Gaussian source with unit variance
%d - excess distortion
%n - block length (scalar)
%e - excess probability
%fun - which function to use for calculation

%
%   Created in 2012 by Victoria Kostina (vkostina@caltech.edu)
%

%starting points - 'persistent' to make optimization faster
persistent x0;
%rate distortion
Rd = -1/2*log2(d);
%precision options
tol = 1e-10;
options = optimset('TolX',tol, 'MaxFunEvals', 500, 'Algorithm', 'active-set', 'Display', 'off');

switch lower(fun)
    case {'shannon', 'spherecoveringa', 'generalc', 'generalcopt'}
        if isempty(x0)
            x0 = [Rd 3*Rd]; %rate
        end
        out = Generic();
    case 'normal'
        %normal approximation:
        out = Normal();
    case 'spherecoveringc'
        %converse via sphere covering
        out = SphereCoveringC();
    case 'rogersa'
        out = RogersA();
    otherwise
        disp('Unknown type.')
end

%--------------------------------------------------------------------------
    function out = Generic()
        out = fzero(@(x)goal(x) - e, x0, options);
        function out = goal(x)
            out = Pexcess(x, d, n, fun);
        end
    end

%--------------------------------------------------------------------------
    function out = Normal()
        %Normal approximation
        out = - .5*log2(d) + 1/sqrt(2*n)*Qinv(e)*log2(exp(1));
    end

%--------------------------------------------------------------------------
    function out = SphereCoveringC()
        %converse via sphere covering

        %find rstar:
        factor = 3;
        rstar = fzero(@(x)chi2cdf(x^2*n, n) - chi2cdf(factor*n, n) + e, 1);
        out = rate(d);

        function out = rate(x)
            out = 1/2*log2(rstar^2/x);
        end
    end

%--------------------------------------------------------------------------
    function out = RogersA()
        rstar = fzero(@(x)1 - chi2cdf(x^2*n, n) - e, 1);
        out = rogers(rstar/sqrt(d));
        function out = rogers(r)
            if r >= n
                out = log2(exp(1))+ log2((n*log(n)+n*log(log(n)) + 5*n)) + n*log2(r);
            elseif r >= n/log(n)
                out = log2(n*(n*log(n)+n*log(log(n)) + 5*n)) + n*log2(r);
            elseif r >2
                out = log2(7^(4*log(7)/7)/4*sqrt(2*pi)) + log2(n^(3/2)*((n-1)*log(r*n)+(n-1)*log(log(n)+log(n)/2+log((pi*sqrt(2*n))/(sqrt(pi*n)-2))))) - log2(r*(1-2/log(n))*(1 - 2/sqrt(pi*n))*(log(n))^2) + n*log2(r);
            else
                out = log2(sqrt(2*pi)) + log2(sqrt(n)*((n-1)*log(r*n)+(n-1)*log(log(n)+log(n)/2+log((pi*sqrt(2*n))/(sqrt(pi*n)-2))))) - log2(r*(1-2/log(n))*(1 - 2/sqrt(pi*n))) + n*log2(r);
            end
            out = out/n;
        end
    end

end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
% std_topomovie - make movie in the frequency domain
%
% Usage:
%  >> [STUDY] = std_movie(STUDY, ALLEEG, key1, val1, key2, val2, ...);
%
% Inputs:
%   STUDY      - STUDY structure comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - vector of EEG dataset structures for the dataset(s) in the STUDY,
%                typically created using load_ALLEEG().
%   'channels' - [numeric vector]  specific channel group to plot. By
%                default, the grand mean channel spectrum is plotted (using the
%                same format as for the cluster component means described above)
%   'moviemode' - ['erp'|'spec'|'ersptime'] movie mode. Currently only
%                 'spec' is implemented.
%   'erspfreq'  - [min max] frequency range when making movie of ERSP. The
%                 ERSP values are averaged over the selected frequency
%                 range. Not implemented
%   'limitbeg'  - [min max] limits at the beginning of the movie
%   'limitend'  - [min max] limits at the end of the movie
%   'freqslim'  - [freqs] array of landmark frequencies to set color
%                 limits.
%   'movieparams' - [low inc high] lower limit, higher limit and increment. If
%                increment is omited, movie is generate at every possible
%                increment (max freq resolution or max time resolution)
%
% Authors: Arnaud Delorme, CERCO, August, 2006
%
% See also: std_specplot(), std_erppplot(), std_erspplot()

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [STUDY, M] = std_movie(STUDY, ALLEEG, varargin);

if nargin < 2
    help std_specplot;
    return;
end;

[opt addparams ] = finputcheck( varargin, ...
                             { 'erspfreq'    'real'    [] [];
                               'movieparams' 'real'    [] [];
                               'channels'    'cell'    [] {};
                               'freqslim'    'real'    [] [];
                               'limitbeg'    'real'    [] [];
                               'limitend'    'real'    [] [];
                               'moviemode'   'string'  { 'ersptime' 'erp' 'spec' } 'spec' }, 'std_movie', 'ignore');

if isstr(opt), error(opt); end;
tmpchanlocs =  ALLEEG(1).chanlocs;
if isempty(opt.channels), opt.channels = { tmpchanlocs.labels }; end;
if ~strcmpi(opt.moviemode, 'spec'), error('Only spec has been implemented so far'); end;

% read data once
% --------------
STUDY = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', [10 11]); close;

% find first data channel with info
% ---------------------------------
for cind = 1:length(STUDY.changrp)
    if ~isempty(STUDY.changrp(cind).specdata), break; end;
end;

% generate movie
% --------------
if isempty(opt.movieparams),
    opt.movieparams(1) = STUDY.changrp(cind).specfreqs(1);
    opt.movieparams(2) = STUDY.changrp(cind).specfreqs(end);
end;
if length(opt.movieparams) == 2
    opt.movieparams(3) = opt.movieparams(2);
    opt.movieparams(2) = STUDY.changrp(cind).specfreqs(2)-STUDY.changrp(cind).specfreqs(1);
end;
if length(opt.movieparams) == 3
    opt.movieparams = [opt.movieparams(1):opt.movieparams(2):opt.movieparams(3)];
end;

% find limits
% -----------
if isempty(opt.limitbeg)
    [STUDY specdata] = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(1)); close;
    opt.limitbeg = [ min(specdata{1}) max(specdata{1}) ];
    [STUDY specdata] = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(end)); close;
    opt.limitend = [ min(specdata{1}) max(specdata{1}) ];
end;
lowlims  = linspace(opt.limitbeg(1), opt.limitend(1), length(opt.movieparams));
highlims = linspace(opt.limitbeg(2), opt.limitend(2), length(opt.movieparams));

% limits at specific frequencies
% ------------------------------
if ~isempty(opt.freqslim)
    if opt.freqslim(1)   ~= opt.movieparams(1)  , opt.freqslim = [ opt.movieparams(1) opt.freqslim ]; end;
    if opt.freqslim(end) ~= opt.movieparams(end), opt.freqslim = [ opt.freqslim opt.movieparams(end) ]; end;

    for ind = 1:length(opt.freqslim)
        [tmp indf(ind)] = min(abs(opt.freqslim(ind) - opt.movieparams));
        [STUDY specdata] = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(indf(ind))); close;
        minimum(ind) = min(specdata{1});
        maximum(ind) = max(specdata{1});
    end;
    indf(1)  = 0;
    lowlims  = [ ];
    highlims = [ ];
    for ind = 2:length(opt.freqslim)
        lowlims  = [lowlims linspace(minimum(ind-1), minimum(ind), indf(ind)-indf(ind-1)) ];
        highlims = [highlims linspace(maximum(ind-1), maximum(ind), indf(ind)-indf(ind-1)) ];
    end;
end;

% make movie
% ----------
for ind = 1:length(opt.movieparams)
    STUDY = std_specplot(STUDY, ALLEEG, 'channels', opt.channels, 'topofreq', opt.movieparams(ind), 'caxis', [lowlims(ind) highlims(ind)]);
    pos = get(gcf, 'position');
    set(gcf, 'position', [pos(1) pos(2) pos(3)*2 pos(4)*2]);
    M(ind) = getframe(gcf);
    close;
end;
figure; axis off; movie(M);

author: PatternRecognition

repo: OpenBMI


text:
filename='Throne_Tetrahedra_SYM';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';

nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target=3.5;
optimality_final =1e-3;
constr_final =1e-3;

BCscale_factor = 0.3;
rotation_per_it = 20;
HJiter0 = 1;
e2 = 1;
N_holes = [6 3 3];
R_holes = 0.4;
phase_holes = [0 0 0];

Vfrac_initial = 0.7;
optimality_initial = 5e-2;
constr_initial = 5e-2;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;

% maxiter = 1;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function D = dijkstra( G , S )
% --------------------------------------------------------------------
%      Mark Steyvers, Stanford University, 12/19/00
% --------------------------------------------------------------------
%
% DIJKSTRA  Find shortest paths in graphs
% 	D = dijkstra( G , S ) use the full or sparse matrix G in which
% 	an entry (i,j) represents the arc length between nodes i and j in a
%	graph. In a full matrix, the value INF represents the absence of an arc;
%	in a sparse matrix, no entry at (i,j) naturally represents no arc.
%
%	S is the one-dimensional matrix of source nodes for which the shortest
%	to ALL other nodes in the graphs will be calculated. The output matrices
%	D and P contain the shortest distances and predecessor indices respectively.
%	An infinite distance is represented by INF. The predecessor indices contain
%	the node indices of the node along the shortest path before the destination
%	is reached. These indices are useful to construct the shortest path with the
%	function pred2path (by Michael G. Kay).
%
%	This function was implemented in C++. The source code can be compiled to a
%	Matlab compatible mex file by the command "mex -O dijkstra.cpp" at the Matlab
%	prompt. In this package, we provide a compiled .dll version that is
%       compatible all Windows based machines.  If you are not working on a
%       Windows platform, delete the .dll version provided and recompile from
%       the .cpp source file.  If you do not have the Matlab compiler or a Windows
%       platform, delete the .dll version and dijkstra will then call the
%	Matlab function dijk.m (by Michael G. Kay).  Note that this Matlab
%	code is several orders of magnitude slower than the C based mex file.

N = size( G , 1 );
D = dijk( G , S , 1:N );

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function name = p12_name ( )

%*****************************************************************************80
%
%% P12_NAME returns the name of problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'BFN02';

  return
end


text:
function lik = lik_lgp(varargin)
%LIK_LGP  Create a logistic Gaussian process likelihood structure
%
%  Description
%    LIK = LIK_LGP creates a logistic Gaussian process likelihood structure
%
%    The likelihood is defined as follows:
%               __ n
%      p(y|f) = || i=1 exp(f_i) / Sum_{j=1}^n exp(f_j),
%
%      where f contains latent values.
%
%  Reference
%
%    Jaakko Riihimäki and Aki Vehtari (2014). Laplace approximation
%    for logistic Gaussian process density estimation and
%    regression. Bayesian analysis, in press.
%
%  See also
%    LGPDENS, GP_SET, LIK_*
%
% Copyright (c) 2011 Jaakko Riihimäki and Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'LIK_LGP';
  ip.addOptional('lik', [], @isstruct);
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'LGP';
    lik.nondiagW = true;
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'LGP')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_lgp_pak;
    lik.fh.unpak = @lik_lgp_unpak;
    lik.fh.ll = @lik_lgp_ll;
    lik.fh.llg = @lik_lgp_llg;
    lik.fh.llg2 = @lik_lgp_llg2;
    lik.fh.llg3 = @lik_lgp_llg3;
    lik.fh.tiltedMoments = @lik_lgp_tiltedMoments;
    lik.fh.predy = @lik_lgp_predy;
    lik.fh.invlink = @lik_lgp_invlink;
    lik.fh.recappend = @lik_lgp_recappend;
  end

end

function [w,s,h] = lik_lgp_pak(lik)
%LIK_LGP_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_LGP_PAK(LIK) takes a likelihood structure LIK
%    and returns an empty verctor W. If LGP likelihood had
%    parameters this would combine them into a single row vector
%    W (see e.g. lik_negbin). This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%  See also
%    LIK_LGP_UNPAK, GP_PAK

  w = []; s = {}; h=[];
end

function [lik, w] = lik_lgp_unpak(lik, w)
%LIK_LGP_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_LGP_UNPAK(W, LIK) Doesn't do anything.
%
%    If LGP likelihood had parameters this would extract them
%    parameters from the vector W to the LIK structure. This
%    is a mandatory subfunction used for example in energy
%    and gradient computations.
%
%
%  See also
%    LIK_LGP_PAK, GP_UNPAK

  lik=lik;
  w=w;

end

function logLik = lik_lgp_ll(lik, y, f, z)
%LIK_LGP_LL    Log likelihood
%
%  Description
%    E = LIK_LGP_LL(LIK, Y, F, Z) takes a likelihood data
%    structure LIK, incedence counts Y, expected counts Z, and
%    latent values F. Returns the log likelihood, log p(y|f,z).
%    This subfunction is needed when using Laplace approximation
%    or MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria (DIC, WAIC)
%    computations.
%
%  See also
%    LIK_LGP_LLG, LIK_LGP_LLG3, LIK_LGP_LLG2, GPLA_E

  n=sum(y);
  qj=exp(f);
  logLik = sum(f.*y)-n*log(sum(qj));
end

function deriv = lik_lgp_llg(lik, y, f, param, z)
%LIK_LGP_LLG    Gradient of the log likelihood
%
%  Description
%    G = LIK_LGP_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z
%    and latent values F. Returns the gradient of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    'param' or 'latent'. This subfunction is needed when using Laplace
%    approximation or MCMC for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LGP_LL, LIK_LGP_LLG2, LIK_LGP_LLG3, GPLA_E

  switch param
    case 'latent'
      n=sum(y);
      qj=exp(f);
      pj=qj./sum(qj);
      deriv=y-n*pj;
  end
end

function g2 = lik_lgp_llg2(lik, y, f, param, z)
%function g2 = lik_lgp_llg2(lik, y, f, param, z)
%LIK_LGP_LLG2  Second gradients of the log likelihood
%
%  Description
%    G2 = LIK_LGP_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z,
%    and latent values F. Returns the Hessian of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. G2 is a vector with diagonal elements of the
%    Hessian matrix (off diagonals are zero). This subfunction
%    is needed when using Laplace approximation or EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LGP_LL, LIK_LGP_LLG, LIK_LGP_LLG3, GPLA_E

  switch param
    case 'latent'
      qj=exp(f);

      % g2 is not the second gradient of the log likelihood but only a
      % vector to form the exact gradient term in gpla_nd_e, gpla_nd_g and
      % gpla_nd_pred functions
      g2=qj./sum(qj);
  end
end

function g3 = lik_lgp_llg3(lik, y, f, param, z)
%LIK_LGP_LLG3  Third gradients of the log likelihood
%
%  Description
%    G3 = LIK_LGP_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z
%    and latent values F and returns the third gradients of the
%    log likelihood with respect to PARAM. At the moment PARAM
%    can be only 'latent'. G3 is a vector with third gradients.
%    This subfunction is needed when using Laplace approximation
%    for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_LGP_LL, LIK_LGP_LLG, LIK_LGP_LLG2, GPLA_E, GPLA_G

  switch param
    case 'latent'
      qj=exp(f);

      % g3 is not the third gradient of the log likelihood but only a
      % vector to form the exact gradient term in gpla_nd_e, gpla_nd_g and
      % gpla_nd_pred functions
      g3=qj./sum(qj);

      %n=sum(y);
      %nf=size(f,1);
      %g3d=zeros(nf,nf);
      %for i1=1:nf
      %  g3dtmp=-g3*g3(i1);
      %  g3dtmp(i1)=g3dtmp(i1)+g3(i1);
      %  g3d(:,i1)=g3dtmp;
      %  %g3i1= n*(-diag(g3d(:,i1)) + bsxfun(@times,g3,g3d(:,i1)') + bsxfun(@times,g3d(:,i1),g3'));
      %end
  end
end

function [logM_0, m_1, sigm2hati1] = lik_lgp_tiltedMoments(lik, y, i1, sigm2_i, myy_i, z)
%LIK_LGP_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_LGP_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, incedence counts
%    Y, expected counts Z, index I and cavity variance S2 and
%    mean MYY. Returns the zeroth moment M_0, mean M_1 and
%    variance M_2 of the posterior marginal (see Rasmussen and
%    Williams (2006): Gaussian processes for Machine Learning,
%    page 55). This subfunction is needed when using EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    GPEP_E

  error('Not implemented')

end

function [lpy, Ey, Vary] = lik_lgp_predy(lik, Ef, Varf, yt, zt)
%LIK_LGP_PREDY    Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_LGP_PREDY(LIK, EF, VARF YT, ZT)
%    Returns also the predictive density of YT, that is
%        p(yt | y,zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires also the incedence counts YT, expected counts ZT.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_LGP_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This
%    subfunction is needed when computing posterior predictive
%    distributions for future observations.
%

%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  error('Not implemented')

end

function [df,minf,maxf] = init_lgp_norm(yy,myy_i,sigm2_i,myy)
%INIT_LGP_NORM
%
%  Description
%    Return function handle to a function evaluating LGP *
%    Gaussian which is used for evaluating (likelihood * cavity)
%    or (likelihood * posterior) Return also useful limits for
%    integration. This is private function for lik_lgp. This
%    subfunction is needed by sufunctions tiltedMoments, siteDeriv
%    and predy.
%
%  See also
%    LIK_LGP_TILTEDMOMENTS, LIK_LGP_PREDY

% Not applicable

end

function mu = lik_lgp_invlink(lik, f, z)
%LIK_LGP_INVLINK  Returns values of inverse link function
%
%  Description
%    P = LIK_LGP_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values MU of inverse link function.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_LGP_LL, LIK_LGP_PREDY

  mu = exp(f);
  mu = mu./sum(mu);

end

function reclik = lik_lgp_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = LIK_LGP_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK. This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    reclik.type = 'LGP';

    % Set the function handles
    reclik.fh.pak = @lik_lgp_pak;
    reclik.fh.unpak = @lik_lgp_unpak;
    reclik.fh.ll = @lik_lgp_ll;
    reclik.fh.llg = @lik_lgp_llg;
    reclik.fh.llg2 = @lik_lgp_llg2;
    reclik.fh.llg3 = @lik_lgp_llg3;
    reclik.fh.tiltedMoments = @lik_lgp_tiltedMoments;
    reclik.fh.predy = @lik_lgp_predy;
    reclik.fh.invlink = @lik_lgp_invlink;
    reclik.fh.recappend = @lik_lgp_recappend;
    return
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
function Ne = som_neighborhood(Ne1,n)

%SOM_NEIGHBORHOOD Calculate neighborhood matrix.
%
% Ne = som_neighborhood(Ne1,n)
%
%  Ne = som_neighborhood(Ne1);
%  Ne = som_neighborhood(som_unit_neighs(topol),2);
%
%  Input and output arguments ([]'s are optional):
%   Ne1       (matrix, size [munits m]) a sparse matrix indicating
%                      the units in 1-neighborhood for each map unit
%   [n]       (scalar) maximum neighborhood which is calculated, default=Inf
%
%   Ne        (matrix, size [munits munits]) neighborhood matrix,
%                      each row (and column) contains neighborhood
%                      values from the specific map unit to all other
%                      map units, or Inf if the value is unknown.
%
% For more help, try 'type som_neighborhood' or check out online documentation.
% See also SOM_UNIT_NEIGHS, SOM_UNIT_DISTS, SOM_UNIT_COORDS, SOM_CONNECTION.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_neighborhood
%
% PURPOSE
%
% Calculate to which neighborhood each map unit belongs to relative to
% each other map unit, given the units in 1-neighborhood of each unit.
%
% SYNTAX
%
%  Ne = som_neighborhood(Ne1);
%  Ne = som_neighborhood(Ne1,n);
%
% DESCRIPTION
%
% For each map unit, finds the minimum neighborhood to which it belongs
% to relative to each other map unit. Or, equivalently, for each map
% unit, finds which units form its k-neighborhood, where k goes from
% 0 to n.
%
% The neighborhood is calculated iteratively using the reflexivity of
% neighborhood.
%     let  N1i  be the 1-neighborhood set a unit i
% and let  N11i be the set of units in the 1-neighborhood of any unit j in N1i
%     then N2i  (the 2-neighborhood set of unit i) is N11i \ N1i
%
% Consider, for example, the case of a 5x5 map. The neighborhood in case of
% 'rect' and 'hexa' lattices (and 'sheet' shape) for the unit at the
% center of the map are depicted below:
%
%   'rect' lattice           'hexa' lattice
%   --------------           --------------
%   4  3  2  3  4            3  2  2  2  3
%   3  2  1  2  3             2  1  1  2  3
%   2  1  0  1  2            2  1  0  1  2
%   3  2  1  2  3             2  1  1  2  3
%   4  3  2  3  4            3  2  2  2  3
%
% Because the iterative procedure is rather slow, the neighborhoods
% are calculated upto given maximal value. The uncalculated values
% in the returned matrix are Inf:s.
%
% REQUIRED INPUT ARGUMENTS
%
%  Ne1   (matrix) Each row contains 1, if the corresponding unit is adjacent
%                 for that map unit, 0 otherwise. This can be calculated
%                 using SOM_UNIT_NEIGHS. The matrix can be sparse.
%                 Size munits x munits.
%
% OPTIONAL INPUT ARGUMENTS
%
%  n     (scalar) Maximal neighborhood value which is calculated,
%                 Inf by default (all neighborhoods).
%
% OUTPUT ARGUMENTS
%
%  Ne    (matrix) neighborhood values for each map unit, size is
%                 [munits, munits]. The matrix contains the minimum
%                 neighborhood of unit i, to which unit j belongs,
%                 or Inf, if the neighborhood was bigger than n.
%
% EXAMPLES
%
%  Ne = som_neighborhood(Ne1,1);    % upto 1-neighborhood
%  Ne = som_neighborhood(Ne1,Inf);  % all neighborhoods
%  Ne = som_neighborhood(som_unit_neighs(topol),4);
%
% SEE ALSO
%
%  som_unit_neighs   Calculate units in 1-neighborhood for each map unit.
%  som_unit_coords   Calculate grid coordinates.
%  som_unit_dists    Calculate interunit distances.
%  som_connection    Connection matrix.

% Copyright (c) 1999-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 1.0beta juuso 141097
% Version 2.0beta juuso 101199

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Check arguments

error(nargchk(1, 2, nargin));

if nargin<2, n=Inf; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Action

% initialize
if issparse(Ne1), Ne = full(Ne1); else Ne = Ne1; end
clear Ne1
[munits dummy] = size(Ne);
Ne(find(Ne==0)) = NaN;
for i=1:munits, Ne(i,i)=0; end

% Calculate neighborhood distance for each unit using reflexsivity
% of neighborhood:
%   let  N1i be the 1-neighborhood set a unit i
%   then N2i is the union of all map units, belonging to the
%        1-neighborhood of any unit j in N1i, not already in N1i
k=1;
if n>1,
  fprintf(1,'Calculating neighborhood: 1 ');
  N1 = Ne;
  N1(find(N1~=1)) = 0;
end
while k<n && any(isnan(Ne(:))),
  k=k+1;
  fprintf(1,'%d ',k);
  for i=1:munits,
    candidates = isnan(Ne(i,:));              % units not in any neighborhood yet
    if any(candidates),
      prevneigh = find(Ne(i,:)==k-1);         % neighborhood (k-1)
      N1_of_prevneigh = any(N1(prevneigh,:)); % union of their N1:s
      Nn = find(N1_of_prevneigh & candidates);
      if length(Nn), Ne(i,Nn) = k; Ne(Nn,i) = k; end
    end
  end
end
if n>1, fprintf(1,'\n'); end

% finally replace all uncalculated distance values with Inf
Ne(find(isnan(Ne))) = Inf;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% faster version?

l = size(Ne1,1); Ne1([0:l-1]*(l+1)+1) = 1; Ne = full(Ne1); M0 = Ne1; k = 2;
while any(Ne(:)==0), M1=(M0*Ne1>0); Ne(find(M1-M0))=k; M0=M1; k=k+1; end
Ne([0:l-1]*(l+1)+1) = 0;

author: ilarinieminen

repo: SOM-Toolbox


text:
function Centers = NSGA2ESelection(tr_x, tr_y, models, str, Problem, Ke)
% NSGA-II based selection

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    nn=0;
    while (nn==100||nn==0)
        population = rand(Problem.N,Problem.D); %N*D
        population = population.*repmat(Problem.upper,Problem.N,1)+(1-population).*repmat(Problem.lower,Problem.N,1);%N*D
        population(:,Problem.D+1:Problem.D+Problem.M)=Estimate(population(:,1:Problem.D), tr_x, tr_y, models, str, Problem.M);%The value of objective by model

        [~,FrontNo,CrowdDis] = EnvironmentalSelection(population,Problem.N,Problem.D,Problem.M);%N*(D+M)

        %% Optimization of NSGA-II
        for i=1:floor(Problem.maxFE/Problem.N)
            MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
            parent     = population(MatingPool,:);
            offspring  = OperatorGA(Problem,parent(:,1:Problem.D),{1,20,1,20});
            offspring(:,Problem.D+1:Problem.D+Problem.M)=Estimate(offspring(:,1:Problem.D), tr_x, tr_y, models, str, Problem.M);%The value of objective by model
            [population,FrontNo,CrowdDis] = EnvironmentalSelection([population;offspring],Problem.N,Problem.D,Problem.M);
        end
        [Centers, nn]=Kmean([population,FrontNo',CrowdDis'], tr_x, Problem.N, Problem.D, Ke);
    end
end

function objv_LCB=Estimate(x, tr_x, tr_y, models, str,m)
    alpha=2;
    TestSamNum=size(x,1);n=length(models);
    sum=zeros(TestSamNum, size(tr_y,2));
    for i=1:n
        clear y;
        stri=str(i,:);
        M=models(i).M;p=models(i).p;
        if strcmp(stri(1), 'FE')%PCA
            x_pca=x*M;
            if strcmp(stri(2), 'RBF1')%RBF
                y= sim(p,x_pca');y=y';
            elseif strcmp(stri(2), 'SVM')
                y=SVMtest(x_pca, p, TestSamNum,m);
            elseif strcmp(stri(2), 'RBF2')%RBF
                y=RBF2test(x_pca, p, TestSamNum);
            end

        elseif strcmp(stri(1), 'NONE')%NONE
            if strcmp(stri(2), 'RBF1')%RBF
                y= sim(p,x');y=y';
            elseif strcmp(stri(2), 'SVM')
                y=SVMtest(x, p, TestSamNum, m);
            elseif strcmp(stri(2), 'RBF2')%RBF
                y=RBF2test(x, p, TestSamNum);
            end

        elseif strcmp(stri(1), 'FS')%CSO
            x_cso=x(:,M);
            if strcmp(stri(2), 'RBF1')%RBF
                y= sim(p,x_cso');y=y';
            elseif strcmp(stri(2), 'SVM')
                y=SVMtest(x_cso, p, TestSamNum,m);
            elseif strcmp(stri(2), 'RBF2')%RBF
                y=RBF2test(x_cso, p, TestSamNum);
            end
        end
        result(i).y=y;
        sum=sum+y;
    end
    me=sum/n;
    sum=zeros(TestSamNum, size(tr_y,2));
    for i=1:n
        y=result(i).y;
        sum=sum+(y-me).^2;
    end
    s2=sum/(n-1);
    %fmin=min(tr_y);
    %fmin_matrix=repmat(fmin,TestSamNum,1);
    objv_LCB=me-alpha*sqrt(s2);
end

function y=SVMtest(x, p, N, M)
    ps1=p{1};ps2=p{2};
    xtest0=mapminmax('apply',x',ps1);xtest0=xtest0';
%     py = zeros(N,M);
    for j=1:M
        for k=1:N
            py(k,j) = svmpredict(0,xtest0(k,:),p{j+2}, '-q');
        end
    end
    y=mapminmax('reverse',py',ps2);y=y';
end

function y=RBF2test(x, p, N)
    Centers=p.Centers;Spreads=p.Spreads;
    W2=p.W2;B2=p.B2;
    TestDistance = dist(Centers,x');
    TestSpreadsMat = repmat(Spreads,1,N);
    TestHiddenUnitOut = radbas(TestDistance./TestSpreadsMat);
    y= (W2*TestHiddenUnitOut+repmat(B2,1,N))';
end

function [population,FrontNo,CrowdDis] = EnvironmentalSelection(population,N,D,M)
% The environmental selection of NSGA-II

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(population(:,D+1:D+M),[],N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(population(:,D+1:D+M),FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    population = population(Next,:);
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
end

function [Centers, nn]=Kmean(chromosome, tr_x, N, D, Ke)
    % Delete the multiple solutions between chromosome and the train data
    Q=[];P=[];
    for i=1:N
        PandQ=[tr_x;Q];
        matrix=dist(chromosome(i,1:D),PandQ');%1*size(PandQ,1) array
        if isempty(find(matrix<=10^-06))
            Q=[Q;chromosome(i,1:D)];
            P=[P;chromosome(i,:)];
        end
    end
    if size(Q,1)<Ke
        Centers=[];nn=100;
    else
        index=randperm(size(Q,1));
        Centers=Q(index(1:Ke),1:D);n=1;si=size(Q,1);
        % Cluster by the distance of the decision space
        while n<100

            NumberInClusters = zeros(Ke,1); % Number of samples in each class ,default is 0
            IndexInClusters = zeros(Ke,N); % Index of samples in each class
            % Classify all samples by the least distance principle
            for i = 1:si
                AllDistance = dist(Centers,Q(i,1:D)');% Calculate the distance between the i-th solution and each clustering center
                [~,Pos] = min(AllDistance);   % Minimum distance,training input is the index of clustering center
                NumberInClusters(Pos) = NumberInClusters(Pos) + 1;
                IndexInClusters(Pos,NumberInClusters(Pos)) = i;% Stores the training indexes belonging to this class in turn
            end
            % Store the old clustering centers
            OldCenters = Centers;
            % Recalculate the clustering centers
            for i = 1:Ke
                Index = IndexInClusters(i,1:NumberInClusters(i));% Extract the training input index belonging to this class
                Centers(i,:) = mean(Q(Index,1:D),1);    %Take the average of each class as the new clustering center
            end
            % Judge whether the old and new clustering centers are consistent
            EqualNum = sum(sum(Centers==OldCenters));% Centers and OldCenters are subtracted from each other to sum up all corresponding bits
            if EqualNum ==D*Ke % The old and new clustering centers are consistent
                break,
            end
            n=n+1;
        end
        nn=n;fprintf('k-means clustering %d\n',n);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO