text:
function com = pop_topohistplot( EEG, typeplot, components, topotitle, rowcols, varargin);

com = '';
if nargin < 1
    help pop_topoplot;
    return;
end;
if nargin < 2
    typeplot = 1;
end;
if typeplot == 0 & isempty(EEG.icasphere)
    disp('Error: no ICA data for this set, first run ICA or Load AMICA components'); return;
end;
if isempty(EEG.chanlocs)
    disp('Error: cannot plot topography without channel location file'); return;
end;
if ~isfield(EEG.etc,'amica') || isempty(EEG.etc.amica)
    error('No AMICA solution found. You should first load AMICA components');return;
end

if nargin < 3
    % which set to save
    % -----------------
    if typeplot
        txtwhat2plot1 = 'Plotting ERP scalp maps at these latencies';
        txtwhat2plot2 = sprintf('(range: %d to %d ms, NaN -> empty):', ...
            round(EEG.xmin*1000), round(EEG.xmax*1000));
        editwhat2plot = [''];
    else
        txtwhat2plot1 = 'Component numbers';
        txtwhat2plot2 = '(negate index to invert component polarity; NaN -> empty subplot; Ex: -1 NaN 3)';
        editwhat2plot = ['1:' int2str(size(EEG.icaweights,1))];
    end;
    allmodels = 1:EEG.etc.amica.num_models;

    if ~isfield(EEG.etc.amica,'modnames')
        EEG.etc.amica.modnames{length(allmodels)} = '';
    end

    str = [num2str(allmodels(1)) fastif(~isempty(EEG.etc.amica.modnames{1}) ,[' - ' EEG.etc.amica.modnames{1}],'')];
    for i = 2:EEG.etc.amica.num_models
        str = [str '|' num2str(allmodels(i)) fastif(~isempty(EEG.etc.amica.modnames{i}) ,[' - ' EEG.etc.amica.modnames{i}],'')];
    end
    editwhat2select = str;

    if EEG.nbchan > 64,
        elecdef = ['''electrodes'', ''off'''];
    else,
        elecdef = ['''electrodes'', ''on'''];
    end;
    cb_enable = ['set(findobj(''parent'', gcbf, ''tag'', ''block'')    , ''enable'', ''on'');'];
    cb_disable = ['set(findobj(''parent'', gcbf, ''tag'', ''block'')    , ''enable'', ''off'');'];
    cb_block = ['if get(gcbo,''value''),' cb_enable ...
        'else,' cb_disable ...
        'end;'];
    cb_chosen = ['m = get(gcbo,''value''); set(findobj(''parent'',gcbf,''tag'',''title''),''string'',[EEG.setname '', Model '' num2str(m)])'];

    %EEG.etc.amica
    uilist = { { 'style'   'text'     'string'    txtwhat2plot1 } ...
        { 'style'   'edit'     'string'    editwhat2plot } ...
        { 'style'   'text'     'string'    txtwhat2plot2 } ...
        { } ...
        {'style' 'text' 'string' 'Select AMICA model'} ...
        {'style' 'popupmenu' 'string' editwhat2select 'value' 1 'callback' cb_chosen} ...
        { 'style'   'text'     'string'    'Plot title' } ...
        { 'style'   'edit'     'tag' 'title' 'string'    fastif(~isempty(EEG.setname), [EEG.setname ', Model 1'], '') } ...
        { 'style'   'text'     'string'    'Plot geometry (rows,col.); [] -> near square' } ...
        { 'style'   'edit'     'string'    '[]' } ...
        { 'style'   'text'     'string'    'Plot associated dipole(s) (if present)' } ...
        { 'style'   'checkbox' 'string'    '' } { } ...
        { } ...
        { 'style'   'text'     'string'    [ '-> Additional topoplot()' fastif(typeplot,'',' (and dipole)') ...
        ' options (see Help)' ] } ...
        { 'style'   'edit'     'string'    elecdef } ...
        {'style' 'checkbox' 'string' 'Show IC histograms' 'value' 1 'callback' cb_block} ...
        {'style' 'checkbox' 'tag' 'block' 'string' 'Use block analysis when determining histogram range' 'enable' 'on'}};
    uigeom = { [1.5 1] [1] [1] [1 1] [1.5 1] [1.5 1] [1.55 0.2 0.8] [1] [1] [1] [1] [1]};
    if typeplot
        uilist(9:11) = [];
        uigeom(6) = [];
    end;
    guititle = fastif( typeplot, 'Plot ERP scalp maps in 2-D -- pop_topoplot()', ...
        'Plot AMICA component scalp maps in 2-D -- pop_topohistplot()');

    result = inputgui( uigeom, uilist, 'pophelp(''pop_topoplot'')', guititle, [], 'normal');
    if isempty(result)
        return;
    end

    % reading first param
    % -------------------
    components   	     = eval( [ '[' result{1} ']' ] );
    if length(components) > EEG.nbchan
        tmpbut = questdlg2(...
            ['This involves drawing ' int2str(length(components)) ' plots. Continue ?'], ...
            '', 'Cancel', 'Yes', 'Yes');
        if strcmp(tmpbut, 'Cancel'), return; end;
    end;
    if isempty(components), error('Nothing to plot; enter parameter in first edit box'); end;

    % reading other params
    % --------------------
    h = result{2}; % chosen AMICA model
    topotitle   = result{3};
    rowcols     = eval( [ '[' result{4} ']' ] );
    options = [];
    if typeplot
        plotdip = 0;
        try, options      = eval( [ '{ ' result{5} ' }' ]);
        catch, error('Invalid scalp map options'); end;
    else
        plotdip     = result{5};
        try, options      = eval( [ '{ ' result{6} ' }' ]);
        catch, error('Invalid scalp map options'); end;
    end;
    if length(components) == 1,
        figure('paperpositionmode', 'auto'); curfig=gcf;
        try, icadefs;
            set(curfig, 'color', BACKCOLOR);
        catch, end;
    end;
    showhist = result{7};
    use_block = result{8};
    options = {options{:} 'showhist' showhist 'use_block' use_block 'model' h};
else
    if ~isempty(varargin) & isnumeric(varargin{1})
        plotdip = varargin{1};
        varargin = varargin(2:end);
    else
        plotdip = 0;
    end;
    options = varargin;
end;

% additional options
% ------------------
outoptions = { options{:} }; % for command
options    = { options{:} 'masksurf' 'on' };

% find maplimits
% --------------
maplimits = [];
for i=1:2:length(options)
    if isstr(options{i})
        if strcmpi(options{i}, 'maplimits')
            maplimits = options{i+1};
            options(i:i+1) = [];
            break;
        end;
    end;
end;

%--------------Added for IC histogram plotting--------------------------
for i=1:2:length(options)
    if isstr(options{i})
        if strcmpi(options{i}, 'showhist')
            showhist = options{i+1};
            options(i:i+1) = [];
            break;
        end;
    end;
end;

for i=1:2:length(options)
    if isstr(options{i})
        if strcmpi(options{i}, 'use_block')
            use_block = options{i+1};
            options(i:i+1) = [];
            break;
        end;
    end;
end;

for i=1:2:length(options)
    if isstr(options{i})
        if strcmpi(options{i}, 'model')
            h = options{i+1};
            options(i:i+1) = [];
            break;
        end;
    end;
end;
if ~exist('h','var')
    h = 1;
end
% ----------------------------------------------------------------------

nbgraph = size(components(:),1);
if ~exist('topotitle')
    topotitle = '';
end;
if ~exist('rowcols') | isempty(rowcols) | rowcols == 0
    rowcols(2) = ceil(sqrt(nbgraph));
    rowcols(1) = ceil(nbgraph/rowcols(2));
end;

SIZEBOX = 150;

fprintf('Plotting...\n');
if isempty(EEG.chanlocs)
    fprintf('Error: set has no channel location file\n');
    return;
end;

% Check if pop_topoplot input 'colorbar' was called, and don't send it to topoplot
loc = strmatch('colorbar', options(1:2:end), 'exact');
loc = loc*2-1;
if ~isempty(loc)
    colorbar_switch = strcmp('on',options{ loc+1 });
    options(loc:loc+1) = [];
else
    colorbar_switch = 1;
end

% determine the scale for plot of different times (same scales)
% -------------------------------------------------------------
if typeplot
    SIGTMP = reshape(EEG.data, EEG.nbchan, EEG.pnts, EEG.trials);
    pos = round( (components/1000-EEG.xmin)/(EEG.xmax-EEG.xmin) * (EEG.pnts-1))+1;
    nanpos = find(isnan(pos));
    pos(nanpos) = 1;
    SIGTMPAVG = mean(SIGTMP(:,pos,:),3);
    SIGTMPAVG(:, nanpos) = NaN;
    if isempty(maplimits)
        maxlim = max(SIGTMPAVG(:));
        minlim = min(SIGTMPAVG(:));
        maplimits = [ -max(maxlim, -minlim) max(maxlim, -minlim)];
    end;
else
    if isempty(maplimits)
        maplimits = 'absmax';
    end;
end;

if plotdip
    if strcmpi(EEG.dipfit.coordformat, 'CTF')
        disp('Cannot plot dipole on scalp map for CTF MEG data');
    end;
end;

%Calculate the winv matrix for the model chosen
A = pinv(EEG.etc.amica.W(:,:,h)*EEG.etc.amica.S);
%------------------------------------------------

% plot the graphs
% ---------------
counter = 1;
countobj = 1;
allobj = zeros(1,1000);
curfig = get(0, 'currentfigure');
if isfield(EEG, 'chaninfo'), options = { options{:} 'chaninfo' EEG.chaninfo }; end

if showhist

    alpha = EEG.etc.amica.alpha;
    sbeta = EEG.etc.amica.sbeta;
    rho = EEG.etc.amica.rho;
    mu = EEG.etc.amica.mu;
    nbgraph = length(components);
    mn = EEG.etc.amica.data_mean;
    x = EEG.data;

    models = EEG.etc.amica;
    if EEG.trials>1
        x = reshape(x,EEG.nbchan,EEG.trials*EEG.pnts);
        models.LLt = reshape(models.LLt,EEG.etc.amica.num_models,size(models.LLt,2)*size(models.LLt,3));
    end
    if max(size(models.LLt)) > 1
        indx = find(models.LLt(1,:)~=0);
        lim = find(indx>size(x,2));
        if ~isempty(lim)
            indx = indx(1:lim(1)-1);
        end
        x = x(:,indx);
        models.LLt = models.LLt(:,indx);
    end
    [nx,N] = size(x);
    n = size(models.c,1);
    B = 1000;
    num_blocks = floor(N/B);
    nbins = min(100,sqrt(N));
    [mval,mind] = max(models.LLt);
    n = size(models.c,1);
    s = zeros(n,N);
    disp('Sphering the data...')
    if use_block

        % sphere the data

        for k = 1:num_blocks
            xstrt = (k-1)*B + 1;
            if k < num_blocks
                xstp = k*B;
            else
                xstp = N;
            end
            for i = 1:nx
                x(i,xstrt:xstp) = x(i,xstrt:xstp) - mn(i);
                %x(i,xstrt:xstp) = s(1,xstrt:xstp) - mn(i);
            end
            x(1:n,xstrt:xstp) = models.S(1:n,:) * x(:,xstrt:xstp);
        end

        for k = 1:num_blocks
            xstrt = (k-1)*B + 1;
            if k < num_blocks
                xstp = k*B;
            else
                xstp = N;
            end
            s(components,xstrt:xstp) = models.W(components,:,h) * x(1:n,xstrt:xstp);
            for i = components
                s(i,xstrt:xstp) = s(i,xstrt:xstp) - models.c(i,h);
            end
        end

    else
        for i = 1:nx
            x(i,:) = x(i,:) - mn(i);
            %x(i,xstrt:xstp) = s(1,xstrt:xstp) - mn(i);
        end

        x(1:n,:) = models.S(1:n,:) * x(:,:);

        s(components,:) = models.W(components,:,h) * x(1:n,:);
        for i = components
            s(i,:) = s(i,:) - models.c(i,h);
        end

    end

