text:
function out=add_circular_mask(in,th,norm)
if nargin < 3
norm=0;
end

	thickness=4;
	[N,M]=size(in);
	[ix,iy]=meshgrid([-N/2+1/2:N/2-1/2],[-N/2+1/2:N/2-1/2]);
	rad=ix.^2+iy.^2;
	mask=max(0,thickness-abs(sqrt(rad)-(th*N/2)));
	outermask=(rad < (th*N/2)^2);
	if norm==0
	out=in.*outermask+min(0,min(in(:)))*(1-outermask);
	mask=mask*(max(64,0*max(out(:)))/max(mask(:)));
	else
	out=in.*outermask+min(in(:))*(1-outermask);
	mask=mask*((max(out(:)))/max(mask(:)));
	end
	tope=max(out(:));
	%out=min(tope,out+mask);
	out=(mask/tope).*mask+(1-mask/tope).*out;

if 0
	cjet=colormap(gray);
	cjet=1-cjet;
	inr=in(:);
	if norm==1
	inr=(inr-min(inr))*(55/(max(inr)-min(inr)));
	end
	tempo=cjet(max(1,min(64,round(inr))),:);
	out=reshape(tempo,N,M,3);
	mtempo(:,:,1)=outermask;
	mtempo(:,:,2)=outermask;
	mtempo(:,:,3)=outermask;
	out=out.*mtempo+ (1-mtempo);
end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function []=panel2fdisp(N,n,T,Units,endo,Ymat,stringdates2,decimaldates2,Fstartlocation,Fendlocation,forecast_estimates,pref)