end

for index = 1:size(components(:),1)
    if nbgraph > 1
        if mod(index, rowcols(1)*rowcols(2)) == 1

            if index> 1, figure(curfig); a = textsc(0.5, 0.05, topotitle); set(a, 'fontweight', 'bold'); end;
            curfig = figure('paperpositionmode', 'auto');
            pos = get(curfig,'Position');

            posx = max(0, pos(1)+(pos(3)-SIZEBOX*rowcols(2))/2);

            posy = pos(2)+pos(4)-SIZEBOX*rowcols(1);
            if showhist

                set(curfig,'Position', [posx posy  2*SIZEBOX*rowcols(2)  SIZEBOX*rowcols(1)]);
            else
                set(curfig,'Position', [posx posy  SIZEBOX*rowcols(2)  SIZEBOX*rowcols(1)]);
            end
            try, icadefs; set(curfig, 'color', BACKCOLOR); catch, end;
        end;

        if showhist
            curax = subplot( rowcols(1), 2*rowcols(2), 2*mod(index-1, rowcols(1)*rowcols(2))+1,'parent',curfig);
        else
            curax = subplot( rowcols(1), rowcols(2), mod(index-1, rowcols(1)*rowcols(2))+1,'parent',curfig);
        end
        set(curax, 'visible', 'off')
        if index == 1
            constantPosition = get(curax,'position');
        end

    else
        if showhist
            curax = subplot(1,2,1);
            set(curax,'visible','off')
        end

    end

    % add dipole location if present
    % ------------------------------
    dipoleplotted = 0;
    if plotdip && typeplot == 0
        if isfield(EEG, 'dipfit') & isfield(EEG.dipfit, 'model')
            if length(EEG.dipfit.model) >= index & ~strcmpi(EEG.dipfit.coordformat, 'CTF')
                %curpos = EEG.dipfit.model(components(index)).posxyz/EEG.dipfit.vol.r(end);
                curpos = EEG.dipfit.model(components(index)).posxyz;
                curmom = EEG.dipfit.model(components(index)).momxyz;
                try,
                    select = EEG.dipfit.model(components(index)).select;
                catch ,
                    select = 0;
                end;
                if ~isempty(curpos)
                    if strcmpi(EEG.dipfit.coordformat, 'MNI') % from MNI to sperical coordinates
                        transform = pinv( sph2spm );
                        tmpres = transform * [ curpos(1,:) 1 ]'; curpos(1,:) = tmpres(1:3);
                        tmpres = transform * [ curmom(1,:) 1 ]'; curmom(1,:) = tmpres(1:3);
                        try, tmpres = transform * [ curpos(2,:) 1 ]'; curpos(2,:) = tmpres(1:3); catch, end;
                        try, tmpres = transform * [ curmom(2,:) 1 ]'; curmom(2,:) = tmpres(1:3); catch, end;
                    end;
                    curpos = curpos / 85;
                    if size(curpos,1) > 1 && length(select) == 2

                        dipole_index = find(strcmpi('dipole',options),1);
                        if  ~isempty(dipole_index) % if 'dipoles' is already defined in options{:}
                            options{dipole_index+1} = [ curpos(:,1:2) curmom(:,1:3) ];
                        else
                            options = { options{:} 'dipole' [ curpos(:,1:2) curmom(:,1:3) ] };
                        end
                        dipoleplotted = 1;
                    else
                        if any(curpos(1,:) ~= 0)
                            dipole_index = find(strcmpi('dipole',options),1);
                            if  ~isempty(dipole_index) % if 'dipoles' is already defined in options{:}
                                options{dipole_index+1} = [ curpos(1,1:2) curmom(1,1:3) ];
                            else
                                options = { options{:} 'dipole' [ curpos(1,1:2) curmom(1,1:3) ] };
                            end
                            dipoleplotted = 1;
                        end
                    end

                end
                if nbgraph ~= 1
                    dipscale_index = find(strcmpi('dipscale',options),1);
                    if ~isempty(dipscale_index) % if 'dipscale' is already defined in options{:}
                        options{dipscale_index+1} = 0.6;
                    else
                        options = {  options{:} 'dipscale' 0.6 };
                    end
                end
                %options = { options{:} 'dipsphere' max(EEG.dipfit.vol.r) };
            end
        end
    end
    %curax = subplot( rowcols(1), 2*rowcols(2), 2*mod(index-1, rowcols(1)*rowcols(2))+1);
    % plot scalp map
    % --------------
    if index == 1
        addopt = { 'verbose', 'on' };
    else
        addopt = { 'verbose', 'off' };
    end;
    %fprintf('Printing to figure %d.\n',curfig);
    options = {  'maplimits' maplimits options{:} addopt{:} };
    if ~isnan(components(index))
        if typeplot
            if nbgraph > 1, axes(curax); end;
            tmpobj = topoplot( SIGTMPAVG(:,index), EEG.chanlocs, options{:});
            if nbgraph == 1,
                figure(curfig); if nbgraph > 1, axes(curax); end;
                title( [ 'Latency ' int2str(components(index)) ' ms from ' topotitle]);
            else
                figure(curfig); if nbgraph > 1, axes(curax); end;
                title([int2str(components(index)) ' ms'] );
            end;
        else
            if components(index) < 0
                figure(curfig);  if nbgraph > 1, axes(curax); end;
                tmpobj = topoplot( -A(:, -components(index)), EEG.chanlocs, options{:} );
            else
                figure(curfig);  if nbgraph > 1, axes(curax); end;
                tmpobj = topoplot( A(:, components(index)), EEG.chanlocs, options{:} );
            end;
            if nbgraph == 1, texttitle = [ 'IC ' int2str(components(index)) ' from ' topotitle];
            else             texttitle = ['' int2str(components(index))];
            end;
            if dipoleplotted, texttitle = [ texttitle ' (' num2str(EEG.dipfit.model(components(index)).rv*100,2) '%)']; end;
            figure(curfig);  if nbgraph > 1, axes(curax); end; title(texttitle);
        end;
        allobj(countobj:countobj+length(tmpobj)-1) = tmpobj;
        countobj = countobj+length(tmpobj);
        pos = get(curax,'position');

        set(curax,'position',[pos(1) pos(2) constantPosition(3) constantPosition(4)]);
        drawnow;
        axis square;
    else
        axis off
    end;

    % plot IC histogram
    % --------------------------

    if showhist

        if EEG.etc.amica.num_models > 1
            inds = find(mind == h);
            if ~isempty(inds)
                smn = mean(s(components(index),inds));
                [~,b] = hist(s(components(index),inds)-smn,nbins);

            else
                b = -10:0.01:10;

            end
        else
            [~,b] = hist(s(components(index),:),nbins);
        end

        if nbgraph > 1
            subplot(rowcols(1),2*rowcols(2),2*index,'parent',curfig);
        else
            subplot(1,2,2,'parent',curfig);
        end
        hold on
        sum = 0;
        for j = 1:size(alpha,1)
            mixtpdf = alpha(j,components(index),h)*sbeta(j,components(index),h)*(1/(2*gamma(1+1/rho(j,components(index),h))))*exp(-abs(sbeta(j,components(index),h)*(b-mu(j,components(index),h))).^rho(j,components(index),h));
            if ~isnan(mixtpdf)
                sum = mixtpdf + sum;
                plot(b,mixtpdf,'g');
            end
        end

        plot(b,sum,'b');

        drawnow;
    end
    % ---------------------------------------------
end

% Draw colorbar
if colorbar_switch
    if nbgraph == 1
        if ~isstr(maplimits)
            ColorbarHandle = cbar(0,0,[maplimits(1) maplimits(2)]);
        else
            ColorbarHandle = cbar(0,0,get(gca, 'clim'));
        end;
        pos = get(ColorbarHandle,'position');  % move left & shrink to match head size
        if showhist %if IC histogram is also plotted
            set(ColorbarHandle,'position',[pos(1)+.02 pos(2)+0.13 pos(3)*0.7 pos(4)-0.26]);
        else
            set(ColorbarHandle,'position',[pos(1)-.05 pos(2)+0.13 pos(3)*0.7 pos(4)-0.26]);
        end
    elseif ~isstr(maplimits)
        cbar('vert',0,[maplimits(1) maplimits(2)]);
    else cbar('vert',0,get(gca, 'clim'));
    end
    if ~typeplot    % Draw '+' and '-' instead of numbers for colorbar tick labels
        tmp = get(gca, 'ytick');
        set(gca, 'ytickmode', 'manual', 'yticklabelmode', 'manual', 'ytick', [tmp(1) tmp(end)], 'yticklabel', { '-' '+' });
    end
end

if nbgraph> 1,
    figure(curfig); a = textsc(0.5, 0.05, topotitle);
    set(a, 'fontweight', 'bold');
end;
if nbgraph== 1,
    com = 'figure;';
end;
set(allobj(1:countobj-1), 'visible', 'on');

figure(curfig);
axcopy(curfig, 'set(gcf, ''''units'''', ''''pixels''''); postmp = get(gcf, ''''position''''); set(gcf, ''''position'''', [postmp(1) postmp(2) 560 420]); clear postmp;');

com = [com sprintf('pop_topohistplot(%s,%d, %s);', ...
    inputname(1), typeplot, vararg2str({components topotitle rowcols plotdip outoptions{:} }))];
disp('Done.');
return;

author: PatternRecognition

repo: OpenBMI


text:
function [g,info] = wilwin(g,M,L,callfun);
%WILWIN  Compute a Wilson/WMDCT window from text or cell array
%   Usage: [g,info] = wilwin(g,M,L);
%
%   `[g,info]=wilwin(g,M,L)` computes a window that fits well with the
%   specified number of channels *M* and transform length *L*. The window itself
%   is specified by a text description or a cell array containing additional
%   parameters.
%
%   The window can be specified directly as a vector of numerical
%   values. In this case, `wilwin` only checks assumptions about transform
%   sizes etc.
%
%   `[g,info]=wilwin(g,M)` does the same, but the window must be a FIR
%   window, as the transform length is unspecified.
%
%   The window can be specified as one of the following text strings:
%
%     'gauss'      Gaussian window with optimal concentration
%
%     'dualgauss'  Riesz dual of Gaussian window with optimal concentration.
%
%     'tight'      Window generating an orthonormal basis
%
%   In these cases, a long window is generated with a length of *L*.
%
%   It is also possible to specify one of the window names from |firwin|. In
%   such a case, `wilwin` generates the specified FIR window with a length
%   of *M*.
%
%   The window can also be specified as cell array. The possibilities are:
%
%     `{'gauss',...}`
%       Additional parameters are passed to PGAUSS
%
%     `{'dual',...}`
%       Dual window of whatever follows. See the examples below.
%
%     `{'tight',...}`
%       Orthonormal window of whatever follows.
%
%   It is also possible to specify one of the window names from |firwin| as
%   the first field in the cell array. In this case, the remaining
%   entries of the cell array are passed directly to |firwin|.
%
%   Some examples::
%
%     g=wilwin('gauss',M,L);
%
%   This computes a Gaussian window of length L fitted for a system with
%   *M* channels. ::
%
%     g=wilwin({'gauss',1},M,L);
%
%   This computes a Gaussian window with equal time and frequency support
%   irrespective of *M*. ::
%
%     gd=wilwin('gaussdual',M,L);
%
%   This computes the dual of a Gaussian window fitted for a system with *M*
%   channels. ::
%
%     gd=wilwin({'tight','gauss'},M,L);
%
%   This computes the orthonormal window of the Gaussian window fitted for
%   the system. ::
%
%     g=wilwin({'dual',{'hann',20}},M,L);
%
%   This computes the dual of a Hann window of length 20.
%
%   The structure info provides some information about the computed
%   window:
%
%     `info.gauss`
%        True if the window is a Gaussian.
%
%     `info.tfr`
%        Time/frequency support ratio of the window. Set whenever it makes sense.
%
%     `info.wasrow`
%        Input was a row window
%
%     `info.isfir`
%        Input is an FIR window
%
%     `info.isdual`
%        Output is the dual window of the auxiliary window.
%
%     `info.istight`
%        Output is known to be a tight window.
%
%     `info.auxinfo`
%        Info about auxiliary window.
%
%     `info.gl`
%        Length of window.
%
%   See also: pgauss, firwin, gabwin