% preliminary task: gather in a cell the values to be plotted
% initiate the cell
plotdata={};
% because forecasts have to be computed for each unit, loop over units
for ii=1:N
% each cell entry is a matrix of actual and forecast values
% this matrix comprises 4 rows: the first row is actual data, while the three other rows are the estimates (point estimates and confidence bands) for the forecasts
% also, this matrix has a number of rows equal to the dimension of decimaldates2, which comprises the total period sample+forecasts
   % loop over variables
   for jj=1:n
      plotdata{jj,1,ii}=nan(4,size(decimaldates2,1));
      % record actual sample values
      plotdata{jj,1,ii}(1,1:T)=Ymat(:,jj,ii)';
      % copy the last point of the actual sample for the forecast part of the matrix (required to have a clean plot)
      plotdata{jj,1,ii}(:,Fstartlocation-1)=repmat(Ymat(Fstartlocation-1,jj,ii),4,1);
      % record forecast, lower bound
      plotdata{jj,1,ii}(2,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(1,:);
      % record forecast, point estimate
      plotdata{jj,1,ii}(3,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(2,:);
      % record forecast, upper bound
      plotdata{jj,1,ii}(4,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(3,:);
   end
end
% then plot the figure
if pref.plot
forecast=figure('Tag','BEARresults');
set(forecast,'Color',[0.9 0.9 0.9]);
set(forecast,'name','unconditional forecasts');
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(N,n,count)
   hold on
   Xpatch=[decimaldates2(Fstartlocation-1:Fendlocation,1)' fliplr((decimaldates2(Fstartlocation-1:Fendlocation,1))')];
   Ypatch=[plotdata{jj,1,ii}(2,Fstartlocation-1:Fendlocation) fliplr(plotdata{jj,1,ii}(4,Fstartlocation-1:Fendlocation))];
   Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(Fpatch,'facealpha',0.6);
   set(Fpatch,'edgecolor','none');
   plot(decimaldates2,plotdata{jj,1,ii}(3,:),'Color',[0.4 0.4 1],'LineWidth',2);
   plot(decimaldates2,plotdata{jj,1,ii}(1,:),'Color',[0 0 0],'LineWidth',2);
   hold off
   set(gca,'XLim',[decimaldates2(1,1) decimaldates2(end,1)],'FontName','Times New Roman');
   set(gca,'XGrid','on');
   set(gca,'YGrid','on');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(Units{ii,1},'FontWeight','normal');
      end
   end
end
end

% save on Excel
% create the cell that will be saved on excel
forecastcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},size(stringdates2,1)+3,1);
horzspace=repmat({''},3,6*n);
% loop over units
for ii=1:N
% initiate the cell of results
unitcell={};
   % loop over endogenous variables (horizontal dimension)
   for jj=1:n
   % create a header
   header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''} {''} {''};{''} {''} {''} {''} {''};{''} {'actual'} {'lower bound'} {'median'} {'upper bound'}];
   % complete the cell
   endocell=[[header;stringdates2 num2cell((plotdata{jj,1,ii})')] vertspace];
   % concatenate to the previous parts of unitcell
   unitcell=[unitcell endocell];
   end
% concatenate to the previous parts of afcell
forecastcell=[forecastcell;horzspace;unitcell];
end
% trim
forecastcell=forecastcell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),forecastcell,'forecasts','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function varargout = MDoc(varargin)
    %MDOC  mexopencv documentation utility
    %
    %     MDoc()
    %     MDoc(action)
    %     MDoc(topic)
    %     MDoc(..., 'OptionName',optionValue, ...)
    %     filename = MDoc(...)
    %
    % ## Input
    % * __action__ Special targets (default `-all`). One of:
    %   * `-clean`: delete all generated docs.
    %   * `-wiki`: generate HTML wiki files.
    %   * `-all`: generate full documentation for all mexopencv functions.
    %   * `-index`: only generate table of contents.
    %   * `-indexcat`: only generate table of contents (grouped by modules).
    %   * `-contents`: generate Contents.m file.
    %   * `-helptoc`: generate helpfuncbycat.xml file (TOC function reference).
    % * __topic__ generate docs for specified mexopencv function.
    %
    % ## Output
    % * __filename__ optional, generated HTML file path.
    %
    % ## Options
    % * __Force__ Unconditionally re-create HTML output files. default true
    % * __MarkdownRender__ Process markdown. default true
    % * __MarkdownBackend__ Markdown parser implementation to use. The
    %   following backends are available to choose from:
    %   * 'CommonMark-Java' (https://github.com/atlassian/commonmark-java/)
    %   * 'MarkdownPapers' (https://github.com/lruiz/MarkdownPapers/)
    %   * 'pegdown' (https://github.com/sirthias/pegdown/)
    %   * 'MarkdownJ' (https://github.com/myabc/markdownj/)
    %   * 'Markdown4j' (https://github.com/jdcasey/markdown4j/)
    %   * 'Txtmark' (https://github.com/rjeschke/txtmark/)
    % * __AutoLinking__ search for "http://" in text and convert them into
    %   hyperlinks. default false
    % * __AutoLinkingDeep__ search for "www." in text and convert them into
    %   hyperlinks. default false
    % * __AutoLinkingCV__ search for "cv.*" functions in text, and convert
    %   them into hyperlinks. default true
    % * __SyntaxHighlighting__ syntax highlighted code. default true
    % * __PrettifyLang__ language of syntax highlighting. default 'lang-matlab'
    % * __TableSorting__ sortable HTML tables (used in index.html).
    %   default true
    % * __NoExternals__ download all CSS/JS external resources locally.
    %   default true
    %
    % ## Options (batch mode)
    % * __IndexGroups__ group by OpenCV module in index. default true
    % * __IgnoreHandleClass__ skip generating redundant docs for methods
    %   inherited from "handle" base class. default true
    % * __Verbose__ Display info messages. default true
    % * __Progress__ Show progress bar. default true
    %
    % This function generates simple MATLAB HTML documentation using
    % MATLAB's internal utility help2html.
    %
    % See also: doc, help, help2html, publish
    %

    % validate input arguments
    nargoutchk(0, 1);
    opts = parse_options(varargin{:});

    % initialize Java dependencies
    init_jsoup(opts);
    init_markdown(opts);

    if strcmp(opts.topic, '-clean')
        % clean all docs
        if isdir(opts.DIR)
            if opts.Verbose, fprintf('Cleaning %s\n', opts.DIR); end
            rmdir(opts.DIR, 's');
        end
        return;
    elseif ~any(strcmp(opts.topic, {'-wiki', '-contents', '-helptoc'}))
        if ~isdir(opts.DIR)
            % prepare output directory and CSS stylesheets
            if opts.Verbose, fprintf('Creating %s\n', opts.DIR); end
            mkdir(opts.DIR);
            copy_stylesheets(opts);
        end
    end

    % generate docs
    [filename, isbatch, needed] = get_output_filename(opts, opts.topic);
    if needed
        switch opts.topic
            case '-wiki'
                download_wiki(opts);
                generate_wiki(opts);
            case '-all'
                generate_all_docs(opts);
            case {'-index', '-indexcat'}
                txt = generate_index(opts, opts.topic);
            case '-contents'
                txt = generate_contents_m(opts);
            case '-helptoc'
                txt = generate_helptoc_xml(opts);
            otherwise
                txt = generate_doc(opts, opts.topic);
        end
        if ~isbatch
            % write to file
            assert(~isempty(txt), 'No doc generated');
            filewrite(filename, txt);
        end
    end

    % return/show result
    if nargout > 0
        % return output filename
        varargout{1} = filename;
    else
        if ~isbatch
            % open documentation in embedded browser or editor
            [~,~,ext] = fileparts(filename);
            if strcmp(ext, '.html')
                web(filename, '-new');
            else
                edit(filename);
            end
        else
            % show output directory
            if ispc()
                winopen(filename);
            else
                disp(filename);
            end
        end
    end
end

function opts = parse_options(varargin)
    %PARSE_OPTIONS  Parse function input argument

    % helper function to validate true/false arguments
    isbool = @(x) isscalar(x) && (islogical(x) || isnumeric(x));

    % parse input arguments
    p = inputParser();
    p.addOptional('topic', '-all', @ischar);
    p.addParameter('Force', true, isbool);
    p.addParameter('MarkdownRender', true, isbool);
    p.addParameter('MarkdownBackend', 'commonmark-java', @ischar);
    p.addParameter('AutoLinking', false, isbool);
    p.addParameter('AutoLinkingDeep', false, isbool);
    p.addParameter('AutoLinkingCV', true, isbool);
    p.addParameter('SyntaxHighlighting', true, isbool);
    p.addParameter('PrettifyLang', 'lang-matlab', @ischar);
    p.addParameter('TableSorting', true, isbool);
    p.addParameter('NoExternals', true, isbool);
    p.addParameter('IndexGroups', true, isbool);
    p.addParameter('IgnoreHandleClass', true, isbool);
    p.addParameter('Verbose', true, isbool);
    p.addParameter('Progress', true, isbool);

    % options struct
    p.parse(varargin{:});
    opts = p.Results;

    % documentation output directory
    opts.DIR = fullfile(mexopencv.root(), 'doc', 'matlab');

    % validate topic/action
    if strncmp(opts.topic, '-', 1)
        opts.topic = validatestring(opts.topic, ...
            {'-clean', '-wiki', '-all', '-index', '-indexcat', '-contents', '-helptoc'});
    elseif ~strncmp(opts.topic, 'cv.', 3)
        % prepend "cv." prefix to topic name if necessary
        opts.topic = ['cv.' opts.topic];
    end

    % list of recognized backends
    backends = {
        'commonmark-java'
        'MarkdownPapers'
        'pegdown'
        'markdownj'
        'markdown4j'
        'txtmark'
    };
    opts.MarkdownBackend = validatestring(opts.MarkdownBackend, backends);

    opts.Ignore = {};
    if opts.IgnoreHandleClass
        % methods/properties/events inherited from "handle"
        opts.Ignore = union(opts.Ignore, methods('handle'));
        opts.Ignore = union(opts.Ignore, properties('handle'));
        opts.Ignore = union(opts.Ignore, events('handle'));
        opts.Ignore = setdiff(opts.Ignore, 'delete');  % exclude dtor
    end
end

function init_jsoup(opts)
    %INIT_JSOUP  Download and setup the jsoup Java library

    persistent initialized
    if isempty(initialized), initialized = false; end
    if ~initialized
        % download if it doesn't exist
        dname = fullfile(mexopencv.root(), 'utils', 'jars', 'jsoup');
        if ~isdir(dname)
            MavenDownload('org.jsoup', 'jsoup', 'OutDir',dname, ...
                'Verbose',opts.Verbose);
        end

        % add all JAR files to classpath
        if opts.Verbose, disp('Initializing jsoup'); end
        JavaAddJarDir(dname);

        % mark as initialized
        initialized = true;
    end
end

function init_markdown(opts)
    %INIT_MARKDOWN  Download and setup backend Markdown Java libraries

    persistent initialized
    if isempty(initialized), initialized = {}; end
    if ~ismember(opts.MarkdownBackend, initialized)
        % download if it doesn't exist
        dname = fullfile(mexopencv.root(), 'utils', 'jars', opts.MarkdownBackend);
        if ~isdir(dname)
            % group and artifact ids
            switch opts.MarkdownBackend
                case 'commonmark-java'
                    groupId = 'com.atlassian.commonmark';
                    artifactId = 'commonmark';
                case 'MarkdownPapers'
                    groupId = 'org.tautua.markdownpapers';
                    artifactId = 'markdownpapers-core';
                case 'pegdown'
                    groupId = 'org.pegdown';
                    artifactId = 'pegdown';
                case 'markdownj'
                    groupId = 'org.markdownj';
                    artifactId = 'markdownj-core';
                case 'markdown4j'
                    groupId = 'org.commonjava.googlecode.markdown4j';
                    artifactId = 'markdown4j';
                case 'txtmark'
                    groupId = 'com.github.rjeschke';
                    artifactId = 'txtmark';
                otherwise
                    error('Unrecognized Markdown backend')
            end
            MavenDownload(groupId, artifactId, 'OutDir',dname, ...
                'Verbose',opts.Verbose);
        end

        % add all JAR files to classpath
        if opts.Verbose, fprintf('Initializing %s\n', opts.MarkdownBackend); end
        JavaAddJarDir(dname);

        % mark as initialized
        initialized = union(initialized, opts.MarkdownBackend);
    end
end

function JavaAddJarDir(dname)
    %JAVAADDJARDIR  Add all jar files inside a directory to Java class path

    % list of JAR files inside directory
    jars = dir(fullfile(dname, '*.jar'));
    jars = {jars.name};

    % only keep those not already on the Java class path
    jcp = javaclasspath();
    idx = cellfun(@(j) all(cellfun(@isempty, strfind(jcp, j))), jars);
    jars = jars(idx);

    % put remaining jars on classpath
    if ~isempty(jars)
        jars = cellfun(@(j) fullfile(dname, j), jars, 'UniformOutput',false);
        javaaddpath(jars);
    end
end

function copy_stylesheets(opts)
    %COPY_STYLESHEETS  Copy CSS stylesheets to output directory

    % copy default stylesheet locally
    CSS = fullfile(toolboxdir('matlab'), 'helptools', 'private', 'helpwin.css');
    copyfile(CSS, opts.DIR);

    % copy second stylesheet with customization
    CSS = fullfile(fileparts(mfilename('fullpath')), 'helpwin_custom.css');
    copyfile(CSS, opts.DIR);
end

function [filename, isbatch, needed] = get_output_filename(opts, topic, fname)
    %GET_OUTPUT_FILENAME  Output filename according to action

    isbatch = false;
    switch topic
        case '-wiki'
            if nargin > 2
                [~,fname,~] = fileparts(fname);
                filename = [strrep(lower(fname), '-', '_') '.html'];
                filename = fullfile(mexopencv.root(), 'doc', 'wiki', filename);
            else
                filename = fullfile(mexopencv.root(), 'doc', 'wiki');
                isbatch = true;
            end
        case '-all'
            filename = opts.DIR;
            isbatch = true;
        case '-index'
            filename = fullfile(opts.DIR, 'index.html');
        case '-indexcat'
            filename = fullfile(opts.DIR, 'index_cat.html');
        case '-contents'
            filename = fullfile(mexopencv.root(), 'Contents.m');
        case '-helptoc'
            filename = fullfile(mexopencv.root(), 'doc', 'helpfuncbycat.xml');
        otherwise
            filename = fullfile(opts.DIR, [topic '.html']);
    end

    % checks if file already exists
    needed = isbatch || opts.Force || exist(filename, 'file') ~= 2;
end

function files = download_wiki(opts)
    %DOWNLOAD_WIKI  Download wiki files from GitHub

    % list of Wiki files
    dname = get_output_filename(opts, '-wiki');
    files = {
        'README.md'
        'Home.md'
        'Getting-Started.md'
        'Developing-a-new-MEX-function.md'
        'Gotchas.md'
        'Installation-(Windows,-MATLAB,-OpenCV-3).md'
        'Installation-(Windows,-Octave,-OpenCV-3).md'
        'Installation-(Linux,-MATLAB,-OpenCV-3).md'
        'Installation-(Linux,-Octave,-OpenCV-3).md'
        'Installation-(Windows,-MATLAB,-OpenCV-2).md'
        'Troubleshooting-(Windows).md'
        'Troubleshooting-(UNIX).md'
    };

    % download from GitHub
    % (git clone https://github.com/kyamagu/mexopencv.wiki.git)
    url = 'https://rawgit.com/wiki/kyamagu/mexopencv/';
    if ~isdir(dname), mkdir(dname); end
    copyfile(fullfile(mexopencv.root(), 'README.markdown'), fullfile(dname, files{1}));
    for i=2:numel(files)
        fname = fullfile(dname, files{i});
        if exist(fname, 'file') ~= 2 % || opts.Force
            if opts.Verbose, fprintf('Downloading %s...\n', files{i}); end
            urlwrite([url files{i}], fname);
        end
    end
end

function generate_wiki(opts)
    %GENERATE_WIKI  Generate HTML docs from Wiki files

    % process markdown files
    dname = get_output_filename(opts, '-wiki');
    files = dir(fullfile(dname, '*.md'));
    for i=1:numel(files)
        % convert and save as HTML
        fname = fullfile(dname, files(i).name);
        [filename, ~, needed] = get_output_filename(opts, '-wiki', fname);
        if needed
            txt = wiki_md2html(opts, fname);
            assert(~isempty(txt), 'Failed to convert Markdown');
            filewrite(filename, txt);
        end
    end
end

function generate_all_docs(opts)
    %GENERATE_ALL_DOCS  Generate full docs for all functions

    % create index.html
    [filename, ~, needed] = get_output_filename(opts, '-index');
    if needed
        if opts.Verbose, disp('Creating index...'); end
        txt = generate_index(opts, '-index');
        assert(~isempty(txt), 'Failed to generate index');
        filewrite(filename, txt);
    end

    % create index_cat.html
    [filename, ~, needed] = get_output_filename(opts, '-indexcat');
    if needed
        if opts.Verbose, disp('Creating index...'); end
        txt = generate_index(opts, '-indexcat');
        assert(~isempty(txt), 'Failed to generate index');
        filewrite(filename, txt);
    end

    % get a list of functions/classes
    list = enumerate_mexopencv_members(opts);

    % show progress
    if opts.Progress
        hWait = waitbar(0, 'Generating docs...', 'Name',mfilename(), ...
            'CreateCancelBtn','setappdata(gcbf,''cancel'',true)');
        setappdata(hWait, 'cancel',false);
        wbCleanObj = onCleanup(@() delete(hWait));
    end

    % process all functions
    if opts.Verbose, disp('Creating docs...'); end
    for i=1:numel(list)
        % next topic to process
        topic = list{i};
        if opts.Progress
            waitbar(i/numel(list), hWait, strrep(topic, '_', '\_'));
            if getappdata(hWait, 'cancel'), break; end
        end

        % generate HTML doc file
        [filename, ~, needed] = get_output_filename(opts, topic);
        if needed
            txt = generate_doc(opts, topic);
            if ~isempty(txt)
                filewrite(filename, txt);
                status = 'DONE';
            else
                % non-existant cv function
                % (doc/help/help2html cannot process hidden and private methods)
                %TODO: it seems thats no longer the case in recent versions,
                % we need to filter them in enumerate_mexopencv_members
                status = 'MISSING';
            end
            if opts.Verbose, fprintf('[%s] %s\n', status, topic); end
        end
    end
end

function txt = generate_doc(opts, topic)
    %GENERATE_DOC  Generate HTML docs for a function, class method/member, etc.

    % convert M-help into parsed HTML document
    jdoc = parse_help2html(opts, topic);
    if isempty(jdoc)
        txt = '';
        return;
    end

    % convert document back to HTML string
    txt = char(jdoc.toString());
end

function txt = generate_index(opts, topic)
    %GENERATE_INDEX  Create index.html/index_cat.html containing list of all functions

    % HTML doc of cv package function/class names and their H1-lines
    opts.MarkdownRender = false;  % turn off markdown processing for index
    jdoc = parse_help2html(opts, 'cv');
    if isempty(jdoc)
        txt = '';
        return;
    end

    % replace title near the top
    div = jdoc.select('div.title').empty().first();
    assert(~isempty(div), '
 not found');
    a = div.appendElement('a');
    a.text('mexopencv').attr('href', 'https://github.com/kyamagu/mexopencv');
    % append description after it
    p = jdoc.createElement('p').addClass('h1line');
    p.appendText('Collection and development kit of MATLAB MEX functions for ');
    p.appendElement('a').text('OpenCV').attr('href', 'https://opencv.org/');
    p.appendText(' library.');
    div.after(p);

    % extract listing of cv package function/class names and their H1-lines
    pre = jdoc.select('div.helptext > pre').first();
    assert(~isempty(pre), '
 not found');
    txt = char(pre.text());
    t = create_files_table(opts, txt);

    % format it as HTML table(s)
    switch topic
        case '-index'
            % one big table
            txt = table_to_index_html(opts, t);
        case '-indexcat'
            % one table per module
            txt = table_to_index_cat_html(opts, t);
        otherwise
            txt = '';
    end

    % insert table(s)
    div = pre.parent();
    assert(~isempty(div), '
 not found');
    div.removeClass('helptext').addClass('helpcontent');
    div.html(txt);

    % style table(s)
    tabl = div.select('table');
    assert(~isempty(tabl), '
 not found');
    tabl.addClass('table');
    if opts.TableSorting
        tabl.addClass('sortable');  % enable_sortable_tables
    end

    % convert document back to HTML string
    txt = char(jdoc.toString());
end

function txt = generate_contents_m(opts)
    %GENERATE_CONTENTS_M  Create Contents.m containing list of all functions

    % cv package function/class names and their H1-lines
    txt = help('cv');
    if isempty(txt)
        return;
    end

    % parse text into table
    t = create_files_table(opts, txt);

    % format table as contents.m file
    txt = table_to_contents_m(opts, t);
end

function txt = generate_helptoc_xml(opts)
    %GENERATE_HELPTOC_XML  Create helpfuncbycat.xml containing list of all functions

    % cv package function/class names and their H1-lines
    txt = help('cv');
    if isempty(txt)
        return;
    end

    % parse text into table
    t = create_files_table(opts, txt);

    % format table as helpfuncbycat.xml file
    txt = table_to_helptoc_xml(opts, t);
end

function t = create_files_table(opts, txt)
    %CREATE_FILES_TABLE  Build table containing list of all mexopencv files

    % split by lines
    C = textscan(txt, '%s', 'Delimiter','\n');
    C = C{1};
    C(cellfun(@isempty, C)) = [];

    % remove section headers
    str = 'Contents of ';
    C(strncmp(C, str, numel(str))) = [];

    % remove simcoverage entry cv.schema
    str = 'schema ';
    C(strncmp(C, str, numel(str))) = [];

    % split each line: "name - description"
    C = regexp(C, ' - ', 'split', 'once');
    C = strtrim(cat(1, C{:}));

    % add cv. prefix to name
    C(:,1) = strcat('cv.', C(:,1));

    % create table
    t = cell2table(C, 'VariableNames',{'names','descriptions'});

    if opts.IndexGroups
        % extract OpenCV module names for each mexopencv function
        tt = enumerate_mexopencv_files(opts);

        % join the two tables by name
        t = innerjoin(t, tt, 'Keys','names');

        % sort table
        t = sortrows(t, {'repos', 'modules', 'names'});
    end
end

function txt = table_to_index_html(opts, t)
    %TABLE_TO_INDEX_HTML  Format table of mexopencv files as HTML table

    if opts.IndexGroups
        linkify = @(p) strrep(strrep(p, mexopencv.root(), ...
            'https://github.com/kyamagu/mexopencv/blob/master'), '\', '/');
        txt = strcat('
', ...
            '', t.names, '	', ...
            '', t.modules, '	', ...
            '', t.repos, '	', ...
            '', t.descriptions, '	', ...
            'C ', ...
            'M ', ...
            'T	', ...
            '
');
        txt = strjoin({
            '
'
            '

Name	Module	Repo	'
            'Description	Source Files

'
            '
'
            strjoin(txt, '\n')
            '
'
            '
'
        }, '\n');
    else
        txt = strcat(...
            '
', t.names, '	', ...
            '', t.descriptions, '
');
        txt = strjoin({
            '
'
            '

Name	Description

'
            '
'
            strjoin(txt, '\n')
            '
'
            '
'
        }, '\n');
    end
end

function txt = table_to_index_cat_html(opts, t)
    %TABLE_TO_INDEX_CAT_HTML  Format table of mexopencv files as HTML tables

    opts.IndexGroups = false; % for the table_to_index_html calls below

    % format list of files as HTML tables, one per module
    txt = {};
    [repos,~,rid] = unique(t.repos);
    for r=1:max(rid)
        tt = t(rid == r,:);
        [modules,~,mid] = unique(tt.modules);
        txt{end+1} = sprintf('
%s
', repos{r}, repos{r});
        for m=1:max(mid)
            ttt = tt(mid == m,:);
            txt{end+1} = sprintf('
%s
', modules{m}, modules{m});
            txt{end+1} = table_to_index_html(opts, ttt);
        end
    end
    txt = strjoin(txt, '\n');
end

function txt = table_to_contents_m(opts, t)
    %TABLE_TO_CONTENTS_M  Format table of mexopencv files as contents.m file

    % formatting string for each entry
    maxlen = max(cellfun(@length, t.names));
    frmt = ['%%   %-' int2str(maxlen) 's - %s\n'];

    % format list of files as Contents.m text
    if opts.IndexGroups
        txt = {};
        [repos,~,rid] = unique(t.repos);
        for r=1:max(rid)
            tt = t(rid == r,:);
            [modules,~,mid] = unique(tt.modules);
            txt{end+1} = ['%% ' repos{r} ':'];
            txt{end+1} = '%';
            for m=1:max(mid)
                ttt = tt(mid == m,:);
                C = [ttt.names, ttt.descriptions]';
                txt{end+1} = ['% ' modules{m} ':'];
                txt{end+1} = [sprintf(frmt, C{:}) '%'];
            end
        end
        txt = strjoin(txt, '\n');
    else
        C = [t.names, t.descriptions]';
        txt = [sprintf(frmt, C{:}) '%'];
    end

    % add header
    txt = strjoin({
        '% mexopencv'
        sprintf('%% Version %s (R%s) %s', cv.Utils.version(), ...
            version('-release'), datestr(now(), 'dd-mmmm-yyyy'))
        '%'
        txt
        ''
    }, '\n');
end

function txt = table_to_helptoc_xml(opts, t)
    %TABLE_TO_HELPTOC_XML  Format table of mexopencv files as helpfuncbycat.xml file

    %HACK: description should not be included in R2014b and up
    % (only supported in older help browser "doc -classic")
    withPurpose = verLessThan('matlab','8.4');

    % format list of files as functions reference for helptoc.xml/helpfuncbycat.xml
    txt = {};
    txt{end+1} = '<?xml version="1.0" encoding="utf-8"?>';
    txt{end+1} = '<toc version="2.0">';
    txt{end+1} = '  <tocitem image="HelpIcon.FUNCTION" target="matlab/index.html">Functions';
    if opts.IndexGroups
        [repos,~,rid] = unique(t.repos);
        for r=1:max(rid)
            tt = t(rid == r,:);
            [modules,~,mid] = unique(tt.modules);
            txt{end+1} = sprintf(...
                '    <tocitem target="matlab/index_cat.html#%s">%s', ...
                repos{r}, repos{r});
            for m=1:max(mid)
                ttt = tt(mid == m,:);
                if withPurpose
                    C = strcat(...
                        '        <tocitem image="HelpIcon.FUNCTION" target="matlab/', ...
                        ttt.names, '.html"><name>', ttt.names, '</name><purpose>', ...
                        xml_escape(ttt.descriptions), '</purpose></tocitem>');
                else
                    C = strcat(...
                        '        <tocitem image="HelpIcon.FUNCTION" target="matlab/', ...
                        ttt.names, '.html">', ttt.names, '</tocitem>');
                end
                txt{end+1} = sprintf(...
                    '      <tocitem target="matlab/index_cat.html#%s">%s', ...
                    modules{m}, modules{m});
                txt{end+1} = strjoin(C, '\n');
                txt{end+1} = '      </tocitem>';
            end
            txt{end+1} = '    </tocitem>';
        end
    else
        if withPurpose
            C = strcat(...
                '    <tocitem image="HelpIcon.FUNCTION" target="matlab/', ...
                t.names, '.html"><name>', t.names, '</name><purpose>', ...
                xml_escape(t.descriptions), '</purpose></tocitem>');
        else
            C = strcat(...
                '    <tocitem image="HelpIcon.FUNCTION" target="matlab/', ...
                t.names, '.html">', t.names, '</tocitem>');
        end
        txt{end+1} = strjoin(C, '\n');
    end
    txt{end+1} = '  </tocitem>';
    txt{end+1} = '</toc>';
    txt{end+1} = '';
    txt = strjoin(txt, '\n');
end

function txt = wiki_md2html(opts, mdFile)
    %WIKI_MD2HTML  Convert Wiki Markdown file to HTML

    % read and convert Markdown to HTML
    txt = fileread(mdFile);
    txt = MarkdownProcess(txt, opts.MarkdownBackend);

    % full HMTL page from fragment
    [~,name] = fileparts(mdFile);
    C = {
        '<!DOCTYPE html>'
        ''
        ''
        '<meta charset="utf-8">'
        ['<title>' strrep(name, '-', ' ') '</title>']
        '<link rel="stylesheet" type="text/css" href="../matlab/helpwin.css">'
        '<link rel="stylesheet" type="text/css" href="../matlab/helpwin_custom.css">'
    };
    if opts.SyntaxHighlighting
        C = [
            C
            '<link rel="stylesheet" type="text/css" href="../matlab/matlab.min.css">'
            ''
            ''
            txt
            ''
            ''
            ''
            ''
            ''
        ];
    else
        C = [
            C
            ''
            ''
            txt
            ''
            ''
        ];
    end
    C{end+1} = '';
    txt = strjoin(C, '\n');
end

function jdoc = parse_help2html(opts, topic)
    %PARSE_HELP2HTML  Parse HELP2HTML output and perform common filtering

    % convert M-help into HTML document
    [txt, found] = help2html(topic);
    if ~found
        jdoc = [];
        return;
    end

    % parse HTML string
    jdoc = javaMethod('parse', 'org.jsoup.Jsoup', txt);
    jdoc.outputSettings.prettyPrint(false);  % for proper markdown parsing

    % replace CSS stylesheet with our local version (relative URI)
    jdoc.select('link[href$=helpwin.css]').attr('href', 'helpwin.css');

    % add another sytlesheet for customization (relative URI)
    inject_css_file(jdoc, 'helpwin_custom.css');

    % change title
    txt = char(jdoc.title());
    txt = strrep(txt, 'MATLAB File Help: ', '');
    txt = strrep(txt, ' - MATLAB File Help', '');
    jdoc.title([txt ' - mexopencv']);

    % add generator meta tag
    mt = append_element(jdoc.head(), 'meta');
    mt.attr('name', 'generator').attr('content', ['MATLAB ' version()]);

    % remove highlighting of function names
    % (turn func into func)
    jdoc.select('span.helptopic').unwrap();

    % rewrite .. hyperlinks (mainly in header/footer)
    rewrite_hyperlinks(opts, jdoc);

    if opts.MarkdownRender
        % render Markdown
        render_markdown(opts, jdoc);

        % mark the first paragraph which corresponds to H1-line
        p = jdoc.select('div.helpcontent > p:first-child').first();
        if ~isempty(p)
            p.addClass('h1line');

            % add description meta tag (using H1-line description)
            desc = [topic ' - ' char(p.text())];
            mt = append_element(jdoc.head(), 'meta');
            mt.attr('name', 'description').attr('content', desc);
        end
    end

    % set external links to open in new tabs
    jdoc.select('a[href^=http]').attr('target', '_blank');

    % enable syntax highlighting of code blocks using google-code-prettify
    if opts.SyntaxHighlighting && ~strcmp(topic, 'cv')
        enable_syntax_highlighting(opts, jdoc);
    end

    % enable sorting tables (only used in index.html)
    if opts.TableSorting && strcmp(topic, 'cv')
        enable_sortable_tables(opts, jdoc);
    end
end

function rewrite_hyperlinks(opts, jdoc)
    %REWRITE_HYPERLINKS  Rewrite hyperlinks for offline use

    % replace "matlab:helpwin" (i.e "Default Topics") with index link
    jdoc.select('a[href=matlab:helpwin]').html('Index').attr('href', 'index.html');

    % process rest of "matlab:(helpwin|doc) cv.fcn" links
    links = jdoc.select('a[href^=matlab:helpwin], a[href^=matlab:doc]');
    it = links.iterator();  % iterate over links as ArrayList<Element>
    while it.hasNext()
        % get href
        a = it.next();
        topic = char(a.attr('href'));

        % remove "matlab:(helpwin|doc)" prefix
        topic = regexprep(topic, 'matlab:(helpwin|doc)\s*', '', 'once');

        % remove enclosing parens/quotes in case of function syntax
        % "matlab:helpwin('x')" vs. command syntax "matlab:helpwin x"
        topic = regexprep(topic, '^\(''', '', 'once');
        topic = regexprep(topic, '''\);?$', '', 'once');

        % fix class methods (i.e: cv.Class/method -> cv.Class.method)
        topic = strrep(topic, '/', '.');

        % check if link to another cv function
        if strncmp(topic, 'cv.', 3)
            if is_handle_member(opts, topic)
                % replace link to handle.* methods
                url = 'https://www.mathworks.com/help/matlab/ref/handle-class.html';
                a.attr('href', url);
            else
                % replace link with a relative "fcn.html" link
                a.attr('href', [topic '.html']);
            end
        else
            % check if link to an official TMW function
            url = is_mathworks_function(opts, topic);
            if ~isempty(url)
                a.attr('href', url);
            end
        end
    end

    % remove some raw-code links ("matlab:open .." and "matlab:edit ..")
    jdoc.select('a[href~=matlab:(open|edit)]').remove();

    % turn any remaning "matlab:" links to plain text
    jdoc.select('a[href^=matlab:]').unwrap();
end

function render_markdown(opts, jdoc)
    %RENDER_MARKDOWN  Convert help text from Markdown to HTML + autolinking

    % extract markdown text
    pre = jdoc.select('div.helptext > pre').first();
    assert(~isempty(pre), '
 not found');
    %NOTE: this also undoes http:// autolinking done by help2html
    txt = char(pre.text());

    %HACK: for correct Markdown indentation, we remove two leading spaces:
    % - first one inserted by help2html for some reason
    % - second one is part of when we write a comment (e.g "% some text")
    %   being the space between the % symbol and the first character
    txt = regexprep(txt, '^  ', '', 'lineanchors');

    % remove redundant function name at the beginning (exclude property docs)
    if jdoc.select('div.sectiontitle:containsOwn(Property Details)').isEmpty()
        txt = regexprep(txt, '^\w+\s+', '', 'once');
    end

    % autolink "http://" and "www." text
    % (before rendering MD so that we dont later double-link  links)
    if opts.AutoLinking
        txt = auto_link_http(opts, txt);
    end
    if opts.AutoLinkingDeep
        txt = auto_link_www(opts, txt);
    end

    % convert Markdown to HTML
    txt = MarkdownProcess(txt, opts.MarkdownBackend);

    % autolink "cv." functions text
    % (after rendering MD so that we dont break indented code with cv. calls)
    if opts.AutoLinkingCV
        txt = auto_link_cv(opts, txt);
    end

    % replace the "helptext" part with the rendered Markdown
    div = jdoc.select('div.helptext').first();
    assert(~isempty(div), '
 not found');
    div.removeClass('helptext').addClass('helpcontent');
    div.html(txt);
end

function klass = enable_syntax_highlighting(opts, jdoc)
    %ENABLE_SYNTAX_HIGHLIGHTING  Add syntax highlighting using google-code-prettify

    % add prettyprint CSS class to pre+code blocks
    klass = 'prettyprint';
    jdoc.select('pre:has(code)').addClass(klass);

    % set language hint (only if one is not already present, which could be
    % specified in case MD backend supports fenced code blocks "```lang")
    if ~isempty(opts.PrettifyLang)
        jdoc.select('pre:has(code:not([class]))').addClass(opts.PrettifyLang);
    end

    % append code-prettify matlab theme CSS stylesheet
    url = 'https://cdn.rawgit.com/amroamroamro/prettify-matlab/master/dist/css/matlab.min.css';
    inject_css_file(jdoc, get_resource(opts, url));

    % append code-prettify JavaScript
    url = 'https://cdn.rawgit.com/google/code-prettify/master/loader/prettify.js';
    inject_js_file(jdoc, get_resource(opts, url));

    % append lang-matlab extension for code-prettify JavaScript
    url = 'https://cdn.rawgit.com/amroamroamro/prettify-matlab/master/dist/js/full/lang-matlab.min.js';
    inject_js_file(jdoc, get_resource(opts, url));

    % execute JavaScript when page is ready
    jdoc.body().attr('onload', 'PR.prettyPrint();');
end

function klass = enable_sortable_tables(opts, jdoc)
    %ENABLE_SORTABLE_TABLES  Adds ability to sort tables by clicking on headers

    % set CSS class to tables with header
    klass = 'sortable';
    jdoc.select('table:has(thead)').addClass(klass);

    % append sorttable JavaScript library
    url = 'https://cdn.jsdelivr.net/sorttable/2/sorttable.min.js';
    inject_js_file(jdoc, get_resource(opts, url));
end

function fname = get_resource(opts, url)
    %GET_RESOURCE  Return resource URI, optionally downloaded locally

    if opts.NoExternals
        % download resource
        [~, f, ext] = fileparts(url);
        fname = [f ext];
        filename = fullfile(opts.DIR, fname);
        if exist(filename, 'file') ~= 2
            urlwrite(url, filename);
        end
    else
        fname = url;
    end
end

function skip = is_handle_member(opts, topic)
    %IS_HANDLE_MEMBER  Whether to ignore processing a function

    % decide to skip if it matches the pattern "cv.Class.method",
    % and the "method" is inherited from the HANDLE class
    C = strsplit(topic, {'.', '/'});
    skip = (numel(C) >= 2 && ismember(C{end}, opts.Ignore));
end

function url = is_mathworks_function(~, topic)
    %IS_MATHWORKS_FUNCTION  Get online address for official MathWorks function

    % template for URL of online function help
    tmpl = 'https://www.mathworks.com/help/%s/ref/%s.html';

    % check if it is a builtin function
    if exist(topic, 'builtin') == 5
        url = sprintf(tmpl, 'matlab', lower(topic));
        return;
    end

    % locate function and test if it is a toolbox function
    str = which(topic);
    base = toolboxdir('');
    if ~isempty(str) && strncmp(str, base, numel(base))
        tbx = strsplit(strrep(str, base, ''), filesep());
        url = sprintf(tmpl, tbx{2}, lower(topic));
        return;
    end

    % function not found
    url = '';
end

function txt = auto_link_http(~, txt)
    %AUTO_LINK_HTTP  Autolinks "http://" text

    % http:// or https://
    % (neg lookbehind to avoid breaking existing [..](http://..) MD links/imgs)
    re = '(?<!\]\()https?://[\w&@#\/%=~|!,;:.?+-]*[\w&@#\/%=~|+-]';
    rep = '[$0]($0)';
    txt = regexprep(txt, re, rep, 'all', 'preservecase', 'lineanchors');
end

function txt = auto_link_www(~, txt)
    %AUTO_LINK_WWW  Autolinks "www." text

    % www. without http:// or https://
    re = '(^|[^/-])(www\.\S+(?:\>|$))';
    rep = '$1[$2](http://$2)';
    txt = regexprep(txt, re, rep, 'all', 'preservecase', 'lineanchors');
end

function txt = auto_link_cv(~, txt)
    %AUTO_LINK_HTTP  Autolinks cv.* functions text

    % linkify cv.func (relative URLs)
    re = '\<cv\.[\w\.]*\w+\>';
    rep = '$0';
    txt = regexprep(txt, re, rep, 'all', 'preservecase', 'lineanchors');
end

function names = enumerate_mexopencv_members(opts)
    %ENUMERATE_MEXOPENCV_MEMBERS  Enumerate all members in mexopencv

    % names of ignored methods/properties/events
    if opts.IgnoreHandleClass
        % those inherited from handle class (exclude dtor and static empty)
        mt0 = meta.class.fromName('handle');
        m0 = setdiff({mt0.MethodList.Name}, {'delete', 'empty'});
        p0 = {mt0.PropertyList.Name};
        e0 = {mt0.EventList.Name};
    else
        m0 = {};
        p0 = {};
        e0 = {};
    end

    %NOTE: we're taking all methods/props/events, including hidden and private
    %TODO: filter out hidden and private members

    mt = meta.package.fromName('cv');

    % list of functions
    f = strcat('cv.', {mt.FunctionList.Name});

    % list of classes
    c = {mt.ClassList.Name};

    % list of class methods
    m = arrayfun(@(k) strcat(k.Name, '.', setdiff({k.MethodList.Name}, m0)), ...
        mt.ClassList, 'UniformOutput',false);
    m = m(~cellfun(@isempty, m));
    m = cat(2, m{:});

    % conditionaly exclude .empty() method if inherited from handle class
    ind = arrayfun(@(k) find(strcmp({k.MethodList.Name}, 'empty'), ...
        1, 'first'), mt.ClassList);
    mask = arrayfun(@(k,i) k.MethodList(i).Static, mt.ClassList, ind);
    m1 = strcat({mt.ClassList(mask).Name}, '.empty');
    m = setdiff(m, m1);

    % list of class properties
    p = arrayfun(@(k) strcat(k.Name, '.', setdiff({k.PropertyList.Name}, p0)), ...
        mt.ClassList, 'UniformOutput',false);
    p = p(~cellfun(@isempty, p));
    p = cat(2, p{:});

    % list of class events
    e = arrayfun(@(k) strcat(k.Name, '.', setdiff({k.EventList.Name}, e0)), ...
        mt.ClassList, 'UniformOutput',false);
    e = e(~cellfun(@isempty, e));
    e = cat(2, e{:});

    % combined list
    names = unique([f(:); c(:); m(:); p(:); e(:)]);
end

function t = enumerate_mexopencv_files(~)
    %ENUMERATE_MEXOPENCV_FILES  Return a table of all mexopencv files

    % get list of all C++ sources
    paths = {
        fullfile(mexopencv.root(), 'src', '+cv', '*.cpp')
        fullfile(mexopencv.root(), 'src', '+cv', 'private', '*.cpp')
        fullfile(mexopencv.root(), 'opencv_contrib', 'src', '+cv', '*.cpp')
        fullfile(mexopencv.root(), 'opencv_contrib', 'src', '+cv', 'private', '*.cpp')
    };
    cppfiles = cellfun(@cv.glob, paths, 'UniformOutput',false);
    cppfiles = [cppfiles{:}];
    cppfiles = cppfiles(:);

    % determine which ones are from "opencv_contrib"
    iscontrib = regexp(cppfiles, '\<opencv_contrib\>', 'once');
    iscontrib = ~cellfun(@isempty, iscontrib);

    repos = {'opencv'; 'opencv_contrib'};
    repos = repos(iscontrib + 1);

    % get module name of each file
    modules = cellfun(@get_opencv_module, cppfiles, 'UniformOutput',false);

    % get base file name
    [~,names] = cellfun(@fileparts, cppfiles, 'UniformOutput',false);
    names = regexprep(names, '_$', '', 'once');

    % corresponding M-files
    paths = {
        fullfile(mexopencv.root(), '+cv')
        fullfile(mexopencv.root(), 'opencv_contrib', '+cv')
    };
    mfiles = strcat(paths(iscontrib + 1), filesep(), names, '.m');

    % corresponding test files
    paths = {
        fullfile(mexopencv.root(), 'test', 'unit_tests')
        fullfile(mexopencv.root(), 'opencv_contrib', 'test', 'unit_tests')
    };
    tfiles = regexprep(names, '^[a-z]', '${upper($0)}', 'once');
    tfiles = regexprep(tfiles, '_(\w)', '${upper($1)}', 'once');
    tfiles = strcat(paths(iscontrib + 1), filesep(), 'Test', tfiles, '.m');

    % add cv. prefix to name
    names = strcat('cv.', names);

    % create table
    t = table(names, cppfiles, mfiles, tfiles, modules, repos);
end

function module = get_opencv_module(cppfile)
    %GET_OPENCV_MODULE  Extract OpenCV module name for a mexopencv C++ source file

    module = '';
    fid = fopen(cppfile, 'rt');
    tline = fgetl(fid);
    while ischar(tline)
        if ~isempty(strfind(tline, '@ingroup'))
            module = regexp(tline, '@ingroup (\w+)$', 'tokens', 'once');
            module = module{1};
            break;
        end
        tline = fgetl(fid);
    end
    fclose(fid);
end

function inject_css_file(jdoc, url)
    %INJECT_CSS_FILE  Inject CSS stylesheet

    link = append_element(jdoc.head(), 'link');
    link.attr('rel', 'stylesheet');
    link.attr('type', 'text/css');
    link.attr('href', url);
end

function inject_css(jdoc, txt)
    %INJECT_CSS  Inject CSS code

    style = append_element(jdoc.head(), 'style');
    style.attr('type', 'text/css');
    append_data_node(style, txt);
end

function inject_js_file(jdoc, url)
    %INJECT_JS_FILE  Inject JavaScript file

    script = append_element(jdoc.body(), 'script');
    script.attr('type', 'text/javascript');
    script.attr('src', url);
end

function inject_js(jdoc, txt)
    %INJECT_JS  Inject JavaScript code

    script = append_element(jdoc.body(), 'script');
    script.attr('type', 'text/javascript');
    append_data_node(script, txt);
end

function append_data_node(node, txt)
    %APPEND_DATA_NODE  Used to set a style/script tag CSS/JavaScript text without escaping

    % append a child DataNode to specified node
    if ~isempty(node)
        node.appendChild(...
            javaObject('org.jsoup.nodes.DataNode', txt, node.baseUri()));
    end
end

function el = append_element(node, tagName)
    %APPEND_ELEMENT  Create new element by tag name and append it to node

    % create by tag name, append to node, and return newly created element
    assert(~isempty(node), 'empty node');
    el = node.appendElement(tagName);

    %HACK: insert a newline to avoid long lines (assuming prettyPrint=false)
    if true
        el.after(sprintf('\n'));
    else
        node.append(sprintf('\n'));
    end
end

function txt = xml_escape(txt)
    %XML_ESCAPE  Escape XML special characters

    txt = strrep(txt, '&', '&');
    txt = strrep(txt, '<', '<');
    txt = strrep(txt, '>', '>');
    txt = strrep(txt, '"', '"');
    txt = strrep(txt, '''', ''');
end

function filewrite(filename, str)
    %FILEWRITE  Write text to file

    fid = fopen(filename, 'wt', 'n', 'UTF-8');
    fprintf(fid, '%s', str);
    fclose(fid);
end

function outHTML = MarkdownProcess(inMD, backend)
    %MARKDOWNPROCESS  Converts Markdown text to HTML
    %
    %     outHTML = MarkdownProcess(inMD)
    %     outHTML = MarkdownProcess(inMD, backend)
    %
    % ## Input
    % * __inMD__ Input string containing Markdown source.
    %
    % ## Output
    % * __outHTML__ Output HTML formatted string.
    %
    % This function takes a a string of Markdown source and converts it
    % into a formatted HTML output string.
    %
    % [Markdown language](https://daringfireball.net/projects/markdown/) was
    % created by John Gruber and provides an easy-to-read, easy-to-write plain
    % text format that takes many cues from existing conventions for marking
    % up plain text in email.
    %
    % ## Example
    %
    %     % benchmark the different backends
    %     str = repmat(sprintf('*Hello* __world__!\n\n'), 1, 100);
    %     %web(['text://' MarkdownProcess(str)])
    %     backends = {..};
    %     t = zeros(size(backends));
    %     for i=1:numel(backends)
    %         t(i) = timeit(@() MarkdownProcess(str, backends{i}));
    %     end
    %     disp([backends num2cell(t*1000)])
    %

    % process Markdown
    switch backend
        case 'commonmark-java'
            builder = javaMethod('builder', 'org.commonmark.parser.Parser');
            parser = builder.build();
            document = parser.parse(inMD);
            builder = javaMethod('builder', ...
                'org.commonmark.renderer.html.HtmlRenderer');
            renderer = builder.build();
            outHTML = char(renderer.render(document));

        case 'MarkdownPapers'
            parser = javaObject('org.tautua.markdownpapers.Markdown');
            reader = javaObject('java.io.StringReader', inMD);
            writer = javaObject('java.io.StringWriter');
            parser.transform(reader, writer);
            outHTML = char(writer.toString());

        case 'pegdown'
            parser = javaObject('org.pegdown.PegDownProcessor');
            outHTML = char(parser.markdownToHtml(inMD));

        case 'markdownj'
            parser = javaObject('org.markdownj.MarkdownProcessor');
            outHTML = char(parser.markdown(inMD));

        case 'markdown4j'
            %NOTE: it treats newlines as real line breaks!
            parser = javaObject('org.markdown4j.Markdown4jProcessor');
            outHTML = char(parser.process(inMD));

        case 'txtmark'
            outHTML = char(javaMethod('process', ...
                'com.github.rjeschke.txtmark.Processor', inMD));

        case 'none'
            outHTML = inMD;

        otherwise
            error('Unrecognized Markdown backend: %s', backend);
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function predictions = predictOneVsAll(all_theta, X)
    %% PREDICT Predict the label for a trained one-vs-all classifier. The labels
    %are in the range 1..K, where K = size(all_theta, 1).
    %  predictions = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
    %  for each example in the matrix X. Note that X contains the examples in
    %  rows. all_theta is a matrix where the i-th row is a trained logistic
    %  regression theta vector for the i-th class. You should set p to a vector
    %  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
    %  for 4 examples)

    % Add ones to the X data matrix
    n = size(X, 1);
    X = [ones(n, 1) X];

    % ====================== YOUR CODE HERE ======================
    % Instructions: Complete the following code to make predictions using
    %               your learned logistic regression parameters (one-vs-all).
    %               You should set p to a vector of predictions (from 1 to
    %               num_labels).
    %
    % Hint: This code can be done all vectorized using the max function.
    %       In particular, the max function can also return the index of the
    %       max element, for more information see 'help max'. If your examples
    %       are in rows, then, you can use max(A, [], 2) to obtain the max
    %       for each row.
    %

    % try every example against every classifier
    activations = sigmoid(X * all_theta');

    % return both probabilities, predictions
    % but throw away probabilities for now
    [~, predictions] = max(activations, [], 2);

end

author: worldveil

repo: coursera-ml


text:
function se_hdr = GE_readHeaderSeries(fid, byte_align)
%
% se_hdr = read_series_header(fid, byte_align)
%
% Loads the series header from a file with filed id fid
% and returns it as a structure.
% if byte_align = 1 then 32-bit alignment (SGI, LX2 format)
% if byte_align = 0 then 16-bit alignment (Sun, 5.X format)
%
%
% Souheil J. Inati
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%

% define the structure and read in the data
% to overcome the byte alignment problems
% break up the assignment into pieces using the setfield function
se_hdr = struct('se_suid', fread(fid,4,'uchar'));                       %Suite ID for this Series%
se_hdr = setfield(se_hdr, 'se_uniq', fread(fid,1,'int16'));            %The Make-Unique Flag%
se_hdr = setfield(se_hdr, 'se_diskid', fread(fid,1,'uchar'));          %Disk ID for this Series%
fseek(fid,1,0); % 16-bit alignment
se_hdr = setfield(se_hdr, 'se_exno', fread(fid,1,'uint16'));            %Exam Number%
se_hdr = setfield(se_hdr, 'se_no ', fread(fid,1,'int16'));              %Series Number%
se_hdr = setfield(se_hdr, 'se_datetime', fread(fid,1,'int32'));        %Allocation Series Data/Time stamp%
se_hdr = setfield(se_hdr, 'se_actual_dt', fread(fid,1,'int32'));       %Actual Series Data/Time stamp%
se_hdr = setfield(se_hdr, 'se_desc', fread(fid,30,'uchar'));       %Series Description%
se_hdr = setfield(se_hdr, 'pr_sysid', fread(fid,9,'uchar'));       %Primary Receiver Suite and Host%
se_hdr = setfield(se_hdr, 'pansysid', fread(fid,9,'uchar'));       %Archiver Suite and Host%
se_hdr = setfield(se_hdr, 'se_typ ', fread(fid,1,'int16'));             %Series Type%
se_hdr = setfield(se_hdr, 'se_source ', fread(fid,1,'int16'));          %Series from which prescribed%
se_hdr = setfield(se_hdr, 'se_plane ', fread(fid,1,'int16'));           %Most-like Plane (for L/S)%
se_hdr = setfield(se_hdr, 'scan_type ', fread(fid,1,'int16'));          %Scout or Axial (for CT)%
se_hdr = setfield(se_hdr, 'position', fread(fid,1,'int32'));           %Patient Position%
se_hdr = setfield(se_hdr, 'entry', fread(fid,1,'int32'));              %Patient Entry%
se_hdr = setfield(se_hdr, 'anref', fread(fid,3,'uchar'));          %Anatomical reference%
fseek(fid,1,0); % 16-bit alignment
se_hdr = setfield(se_hdr, 'lmhor', fread(fid,1,'float32'));              %Horizontal Landmark%
se_hdr = setfield(se_hdr, 'prtcl', fread(fid,25,'uchar'));         %Scan Protocol Name%
fseek(fid,1,0); % 16-bit alignment
se_hdr = setfield(se_hdr, 'se_contrast ', fread(fid,1,'int16'));        %Non-zero if > 0 image used contrast(L/S)%
se_hdr = setfield(se_hdr, 'start_ras', fread(fid,1,'uchar'));          %RAS letter for first scan location (L/S)%
fseek(fid,1,0); % 16-bit alignment
if byte_align; fseek(fid,2,0); end % 32-bit alignment
se_hdr = setfield(se_hdr, 'start_loc', fread(fid,1,'float32'));          %First scan location (L/S)%
se_hdr = setfield(se_hdr, 'end_ras', fread(fid,1,'uchar'));            %RAS letter for last scan location (L/S)%
fseek(fid,1,0); % 16-bit alignment
if byte_align; fseek(fid,2,0); end % 32-bit alignment
se_hdr = setfield(se_hdr, 'end_loc', fread(fid,1,'float32'));            %Last scan location (L/S)%
se_hdr = setfield(se_hdr, 'se_pseq ', fread(fid,1,'int16'));            %Last Pulse Sequence Used (L/S)%
se_hdr = setfield(se_hdr, 'se_sortorder ', fread(fid,1,'int16'));       %Image Sort Order (L/S)%
se_hdr = setfield(se_hdr, 'se_lndmrkcnt', fread(fid,1,'int32'));       %Landmark Counter%
se_hdr = setfield(se_hdr, 'se_nacq ', fread(fid,1,'int16'));            %Number of Acquisitions%
se_hdr = setfield(se_hdr, 'xbasest ', fread(fid,1,'int16'));            %Starting number for baselines%
se_hdr = setfield(se_hdr, 'xbaseend', fread(fid,1,'int16'));            %Ending number for baselines%
se_hdr = setfield(se_hdr, 'xenhst', fread(fid,1,'int16'));             %Starting number for enhanced scans%
se_hdr = setfield(se_hdr, 'xenhend', fread(fid,1,'int16'));            %Ending number for enhanced scans%
if byte_align; fseek(fid,2,0); end % 32-bit alignment
se_hdr = setfield(se_hdr, 'se_lastmod', fread(fid,1,'int32'));         %Date/Time of Last Change%
se_hdr = setfield(se_hdr, 'se_alloc_key', fread(fid,13,'uchar'));  %Process that allocated this record%
fseek(fid,1,0); % 16-bit alignment
if byte_align; fseek(fid,2,0); end % 32-bit alignment
se_hdr = setfield(se_hdr, 'se_delta_cnt', fread(fid,1,'int32'));       %Indicates number of updates to header%
se_hdr = setfield(se_hdr, 'se_verscre', fread(fid,2,'uchar'));     %Genesis Version - Created%
se_hdr = setfield(se_hdr, 'se_verscur', fread(fid,2,'uchar'));     %Genesis Version - Now%
se_hdr = setfield(se_hdr, 'se_pds_a', fread(fid,1,'float32'));           %PixelData size - as stored%
se_hdr = setfield(se_hdr, 'se_pds_c', fread(fid,1,'float32'));           %PixelData size - Compressed%
se_hdr = setfield(se_hdr, 'se_pds_u', fread(fid,1,'float32'));           %PixelData size - UnCompressed%
se_hdr = setfield(se_hdr, 'se_checksum', fread(fid,1,'uint32'));        %Series Record checksum%
se_hdr = setfield(se_hdr, 'se_complete', fread(fid,1,'int32'));        %Series Complete Flag%
se_hdr = setfield(se_hdr, 'se_numarch', fread(fid,1,'int32'));         %Number of Images Archived%
se_hdr = setfield(se_hdr, 'se_imagect', fread(fid,1,'int32'));         %Last Image Number Used%
se_hdr = setfield(se_hdr, 'se_numimages', fread(fid,1,'int32'));       %Number of Images Existing%
se_hdr = setfield(se_hdr, 'se_images', struct('length', fread(fid,1,'uint32'), ...
                                               'data', fread(fid,1,'uint32'))); %Image Keys for this Series%
se_hdr = setfield(se_hdr, 'se_numunimg', fread(fid,1,'int32'));        %Number of Unstored Images%
se_hdr = setfield(se_hdr, 'se_unimages', struct('length', fread(fid,1,'uint32'), ...
                                               'data', fread(fid,1,'uint32'))); %Unstored Image Keys for this Series%
se_hdr = setfield(se_hdr, 'se_toarchcnt', fread(fid,1,'int32'));       %Number of Unarchived Images%
se_hdr = setfield(se_hdr, 'se_toarchive', struct('length', fread(fid,1,'uint32'), ...
                                               'data', fread(fid,1,'uint32'))); %Unarchived Image Keys for this Series%
se_hdr = setfield(se_hdr, 'echo1_alpha', fread(fid,1,'float32'));        %Echo 1 Alpha Value%
se_hdr = setfield(se_hdr, 'echo1_beta', fread(fid,1,'float32'));         %Echo 1 Beta Value%
se_hdr = setfield(se_hdr, 'echo1_window', fread(fid,1,'uint16'));       %Echo 1 Window Value%
se_hdr = setfield(se_hdr, 'echo1_level', fread(fid,1,'int16'));        %Echo 1 Level Value%
se_hdr = setfield(se_hdr, 'echo2_alpha', fread(fid,1,'float32'));        %Echo 2 Alpha Value%
se_hdr = setfield(se_hdr, 'echo2_beta', fread(fid,1,'float32'));         %Echo 2 Beta Value%
se_hdr = setfield(se_hdr, 'echo2_window', fread(fid,1,'uint16'));       %Echo 2 Window Value%
se_hdr = setfield(se_hdr, 'echo2_level', fread(fid,1,'int16'));        %Echo 2 Level Value%
se_hdr = setfield(se_hdr, 'echo3_alpha', fread(fid,1,'float32'));        %Echo 3 Alpha Value%
se_hdr = setfield(se_hdr, 'echo3_beta', fread(fid,1,'float32'));         %Echo 3 Beta Value%
se_hdr = setfield(se_hdr, 'echo3_window', fread(fid,1,'uint16'));       %Echo 3 Window Value%
se_hdr = setfield(se_hdr, 'echo3_level', fread(fid,1,'int16'));        %Echo 3 Level Value%
se_hdr = setfield(se_hdr, 'echo4_alpha', fread(fid,1,'float32'));        %Echo 4 Alpha Value%
se_hdr = setfield(se_hdr, 'echo4_beta', fread(fid,1,'float32'));         %Echo 4 Beta Value%
se_hdr = setfield(se_hdr, 'echo4_window', fread(fid,1,'uint16'));       %Echo 4 Window Value%
se_hdr = setfield(se_hdr, 'echo4_level', fread(fid,1,'int16'));        %Echo 4 Level Value%
se_hdr = setfield(se_hdr, 'echo5_alpha', fread(fid,1,'float32'));        %Echo 5 Alpha Value%
se_hdr = setfield(se_hdr, 'echo5_beta', fread(fid,1,'float32'));         %Echo 5 Beta Value%
se_hdr = setfield(se_hdr, 'echo5_window', fread(fid,1,'uint16'));       %Echo 5 Window Value%
se_hdr = setfield(se_hdr, 'echo5_level', fread(fid,1,'int16'));        %Echo 5 Level Value%
se_hdr = setfield(se_hdr, 'echo6_alpha', fread(fid,1,'float32'));        %Echo 6 Alpha Value%
se_hdr = setfield(se_hdr, 'echo6_beta', fread(fid,1,'float32'));         %Echo 6 Beta Value%
se_hdr = setfield(se_hdr, 'echo6_window', fread(fid,1,'uint16'));       %Echo 6 Window Value%
se_hdr = setfield(se_hdr, 'echo6_level', fread(fid,1,'int16'));        %Echo 6 Level Value%
se_hdr = setfield(se_hdr, 'echo7_alpha', fread(fid,1,'float32'));        %Echo 7 Alpha Value%
se_hdr = setfield(se_hdr, 'echo7_beta', fread(fid,1,'float32'));         %Echo 7 Beta Value%
se_hdr = setfield(se_hdr, 'echo7_window', fread(fid,1,'uint16'));       %Echo 7 Window Value%
se_hdr = setfield(se_hdr, 'echo7_level', fread(fid,1,'int16'));        %Echo 7 Level Value%
se_hdr = setfield(se_hdr, 'echo8_alpha', fread(fid,1,'float32'));        %Echo 8 Alpha Value%
se_hdr = setfield(se_hdr, 'echo8_beta', fread(fid,1,'float32'));         %Echo 8 Beta Value%
se_hdr = setfield(se_hdr, 'echo8_window', fread(fid,1,'uint16'));       %Echo 8 Window Value%
se_hdr = setfield(se_hdr, 'echo8_level', fread(fid,1,'int16'));        %Echo 8 Level Value%
se_hdr = setfield(se_hdr, 'series_uid', fread(fid,32,'uchar'));    %Series Entity Unique ID%
se_hdr = setfield(se_hdr, 'landmark_uid', fread(fid,32,'uchar'));  %Landmark Unique ID%
se_hdr = setfield(se_hdr, 'equipmnt_uid', fread(fid,32,'uchar'));  %Equipment Unique ID%
se_hdr = setfield(se_hdr, 'se_padding', fread(fid,588,'uchar'));   %Spare Space%

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [eigvector, eigvalue] = IsoP(options, data)
% IsoP: Isometric Projection
%
%       [eigvector, eigvalue] = IsoP(options, data)
%
%             Input:
%               data    - Data matrix. Each row vector of data is a data point.
%
%               options - Struct value in Matlab. The fields in options
%                         that can be set:
%
%                      NeighborMode -  Indicates how to construct the graph. Choices
%                           are:
%                           'KNN'     -  Put an edge between two nodes if and
%                                        only if they are among the k nearst
%                                        neighbors of each other. Default
%                                        option.
%                       'Supervised'  -  Two variations:
%                                       1. k=0, Put an edge between two nodes
%                                          if and only if they belong to
%                                          same class.
%                                       2. k>0, The distance between two nodes
%                                          in the same class will be smaller than
%                                          two nodes have diff. labels
%                                          The label information 'gnd' should be
%                                          provided.
%
%                       k           -   The number of neighbors.
%                                       Default k = 5;
%                       gnd         -   The parameter needed under 'Supervised'
%                                       NeighborMode.  Colunm vector of the label
%                                       information for each data point.
%
%                         Please see LGE.m for other options.
%
%
%             Output:
%               eigvector - Each column is an embedding function, for a new
%                           data point (row vector) x,  y = x*eigvector
%                           will be the embedding result of x.
%               eigvalue  - The eigvalue of LPP eigen-problem. sorted from
%                           smallest to largest.
%               elapse    - Time spent on different steps
%
%
%    Examples:
%
%
%
%       fea = rand(50,70);
%       gnd = [ones(10,1);ones(15,1)*2;ones(10,1)*3;ones(15,1)*4];
%       options = [];
%       options.k = 0;
%       options.NeighborMode = 'Supervised';
%       options.gnd = gnd;
%       [eigvector, eigvalue] = IsoP(options, fea);
%       Y = fea*eigvector;
%
%
%
% See also LPP, LGE
%
%Reference:
%
%   Deng Cai, Xiaofei He, and Jiawei Han, "Isometric Projection",
%   Twenty-Second Conference on Artificial Intelligence (AAAI-07), 2007
%
%   Deng Cai, Xiaofei He and Jiawei Han, "Isometric Projection", Technical
%   report, Computer Science Department, UIUC, UIUCDCS-R-2006-2747, July 2006
%
%   Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. "A Global
%   Geometric Framework for Nonlinear Dimensionality Reduction", Science,
%   v.290 no.5500 , Dec.22, 2000. pp.2319-2323.
%
%
%   version 2.1 --June/2007
%   version 2.0 --May/2007
%   version 1.1 --May/2006
%   version 1.0 --Nov/2005
%
%   Written by Deng Cai (dengcai2 AT cs.uiuc.edu)

INFratio = 1000;

if (~exist('options','var'))
   options = [];
end

if ~isfield(options,'NeighborMode')
    options.NeighborMode = 'KNN';
end

if ~isfield(options,'k')
    options.k = 5;
end

nSmp = size(data,1);

if options.k >= nSmp
    error('k is too large!');
end

if options.k <= 0  % Always supervised!
    if ~isfield(options,'gnd')
        error('gnd should be provided!');
    end
    if length(options.gnd) ~= nSmp
        error('gnd and data mismatch!');
    end

    Label = unique(options.gnd);
    nLabel = length(Label);

    G = zeros(nSmp,nSmp);
    for i=1:nLabel
        classIdx = find(options.gnd==Label(i));
        D = EuDist2(data(classIdx,:),[],1);
        G(classIdx,classIdx) = D;
    end
    maxD = max(max(G));
    INF = maxD*INFratio;  % effectively infinite distance

    D = INF*ones(nSmp,nSmp);
    for i=1:nLabel
        classIdx = find(options.gnd==Label(i));
        D(classIdx,classIdx) = G(classIdx,classIdx);
    end

    clear G
else
    switch lower(options.NeighborMode)
        case {lower('KNN')}
            D = EuDist2(data);
            maxD = max(max(D));
            INF = maxD*INFratio;  % effectively infinite distance

            [dump,iidx] = sort(D,2);
            iidx = iidx(:,(2+options.k):end);
            for i=1:nSmp
                D(i,iidx(i,:)) = 0;
            end
            D = max(D,D');

            D = sparse(D);
            D = dijkstra(D, 1:nSmp);

            D = reshape(D,nSmp*nSmp,1);
            infIdx = find(D==inf);
            if ~isempty(infIdx)
                D(infIdx) = INF;
            end
            D = reshape(D,nSmp,nSmp);

        case {lower('Supervised')}
            if ~isfield(options,'gnd')
                error('gnd should be provided!');
            end
            if length(options.gnd) ~= nSmp
                error('gnd and data mismatch!');
            end

            Label = unique(options.gnd);
            nLabel = length(Label);

            G = zeros(nSmp,nSmp);
            maxD = 0;
            for idx=1:nLabel
                classIdx = find(options.gnd==Label(idx));
                nSmpClass = length(classIdx);
                D = EuDist2(data(classIdx,:),[],1);
                if maxD < max(max(D))
                    maxD = max(max(D));
                end
                if options.k >= nSmpClass
                    G(classIdx,classIdx) = D;
                else
                    [dump,iidx] = sort(D,2);
                    iidx = iidx(:,(2+options.k):end);
                    for i=1:nSmpClass
                        D(i,iidx(i,:)) = 0;
                    end
                    D = max(D,D');
                    D = sparse(D);
                    D = dijkstra(D, 1:nSmpClass);
                    G(classIdx,classIdx) = D;
                end
            end

            INF = maxD*INFratio;  % effectively infinite distance

            D = INF*ones(nSmp,nSmp);
            for i=1:nLabel
                classIdx = find(options.gnd==Label(i));
                D(classIdx,classIdx) = G(classIdx,classIdx);
            end
            clear G

        otherwise
            error('NeighborMode does not exist!');
    end
end

S = D.^2;
sumS = sum(S);
H = sumS'*ones(1,nSmp)/nSmp;
TauDg = -.5*(S - H - H' + sum(sumS)/(nSmp^2));

TauDg = max(TauDg,TauDg');

%==========================
% If data is too large, the following centering codes can be commented
%==========================
if isfield(options,'keepMean') && options.keepMean
else
    if issparse(data)
        data = full(data);
    end
    sampleMean = mean(data);
    data = (data - repmat(sampleMean,nSmp,1));
end
%==========================

[eigvector, eigvalue] = LGE(TauDg, [], options, data);

eigIdx = find(eigvalue < 1e-3);
eigvalue (eigIdx) = [];
eigvector(:,eigIdx) = [];

author: ZJULearning

repo: MatlabFunc


text:
function d = mahalanobis(varargin)
%MAHALANOBIS Computes the Mahalanobis distance.
%   D = MAHALANOBIS(Y, X) computes the Mahalanobis distance between
%   each vector in Y to the mean (centroid) of the vectors in X, and
%   outputs the result in vector D, whose length is size(Y, 1).  The
%   vectors in X and Y are assumed to be organized as rows.  The
%   input data can be real or complex. The outputs are real
%   quantities.
%
%   D = MAHALANOBIS(Y, CX, MX) computes the Mahalanobis distance
%   between each vector in Y and the given mean vector, MX. The
%   results are output in vector D, whose length is size(Y, 1).  The
%   vectors in Y are assumed to be organized as the rows of this
%   array. The input data can be real or complex. The outputs are
%   real quantities. In addition to the mean vector MX, the
%   covariance matrix CX of a population of vectors X also must be
%   provided. Use function COVMATRIX (Section 11.5) to compute MX and
%   CX.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2005/01/18 13:44:47 $

% Reference: Acklam, P. J. [2002]. "MATLAB Array Manipulation Tips
% and Tricks." Available at
%     home.online.no/~pjacklam/matlab/doc/mtt/index.html
% or at
%     www.prenhall.com/gonzalezwoodseddins

param = varargin; % Keep in mind that param is a cell array.
Y = param{1};
ny = size(Y, 1); % Number of vectors in Y.

if length(param) == 2
   X = param{2};
   % Compute the mean vector and covariance matrix of the vectors
   % in X.
   [Cx, mx] = covmatrix(X);
elseif length(param) == 3 % Cov. matrix and mean vector provided.
   Cx = param{2};
   mx = param{3};
else
   error('Wrong number of inputs.')
end
mx = mx(:)'; % Make sure that mx is a row vector.

% Subtract the mean vector from each vector in Y.
Yc = Y - mx(ones(ny, 1), :);

% Compute the Mahalanobis distances.
d = real(sum(Yc/Cx.*conj(Yc), 2));

author: Ultrasty

repo: Digital-Image-Processing


text:
clear all, close all, clc

xlim([-2 2])
ylim([-2 2])
zlim([-2 2])
% view([-72 18])
view([45 26])
% axis off

set(gca,'FontSize',13);
set(gcf,'Position',[100 100 300 300])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', 'CH01_EX09_1_axis.eps')

%% Define rotation matrix
t1 = pi/15;
t2 = -pi/9;
t3 = -pi/20;

Rx = [1 0 0;
    0 cos(t1) -sin(t1);
    0 sin(t1) cos(t1)];

Ry = [cos(t2) 0 sin(t2);
    0 1 0;
    -sin(t2) 0 cos(t2)];

Rz = [cos(t3) -sin(t3) 0;
    sin(t3) cos(t3) 0;
    0 0 1];

Sigma = diag([3; 1; 0.5]);

R = Rz*Ry*Rx*Sigma;

%% Plot sphere and great circles
[x,y,z] = sphere(250);
% subplot(1,2,1)
plot3([-2 2],[0 0],[0 0],'k','LineWidth',6);
hold on
plot3([0 0],[-2 2],[0 0],'k','LineWidth',6);
plot3([0 0],[0 0],[-2 2],'k','LineWidth',6);

h1=surf(x,y,z);
hold on
set(h1,'EdgeColor','none','FaceColor',[.5 .5 .5],'FaceAlpha',.75)
% set(h1,'FaceColor','none')
lighting phong
shading interp
axis equal

theta = (0:.001:1)*2*pi;
x1 = cos(theta);
y1 = sin(theta);
z1 = 0*theta;
x2 = 0*theta;
y2 = cos(theta);
z2 = sin(theta);
x3 = cos(theta);
y3 = 0*theta;
z3 = sin(theta);

plot3(x1,y1,z1,'k','LineWidth',6);
plot3(x2,y2,z2,'k','LineWidth',6);
plot3(x3,y3,z3,'k','LineWidth',6);

xlim([-2 2])
ylim([-2 2])
zlim([-2 2])
% view([-72 18])
view([45 26])
% axis off

set(gca,'FontSize',13);
set(gcf,'Position',[100 100 100 100])
set(gcf,'PaperPositionMode','auto')
set(gcf,'renderer','opengl')

axis off
%   save2pdf('CH01_EX09_1b.pdf',gcf,600)
% print('-depsc2', '-loose', 'CH01_EX09_1b.eps')
print('-dpng','-r1800' ,'-loose', 'CH01_EX09_1a.png')

%%
xR = 0*x;
yR = 0*y;
zR = 0*z;

for i=1:size(x,1)
    for j=1:size(x,2)
        vec = [x(i,j); y(i,j); z(i,j)];
        vecR = R*vec;
        xR(i,j) = vecR(1);
        yR(i,j) = vecR(2);
        zR(i,j) = vecR(3);
    end
end

vec1 = [x1; y1; z1];
vec2 = [x2; y2; z2];
vec3 = [x3; y3; z3];

vec1R = R*vec1;
vec2R = R*vec2;
vec3R = R*vec3;

eX = [2; 0; 0];
eY = [0; 2; 0];
eZ = [0; 0; 2];
eXR = R*eX;
eYR = R*eY;
eZR = R*eZ;

x1R = vec1R(1,:);
y1R = vec1R(2,:);
z1R = vec1R(3,:);
x2R = vec2R(1,:);
y2R = vec2R(2,:);
z2R = vec2R(3,:);
x3R = vec3R(1,:);
y3R = vec3R(2,:);
z3R = vec3R(3,:);

figure
% plot3([0 -2],[0 0],[0 0],'Color',[.3 .3 .3],'LineWidth',6);
% hold on
% plot3([0 0],[0 -2],[0 0],'Color',[.3 .3 .3],'LineWidth',6);
% plot3([0 0],[0 0],[0 2],'Color',[.3 .3 .3],'LineWidth',6);
h1=surf(xR,yR,zR,z);
hold on
set(h1,'EdgeColor','none','FaceColor',[.5 .5 .5],'FaceAlpha',.75)
% set(h1,'FaceColor','none')
axis equal
shading interp
lighting phong

plot3(x1R,y1R,z1R,'k','LineWidth',6);
plot3(x2R,y2R,z2R,'k','LineWidth',6);
plot3(x3R,y3R,z3R,'k','LineWidth',6);

plot3([-eXR(1) eXR(1)],[-eXR(2) eXR(2)],[-eXR(3) eXR(3)],'k','LineWidth',6);
plot3([-eYR(1) eYR(1)],[-eYR(2) eYR(2)],[-eYR(3) eYR(3)],'k','LineWidth',6);
plot3([-eZR(1) eZR(1)],[-eZR(2) eZR(2)],[-eZR(3) eZR(3)],'k','LineWidth',6);

xlim([-2 2])
ylim([-2 2])
zlim([-2 2])
grid off

% view([-72 18])
view([45 26])
% axis off

set(gca,'FontSize',13);
set(gcf,'Position',[100 100 100 100])
set(gcf,'PaperPositionMode','auto')
set(gcf,'renderer','opengl')

axis off
%   save2pdf('CH01_EX09_1b.pdf',gcf,600)
% print('-depsc2', '-loose', 'CH01_EX09_1b.eps')
print('-dpng','-r1800' ,'-loose', 'CH01_EX09_1b.png')

author: dynamicslab


text:
function varargout = log10(varargin)
%log10 (overloaded)

varargout{1} = log(varargin{1})/log(10);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function I = intersectionHull(varargin)
%intersectionHull - computes a bounded convex polyhedron resulting from the
%intersection of other convex polyhedra. The other polyhedra need not individually
%be bounded. Only the final intersection resulting from them must be bounded.
%Any number of polyhedra can be input to the intersection operation. They can
%be expressed using either vertices or linear in/equalities, according to
%the input scheme,
%
%
%   I = intersectionHull('vert', V1, 'lcon', A2,b2, 'lcon', A3,b3,Aeq3,,beq3,...)
%
%The arguments specifying different polyhedra are separated using labels
%'vert' and 'lcon'. The label 'vert' signifies that an input polyhedron will
%be expressed using vertices and is to be followed by any string of input arguments
%accepted by vert2lcon(). The label 'lcon' signifies that an input polyhedron will
%be expressed using linear constraints and is to be followed by any string of
%input arguments accepted by lcon2vert().
%
%The output, I, is a struct containing fields
%
%   I.vert: A matrix whose rows are the vertices of the polyhedron formed from
%           the intersection.
%   I.lcon: The quadruplet of linear constraint data {A,b,Aeq,beq}
%           describing the polyhedral intersection.
%
%EXAMPLE 1: This example computes the intersection of a unit square and an
%oblique 2D line segment, both expressed in terms of their vertices.
%
%     V1=dec2bin(0:2^2-1,2)-'0';   %vertices of unit square
%
%     V2=[1,1;0,-1];               %vertices of 2D line segment
%
%     I=intersectionHull('vert',V1,'vert',V2);  %compute intersection
%
%The intersection is another line segment with vertices
%
%        >> I.vert
%
%         ans =
%
%             0.5000         0
%             1.0000    1.0000
%
%EXAMPLE 2: This example computes the intersection of a unit cube, expressed in
%terms of its vertices, and an infinite oblique 3D line, expressed in terms of linear equalities.
%Note that the line is an unbounded polyhedron. This is okay, since we know in advance
%that the final polyhedron formed from the intersection is bounded.
%
%     V=dec2bin(0:2^3-1,3)-'0';    %vertices of unit cube
%
%     Aeq=[1 -1 0; 0 1 -1]; beq=[0;0]; %oblique line in 3D
%
%     I=intersectionHull('vert',V,'lcon',[],[],Aeq,beq);  %compute intersection
%
%Once again, the intersection is a line segment. Its vertices are
%
%      >> I.vert   %vertices of line segment of intersection
%
%         ans =
%
%             0.0000    0.0000    0.0000
%             1.0000    1.0000    1.0000

%%%%begin parsing

if isnumeric(varargin{1})
   TOL=varargin{1};
   varargin(1)=[];
else
    TOL=[];
end

N=length(varargin);
idxType = [find(cellfun(@ischar,varargin)),N+1];

L=length(idxType)-1;
S(L).type=[];
S(L).args={};
S(L).A=[];
S(L).b=[];
S(L).Aeq=[];
S(L).beq=[];

for i=1:L

    j=idxType(i);
    k=idxType(i+1);
    S(i).type=varargin{j};
    S(i).args=varargin(j+1:k-1);

    if isempty(S(i).args)
     error 'Syntax error - arguments missing'
    end

    lcon=cell(1,4);

      switch S(i).type

          case 'vert'

              [lcon{1:4}] = vert2lcon(S(i).args{:});

          case 'lcon'

               lcon(1:k-j-1) = S(i).args;

          case 'qlcon' %deliberately undocumented - no point in using this

               lcon(1:k-j-1) = S(i).args(2:end);

          otherwise

            error(['Unrecognized representation label of polyhedron ' num2str(i)]);

      end

    [S(i).A, S(i).b, S(i).Aeq, S(i).beq] = deal(lcon{:});

end

%%%%end parsing

   A=vertcat(S.A);
   b=vertcat(S.b);
   Aeq=vertcat(S.Aeq);
   beq=vertcat(S.beq);

   [V,nr,nre]=lcon2vert(A,b,Aeq,beq,TOL);

   I.vert=V;
   I.lcon={A(nr,:),b(nr,:), Aeq(nre,:),beq(nre,:)};

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [y,newp,retcode]=sstate_file(obj,y,p,d,id) %#ok<INUSL>
% sstate_model -- shows the way of writing a RISE steady state file
%
% ::
%
%
%   [y,newp,retcode]=sstate_model(obj,y,p,d,id)
%
% Args:
%
%    - **obj** [rise|dsge]: model object (not always needed)
%
%    - **y** [vector]: endo_nbr x 1 vector of initial steady state
%
%    - **p** [struct]: parameter structure
%
%    - **d** [struct]: definitions
%
%    - **id** [vector]: location of the variables to calculate
%
% Returns:
%    :
%
%      CASE 1: one input argument
%
%    - **y** [cell array]: list of the variables for which the steady state
%    will be calculated within the steady state function
%
%    - **newp** [cell array]: List of the parameters calculated inside the
%    steady state function
%
%    - **retcode** [struct]: possible fields are "imposed", "unique", "loop".
%    The default value for all of those is false.
%      - "imposed": This tells RISE not to check that this is actually solves
%          the steady state. Hence, RISE will attempt to approximate the model
%          around the chosen point
%      - "unique": This tells RISE that the steady state is the same across
%          all regimes. RISE will call the function only once but instead of
%          using just any parameter vector, it will use the ergodic mean of
%          the parameter vector (mean across all regimes).
%      - "loop": This tells RISE that if the user was given the steady state
%          values for some of the variables, he would be able to compute the
%          steady state for the remaining variables. RISE will then exploit
%          this information to optimize over the variables that the user needs
%          for computing the steady state.
%
%      CASE 2: More than one input argument
%
%    - **y** []: endo_nbr x 1 vector of updated steady state
%
%    - **newp** [struct]: structure containing updated parameters if any
%
%    - **retcode** [0|number]: return 0 if there are no problems, else return
%      any number different from 0
%
% Note:
%
%    - If the user knows the steady state, it is always an advantage. If the
%    steady state is computed numerically, we don't know whether it is unique
%    or not. Not that it really matters but... some economists have a strong
%    aversion towards models with multiple equilibria.
%
%    - If the user does not know the solution for all the variables in the
%    steady state, it is a good idea to take a log-linear approximation for
%    the variables that potentially have a nonzero steady state. Hence the
%    user should give that information to RISE.
%
%    - One can potentially improve on the above point by explicit bounds on
%    the variables. But this is not implemented.
%
%    - An alternative that potentially avoids taking a loglinearization is to
%    to reset the values proposed by the optimizer whenever they are in a bad
%    region. It is unclear whether this always works.
%
%    - So be on the safe side, i.e. don't do like me: compute your steady
%    state analytically.
%
% Example:
%
%    See also:

retcode=0;
if nargin==1

    % list the variables for which we compute the steady state. Not all
    % variables need to be listed. The variables NOT in this list
    % automatically get a value of 0. RISE is responsible for calculating
    % the steady state of any variable it may have created. The user does
    % not have to worry about those.
    %----------------------------------------------------------------------
    y={'X','A','E','Z','V','MU','PAI','R','Q','LAMBDA','XI','C','M','Y',...
        'K','I','H','W','D','N','LC','LI','LM','LPI','LR'};

	% list of parameters herein calculated
	%-------------------------------------
	newp={};

    % flags on the calculation
    %--------------------------
    retcode=struct('unique',true,'imposed',true);
else
    % no parameters to update
    %-------------------------
    newp=[];

    % create definitions (auxiliary parameters) that are useful for solving
    % the steady state
    %----------------------------------------------------------------------
    z_ss=p.z_ss_trans*10000;
% 	phi_p = 100*abs(p.phi_p_trans);
% 	phi_k = 100*abs(p.phi_k_trans);
	mu_ss = 1 + p.mu_ss_trans;

    X=1;
    A=1;
    E=p.e_ss;
    Z=z_ss;
    V=1;
    % mu determined by policy
    MU=mu_ss;
    PAI=MU/p.g;
    R=MU/p.beta;
    Q=p.g/p.beta-1+p.delta;
    % lambda comes here
    l1ss=1/(p.theta/(p.theta-1)-(p.g-1+p.delta)*p.alpha/Q);
    l2ss=1/(1+E*(R/(R-1))^(p.gam-1));
    l3ss=(1-p.alpha)*Z*((p.theta-1)/p.theta)^(1/(1-p.alpha))*(p.alpha/Q)^(p.alpha/(1-p.alpha));
    LAMBDA=(p.eta+(1-p.alpha)*l1ss*l2ss)/l3ss;
    XI = (p.theta-1)/p.theta*LAMBDA;
    C = 1/(1+E*(R/(R-1))^(p.gam-1))*(1/LAMBDA);
    M = E*(R/(R-1))^p.gam*C;
    Y = 1/(1-(p.g-1+p.delta)*(p.theta-1)/p.theta*p.alpha/Q)*C;
    K = (p.theta-1)/p.theta*p.alpha*Y/Q;
    I = (p.g-1+p.delta)*K;
    H=1/Z*(Y/K^p.alpha)^(1/(1-p.alpha));
    W=(1-p.alpha)*(p.theta-1)/p.theta*Y/H;
    D = Y-W*H-Q*K;
    N=Y/H;
    LC=log(C);
    LI=log(I);
    LM=log(M);
    LPI=log(PAI);
    LR=log(R);
    ys=[X A E Z V MU PAI R Q LAMBDA XI C M Y K I H W D N LC LI LM LPI LR].';

    % check the validity of the calculations
    %----------------------------------------
    if ~utils.error.valid(ys)||any(abs(ys)>100000)
        retcode=1;
    else
        % push the calculations
        %----------------------
        y(id)=ys;
    end
end

end

author: jmaih


text:
%NUMROWS Number of rows in matrix
%
% NR = NUMROWS(M) is the number of rows in the matrix M.
%
% Notes::
% - Readable shorthand for SIZE(M,1);
%
% See also NUMCOLS, SIZE.

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

function r = numrows(m)

	r = size(m, 1);

author: petercorke

repo: spatialmath-matlab


text:
function jed = epoch_to_jed_rd ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_RD returns the epoch of the RD calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1721425.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = huber_circ( x, dim, varargin )

%HUBER_CIRC   Huber penalty function with circular symmetry.
%   For a vector X, HUBER_CIRC(X) computes the Huber penalty function
%
%       HUBER_CIRC(X) =   NORM(X,2)^2 if NORM(X,2)<=1,
%                       2*NORM(X,2)-1 if NORM(X,2)>=1.
%
%   For matrices and N-D arrays, the penalty function is applied to the
%   first dimens
%
%   HUBER_CIRC(X,[],M) computes the penalty function with halfwidth M,
%   M.^2.*HUBER_CIRC(X./M). M must be real and positive.
%
%   HUBER_CIRC(X,[],M,T) computes the penalty function with halfwidth M
%   and concomitant scale T:
%
%       HUBER_CIRC(X,[],M,T) = T.*HUBER_CIRC(X./T,[],M) if T > 0
%                              +Inf                     if T <= 0
%
%   See the help file for HUBER for information about this usage.
%
%   If X is a matrix, the penalty function is applied to the columns of the
%   matrix X, and a row vector is returned. If X is an N-D matrix, the
%   penalties are computed along the first non-singleton dimension.
%
%   HUBER_CIRC(X,DIM), HUBER_CIRC(X,DIM,M), and HUBER_CIRC(X,DIM,M,T)
%   computes the penalty along the dimension DIM.
%
%   Disciplined convex programming information:
%       HUBER_CIRC is jointly convex in X and T. It is nonomonotonic in X
%       and nonincreasing in T. Therefore, when used in CVX specifications,
%       X must be affine and T must be concave (or affine). T must be real.
%       X, on the other hand, may be real or complex.

if nargin < 2, dim = []; end
y = huber_pos( norms( x, 2, dim ), varargin{:} );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function cc_grids_animate ( )

%*****************************************************************************80
%
%% CC_GRIDS_ANIMATE displays a sequence of merged 2D Clenshaw-Curtis grids.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_GRIDS_ANIMATE:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display a sequence of merged 2D Clenshaw-Curtis grids.\n' );

  dim_num = 2;

  while ( 1 )

    fprintf ( 1, '\n' );
    order_max_1d = input ( 'Enter the maximum orders [A,B] or RETURN to exit;' );

    if ( isempty ( order_max_1d ) )
      break
    end

    points_old = [];
    points_old_num = 0;

    points = [];
    order_nd = 0;

    for order = 1 : max ( order_max_1d(1:2) )

      if ( 0 < order_nd )
        points_old(1:2,points_old_num+1:points_old_num+order_nd) = points(1:2,1:order_nd);
        points_old_num = points_old_num + order_nd;
      end

      clf
%
%  We have to name the axes in order to control the grid.
%
      axes_handle = axes;
%
%  Compute the new points and plot them in red.
%
      order_1d(1:2) = min ( order, order_max_1d(1:2) );
      order_nd = prod ( order_1d(1:2) );

      points = cc_grid ( dim_num, order_1d, order_nd );

      handle_new = scatter ( points(1,:), points(2,:), 'r', 'filled' );
%
%  Plot the old points in blue.
%
      hold on

      if ( 0 < points_old_num )
        handle_old = scatter ( points_old(1,:), points_old(2,:), 'b', 'filled' );
      end
%
%  Force the plotting region to be square, not rectangular.
%
      axis square
%
%  Request grid lines.
%
      grid on
%
%  Specify the location of the grid lines, and suppress labeling.
%
      set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
      set ( axes_handle, 'xticklabel', [] );
      set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
      set ( axes_handle, 'yticklabel', [] );
%
%  Make the plotting region slightly bigger than the data.
%
      axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
%
%  Title
%
      s = sprintf ( '%dx%d Clenshaw-Curtis Grid', order_1d(1), order_1d(2) );
      title ( s );

      fprintf ( 1, '+%dx%d Clenshaw-Curtis Grid, Press return\n', ...
        order_1d(1), order_1d(2) );

      pause

      hold off

    end

    fprintf ( 1, 'Press return to CLEAR the grid!\n' );
    pause
    clf

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_GRIDS_ANIMATE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function h = PlotSpectrogram (x, varargin)
% Plot a gray-level spectrogram. The spectrogram is plotted as an image
% with the intensities encoding the levels. The spectrogram has time on the
% abscissa and frequency on the ordinate axis. The spectrogram consists of
% vertical slices displaying the spectral response as intensities in dB.
% The slices are calculated as the DFT of windowed samples. The time window
% for each slice is centred at the corresponding abscissa point. The length
% of the time window determines the time-frequency resolution trade-off.
% The default window length (190 samples) gives relatively good frequency
% resolution.
%
% The intensity values in the plot is calibrated in terms of the decibel
% value relative to a full scale (dBov, see ITU-T G.100.1). A full scale
% sine wave with frequency at one of the DFT values appears  at -6 dBov.
% The default mapping makes the maximum value appear as the darkest value
% and -80 dB below that correspond to the lightest value.
%
%   h = PlotSpectrogram (x, [ts, tf], Fs, Options)
% Simplified forms
%   h = PlotSpectrogram (x, Fs)
%   h = PlotSPectrogram (x)
% Options can appear as the last arguments in any of these forms. Options
% take the form of pairs of arguments, the first being a keyword, the
% second being a value.
%
% h    - handle to the image
% x    - data vector. The data is considered to be extended with zeros at
%        each end in case the windows extend into those regions.
% [ts, tf] - Specifies the start and stop times. The first spectrogram
%        slice is centred at ts and the last spectrogram slice is centred
%        at tf. The default is to create slices spanning the entire data
%        vector, i.e., ts=0, tf=(N-1)/Fs, where N is the number of samples
%        in the data vector.
% Fs   - sampling frequency, default 1.
%
% Options
%  Data window
%    The data window used for each spectrum slice can be specified with
%    the keywords 'Win', 'Lwin' or 'BW'.
%  'Win': the data window vector is explicitly specified.
%  'LWin': the length of the data window (in samples) is specified.
%    A Hamming window of that length is used.
%  'BW': the 3 dB bandwidth (in Hz) of the Hamming window to be used
%    is specified. The length of the Hamming window is chosen using the
%    formula, LWin = 1.18523 * Fs / BW. The default bandwidth is Fs/160.
%    The bandwidth can also be specified with the character strings
%    'NB' (narrowband spectrogram - good frequency resolution, but poor
%    time resolution, same as Fs/160), or 'WB' (wideband spectrogram - poor
%    frequency resolution, but good time resolution, same as 6/160 Fs).
% 'NSlice'
%   The number of spectrogram slices, default 500.
% 'Nfft'
%   The number of rows in the spectrogram is determined by the DFT size,
%   Nfft. The Nfft/2+1 rows correspond to frequencies from 0 to Fs/2. The
%   frequency resolution of the DFT is Fs/Nfft. If not specified, Nfft is
%   chosen to be 1024 or 2048, with the larger value chosen if the window
%   length is greater than 0.9*1024 (0.9 is the Kell factor).
% 'PLimdB'
%   Minimum and maximum levels for the spectrogram on dB. The spectrogram
%   levels are clipped outside of these values. The value for 'PLimdB' can
%   be a two element vector giving the minimum and maximum values. The
%   'PLimdB' value can also be a single value giving the dynamic range in
%   dB. In that case, the maximum is determined automatically and the
%   minimum is set according to the specified dynamic range. The default is
%   to determine the maximum automatically and apply a dynamic range of 80
%   dB.
% 'Amax'
%   The data is expected to take on values between -Amax to +Amax. The
%   default is 1 is suitable for data acquired from sound files. This value
%   affects the scaling of the intensities.
% 'FLim'
%   This parameter is of the form [Fs, Ff] specifying the frequency range
%   of the spectrogram to be plotted. The spectrogram is always calculated
%   over the whole frequency range. This parameter can be used to select
%   a subrange of frequencies for plotting. The default is [0, Fs/2].
% 'preF'
%   The frequency response values can be pre-emphasized to better show
%   lower amplitude high-frequency components. The pre-emphasis response is
%   that of a first order difference filter with parameter preF. The
%   default value of preF is 0. A suitable value for pre-emphasis might be
%   0.97.
%
% Notes:
% - To get a colour spectrogram with a colorbar legend:
%       PlotSpectrogram(...);
%       colormap(SpecColorMap);
%       colorbar;
% - To get an intensity scale in dB SPL (sound pressure level):
%   - Assume that a full scale sine results in PmaxdBSPL (often chosen to
%     be 92 dB SPL).
%   - Then to convert dBov to dB SPL,
%       PmaxdBSPL = 92;  % Max level in dB SPL
%       Amax = 1;        % Assume normalized scaling from a sound file
%       PoffsdB = PmaxdBSPL + 20*log10(2);
%       g = 10^(PoffsdB/20);
%       AmaxN = Amax/g;
%       PlotSpectrogram(..., 'Amax', AmaxN);
%   - This conversion brings the peak sine wave response in the display
%     to be PmaxdBSPL.

% $Id: PlotSpectrogram.m,v 1.10 2009/06/01 18:38:22 pkabal Exp $

% Process arguments
[TLim, Fs, Options] = DecodeArgs(varargin{:});

% Process options
[NSlice, Win, Nfft, FLim, PLimdB, Amax, preF] = ...
                                            PSDecodeOptions(Options, Fs);
% Fill in the default time limits
% Set the times of the slices
if (isempty(TLim))
  TLim = [0 (length(x)-1)/Fs];
end
t = linspace(TLim(1), TLim(2), NSlice);

% Get spectrum in dB
[PdB, f] = SpecSlices(t, x, Fs, Win, Nfft);

% Select a subset of frequencies
I = find(f >= FLim(1) & f <= FLim(2));
f = f(I);
PdB = PdB(I,:);

% Convert to dBov
PdBov = SpecCalib(PdB, Amax, Win);

% Pre-emphasis
pre = (1 + preF^2) - 2 * preF * cos(2 * pi * f / Fs);
predB = 10 * log10(pre');
for (i = 1:NSlice)
  PdBov(:,i) = PdBov(:,i) + predB;
end

% Automatic scaling
if (length(PLimdB) == 1)
  PmaxdB = max(max(PdBov));
  PLimdB = [PmaxdB-PLimdB PmaxdB];
end

% Plot the image
h = imagesc(t, f, PdBov, PLimdB);
axis('xy');

% Grayscale map (dark is more intense)
colormap(flipud(gray));

xlabel('Time (s)');
ylabel('Frequency (Hz)');

return

% ===============================
function PdBov = SpecCalib (PdB, Amax, Win)
% Calibration of a spectrum, The input spectrum in dB is gain modified so
% that a sine wave of maximum amplitude Amax gives a total energy level of
% -3 dBov. This mean that each of the two peaks (one at -fc and the other
% at +fc) will have a level of -6 dBov. A dc level of Amax will give a
% 0 dBov level at zero frequency.
% - Amax, maximum level of the sinusoid (typically 1 for scaled data from
%   sound files). With the true maximum level, the scaling gives an
%   output in dBov. However, Amax can be artificially changed to affect a
%   change in scaling. For instance halving Amax, increases the levels by
%   6 dB.

% Calibration
% - Input:
%   - dc peak amplitude Amax gives 0 dBov
%      - DFT gives sum(Win) at dc
%   - sine peak amplitude Amax give -6dBov at fc and -fc
%     - Assumes sine frequency is coincident with a DFT bin, i.e. fc is of
%       the form m * Fs/Nfft
%     - The peak amplitude of the DFT of the windowed response (assuming no
%       overlap of the components due to the sine at -fc and +fc) is
%       Apeak = abs (Amax / 2 * Wpeak)
Apeak = abs (Amax * sum(Win));
GLdB = -20 * log10(Apeak);

PdBov = PdB + GLdB;

return

% ===============================
function [TLim, Fs, Options] = DecodeArgs(varargin)
% Decode [ts, tf], Fs, Options
% These are distinguished by being a 2 element vector, a scalar, and
% a character string (marking the start of the options)

Options = [];
Fs = 1;
TLim = [];

% Find the first argument with a character string
% Store it and the following arguments in Options
Narg = length(varargin);
for (i = 1:Narg)
  if (ischar(varargin{i}))
    Options = varargin(i:end);
    if (mod(length(Options), 2) ~= 0)
      error('PlotSpectrogram: Invalid Options format');
    end
    varargin(i:end) = [];
    break;
  end
end

% Check for a TLim argument
if (length(varargin) > 1 && length(varargin{1}) == 2)
  TLim = varargin{1};
  varargin(1) = [];
end

% Pick off Fs
if (~ isempty(varargin) && length(varargin{1}) == 1)
  Fs = varargin{1};
  varargin(1) = [];
end

if (~ isempty(varargin))
  error('PlotSpectrogram: Too many input arguments');
end

return

author: Sable

repo: mcbench-benchmarks