% Assert correct input.
complainif_argnonotinrange(nargin,2,4,mfilename);

if nargin==2
  L=[];
end;

% Basic discovery: Some windows depend on L, and some windows help define
% L, so the calculation of L is window dependant.

% Default values.
info.gauss=0;
info.wasrow=0;
info.isfir=0;
info.istight=0;
info.isdual=0;

% Manually get the list of window names
definput=arg_firwin(struct);
firwinnames =  definput.flags.wintype;

% Create window if string was given as input.
if ischar(g)
  winname=lower(g);
  switch(winname)
   case {'pgauss','gauss'}
    complain_L(L,callfun);
    g=comp_pgauss(L,2*M*M/L,0,0);
    info.gauss=1;
    info.tfr=2*M*M/L;
   case {'psech','sech'}
    complain_L(L,callfun);
    g=psech(L,2*M*M/L);
    info.tfr=a*M/L;
   case {'dualgauss','gaussdual'}
    complain_L(L,callfun);
    g=comp_pgauss(L,2*M*M/L,0,0);
    g=wildual(g,M);
    info.isdual=1;
    info.tfr=a*M/L;
   case {'tight'}
    complain_L(L,callfun);
    g=wilorth(M,L);
    info.tfr=2*M*M/L;
    info.istight=1;
   case firwinnames
    [g,firinfo]=firwin(winname,2*M,'2');
    info.isfir=1;
    if firinfo.issqpu
      info.istight=1;
    end;
   otherwise
    error('%s: Unknown window type: %s',callfun,winname);
  end;
end;

if iscell(g)
  if isempty(g) || ~ischar(g{1})
    error('First element of window cell array must be a character string.');
  end;

  winname=lower(g{1});

  switch(winname)
   case {'pgauss','gauss'}
    complain_L(L,callfun);
    [g,info.tfr]=pgauss(L,g{2:end});
    info.gauss=1;
   case {'psech','sech'}
    complain_L(L,callfun);
    [g,info.tfr]=psech(L,g{2:end});
   case {'dual'}
    gorig = g{2};
    [g,info.auxinfo] = wilwin(gorig,M,L,callfun);
    g = wildual(g,M,L);

    % gorig can be string or cell array
    if info.auxinfo.isfir && test_isfir(gorig,M)
        info.isfir = 1;
    end

    info.isdual=1;
   case {'tight'}
    gorig = g{2};
    [g,info.auxinfo] = wilwin(g{2},M,L,callfun);
    g = wilorth(g,M,L);
    % gorig can be string or cell array
    if info.auxinfo.isfir && test_isfir(gorig,M)
        info.isfir = 1;
    end
    info.istight=1;
   case firwinnames
    g=firwin(winname,g{2},'energy',g{3:end});
    info.isfir=1;
   otherwise
    error('Unsupported window type.');
  end;
end;

if isnumeric(g)
  if size(g,2)>1
    if size(g,1)==1
      % g was a row vector.
      g=g(:);
      info.wasrow=1;
    end;
  end;
end;

if rem(length(g),2*M)~=0
  % Zero-extend the window to a multiple of 2*M
  g=fir2long(g,ceil(length(g)/(2*M))*2*M);
end;

% Information to be determined post creation.
info.wasreal = isreal(g);
info.gl      = length(g);

if (~isempty(L) && (info.gl<L))
  info.isfir=1;
end;

function complain_L(L,callfun)

  if isempty(L)
    error(['%s: You must specify a length L if a window is represented as a ' ...
           'text string or cell array.'],callfun);
  end;

  function isfir=test_isfir(gorig,M)
    % Original window is FIR, dual window is FIR if length of the original
    % window is <= M. This is true if the length was not explicitly
    % defined (gorig{2}).
    if iscell(gorig) && numel(gorig)>1 && isnumeric(gorig{2}) && gorig{2}<=2*M...
         || ischar(gorig)
        isfir = 1;
    else
       isfir = 0;
    end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ L ] = SCA_phi_step_para( U,v,N, theta_old )
    %%
    Ut=abs(U);
    vt=abs(v);
    tmp1=0;
    for n0=1:N
        for i0=1:N
            if i0~=n0
                tmp1=tmp1+Ut(n0,i0)^2;
            end
        end
    end
    tmp2=0;
    for n0=1:N
        tmp2n=vt(n0);
        for i0=1:N
            if i0~=n0
                tmp2n=tmp2n+Ut(n0,i0);
            end
        end
        tmp2=tmp2+tmp2n^2;
    end
    tmp=tmp1+tmp2;
    L=2*sqrt(tmp);
end

author: guohuayan


text:
function a = moler4 ( )

%*****************************************************************************80
%
%% MOLER4 returns the MOLER4 matrix.
%
%  Example:
%
%    0  2  0 -1
%    1  0  0  0
%    0  1  0  0
%    0  0  1  0
%
%  Properties:
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is the companion matrix of the polynomial X^4-2X^2+1=0.
%
%    A has eigenvalues -1, -1, +1, +1.
%
%    A can cause problems to a standard QR algorithm, which
%    can fail to converge.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%

%
%  Note that matrix entries are listed by row.
%
  a = [ ...
    0.0,  2.0,  0.0, -1.0; ...
    1.0,  0.0,  0.0,  0.0; ...
    0.0,  1.0,  0.0,  0.0; ...
    0.0,  0.0,  1.0,  0.0 ];

  return
end%

author: johannesgerer

repo: jburkardt-m


text:
% You can create a network and convert it to a format that can be viewed in
% Samiam by running this script.

pedigree = struct('parents', [0,0;1,3;0,0;1,3;2,6;0,0;2,6;4,9;0,0]);
pedigree.names = {'Ira','James','Robin','Eva','Jason','Rene','Benjamin','Sandra','Aaron'};
alleleFreqs = [0.1; 0.9];
alleleList = {'F', 'f'};
alphaList = [0.8; 0.6; 0.1];
phenotypeList = {'CysticFibrosis', 'NoCysticFibrosis'};
positions = [520, 600, 520, 500; 650, 400, 650, 300; 390, 600, 390, 500; 260, 400, 260, 300; 780, 200, 780, 100; 1040, 400, 1040, 300; 910, 200, 910, 100; 130, 200, 130, 100; 0, 400, 0, 300];

% This will construct a Bayesian network and convert it into a file that
% can be viewed in SamIam.

factorList = constructGeneticNetwork(pedigree, alleleFreqs, alphaList);
sendToSamiam(pedigree, factorList, alleleList, phenotypeList, positions, 'cysticFibrosisBayesNet');

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function fcv = Calculate_fcv(Population)
% calculate normalized  constraints violation(CV) measuring feasibility

%------------------------------- 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".
%--------------------------------------------------------------------------

    CV_Original = Population.cons;
    CV_Original(CV_Original<=0) = 0;
    CV = CV_Original./max(CV_Original);
    CV(:,isnan(CV(1,:))) = 0;
    fcv = sum(max(0,CV),2)./size(CV_Original,2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function X = mat(x,n)
% Y = MAT(x,n)   or   Y = MAT(x)     (the 2nd argument is optional)
%   Given a vector of length n^2, this produces the n x n matrix
%   Y such that x = vec(Y).  In other words, x contains the columns of the
%   matrix Y, stacked below each other.
%
% See also vec.

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

if nargin < 2
    n = floor(sqrt(length(x)));
    if (n*n) ~= length(x)
        error('Argument X has to be a square matrix')
    end
end
X = reshape(x,n,n);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [p] = eeg_toolbox_defaults(command,p);

% eeg_toolbox_defaults - Create, read, write eeg_toolbox defaults
%
% Useage: [p] = eeg_toolbox_defaults(command,[p])
%
% command  =  'create'
%             'read'
%             'save'   | 'write'
%             'saveas' | 'write_other'
% p        =  structure returned by create|read command
%             structure to write for write command
%
% All read and write commands are to a matlab .mat file.  The
% 'saveas' or 'write_other' command will write out a data specific archive
% that can be accessed from the recent files list.  It will be located in
% the folder where the eeg data is opened.
%
% Examples:
%[p] = eeg_toolbox_defaults; % create new defaults
%[p] = eeg_toolbox_defaults('read'); % read saved defaults
% eeg_toolbox_defaults('save',p);  % write current p as default
% eeg_toolbox_defaults('saveas',p); % write current p data
%
% The save or write command will write to the eeg_toolbox
% installation folder.  The read command will first
% try to find the paramater file in the eeg_toolbox
% installation and otherwise recreates the defaults.
%
% Notes:    Handles parameter structure for the
%           eeg_toolbox routines.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:53 $

% Licence:  GNU GPL, no express or implied warranties
% Created:  01/2002 Darren.Weber_at_radiology.ucsf.edu
% Modified: 02/2002 Darren.Weber_at_radiology.ucsf.edu
%           - adapted read/write output from .txt format
%             to .mat format so it will contain data structures.
%             Hence, it is no longer editable text.
%           08/2002 Darren.Weber_at_radiology.ucsf.edu
%                   added MRI defaults
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eegversion = '$Revision: 1.1 $';
fprintf('\nEEG_TOOLBOX_DEFAULTS [v %s]\n',eegversion(11:15));

if version('-release') < 10,
    msg = printf('matlab release < version 10, eeg_toolbox GUIs may fail.\n');
    warning(msg);
end

if ~exist('command','var'), command = 'create'; end

% try to locate the installation path
eegPath = eeg_toolbox_path;

% try to locate the mri_toolbox installation path
mriPath = fileparts(which('avw_view'));
if isempty(mriPath),
    msg = sprintf('Cannot find mri_toolbox on the matlab path.\nPlease install and use the addpath command.\n\n');
    warning(msg);
    mriPath = eegPath;
else
    mriPath = strcat(mriPath,filesep);
end

% remove any saved preferences
if ispref('eeg_toolbox'),
    rmpref('eeg_toolbox');
end

switch command
case 'create',

    fprintf('...creating eeg_toolbox defaults\n\n');

   [p] = [];

    % default location for parameter file
    p.file = strcat(eegPath,'eeg_example_data',filesep,'eeg_toolbox_defaults.mat');

    % parameters for 'elec_load' routine, see that for more information
    p.elec.path     = strcat(eegPath,'eeg_example_data',filesep);
    p.elec.file     = 'elec_124_cart.txt';
    p.elec.n        = 129;         % how many electrodes to load
    p.elec.type     = 'Cartesian'; % see elec_load for avail. options
    p.elec.shape    = 'sphere';    % original = '', 'sphere', or 'ellipse'
    p.elec.data     = [];
    p.elec.plot     = 0; % plot electrode positions
    p.elec.plotSurf = 1; % plot electrode surface (from convhull)

    % parameters for 'mesh_open' & 'mesh_plot', see these for more information
    p.mesh.path         = strcat(eegPath,'eeg_example_data',filesep);
    p.mesh.file         = 'BrainStorm_subjecttess.mat';
    p.mesh.type         = 'BrainStorm'; % see mesh_open.m for options
    p.mesh.data         = []; % init for mesh_open.m
    p.mesh.plot.open    = 0; % boolean arg for mesh_open.m
    p.mesh.plot.overlay = 0; % Overlay arg in mesh_plot.m
    p.mesh.plotSurf     = 0; % plot scalp mesh surface (eeg_contours_engine)
    p.mesh.reload       = 1; % boolean, 1 = reload mesh data
    p.mesh.samplePoint  = [];
    p.mesh.sampleTime   = []; % Time of current sample point (msec)
    p.mesh.sampleMsec   = []; % Sample rate (msec)
    p.mesh.sampleHz     = []; % Sample rate (Hz)

    % parameters for 'eeg_open' routine, see that for more information
    p.volt.path                 = strcat(eegPath,'eeg_example_data',filesep);
    p.volt.file                 = 'eeg_124ch.dat';
    p.volt.type                 = 'ascii'; % see eeg_open.m for valid options

    p.volt.data                 = []; % potential (uV)
    p.volt.var                  = []; % variance for ERP (uV^2)
    p.volt.timeArray            = []; % The timing of each sample point (msec)
    p.volt.points               = []; % Number of sample points in epoch
    p.volt.sampleHz             = 400; % Sample rate (Hz)
    p.volt.sampleMsec           = 2.5; % Sample rate (msec)
    p.volt.sampleTime           = 400; % Time of current sample point (msec)
    p.volt.samplePoint          = 250; % single or start time point (row of p.volt.data)
    p.volt.channels             = []; % Number of channels/electrodes
    p.volt.epochStart           = -200; % Time of epoch start (msec)
    p.volt.epochEnd             = 1500; % Time of epoch end (msec)
    p.volt.sweeps               = []; % Number of accepted sweeps
    p.volt.peaks                = []; % Data array of ERP peak values
    p.volt.interpZero           = 1;  % Interpolate time zero ERP potential,
    % mainly for ascii versions of NeuroScan data

    p.cnt.path                  = strcat(eegPath,'eeg_example_data',filesep);
    p.cnt.file                  = 'scan41_short.cnt';
    p.cnt.type                  = 'scan4x_cnt';

    p.eeg.path                  = strcat(eegPath,'eeg_example_data',filesep);
    p.eeg.file                  = 'scan41_short.eeg';
    p.eeg.type                  = 'scan4x_eeg';

    % parameters for 'erf_open' routine, see that for more information
    p.erf.path                 = strcat(eegPath,'eeg_example_data',filesep);
    p.erf.file                 = 'eeg_124ch.dat';
    p.erf.type                 = 'ascii'; % see eeg_open.m for valid options

    p.erf.data                 = [];   % potential (Tesla)
    p.erf.var                  = [];   % variance for ERF (T^2)
    p.erf.timeArray            = [];   % The timing of each sample point (msec)
    p.erf.points               = [];   % Number of sample points in epoch
    p.erf.sampleHz             = 400;  % Sample rate (Hz)
    p.erf.sampleMsec           = 2.5;  % Sample rate (msec)
    p.erf.sampleTime           = 400;  % Time of current sample point (msec)
    p.erf.samplePoint          = 250;  % single or start time point (row of p.erf.data)
    p.erf.channels             = [];   % Number of channels/electrodes
    p.erf.epochStart           = -200; % Time of epoch start (msec)
    p.erf.epochEnd             = 1500; % Time of epoch end (msec)
    p.erf.sweeps               = [];   % Number of accepted sweeps
    p.erf.peaks                = [];   % Data array of ERP peak values
    p.erf.interpZero           = 0;    % Interpolate time zero ERF potential

    % parameters for 'eeg_contours_engine', see that for more information
    p.rangeMethod                 = 'minmaxabs' ; % minmax? = 'abs','one','all','giv'
    p.minimumIntensity            =  -10    ;
    p.maximumIntensity            =   10    ;
    p.interpMethod                = 'cubic' ;   % matlab griddata interpolation method

    p.timeMethod                  =    1    ;   % 1=single time point, 2=range of time points
    p.endTime                     =    2    ;   % finish time point (row of p.volt.data)

    p.contour.stepMethod          =  0; % 0 is step size, 1 is number of steps
    p.contour.Nsteps              = 10;
    p.contour.stepSize            =  1; % 1 uV or similar units
    p.contour.raw2D               =  0; % plot 2D contours
    p.contour.plot2D              =  0; % plot projected 2D contours
    p.contour.plot3D              =  0; % plot 3D surface contours

    % Now defunct (02/2003), as topo maps now use triangulations, rather than meshgrid
    % keeping for now, just in case they prove useful at some later date
    p.grid.method                  =    1    ;   % Use grid size method
    p.grid.size                    =  100    ;
    p.grid.sizeMin                 =   10    ;
    p.grid.sizeMax                 =  210    ;
    p.grid.res                     =    0.1  ;   % 1 mm
    p.grid.resMin                  =    0.01 ;
    p.grid.resMax                  =    2.01 ;

    % Miscellaneous
    p.clickTimePoint              = 0;   % select a time point interactively
    p.saveGraphics                = 'No Save Plots'; % Do not save graphic files by default
    p.topoView                    = '2D'; % Default topography view
    p                             = eeg_colormap(p); % Return defaults for 'eeg_colormap'
    p.hold                        = 0;   % default for GUI hold checkboxes

case 'read'

    % Look for default file
    [path,name,ext] = fileparts(strcat(eegPath,'eeg_toolbox_defaults.mat'));
    p.file = fullfile(path,[name ext]);
    file = exist(p.file);
    if file == 2,
        fprintf('...reading eeg_toolbox defaults from:\n%s\n\n',p.file);
        load(p.file);
    else
        fprintf('...cannot locate eeg_toolbox defaults - recreating defaults.\n\n');
       [p] = eeg_toolbox_defaults;
    end

    % verify that path to default files exists;
    % if not, create defaults again (I hope this will avoid
    % new installation problems)
    voltPathExist = exist(p.volt.path);
    if ~voltPathExist,
        fprintf('...voltage path does not exist - reinitializing defaults.\n');
       [p] = eeg_toolbox_defaults('create');
        eeg_toolbox_defaults('write',p);
    end

case {'write','save'}

    if ~exist('p','var'),
        msg = sprintf('p argument is essential: eeg_toolbox_defaults(''write'',p)');
        error(msg);
    end

    [path,name,ext] = fileparts(strcat(eegPath,'eeg_toolbox_defaults.mat'));
    p.file = fullfile(path,[name ext]);

    save(p.file,'p');

    fprintf('...saved eeg_toolbox defaults to:\n%s\n\n',p.file);

case {'write_other','saveas'}

    if ~exist('p','var'),
        msg = sprintf('p argument is essential: eeg_toolbox_defaults(''write_other'',p)');
        error(msg);
    end

    [path,name,ext] = fileparts(strcat(p.volt.path, filesep, p.volt.file));
    ext = '.mat';
    p.file = fullfile(path,[name ext]);

    [newfile,newpath] = uiputfile(p.file,'Save EEG_TOOLBOX Parameters to File:');
    if (newfile == 0) & (newpath == 0),
        fprintf('...aborting save\n\n');
    else
        [path,name,ext] = fileparts(strcat(newpath, filesep, newfile));
        ext = '.mat';
        p.file = fullfile(path,[name ext]);
        save(p.file,'p');
        fprintf('...saved eeg_toolbox parameters to:\n%s\n\n',p.file);
        recentfiles = eeg_toolbox_recent(p.file);
    end

otherwise
    error('...invalid command\n\n');
end

return


text:
function [Es,Vs]=subTet(E,V,splitMethod)

% function [Es,Vs]=subTet(E,V,splitMethod)
% ------------------------------------------------------------------------
% Sub devides the input tetrahedral mesh into a denser mesh by splitting
% the elements. Two methods are available. If splitMethod==1 then the faces
% are split similar to what is expected for [Fs,Vs]=subtri(F,V,1,1). I.e.
% for side faces mid-edge nodes are introduced so edges are split. This
% yields 4 new corner tetrahedral elements which could inherit shape
% quality from the original tetrahedral element. However the entire
% tetrahedron cannot me split this way and a central octahedron remains.
% This octahedron is split using its Delaunay tesselation representation
% and yields 4 new tetrahedrons which differ in shape quality from the
% original tetrahedron. If splitMethod==2 the tetrahedron is split into 4
% tetrahedra by introducing a single new central node for each element and
% by connecting the side faces
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 2014/03/04
% ------------------------------------------------------------------------

switch splitMethod
    case 1 %Split faces + central delaunay of octahedron
        edgeMat=[E(:,[1 2]); E(:,[2 3]);  E(:,[3 1]); E(:,[1 4]); E(:,[3 4]); E(:,[2 4])]; %Edges matrix
        E_sort=sort(edgeMat,2); %Sorted edges matrix
        [~,ind1,~]=unique(E_sort,'rows');
        edgeMat=edgeMat(ind1,:);

        numPoints = size(V,1);
        numEdges = size(edgeMat,1);

        % Get indices of the three edges associated with each face
        A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
        A = max(A,A'); %Copy symmetric

        %Indices for A matrix
        indA_12=E(:,1)+(E(:,2)-1)*numPoints;
        indA_23=E(:,2)+(E(:,3)-1)*numPoints;
        indA_31=E(:,3)+(E(:,1)-1)*numPoints;
        indA_14=E(:,1)+(E(:,4)-1)*numPoints;
        indA_24=E(:,2)+(E(:,4)-1)*numPoints;
        indA_34=E(:,3)+(E(:,4)-1)*numPoints;

        %Get indices for vertex array
        indV_12=full(A(indA_12));
        indV_23=full(A(indA_23));
        indV_31=full(A(indA_31));
        indV_14=full(A(indA_14));
        indV_24=full(A(indA_24));
        indV_34=full(A(indA_34));

        %Create element array
        Es=[E(:,1)  indV_12 indV_31 indV_14;... %Corner tet 1
            indV_12 indV_23 indV_24 E(:,2);... %Corner tet 2
            E(:,3)  indV_31 indV_23 indV_34;... %Corner tet 3
            indV_14 indV_24 indV_34 E(:,4);... %Corner tet 4
            indV_34 indV_31 indV_23 indV_24;... %Octahedral tet 1
            indV_31 indV_12 indV_23 indV_24;... %Octahedral tet 2
            indV_31 indV_12 indV_24 indV_14;... %Octahedral tet 3
            indV_34 indV_31 indV_24 indV_14;... %Octahedral tet 4
            ];

        %Create vertex array
        Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points
        Vs = [V; Vn]; %Join point sets

    case 2 %Add central point and connect all faces to central point to yield elements
        %Get faces
        [F,faceInd]=element2patch(E,(1:size(E,1))');

        %Create element array
        Es=[F faceInd+size(V,1)];

        %Create vertex array
        X=V(:,1); Y=V(:,2); Z=V(:,3);
        if size(E,1)==1 %Indexing behaviour (e.g. X(E)) differs for a single element
            Vn=[mean(X(E),1) mean(Y(E),1) mean(Z(E),1)]; %new mid-element points
        else
            Vn=[mean(X(E),2) mean(Y(E),2) mean(Z(E),2)]; %new mid-element points
        end
        Vs = [V; Vn]; %Join point sets
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [K] = ku0v0(x, xp, hyp, ubar, vbar, ubarp, vbarp, dt, i)

logsigmau = hyp(1);
logthetau = hyp(2);
logsigmav = hyp(3);
logthetav = hyp(4);

a1 = hyp(5);
a2 = hyp(6);

n_x = size(x,1);
n_xp = size(xp,1);

x = repmat(x,1,n_xp);
xp = repmat(xp',n_x,1);

ubar = repmat(ubar,1,n_xp);
vbar = repmat(vbar,1,n_xp);
ubarp = repmat(ubarp',n_x,1);
vbarp = repmat(vbarp',n_x,1);

switch i

case 0

K=dt.*(a2.*exp(1).^(logsigmav+(-1/2).*exp(1).^((-1).*logthetav).*(x+(-1).* ...
  xp).^2).*(ubar.^2+vbar.^2)+(-1).*a2.*exp(1).^(logsigmau+(-1/2).*exp(1) ...
  .^((-1).*logthetau).*(x+(-1).*xp).^2).*(ubarp.^2+vbarp.^2)+a1.*exp(1).^( ...
  logsigmau+(-2).*logthetau+(-1/2).*exp(1).^((-1).*logthetau).*(x+(-1).* ...
  xp).^2).*(exp(1).^logthetau+(-1).*(x+(-1).*xp).^2)+(-1).*a1.*exp(1).^( ...
  logsigmav+(-2).*logthetav+(-1/2).*exp(1).^((-1).*logthetav).*(x+(-1).* ...
  xp).^2).*(exp(1).^logthetav+(-1).*(x+(-1).*xp).^2));

case 1 % logsigmau

K=(-1).*dt.*exp(1).^(logsigmau+(-2).*logthetau+(-1/2).*exp(1).^((-1).* ...
  logthetau).*(x+(-1).*xp).^2).*(a2.*exp(1).^(2.*logthetau).*(ubarp.^2+ ...
  vbarp.^2)+a1.*((-1).*exp(1).^logthetau+(x+(-1).*xp).^2));

case 2 % logthetau

K=(-1/2).*dt.*exp(1).^(logsigmau+(-3).*logthetau+(-1/2).*exp(1).^((-1).* ...
  logthetau).*(x+(-1).*xp).^2).*(a1.*(2.*exp(1).^(2.*logthetau)+(-5).*exp( ...
  1).^logthetau.*(x+(-1).*xp).^2+(x+(-1).*xp).^4)+a2.*exp(1).^(2.* ...
  logthetau).*(ubarp.^2+vbarp.^2).*(x+(-1).*xp).^2);

case 3 % logsigmav

K=dt.*exp(1).^(logsigmav+(-2).*logthetav+(-1/2).*exp(1).^((-1).*logthetav) ...
  .*(x+(-1).*xp).^2).*(a2.*exp(1).^(2.*logthetav).*(ubar.^2+vbar.^2)+a1.*( ...
  (-1).*exp(1).^logthetav+(x+(-1).*xp).^2));

case 4 % logthetav

K=(1/2).*dt.*exp(1).^(logsigmav+(-3).*logthetav+(-1/2).*exp(1).^((-1).* ...
  logthetav).*(x+(-1).*xp).^2).*(a1.*(2.*exp(1).^(2.*logthetav)+(-5).*exp( ...
  1).^logthetav.*(x+(-1).*xp).^2+(x+(-1).*xp).^4)+a2.*exp(1).^(2.* ...
  logthetav).*(ubar.^2+vbar.^2).*(x+(-1).*xp).^2);

case 5 % a1

K=dt.*(exp(1).^(logsigmau+(-2).*logthetau+(-1/2).*exp(1).^((-1).* ...
  logthetau).*(x+(-1).*xp).^2).*(exp(1).^logthetau+(-1).*(x+(-1).*xp).^2)+ ...
  exp(1).^(logsigmav+(-2).*logthetav+(-1/2).*exp(1).^((-1).*logthetav).*( ...
  x+(-1).*xp).^2).*((-1).*exp(1).^logthetav+(x+(-1).*xp).^2));

case 6 % a2

K=dt.*(exp(1).^(logsigmav+(-1/2).*exp(1).^((-1).*logthetav).*(x+(-1).*xp) ...
  .^2).*(ubar.^2+vbar.^2)+(-1).*exp(1).^(logsigmau+(-1/2).*exp(1).^((-1).* ...
  logthetau).*(x+(-1).*xp).^2).*(ubarp.^2+vbarp.^2));

otherwise

        K = zeros(n_x, n_xp);
end

if K == 0

    K = zeros(n_x, n_xp);

end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
%{
MIT License

Copyright (c) 2018 ZjxRS

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
%}
% This is lifted from https://github.com/ZjxRS/ImageProcess/tree/master/Matlab/ENVI_ReadAndWrite

function info = envihdrread(hdrfile)
% ENVIHDRREAD Reads header of ENVI image.
%   INFO = ENVIHDRREAD('HDR_FILE') reads the ASCII ENVI-generated image
%   header file and returns all the information in a structure of
%   parameters.
%
%   Example:
%   >> info = envihdrread('my_envi_image.hdr')
%   info =
%          description: [1x101 char]
%              samples: 658
%                lines: 749
%                bands: 3
%        header_offset: 0
%            file_type: 'ENVI Standard'
%            data_type: 4
%           interleave: 'bsq'
%          sensor_type: 'Unknown'
%           byte_order: 0
%             map_info: [1x1 struct]
%      projection_info: [1x102 char]
%     wavelength_units: 'Unknown'
%           pixel_size: [1x1 struct]
%           band_names: [1x154 char]
%
%   NOTE: This function is used by ENVIREAD to import data.

% Ian M. Howat, Applied Physics Lab, University of Washington
% ihowat@apl.washington.edu
% Version 1: 19-Jul-2007 00:50:57
% Modified by Felix Totir
%

fid = fopen(hdrfile);
while true
    line = fgetl(fid);
    if line == -1
        break
    else
        eqsn = findstr(line,'=');
        if ~isempty(eqsn)
            param = strtrim(line(1:eqsn-1));
            param(findstr(param,' ')) = '_';
            value = strtrim(line(eqsn+1:end));
            if isempty(str2num(value))
                if ~isempty(findstr(value,'{')) && isempty(findstr(value,'}'))
                    while isempty(findstr(value,'}'))
                        line = fgetl(fid);

                        %for polSARpro
                        if line == -1
                            break;
                        end
                        value = [value,strtrim(line)];
                    end
                end
                eval(['info.',param,' = ''',value,''';'])
            else
                eval(['info.',param,' = ',value,';'])
            end
        end
    end
end
fclose(fid);

if isfield(info,'map_info')
    line = info.map_info;
    line(line == '{' | line == '}') = [];

    %originally: line = strtrim(split(line,','));
    %replaced by
    line=textscan(line,'%s', 'Delimiter',','); %behavior is not quite the same if "line" ends in ','
    line=line{:};
    line=strtrim(line);
    %

    info.map_info = [];
    info.map_info.projection = line{1};
    info.map_info.image_coords = [str2num(line{2}),str2num(line{3})];
    info.map_info.mapx = str2num(line{4});
    info.map_info.mapy = str2num(line{5});
    info.map_info.dx  = str2num(line{6});
    info.map_info.dy  = str2num(line{7});
    if length(line) == 9
        info.map_info.datum  = line{8};
        info.map_info.units  = line{9}(7:end);
    elseif length(line) == 11
        info.map_info.zone  = str2num(line{8});
        info.map_info.hemi  = line{9};
        info.map_info.datum  = line{10};
        info.map_info.units  = line{11}(7:end);
    end

    %part below comes form the original enviread
    %% Make geo-location vectors
    %     xi = info.map_info.image_coords(1);
    %     yi = info.map_info.image_coords(2);
    %     xm = info.map_info.mapx;
    %     ym = info.map_info.mapy;
    %     %adjust points to corner (1.5,1.5)
    %     if yi > 1.5
    %         ym =  ym + ((yi*info.map_info.dy)-info.map_info.dy);
    %     end
    %     if xi > 1.5
    %         xm = xm - ((xi*info.map_info.dy)-info.map_info.dx);
    %     end
    %
    %     info.x= xm + ((0:info.samples-1).*info.map_info.dx);
    %     info.y = ym - ((0:info.lines-1).*info.map_info.dy);
end

if isfield(info,'coordinate_system_string')
    line = info.coordinate_system_string;
    line(line == '{' | line == '}' | line == '"' | line == '[' ...
        | line == ']') = [];

    %originally: line = strtrim(split(line,','));
    %replaced by
    line=textscan(line,'%s', 'Delimiter',','); %behavior is not quite the same if "line" ends in ','
    line=line{:};
    line=strtrim(line);
    %

    info.coordinate_system = [];

    if contains(line{1},'UTM')
        info.coordinate_system.PROJCS = extractAfter(line{1},6);

        info.coordinate_system.GEOGCS = extractAfter(line{2},6);
        info.coordinate_system.DATUM = extractAfter(line{3},5);
        info.coordinate_system.SPHEROID = [];
        info.coordinate_system.SPHEROID.datum = extractAfter(line{4},8);
        info.coordinate_system.SPHEROID.a  = line{5};
        info.coordinate_system.SPHEROID.e  = line{6};
        info.coordinate_system.PRIMEM = extractAfter(line{7},6);
        info.coordinate_system.PMlongToGreenwich = line{8};
        info.coordinate_system.UNITGeo = extractAfter(line{9},4);
        info.coordinate_system.UNITGeo_value = line{10};

        info.coordinate_system.PROJECTION = extractAfter(line{11},10);
        info.coordinate_system.PARAMETER1 = extractAfter(line{12},9);
        info.coordinate_system.PARAMETER1_value = line{13};
        info.coordinate_system.PARAMETER2 = extractAfter(line{14},9);
        info.coordinate_system.PARAMETER2_value = line{15};
        info.coordinate_system.PARAMETER3 = extractAfter(line{16},9);
        info.coordinate_system.PARAMETER3_value = line{17};
        info.coordinate_system.PARAMETER4 = extractAfter(line{18},9);
        info.coordinate_system.PARAMETER4_value = line{19};
        info.coordinate_system.PARAMETER5 = extractAfter(line{20},9);
        info.coordinate_system.PARAMETER5_value = line{21};
        info.coordinate_system.UNITPro = extractAfter(line{22},4);
        info.coordinate_system.UNITPro_value = line{23};

    elseif strcmp(info.map_info.projection, 'Geographic Lat/Lon')
        info.coordinate_system.GEOGCS = extractAfter(line{1},6);
        info.coordinate_system.DATUM = extractAfter(line{2},5);
        info.coordinate_system.SPHEROID = [];
        info.coordinate_system.SPHEROID.datum = extractAfter(line{3},8);
        info.coordinate_system.SPHEROID.a  = line{4};
        info.coordinate_system.SPHEROID.e  = line{5};
        info.coordinate_system.PRIMEM = extractAfter(line{6},6);
        info.coordinate_system.PMlongToGreenwich = line{7};
        info.coordinate_system.UNITGeo = extractAfter(line{8},4);
        info.coordinate_system.UNITGeo_value = line{9};
    end

end

if isfield(info,'pixel_size')
    line = info.pixel_size;
    line(line == '{' | line == '}') = [];

    %originally: line = strtrim(split(line,','));
    %replaced by:
    line=textscan(line,'%s','Delimiter',','); %behavior is not quite the same if "line" ends in ','
    line=line{:};
    line=strtrim(line);

    info.pixel_size = [];
    info.pixel_size.x = str2num(line{1});
    info.pixel_size.y = str2num(line{2});
    info.pixel_size.units = line{3}(7:end);
end

% function split is only used when replacements above do not work
% function A = split(s,d)
%This function by Gerald Dalley (dalleyg@mit.edu), 2004
% A = {};
% while (~isempty(s))
%     [t,s] = strtok(s,d);
%     A = {A{:}, t};
% end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function [V,W] = art1s(p,rho,flag,V,W)
%ART1S    ART1 simulation function.
%
%         [V,W] = ART1S(P,rho,flag)
%           P   - F1xQ matrix of input vectors.
%           rho - the vigilance parameter, 0<= rho <=1.
%           flag  - (Optional) Printing flag. Any value of flag
%                   enables printing of events.
%         Returns:
%           V - the new top-down (T-D) weight matrix F1xF2.
%           W - the new bottom-up (B-U) weight matrix F2xF1.
%               F2 is the number of nodes in layer F2 (max # of categories)
%	  Example:  	P = letno;
%			art1s(P,0.7);
%
%		See also LETNO

%	  Author: Val Ninov, e-mail: valninov@total.net
%		  Grad. Student, Dept. of Electrical Engineering
%	          Concordia University, Montreal, Canada
%		  (c) April, 1997
%	  References:
%	   [1]	Carpenter, G. A. and S. Grossberg, "ART2: self-organization
%		of stable category recognition codes for analog input patterns."
%		Applied Optics, vol. 26, no. 23, Dec. 1987, pp. 4919-4930.
%	   [2]	J. Freeman and D. Skapura, Neural Networks: Algorithms,
%		Applications, and Programming Techniques. Addison Wesley

if nargin<2 | nargin>5  error('Wrong number of input arguments.');  end

% NETWORK PARAMETERS
 [R,Q] = size(p);
 F1 = R;
 F2 = Q;
 L = 2;
 categ = zeros(F2,F2); % category table
 count = ones(1,F2);   % counter for patterns in one category

% INITIALIZE WEIGHTS
if nargin < 4
 W = ones(F2,F1)*(L/(L-1+F1));
 V = ones(F1,F2);
 %fprintf('INITIAL TOP-DOWN MATRIX :');  V
 %fprintf('INITIAL BOTTOM-UP MATRIX :'); W
end

% INITIALIZE RETURN VARIABLES
a1 = zeros(F1,Q);   % output of F1
i = zeros(1,Q);     % winner index
win = 1;	    % First time node 1 in F2 is the winner
nActive = 1;        % The number of active neurons in F2

% PRESENT EACH INPUT VECTOR
for q=1:Q
    Reset = 0;
    B2 = zeros(F2,1);
    resonance = 0;
  while ~resonance    % REPEAT UNTIL: LAYERS F1 & F2 RESONATE
          % Initially a1 = p;
          % Calculate the winning node in F2 (among the active nodes)
	  A2 = compet(W(1:nActive,:)*p(:,q)+B2(1:nActive,1));
          i(q) = find(A2 == 1);
          win = i(q);
      % RECALCULATE a1 WITH FEEDBACK FROM A2
      a1(:,q) = (p(:,q) & V(:,win));

      % RESET if the new a1 is too different from p
	S= sum(a1(:,q));
	X= sum(p(:,q));
      Reset = (S/X) < rho;

       % IF RESET: TAKE WINNING NEURON IN F2 OUT OF COMPETITION
       if Reset
          B2(win) = -100;
             if nargin == 3
              fprintf('              RESET: Pattern %0.f resets F2 neuron %0.f.\n',q,win);
             end

           % IF ALL NEURONS IN A2 OUT OF COMPETITION ADD NEURON TO LAYER 2
           if all(B2(1:nActive,1) == -100)
	     nActive = nActive + 1;
                if nargin == 3
                 fprintf('            A new category %d created\n',nActive);
                end
               win = nActive;
           end

      % ELSE RESET NEURON DOES NOT FIRE: LAYERS 1 & 2 RESONATE
      else
           if nargin == 3
             fprintf('Pattern %d classified in %d category.\n',q,win);
           end
      resonance = 1;
      end

    end  % end of WHILE loop

    	% Update B-U LTM and T-D LTM
    	V(:,win) = a1(:,q)&V(:,win);
    	W(win,:) = (a1(:,q)*L/(L-1+sum(a1(:,q))))';
	%V(:,win) = V(:,win).*p(:,q);
	%W(win,:) = (V(:,win).*p(:,q))'/(0.5 + sum(V(:,win).*p(:,q)));
	 categ(win, count(win)) = q;
         count(win) = count(win)+1;
end	% end of FOR loop

% Display final classification
fprintf('\n  Category      Pattern\n  -----------------------------------------------\n');
 for i = 1:nActive
  fprintf('     %d           ',i);
   for j= 1:F2
	if categ(i,j)
         fprintf('%s, ',(categ(i,j)+64));
        end
   end
   fprintf('\n  ------------------------------------------------\n');
 end

author: Sable

repo: mcbench-benchmarks


text:
function triangulate ( prefix, animate )

%*****************************************************************************80
%
%% TRIANGULATE triangulates a polygon using the ear slicing algorithm.
%
%  Discussion:
%
%    The polygon is described by a file which lists the coordinates of the
%    vertices in consecutive order during a counterclockwise traversal.
%
%    Note that no two successive vertices can be equal, nor can the first
%    and last vertices be equal.
%
%    Information about the polygon is stored as a structure we will call
%    a "vertex array".
%
%    There is one entry in the array for each of the N vertices in the
%    polygon.  The I-th entry is a structure with the following fields:
%
%      vertex(i).index  the index of the node.
%      vertex(i).next   the index of the next node.
%      vertex(i).prev   the index of the previous node.
%      vertex(i).x      the X coordinate.
%      vertex(i).y      the Y coordinate.
%      vertex(i).ear    1 if the vertex represents an "ear" that can
%                       be sliced off.
%
%  Usage:
%
%    triangulate ( 'prefix', 'animate' )
%
%    where
%
%    * 'prefix' is the filename prefix
%    * 'animate' is 'y' if you want to see an animation of the triangulation.
%
%    reads
%
%    * 'prefix'_nodes.txt, the node coordinate file,
%
%    creates
%
%    * 'prefix'_diagonals.txt, the internal diagonals file.
%    * 'prefix'_elements.txt, the triangle file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2013
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, string PREFIX, the filename prefix.
%
%    Input, string ANIMATE, is 'y' if the user wants to see the triangulation
%    process animated.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read a file listing the coordinates of the vertices\n' );
  fprintf ( 1, '  of a polygon, and determine a triangulation.\n' );
%
%  If the user did not supply the filename prefix, request it.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATE:\n' );
    prefix = input ( '  Enter the filename prefix (in ''quotes''!).\n' );
  end
%
%  If the user did not supply the animation option, assume it is NO.
%
  if ( nargin < 2 )
    animate = 'n';
  end
%
%  Create filenames.
%
  node_filename     = strcat ( prefix, '_nodes.txt' );
  diagonal_filename = strcat ( prefix, '_diagonals.txt' );
  triangle_filename = strcat ( prefix, '_elements.txt' );
%
%  Read the node coordinates.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Processing the vertex coordinate file "%s".\n', node_filename )
  direction = +1;
  vertex = read_vertices ( node_filename, direction );

  vertex = ear_init ( vertex );

  area = area_poly2 ( vertex );
%
%  If area is zero, refuse to go on.
%
  if ( area == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATE - Fatal area!\n' );
    fprintf ( 1, '  The area is computed to be exactly zero.\n' );
    return
  end
%
%  If area is NEGATIVE, reverse the ordering.
%
  if ( area < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Polygon area is negative!\n' );
    fprintf ( 1, '  Will try reversing vertex ordering.\n' );
    direction = -1;
    vertex = read_vertices ( node_filename, direction );
    vertex = ear_init ( vertex );
    area = area_poly2 ( vertex );
  end

  n = length ( vertex );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of vertices is %d\n', n );

  print_vertices ( vertex );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Area of polygon is %g\n', area );
%
%  Refuse to triangulate a region in which two successive vertices are equal.
%
  for i = 1 : n
    j = vertex(i).next;
    if ( vertex(i).x == vertex(j).x && vertex(i).y == vertex(j).y )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TRIANGULATE - Fatal error!\n' );
      fprintf ( 1, '  Two neighbor vertices are equal.\n' );
      fprintf ( 1, '  Vertex(%d) = Vertex(%d)\n', i, j );
      error ( 'TRIANGULATE - Fatal error!' );
    end
  end
%
%  Determine a reasonable scale for the plot.
%
  xmin = vertex(1).x;
  xmax = vertex(1).x;
  ymin = vertex(1).y;
  ymax = vertex(1).y;
  for i = 2 : n
    xmin = min ( xmin, vertex(i).x );
    xmax = max ( xmax, vertex(i).x );
    ymin = min ( ymin, vertex(i).y );
    ymax = max ( ymax, vertex(i).y );
  end

  dx = 0.0125 * max ( xmax - xmin, ymax - ymin );
%
%  Draw the (non-triangulated) polygon.
%
%  You might think a single PLOT command would suffice to display
%  the vertices.  However, the proper way of extracting a vector of
%  values from a structure eludes me.
%
  figure ( 1 )

  clf

  grid on

  hold on
%
%  Make the plot a little bigger than the data.
%
  plot ( xmin - 2 * dx, ymin - 2 * dx, 'w.' );
  plot ( xmax + 2 * dx, ymin + 2 * dx, 'w.' );

  for i = 1 : n
    plot ( vertex(i).x, vertex(i).y, 'k.', 'MarkerSize', 25 )
  end

  for i = 1 : n
    text ( vertex(i).x + dx, vertex(i).y + dx, num2str ( i ) );
  end

  for i = 1 : n
    line ( [ vertex(i).x, vertex(vertex(i).next).x ], ...
           [ vertex(i).y, vertex(vertex(i).next).y ], ...
           'LineWidth', 3, 'Color', 'r' )
  end

  xlabel ( '<-- X -->' )
  ylabel ( '<-- Y -->' )
  title ( prefix );
  axis equal
%
%  Triangulate the polygon.
%
  triangles = polygon_triangulate ( vertex );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Internal diagonals for triangulation:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n - 3
    fprintf ( 1, '  %2d  %2d\n', triangles(i,1:2) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Triangles:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n - 2
    fprintf ( 1, '  %2d  %2d  %2d\n', triangles(i,1:3) );
  end
%
%  Draw the triangulated polygon.
%
  figure ( 2 )

  clf

  grid on

  hold on
%
%  Make the plot a little bigger than the data.
%
  plot ( xmin - 2 * dx, ymin - 2 * dx, 'w.' );
  plot ( xmax + 2 * dx, ymin + 2 * dx, 'w.' );

  for i = 1 : n
    plot ( vertex(i).x, vertex(i).y, 'k.', 'MarkerSize', 25 )
  end

  for i = 1 : n
    text ( vertex(i).x + dx, vertex(i).y + dx, num2str ( i ) );
  end

  for i = 1 : n
    line ( [ vertex(i).x, vertex(vertex(i).next).x ], ...
           [ vertex(i).y, vertex(vertex(i).next).y ], 'LineWidth', 2, 'Color', 'r' )
  end

  for i = 1 : n - 3
    i1 = triangles(i,1);
    i2 = triangles(i,2);
    line ( [ vertex(i1).x, vertex(i2).x ], ...
           [ vertex(i1).y, vertex(i2).y ], 'LineWidth', 3, 'Color', 'b' )
  end

  xlabel ( '<-- X -->' )
  ylabel ( '<-- Y -->' )
  title ( sprintf ( 'Triangulation of "%s"', prefix ) );
  axis equal
%
%  Optionally, animate the triangulation.
%
  if ( animate == 'y' )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  ANIMATION OPTION SELECTED!\n' );
    fprintf ( 1, '  Press return to see each diagonal added.\n' );

    figure ( 3 )

    clf

    grid on

    hold on
%
%  Make the plot a little bigger than the data.
%
    plot ( xmin - 2 * dx, ymin - 2 * dx, 'w.' );
    plot ( xmax + 2 * dx, ymin + 2 * dx, 'w.' );

    for i = 1 : n
      plot ( vertex(i).x, vertex(i).y, 'k.', 'MarkerSize', 25 )
    end

    for i = 1 : n
      text ( vertex(i).x + dx, vertex(i).y + dx, num2str ( i ) );
    end

    for i = 1 : n
      line ( [ vertex(i).x, vertex(vertex(i).next).x ], ...
             [ vertex(i).y, vertex(vertex(i).next).y ], 'LineWidth', 2, 'Color', 'r' )
    end

    for i = 1 : n - 3
      i1 = triangles(i,1);
      i2 = triangles(i,2);
      i3 = triangles(i,3);
      axis equal
      pause

      patch ( [ vertex(i1).x, vertex(i2).x, vertex(i3).x ], ...
              [ vertex(i1).y, vertex(i2).y, vertex(i3).y ], 'g' );

      line ( [ vertex(i1).x, vertex(i2).x ], ...
             [ vertex(i1).y, vertex(i2).y ], 'LineWidth', 3, 'Color', 'b' )

      axis equal

    end

    pause
    i = n - 2;
    i1 = triangles(i,1);
    i2 = triangles(i,2);
    i3 = triangles(i,3);
    line ( [ vertex(i1).x, vertex(i2).x ], ...
           [ vertex(i1).y, vertex(i2).y ], 'LineWidth', 3, 'Color', 'b' )

    patch ( [ vertex(i1).x, vertex(i2).x, vertex(i3).x ], ...
            [ vertex(i1).y, vertex(i2).y, vertex(i3).y ], 'g' );

    xlabel ( '<-- X -->' )
    ylabel ( '<-- Y -->' )
    title ( sprintf ( 'Triangulation of "%s"', prefix ) );
    axis equal

  end
%
%  Write the diagonals to a file.
%
  i4mat_write ( diagonal_filename, 2, n-3, triangles(1:n-3,1:2)' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Diagonals written to "%s".\n', diagonal_filename );
%
%  Write the triangles to a file.
%
  i4mat_write ( triangle_filename, 3, n-2, triangles(1:n-2,1:3)' );

  fprintf ( 1, '  Triangles written to "%s".\n', triangle_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATE\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function area = area_poly2 ( vertex )

%*****************************************************************************80
%
%% AREA_POLY2 computes the area of a polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, vertex array VERTEX(*).
%
%    Output, real AREA, the area of the polygon.
%
  area = 0.0;

  first = 1;

  i = first;

  ip1 = vertex(i).next;
  if ( ip1 == first )
    return
  end

  ip2 = vertex(ip1).next;
  if ( ip2 == first )
    return
  end

  while ( 1 )

    area = area + triangle_area ( vertex(i), vertex(ip1), vertex(ip2) );

    i = ip1;
    ip1 = ip2;
    ip2 = vertex(ip1).next;

    if ( ip2 == first )
      break
    end

  end

  return
end
function value = between ( va, vb, vc )

%*****************************************************************************80
%
%% BETWEEN is TRUE if vertex VC is between vertices VA and VB.
%
%  Discussion:
%
%    The points must be (numerically) collinear.
%
%    Given that condition, we take the greater of VA.X - VB.X and VA.Y - VB.Y
%    as a "scale" and check where VC's value lies.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, real VA, VB, VC, the vertices.
%
%    Output, integer VALUE, is 1 if VC is between VA and VB,
%    and 0 otherwise.
%
  if ( ~collinear ( va, vb, vc ) )
    value = 0;
    return
  end

  if ( abs ( va.y - vb.y ) < abs ( va.x - vb.x ) )
    xmax = max ( va.x, vb.x );
    xmin = min ( va.x, vb.x );
    value = ( xmin <= vc.x && vc.x <= xmax );
  else
    ymax = max ( va.y, vb.y );
    ymin = min ( va.y, vb.y );
    value = ( ymin <= vc.y && vc.y <= ymax );
  end

  return
end
function value = collinear ( va, vb, vc )

%*****************************************************************************80
%
%% COLLINEAR returns a measure of collinearity for three points.
%
%  Discussion:
%
%    In order to deal with collinear points whose coordinates are not
%    numerically exact, we compare the area of the largest square
%    that can be created by the line segment between two of the points
%    to (twice) the area of the triangle formed by the points.
%
%    If the points are collinear, their triangle has zero area.
%    If the points are close to collinear, then the area of this triangle
%    will be small relative to the square of the longest segment.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, real VA, VB, VC, the vertices.
%
%    Output, integer VALUE, is 1 if the points are judged to be collinear,
%    and 0 otherwise.
%
  area = 0.5 * ( ( vb.x - va.x ) * ( vc.y - va.y ) ...
               - ( vc.x - va.x ) * ( vb.y - va.y ) );
%
%  MATLAB's MAX function has decided to be unfriendly, and I am
%  too annoyed to try to figure out how to placate it when all
%  I want to say is M = max ( A, B, C )!
%
  area_max =                 ( va.x - vb.x )^2 + ( va.y - vb.y )^2;
  area_max = max ( area_max, ( vb.x - vc.x )^2 + ( vb.y - vc.y )^2 );
  area_max = max ( area_max, ( vc.x - va.x )^2 + ( vc.y - va.y )^2 );

  if ( area_max <= eps )
    value = 1;
  elseif ( 2.0 * abs ( area ) <= eps * area_max )
    value = 1;
  else
    value = 0;
  end

  return
end
function value = diagonal ( im1, ip1, vertex )

%*****************************************************************************80
%
%% DIAGONAL: VERTEX(IM1):VERTEX(IP1) is a proper internal diagonal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 May 2014
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, integer IM1, IP1, the indices of two vertices.
%
%    Input, vertex array VERTEX(*).
%
%    Output, integer VALUE, the value of the test.
%
  value = in_cone ( im1, ip1, vertex ) && ...
          in_cone ( ip1, im1, vertex ) && ...
          diagonalie ( im1, ip1, vertex );

  return
end
function value = diagonalie ( im1, ip1, vertex )

%*****************************************************************************80
%
%% DIAGONALIE: VERTEX(IM1):VERTEX(IP1) is a proper diagonal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, integer IM1, IP1, the indices of two vertices.
%
%    Input, vertex array VERTEX(*).
%
%    Output, int VALUE, the value of the test.
%
  first = im1;

  j = vertex(first).index;
  jp1 = vertex(first).next;

  value = 1;
%
%  For each edge VERTEX(J):VERTEX(JP1) of the polygon:

  while ( 1 )
%
%  Skip any edge that includes vertex IM1 or IP1.
%
    if ( j == im1 || j == ip1 || jp1 == im1 || jp1 == ip1 )

    else
      if ( intersect ( vertex(im1), vertex(ip1), vertex(j), vertex(jp1) ) )
        value = 0;
        return
      end
    end

    j = jp1;
    jp1 = vertex(j).next;

    if ( j == vertex(first).index )
      break
    end

  end

  return
end
function vertex = ear_init ( vertex )

%*****************************************************************************80
%
%% EAR_INIT initializes the "earity" of each vertex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input/output, vertex array VERTEX(*).  On output, the "ear" field
%    of the vertex array has been initialized.
%
  im1 = vertex(1).prev;
  i   = vertex(1).index;
  ip1 = vertex(1).next;

  while ( 1 )

    vertex(i).ear = diagonal ( im1, ip1, vertex );

    im1 = i;
    i   = ip1;
    ip1 = vertex(i).next;

    if ( i == vertex(1).index )
      break
    end

  end

  return
end
function value = i4_modp ( i, j )

%*****************************************************************************80
%
%% I4_MODP returns the nonnegative remainder of I4 division.
%
%  Discussion:
%
%    If
%      NREM = I4_MODP ( I, J )
%      NMULT = ( I - NREM ) / J
%    then
%      I = J * NMULT + NREM
%    where NREM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
%
%  Example:
%
%        I     J     MOD  I4_MODP    Factorization
%
%      107    50       7       7    107 =  2 *  50 + 7
%      107   -50       7       7    107 = -2 * -50 + 7
%     -107    50      -7      43   -107 = -3 *  50 + 43
%     -107   -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be divided.
%
%    Input, integer J, the number that divides I.
%
%    Output, integer VALUE, the nonnegative remainder when I is
%    divided by J.
%
  if ( j == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MODP - Fatal error!\n' );
    fprintf ( 1, '  Illegal divisor J = %d\n', j );
    error ( 'I4_MODP - Fatal error!' );
  end

  value = mod ( i, j );

  if ( value < 0 )
    value = value + abs ( j );
  end

  return
end
function value = i4_wrap ( ival, ilo, ihi )

%*****************************************************************************80
%
%% I4_WRAP forces an integer to lie between given limits by wrapping.
%
%  Example:
%
%    ILO = 4, IHI = 8
%
%    I   Value
%
%    -2     8
%    -1     4
%     0     5
%     1     6
%     2     7
%     3     8
%     4     4
%     5     5
%     6     6
%     7     7
%     8     8
%     9     4
%    10     5
%    11     6
%    12     7
%    13     8
%    14     4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IVAL, an integer value.
%
%    Input, integer ILO, IHI, the desired bounds for the integer value.
%
%    Output, integer I4_WRAP, a "wrapped" version of IVAL.
%
  jlo = min ( ilo, ihi );
  jhi = max ( ilo, ihi );

  wide = jhi - jlo + 1;

  if ( wide == 1 )
    value = jlo;
  else
    value = jlo + i4_modp ( ival - jlo, wide );
  end

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function value = in_cone ( im1, ip1, vertex )

%*****************************************************************************80
%
%% IN_CONE is TRUE if the diagonal VERTEX(IM1):VERTEX(IP1) is strictly internal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, integer IM1, IP1, the indices of two vertices.
%
%    Input, vertex array VERTEX(*).
%
%    Output, int VALUE, the value of the test.
%
  im2 = vertex(im1).prev;
  i = vertex(im1).next;

  if ( left_on ( vertex(im1), vertex(i), vertex(im2) ) )

    value = left ( vertex(im1), vertex(ip1), vertex(im2) ) && ...
            left ( vertex(ip1), vertex(im1), vertex(i) );

  else

    value = ~ ( left_on ( vertex(im1), vertex(ip1), vertex(i) ) && ...
                left_on ( vertex(ip1), vertex(im1), vertex(im2) ) );

  end

  return
end
function value = intersect ( va, vb, vc, vd )

%*****************************************************************************80
%
%% INTERSECT is true if lines VA:VB and VC:VD intersect.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, vertex VA, VB, VC, VD, four vertices.
%
%    Output, integer VALUE, the value of the test.
%
  if ( intersect_prop ( va, vb, vc, vd ) )
    value = 1;
  elseif ( between ( va, vb, vc ) || ...
           between ( va, vb, vd ) || ...
           between ( vc, vd, va ) || ...
           between ( vc, vd, vb ) )
    value = 1;
  else
    value = 0;
  end

  return
end
function value = intersect_prop ( va, vb, vc, vd )

%*****************************************************************************80
%
%% INTERSECT_PROP is TRUE if lines VA:VB and VC:VD have a proper intersection.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, vertex VA, VB, VC, VD, four vertices.
%
%    Output, integer VALUE, the result of the test.
%
  if ( collinear ( va, vb, vc ) || ...
       collinear ( va, vb, vd ) || ...
       collinear ( vc, vd, va ) || ...
       collinear ( vc, vd, vb ) )
    value = 0;
  else
    value = xor ( left ( va, vb, vc ), left ( va, vb, vd ) ) && ...
            xor ( left ( vc, vd, va ), left ( vc, vd, vb ) );
  end

  return
end
function value = left ( va, vb, vc )

%*****************************************************************************80
%
%% LEFT is true if VC is to the left of the line VA:VB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, vertex VA, VB, VC, three vertices.
%
%    Output, integer VALUE, the result of the test.
%
  value = ( 0 < triangle_area ( va, vb, vc ) );

  return
end
function value = left_on ( va, vb, vc )

%*****************************************************************************80
%
%% LEFT_ON is TRUE if VC is to the left of, or on, the line VA:VB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, vertex VA, VB, VC, three vertices.
%
%    Output, integer VALUE, the result of the test.
%
  value = ( 0 <= triangle_area ( va, vb, vc ) );

  return
end
function triangles = polygon_triangulate ( vertex )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE determines a triangulation of a polygon.
%
%  Discussion:
%
%    There are N-3 triangles in the triangulation.
%
%    For the first N-2 triangles, the first edge listed is always an
%    internal diagonal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2011
%
%  Author:
%
%    Original C version by Joseph ORourke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joseph ORourke,
%    Computational Geometry in C,
%    Cambridge, 1998,
%    ISBN: 0521649765,
%    LC: QA448.D38.
%
%  Parameters:
%
%    Input, vertex array VERTEX(N).
%
%    Output, integer TRIANGLES(N-2,3), the triangles of the triangulation.
%
  n = length ( vertex );
  first = 1;

  triangles = zeros ( n - 2, 3 );
  triangle_num = 0;

  i2 = first;

  while ( triangle_num < n - 3 )
%
%  If I2 is an ear, gather information necessary to carry out
%  the slicing operation and subsequent "healing".
%
    if ( vertex(i2).ear )

      i3 = vertex(i2).next;
      i4 = vertex(i3).next;
      i1 = vertex(i2).prev;
      i0 = vertex(i1).prev;
%
%  Make vertex I2 disappear.
%
      vertex(i1).next = i3;
      vertex(i2).index = 0;
      vertex(i3).prev = i1;
%
%  Update the earity of I1 and I3, because I2 disappeared.
%
      vertex(i1).ear = diagonal ( i0, i3, vertex );
      vertex(i3).ear = diagonal ( i1, i4, vertex );
%
%  Add the diagonal [I3, I1, I2] to the list.
%
      triangle_num = triangle_num + 1;
      triangles(triangle_num,1) = i3;
      triangles(triangle_num,2) = i1;
      triangles(triangle_num,3) = i2;

    end
%
%  Try the next vertex.
%
    i2 = vertex(i2).next;

  end
%
%  The last triangle is formed from the three remaining vertices.
%
  i3 = vertex(i2).next;
  i1 = vertex(i2).prev;
  triangle_num = triangle_num + 1;
  triangles(triangle_num,1) = i3;
  triangles(triangle_num,2) = i1;
  triangles(triangle_num,3) = i2;

  return
end
function print_vertices ( vertex )

%*****************************************************************************80
%
%% PRINT_VERTICES prints the vertex information.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, vertex array VERTEX(*).
%
  n = length ( vertex );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   I IM1 IP1      X           Y       Earity\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n

    if ( vertex(i).ear == -1 )
      earity = 'Unknown';
    elseif ( vertex(i).ear == 0 )
      earity = 'No';
    else
      earity = 'Yes';
    end

    fprintf ( 1, '  %2d  %2d  %2d  %10f  %10f  %s\n', ...
      vertex(i).index, vertex(i).prev, vertex(i).next, ...
      vertex(i).x, vertex(i).y, earity );

  end

  return
end
function vertex = read_vertices ( filename, direction )

%*****************************************************************************80
%
%% READ_VERTICES reads vertex coordinates and initializes the vertex structure.
%
%  Discussion:
%
%    If the polygon was found to have negative area on the first pass,
%    we try again, reversing the vertex ordering, and hoping the user simply
%    listed the vertices in clockwise order.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the name of the file to be read.
%
%    Input, integer DIRECTION, the ordering for the vertices.
%    +1, first vertex is 1, and count up.
%    -1, first vertex is N, and count down.
%
%    Output, vertex array VERTEX(*).
%
  xy = load ( filename );

  [ n, dim ] = size ( xy );

  if ( direction == -1 )
    xy(1:n,1) = xy(n:-1:1,1);
    xy(1:n,2) = xy(n:-1:1,2);
  end

  for i = 1 : n
    vertex(i).index = i;
    vertex(i).prev = i4_wrap ( i - 1, 1, n );
    vertex(i).next = i4_wrap ( i + 1, 1, n );
    vertex(i).x = xy(i,1);
    vertex(i).y = xy(i,2);
    vertex(i).ear = -1;
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function area = triangle_area ( v1, v2, v3 )

%*****************************************************************************80
%
%% TRIANGLE_AREA computes the signed area of a triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, vertex V1, V2, V3, three vertices.
%
%    Output, real AREA, the signed area of the triangle.
%
  area = 0.5 * ( ( v2.x - v1.x ) * ( v3.y - v1.y ) ...
               - ( v3.x - v1.x ) * ( v2.y - v1.y ) );

  return
end


text:
cur_dir = pwd;
cd(fileparts(mfilename('fullpath')));

try
    fprintf('Downloading region proposals...\n');
    urlwrite('https://onedrive.live.com/download?resid=F371D9563727B96F!91965&authkey=!AErVqYD6NhjxAfw', ...
        'proposals.zip');

    fprintf('Unzipping...\n');
    unzip('proposals.zip', '..');

    fprintf('Done.\n');
    delete('proposals.zip');
catch
    fprintf('Error in downloading, please try links in README.md https://github.com/daijifeng001/R-FCN');
end

cd(cur_dir);

repo: R-FCN


text:
%% The Class poleFigure
%
%%
%

author: mtex-toolbox

repo: mtex

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


function dataout=modeBSS(nummics, numrefs, datain)
%
% Use Aux-IVA to solve dr, aec and separation.
% nummics:              no. of mic channels
% numrefs:              no. of reference channels
% datain:               input data
% dataout:              output data
%

%% perform stft
addpath('stft2');
% fft size
fftsize=512;
stftshift=fftsize/2;

M=nummics;
R=numrefs;
N=M;

Xtf=cell(M+R, 1);
for m=1:M+R
    Xtf{m}=stft(datain(:, m), stftshift, fftsize, false);
end
[K, T]=size(Xtf{1});

Ytf=cell(M, 1);
for m=1:M
    Ytf{m}=zeros(K, T);
end

%% params
% dr filter length
DR_FLEN=5;
% forgetting factor for bss
BF_FORGET=0.999;
% the shape parameter of the source prior
GAMMA=0.2;
%
% used to keep stable
%
VAR_BIAS=0.01;
STABLE_EPS=1e-3;
BF_DIAGLOAD=1e-6;

%% space for bss
% data buffer
fsize=M+R+M*DR_FLEN;
Bufdata=zeros(K, fsize);
% current mic backup
Miccurrent=zeros(K, M);

% the weighted correlation matrices
C1=cell(K, 1);
C2=cell(K, 1);
for k=1:K
    C1{k}=STABLE_EPS*eye(fsize, fsize);
    C2{k}=STABLE_EPS*eye(fsize, fsize);
end

% demixing matrices
Demix=cell(K, 1);
for k=1:K
    Demix{k}=eye(fsize, fsize);
end

%% perform iteration
for tau=1:T
    Bssout=zeros(K, M);

    %
    % shift in new data
    %
    % shift old mic data back
    Bufdata=circshift(Bufdata, M, 2);
    % shift in delayed mic data
    Bufdata(:, M+R+1:M+R+M)=Miccurrent;

    % mic data backup
    for m=1:M
        Miccurrent(:, m)=Xtf{m}(:, tau);
    end
    % current mic data
    Bufdata(:, 1:M)=Miccurrent;

    % shift in reference data
    for r=1:R
        Bufdata(:, M+r)=Xtf{M+r}(:, tau);
    end

    %
    % calculate nonlinearity
    %
    phi1=0;
    phi2=0;

    for k=1:K
        x=Bufdata(k, :).';
        y=Demix{k}*x;
        % output data
        Bssout(k, :)=y(1:M).';

        phi1=phi1+abs(y(1))^2;
        phi2=phi2+abs(y(2))^2;
    end

    phi1=(1-BF_FORGET)*(phi1+VAR_BIAS)^((GAMMA-2)/2);
    phi2=(1-BF_FORGET)*(phi2+VAR_BIAS)^((GAMMA-2)/2);

    % update the demixing matrices
    for k=1:K
        %
        % accumulate the weighted correlation
        %
        x=Bufdata(k, :).';
        C1{k}=BF_FORGET*C1{k}+phi1*(x*x');
        C2{k}=BF_FORGET*C2{k}+phi2*(x*x');

        %
        % update demixing matrix
        %
        H1=inv(Demix{k}*C1{k}+BF_DIAGLOAD*eye(fsize, fsize));
        w1=H1(:, 1);

        H2=inv(Demix{k}*C2{k}+BF_DIAGLOAD*eye(fsize, fsize));
        w2=H2(:, 2);

        D=eye(fsize, fsize);
        D(1, :)=w1';
        D(2, :)=w2';

        %
        % solve the scaling ambiguity
        %
        A=inv(D);

        if abs(A(1, 1))>=abs(A(2, 1))
            a1=A(1, 1);
        else
            a1=A(2, 1);
        end

        if abs(A(2, 2))>=abs(A(1, 2))
            a2=A(2, 2);
        else
            a2=A(1, 2);
        end

        D=eye(fsize, fsize);
        D(1, :)=a1*w1';
        D(2, :)=a2*w2';

        Demix{k}=D;
    end

    for m=1:M
        Ytf{m}(:, tau)=Bssout(:, m);
    end
end

%% perform istft and output signal
dataout=zeros(dataLength(T, stftshift, fftsize ), N);
for n=1:N
    dataout(:, n)=istft(Ytf{n}, stftshift, false);
end

end


text:
clc;
close all;
clearvars;

n = 800000;
k = 100;
trials = 100;

sort_elapsed = 0;
quickselect_elapsed = 0;
for i=1:trials
    data = -randperm(n);
    tstart = tic;
    sort(data);
    sort_elapsed = sort_elapsed + toc(tstart);
    tstart = tic;
    y = spx.fast.quickselect(data, k);
    quickselect_elapsed = quickselect_elapsed + toc(tstart);
    fprintf('%d  %d\n', y(k), sum (y(1:k) < -k));
end

t1 = sort_elapsed;
t2 = quickselect_elapsed;
gain_x = t1 / t2;

fprintf('Sort: %.2f sec, quickselect: %.2f sec, Gain: %.2f\n', t1, t2, gain_x);

author: indigits

repo: sparse-plex


text:
function V=latticetype2matrix(L,a,M,lt);
%LATTICETYPE2MATRIX  Convert lattice description to matrix form
%   Usage: V=latticetype2matrix(L,a,M,lt);
%
%   `V=latticetype2matrix(L,a,M,lt)` converts a standard description of a
%   lattice using the *a*, *M* and *lt* parameters into a $2\times 2$
%   integer matrix description. The conversion is *only* valid for the
%   specified transform length *L*.
%
%   The output will be in lower triangular Hemite normal form.
%
%   For more information, see
%   `<http://en.wikipedia.org/wiki/Hermite_normal_form>`_.
%
%   An example:::
%
%     V = latticetype2matrix(120,10,12,[1 2])
%
%   See also: matrix2latticetype

L2=dgtlength(L,a,M,lt);

if L~=L2
    error('%s: Invalid transform length.',upper(mfilename));
end;

b=L/M;
s=b/lt(2)*lt(1);
V=[a 0;...
   s b];


text:
%  Copyright (c) 2014, Karen Simonyan
%  All rights reserved.
%  This code is made available under the terms of the BSD license (see COPYING file).

function result = evaluate(config, scores, gt)

    scores = reshape(scores, 1, []);

    switch config
        case 'ap'
            [res, extra] = evaluation.ap.eval(config, scores, gt);
        case 'roc'
            [res, extra] = evaluation.roc.eval(config, scores, gt);
        case 'accuracy'
            [res, extra] = evaluation.accuracy.eval(config, scores, gt);
    end

        % measure name
        result.meas_name = config;

        % measure value (a scalar)
        result.measure = res;

        % extra data in a struct (e.g. optimal thresh), or empty
        result.extra = extra;

end

author: AlfredXiangWu