text:
function annModel = annotateGEM(model,annPath,annType,addMiriams,addFields,overwrite)
% Add reaction, metabolite, and/or gene annotation to a model.
%
% Input:
%
%   model        Model structure.
%
%   annPath      Path to the annotation files, which suppose to be named as
%                'reactions.tsv', 'metabolites.tsv', and 'genes.tsv'。
%
%   annType      String or cell array of strings specifying the type(s) of
%                annotation data to add: 'rxn', 'met', and/or 'gene'. To
%                add all annotation types, use 'all'.
%                For example, annType={'rxn','gene'} will add reaction and
%                gene annotation information, but not metabolite.
%                (Opt, default 'all')
%
%   addMiriams   If TRUE, annotation information will be added to the
%                metMiriams field. If the metMiriams field does not exist,
%                it will be created.
%                (Opt, default TRUE)
%
%   addFields    If TRUE, separate model fields will be added for each
%                added ID type, e.g., reaction KEGG IDs will be added to a
%                rxnKEGGID field.
%                (Opt, default TRUE)
%
%   overwrite    If TRUE, any existing metMiriam or model ID field entries
%                will be overwritten. If FALSE, they will be appended with
%                the new data.
%                (Opt, default TRUE)
%
%
% Output:
%
%   annModel     Model structure with added annotation information.
%
%
% Usage:
%
%   annModel = annotateGEM(model,annPath,annType,addMiriams,addFields,overwrite);
%

%% Inputs and setup

if nargin < 2
    [ST, I] = dbstack('-completenames');
    annPath = strcat(fileparts(ST(I).file),'/../model');
end

if nargin < 3 || isempty(annType) || isequal(annType,'all')
    annType = {'rxn','met','gene'};
elseif ~all(ismember(annType,{'rxn','met','gene','reaction','metabolite'}))
    error('annType input(s) not recognized. Valid options are "rxn", "met", and/or "gene", or "all"');
end

if nargin < 4 || isempty(addMiriams)
    addMiriams = true;
end

if nargin < 5 || isempty(addFields)
    addFields = true;
end

if nargin < 6
    overwrite = true;
end

if (~addMiriams) && (~addFields)
    error('addMiriams and addFields are both FALSE - no annotation will be added.')
end

% map field names to those on identifiers.org (used for building Miriams)
id2miriam = {%reactions
             'rxnKEGGID'        'kegg.reaction'
             'rxnBiGGID'        'bigg.reaction'
             'rxnREACTOMEID'    'reactome'
             'rxnRecon3DID'     'vmhreaction'
             'rxnMetaNetXID'    'metanetx.reaction'
             'rxnTCDBID'        'tcdb'
             'rxnRheaID'        'rhea'
             'rxnRheaMasterID'  'rhea'
             % metabolites
             'metBiGGID'        'bigg.metabolite'
             'metKEGGID'        'kegg.compound'
             'metHMDBID'        'hmdb'
             'metChEBIID'       'chebi'
             'metPubChemID'     'pubchem.compound'
             'metLipidMapsID'   'lipidmaps'
             'metRecon3DID'     'vmhmetabolite'
             'metMetaNetXID'    'metanetx.chemical'
             % genes
             'genes'            'ensembl'
             'geneENSTID'       'ensembl'
             'geneENSPID'       'ensembl'
             'geneUniProtID'    'uniprot'
             'geneSymbols'      'hgnc.symbol'
             'geneEntrezID'     'ncbigene'
             };

%% Load and organize annotation data

% load reaction annotation data
if any(ismember({'rxn','reaction'},lower(annType)))
    tmpfile = fullfile(annPath,'reactions.tsv');
    rxnAssoc = importTsvFile(tmpfile);

    % strip "RHEA:" prefix from Rhea IDs since it should not be included in
    % the identifiers.org URL
    rxnAssoc.rxnRheaID = regexprep(rxnAssoc.rxnRheaID, 'RHEA:', '');
    rxnAssoc.rxnRheaMasterID = regexprep(rxnAssoc.rxnRheaMasterID, 'RHEA:', '');

    rxnAssocArray = struct2cell(rxnAssoc);
    numericFields = find(cellfun(@isnumeric, rxnAssocArray));
    for i = 1:numel(numericFields)
        % numeric fields must be converted to text for miriams
        rxnAssocArray{numericFields(i)} = arrayfun(@num2str, rxnAssocArray{numericFields(i)}, 'UniformOutput', false);
    end
    rxnAssocArray = horzcat(rxnAssocArray{:});
else
    rxnAssoc = [];
end

% load metabolite annotation data
if any(ismember({'met','metabolite'},lower(annType)))
    tmpfile = fullfile(annPath,'metabolites.tsv');
    metAssoc = importTsvFile(tmpfile);

    % ChEBI IDs should be of the form "CHEBI:#####"
    metAssoc.metChEBIID = regexprep(metAssoc.metChEBIID,'(^)(\d+)|(;\s*)(\d+)','$1CHEBI:$2');

    metAssocArray = struct2cell(metAssoc);
    numericFields = find(cellfun(@isnumeric, metAssocArray));
    for i = 1:numel(numericFields)
        % numeric fields need to be converted to text for miriams
        metAssocArray{numericFields(i)} = arrayfun(@num2str, metAssocArray{numericFields(i)}, 'UniformOutput', false);
    end
    metAssocArray = horzcat(metAssocArray{:});
else
    metAssoc = [];
end

% load and organize gene annotation data
if ismember('gene',lower(annType))
    tmpfile = fullfile(annPath,'genes.tsv');
    geneAssoc = importTsvFile(tmpfile);

    % add geneEnsemblID field if missing
    if ~isfield(geneAssoc, 'geneEnsemblID')
        geneAssoc.geneEnsemblID = geneAssoc.genes;
    end

    geneAssocArray = struct2cell(geneAssoc);
    geneAssocArray = horzcat(geneAssocArray{:});
else
    geneAssoc = [];
end

%% Add annotation data to model MIRIAM fields

if ( addMiriams )

    % Reactions
    if ~isempty(rxnAssoc)
        % map model rxns to those in the association structures
        [hasRxn,rxnInd] = ismember(model.rxns,rxnAssoc.rxns);

        % add rxnMiriams field if it does not exist or if overwrite is TRUE
        if ~isfield(model,'rxnMiriams') || overwrite
            model.rxnMiriams = repmat({''},size(model.rxns));
        end

        % map rxn ID fields to the proper identifiers.org terms
        [hasRxnField,rxnFieldInd] = ismember(fieldnames(rxnAssoc),id2miriam(:,1));
        miriamNames = id2miriam(rxnFieldInd(hasRxnField),2);
        miriamValues = rxnAssocArray(:,hasRxnField);

        % assign SBO terms to model reactions
        rxnSBO = repmat({'SBO:0000176'},size(model.rxns));  % default is "biochemical reaction"

        biomass_rxn = contains(lower(model.rxns),'biomass');
        biomass_met = ismember(lower(model.metNames),'biomass') & ismember(model.comps(model.metComps),'c');
        biomass_rxn = biomass_rxn | (model.S(biomass_met,:) > 0)';
        rxnSBO(biomass_rxn) = {'SBO:0000629'};  % "biomass production"

        transport_rxn = getTransportRxns(model);
        rxnSBO(transport_rxn) = {'SBO:0000185'};  % "translocation reaction"

        [~,exchange_rxn] = getExchangeRxns(model);
        rxnSBO(exchange_rxn) = {'SBO:0000627'};  % "exchange reaction"

        % add new data to the rxnMiriams field
        for i = 1:numel(model.rxns)
            if hasRxn(i)
                % add annotation
                model.rxnMiriams{i} = appendMiriamData(model.rxnMiriams{i}, miriamNames, miriamValues(rxnInd(i),:)');
            end
            % add SBO terms
            model.rxnMiriams{i} = appendMiriamData(model.rxnMiriams{i}, {'sbo'}, rxnSBO(i));
        end
    end

    % Metabolites
    if ~isempty(metAssoc)
        % map model mets to those in the association structures
        [hasMet,metInd] = ismember(model.mets,metAssoc.mets);

        % add metMiriams field if it does not exist or if overwrite is TRUE
        if ~isfield(model,'metMiriams') || overwrite
            model.metMiriams = repmat({''},size(model.mets));
        end

        % map met ID fields to the proper identifiers.org terms
        [hasMetField,metFieldInd] = ismember(fieldnames(metAssoc),id2miriam(:,1));
        miriamNames = id2miriam(metFieldInd(hasMetField),2);
        miriamValues = metAssocArray(:,hasMetField);

        % add new data to the metMiriams field
        for i = 1:numel(model.mets)
            if hasMet(i)
                % add annotation
                model.metMiriams{i} = appendMiriamData(model.metMiriams{i}, miriamNames, miriamValues(metInd(i),:)');
            end
            % add SBO term (SBO:0000247, "simple chemical" for all mets)
            model.metMiriams{i} = appendMiriamData(model.metMiriams{i}, {'sbo'}, {'SBO:0000247'});
        end
    end

    % Genes
    if ~isempty(geneAssoc)
        % Note: because geneAssoc was built using geneIDs directly from the
        %       model, we do not need to map model.genes to the geneAssoc
        %       structure - the indexing should be identical.

        % map model genes to those in the association structures
        [hasGene,geneInd] = ismember(model.genes,geneAssoc.genes);

        % add geneMiriams field if it does not exist or if overwrite is TRUE
        if ~isfield(model,'geneMiriams') || overwrite
            model.geneMiriams = repmat({''},size(model.genes));
        end

        % map gene ID fields to the proper identifiers.org terms
        [hasGeneField,geneFieldInd] = ismember(fieldnames(geneAssoc),id2miriam(:,1));
        miriamNames = id2miriam(geneFieldInd(hasGeneField),2);
        miriamValues = geneAssocArray(:,hasGeneField);

        % add new data to the geneMiriams field
        for i = 1:numel(model.genes)
            if hasGene(i)
                % add annotation
                model.geneMiriams{i} = appendMiriamData(model.geneMiriams{i}, miriamNames, miriamValues(geneInd(i),:)');
            end
            % add SBO term (SBO:0000243, "gene" for all genes)
            model.geneMiriams{i} = appendMiriamData(model.geneMiriams{i}, {'sbo'}, {'SBO:0000243'});
        end
    end

end

%% Add annotation data to individual model fields

if ( addFields )

    % merge association structures
    allAssoc = mergeStructures({rxnAssoc,metAssoc,geneAssoc});

    % some fields should not be added to the model
    remFields = {'rxns','mets','metsNoComp','genes'};
    allAssoc = rmfield(allAssoc, intersect(fieldnames(allAssoc),remFields));

    % get fields and their types
    f = fieldnames(allAssoc);

    if ~isempty(rxnAssoc)
        fieldType = repmat({'rxn'}, numel(f), 1);
    end

    if ~isempty(metAssoc)
        fieldType(ismember(f, fieldnames(metAssoc))) = {'met'};
    end

    if ~isempty(geneAssoc)
        fieldType(ismember(f, fieldnames(geneAssoc))) = {'gene'};
    end

    % add individual ID fields to the model
    for i = 1:numel(f)

        switch fieldType{i}
            case 'rxn'
                [hasRxn,rxnInd] = ismember(model.rxns,rxnAssoc.rxns);
                if isnumeric(allAssoc.(f{i}))
                    ids = NaN(size(model.rxns));
                else
                    ids = repmat({''},size(model.rxns));
                end
                ids(hasRxn) = allAssoc.(f{i})(rxnInd(hasRxn));
            case 'met'
                [hasMet,metInd] = ismember(model.mets,metAssoc.mets);
                if isnumeric(allAssoc.(f{i}))
                    ids = NaN(size(model.mets));
                else
                    ids = repmat({''},size(model.mets));
                end
                ids(hasMet) = allAssoc.(f{i})(metInd(hasMet));
            case 'gene'
                [hasGene,geneInd] = ismember(model.genes,geneAssoc.genes);
                if isnumeric(allAssoc.(f{i}))
                    ids = NaN(size(model.genes));
                else
                    ids = repmat({''},size(model.genes));
                end
                ids(hasGene) = allAssoc.(f{i})(geneInd(hasGene));
            otherwise
                continue
        end

        if ~isfield(model,f{i}) || overwrite
            model.(f{i}) = ids;
        else
            % merge new IDs with existing IDs
            merged_ids = join([model.(f{i}), ids],';');
            merged_ids = arrayfun(@(id) strjoin(unique(split(id,';')),';'), merged_ids, 'UniformOutput', false);
            model.(f{i}) = regexprep(merged_ids,'^;\s*','');  % remove any preceding semicolons
        end

    end

end

%%

% assign output
annModel = model;

end

%% Additional functions

function miriam_new = appendMiriamData(miriam,names,values)
% Append new entries to an existing miriam structure.

% ignore empty entries
remove = cellfun(@isempty,values);
names(remove) = [];
values(remove) = [];
if isempty(names)
    miriam_new = miriam;
    return
end

% convert existing miriam entry to cell array
if ~isempty(miriam)
    miriam = [miriam.name(:), miriam.value(:)];
end

% append miriam entry with new data
for i = 1:numel(names)
    if contains(values{i},';')
        % handle multiple values separated by semicolon (;)
        val = strtrim(split(values{i},';'));
        name = repmat(names(i),numel(val),1);
        miriam = [miriam; [name, val]];
    else
        miriam = [miriam; [names(i), values(i)]];
    end
end

% remove any repeated entries
[~,miriam_num] = ismember(miriam,miriam);
[~,uniq_ind] = unique(miriam_num,'rows');
miriam = miriam(uniq_ind,:);

% convert back to structure
miriam_new.name = miriam(:,1);
miriam_new.value = miriam(:,2);

end

function mergedStruct = mergeStructures(structs)
% Merge multiple structures into a single structure.

mergedStruct = {};
for i = 1:numel(structs)

    s = structs{i};

    if isempty(s)
        continue
    end

    f = fieldnames(s);
    for j = 1:numel(f)
        mergedStruct.(f{j}) = s.(f{j});
    end

end

end

author: SysBioChalmers

repo: Human-GEM


text:
function vis_footnoteButtonDownFcn

% VIS_FOOTNOTEBUTTONDOWNFCN Callback set by VIS_FOOTNOTE
%
%  som_showtitleButtonDownFcn
%
% Moves the axis of current callback object using DRAGRECT
% command. This callback is set to all texts added to figures by
% VIS_FOOTNOTE function.
%
% See also DRAGRECT, SOM_SHOWTITLE.

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

% Version 2.0beta Johan 080698

% Action

[txt,fig]=gcbo;                     % Get text and figure handles
ax=get(txt,'parent');               % Get axis handle

memunits_fig=get(fig,'units');      % Get figure size in pixels
set(gcf,'units','pixels');
pos_fig=get(fig,'position');

memunits_txt=get(txt,'units');      % Get text field size in pixels
set(txt,'units','pixels');
text_size=get(txt,'extent');

memunits_ax=get(ax,'units');        % Get axis position in pixels
set(ax,'units','pixels');
pos_ax=get(ax,'position');

%%% Move text

pos_final=dragrect([pos_ax(1:2) text_size(3:4)]);

%%% Keep the text inside the figure!

pos_final(1)=max(pos_final(1),0);
pos_final(2)=max(pos_final(2),0);
pos_final(1)=min(pos_final(1),pos_fig(3)-text_size(3));
pos_final(2)=min(pos_final(2),pos_fig(4)-text_size(4));

%%% Set new position

new_pos=[pos_final(1:2) pos_ax(3:4)];
set(ax,'position', new_pos);

%%% Set the text on the top of the object stack

children=get(gcf,'children');
i=find(ismember(children,ax));
new_i=[i 1:i-1 i+1:length(children)];
set(gcf,'children',children(new_i));

set(txt,'position',[0 0 0]);
set(fig,'units',memunits_fig);
set(ax,'units',memunits_ax);
set(txt,'units',memunits_txt);

author: ilarinieminen

repo: SOM-Toolbox


text:
%tstoolbox/mex/amutual
%   Fast, but crude auto mutual information of a scalar timeseries for the
%   timelags from zero to maxtau. The input time series should be much
%   longer than maximal timelag maxtau. The algorithm uses equidistant
%   histogram boxes, so results are bad in a mathematical sense. However,
%   a fast algorithm based on ternary search trees to store only nonempty
%   boxes is used.
%
%   Syntax:
%
%     * a = amutual(ts, maxtau, partitions)
%
%   Input arguments:
%
%     * ts - vector holding time series data
%     * maxtau - maximal time lag
%     * partitions - number of partitions for the one-dimensional
%       histogram
%
%   Output arguments:
%
%     * a - vector of length maxtau+1, holding auto mutual information
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function checkHeader
%CHECKHEADER checks the header of the m-files of matGeom.
%
%   checkHeader()
%   Checks and corrects:
%   - the beginning of the first header line
%   - the "author" section is checked and corrected
%   - the "see also" section
%   Looks for:
%   - to-dos
%   Removes:
%   - the "history" section
%   Additional checks could be added.
%
%   Todo
%       - Standardize E-mail adress to "david.legland@inrae.fr"?
%       - Add checks for "Example" section.
%       - Add additional checks for the "See also" section: No
%       self-reference, references, indentations of the references.
%
%   See also
%       codingConventions
%

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2022-11-08, using MATLAB 9.13.0.2080170 (R2022b) Update 1
% Copyright 2022

cd(strrep(mfilename('fullpath'),mfilename,''))

mFiles = dir('../matGeom/**/*.m');
% Exclude the Contents.m files
mFiles(contains({mFiles.folder}','deprecated'))=[];

for f=1:length(mFiles)
    % Full path of the m-file
    mFile = fullfile(mFiles(f).folder, mFiles(f).name);
    % Read the content of the m-file
    S = readlines(mFile);
    % Make a backup of S
    S_BU = S;
    % File name of the m-file without the file extension
    mFileName = mFiles(f).name(1:end-2);
    % Length of the file name
    LoN = length(mFileName);

    % Assume that the header is invalid
    mFiles(f).validHeader = false;

    %% First header line
    % If the 1st or 2nd line of the file is empty the header is not valid
    if isempty(S{1}) || isempty(S{2})
        continue
    end

    % Check if the 1st header line starts with the capitalized m-file name
    if startsWith(S{1}, 'function')
        H1_Idx = 2;
    else
        H1_Idx = 1;
    end
    % Try to correct the beginning of the 1st header line
    if ~strcmp(S{H1_Idx}(1:2+LoN),['%' upper(mFileName) ' '])
        if startsWith(S{H1_Idx}, '% ') && regexp(S{H1_Idx}, '% ') == 1
            S{H1_Idx} = strrep(S{H1_Idx}, '% ', ['%' upper(mFileName) ' ']);
        end
    end

    if ~endsWith(S(H1_Idx), '.')
        S(H1_Idx) = S{H1_Idx} + ".";
    end

    %% AUTHOR section
    % First line of the header's author section
    % Find the 1st line of the header's author section
    aSec1Idx = find(contains(S,'-----'));
    if isempty(aSec1Idx)
        continue
    elseif length(aSec1Idx) ~= 1
        display(['No line contains or multiple lines contain: ''-----'' in: ' mFileName])
        continue
    end

    if ~strcmp(S(aSec1Idx),'% ------')
        % Remove more than six "-"
        S(aSec1Idx) = regexprep(S(aSec1Idx),'% *-----+','% ------');
    end

    if regexp(S(aSec1Idx+1), '^% *$') == 1
        S(aSec1Idx+1) = [];
    end

    % "Author" line
    % Fix the beginning of the author line
    if contains(S(aSec1Idx+1),'author') && ~contains(S(aSec1Idx+1),'Author')
        S(aSec1Idx+1) = regexprep(S(aSec1Idx+1),'% +author *:','% Author:');
        S(aSec1Idx+1) = regexprep(S(aSec1Idx+1),'% +authors *:','% Authors:');
    elseif contains(S(aSec1Idx+1),'  Author') && ~startsWith(S(aSec1Idx+1), '% Author')
        S(aSec1Idx+1) = regexprep(S(aSec1Idx+1),'% +','% ');
    end

    % "E-mail" line
    if ~contains(S(aSec1Idx+2),'mail','IgnoreCase',1)
        if contains(S(aSec1Idx+2),'created','IgnoreCase',1)
            if contains(S(aSec1Idx+1),'David Legland')
                S = [S(1:aSec1Idx+1); "% E-mail: david.legland@inrae.fr"; S(aSec1Idx+2:end)];
            else
                S = [S(1:aSec1Idx+1); "% E-mail: N/A"; S(aSec1Idx+2:end)];
            end
        elseif contains(S(aSec1Idx+3),'created','IgnoreCase',1)
            coIdx = regexp(S(aSec1Idx+2), 'INRA');
            if any(coIdx)
                % Extract the year of creation for the copyright
                creationYear = regexp(S(aSec1Idx+3), '20\d\d','match');
                if isempty(creationYear); continue; end
                if regexp(S(aSec1Idx+4), '^% *$') == 1
                        S(aSec1Idx+4) = regexprep(S(aSec1Idx+2),'% +INRA','% Copyright ' + creationYear + ' INRA');
                elseif isempty(S{aSec1Idx+4})
                    S(aSec1Idx+4) = "% Copyright " + creationYear + " " + S{aSec1Idx+2}(coIdx:end);
                    S = [S(1:aSec1Idx+4); ''; S(aSec1Idx+5:end)];
                else
                    continue
                end
                if contains(S(aSec1Idx+1),'David Legland')
                    S(aSec1Idx+2) = "% E-mail: david.legland@inrae.fr";
                else
                    S(aSec1Idx+2) = "% E-mail: N/A";
                end
            else
                continue
            end
        end
    end
    % Change e-mail to E-Mail
    S(aSec1Idx+2) = regexprep(S(aSec1Idx+2), '% e-mail:', '% E-mail:');

    % Move affiliation to the copyright line if necesarry
    if contains(S(aSec1Idx+3),'INRA') && contains(S(aSec1Idx+4),'created','IgnoreCase',1)
        coIdx = regexp(S(aSec1Idx+3), 'INRA');
        if contains(S(aSec1Idx+5), 'Copyright')
            S(aSec1Idx+5) = regexprep(S(aSec1Idx+5),'INRA.*',S{aSec1Idx+3}(coIdx:end));
            S(aSec1Idx+3) = [];
        elseif regexp(S(aSec1Idx+5), '^% *$') == 1
            creationYear = regexp(S(aSec1Idx+4), '20\d\d','match');
            if ~isempty(creationYear)
                S(aSec1Idx+5) = "% Copyright " + creationYear + " " + S{aSec1Idx+3}(coIdx:end);
                S(aSec1Idx+3) = [];
            else
                continue
            end
        elseif isempty(S{aSec1Idx+5})
            creationYear = regexp(S(aSec1Idx+4), '20\d\d','match');
            if ~isempty(creationYear)
                S(aSec1Idx+5) = "% Copyright " + creationYear + " " + S{aSec1Idx+3}(coIdx:end);
                S = [S(1:aSec1Idx+5); ''; S(aSec1Idx+6:end)];
                S(aSec1Idx+3) = [];
            else
                continue
            end
        end
    end

    % "Created" line
    if contains(S(aSec1Idx+3),'created the') && ~contains(S(aSec1Idx+3),'Created')
        S(aSec1Idx+3) = regexprep(S(aSec1Idx+3),'% +created the +','% Created: ');
    elseif contains(S(aSec1Idx+3),'created ')
        S(aSec1Idx+3) = regexprep(S(aSec1Idx+3),'% +created +','% Created: ');
    elseif contains(S(aSec1Idx+3),'  Created') && ~startsWith(S(aSec1Idx+3), '% Created')
        S(aSec1Idx+3) = regexprep(S(aSec1Idx+3),'% +','% ');
    elseif contains(S(aSec1Idx+3),' created:')
        S(aSec1Idx+3) = regexprep(S(aSec1Idx+3),'% +created:','% Created:');
    end

    % Standardize creation date
    dateIdxS = regexp(S(aSec1Idx+3),'% Created: ','end') + 1;
    if any(dateIdxS) && length(dateIdxS) == 1
        dateIdxE = regexp(S(aSec1Idx+3), '( by | from|, by |, from|, +using|\.$)') - 1;
        if isempty(dateIdxE)
            dateIdxE = length(S{aSec1Idx+3});
            if length(S{aSec1Idx+3}(dateIdxS:dateIdxE)) >= 10
                if regexp(S{aSec1Idx+3}(dateIdxS:dateIdxS+9),'\d\d\/\d\d\/\d{4}') == 1
                    dateStr = string(datetime(S{aSec1Idx+3}(dateIdxS:dateIdxE),...
                        'InputFormat','dd/MM/uuuu','format','uuuu-MM-dd'));
                else
                    dateStr = string(datetime(S{aSec1Idx+3}(dateIdxS:dateIdxE),'format','uuuu-MM-dd'));
                end
                S(aSec1Idx+3) = "% Created: " + dateStr;
            end
        else
            if length(dateIdxE)~=1; dateIdxE = dateIdxE(1); end
            if length(S{aSec1Idx+3}(dateIdxS:dateIdxE)) >= 10
                if regexp(S{aSec1Idx+3}(dateIdxS:dateIdxS+9),'\d\d\/\d\d\/\d{4}') == 1
                    dateStr = string(datetime(S{aSec1Idx+3}(dateIdxS:dateIdxE),...
                        'InputFormat','dd/MM/uuuu','format','uuuu-MM-dd'));
                else
                    dateStr = string(datetime(S{aSec1Idx+3}(dateIdxS:dateIdxE),'format','uuuu-MM-dd'));
                end
                S(aSec1Idx+3) = "% Created: " + dateStr + S{aSec1Idx+3}(dateIdxE+1:end);
            end
        end
    end

    % "Copyright" line
    % Extract the year of creation for the copyright
    creationYear = regexp(S(aSec1Idx+3), '20\d\d','match');
    if ~isempty(creationYear)
        creationYear = creationYear(1);
    else
        continue
    end
    if isempty(S{aSec1Idx+4})
        S = [S(1:aSec1Idx+3); "% Copyright " + creationYear; S(aSec1Idx+4:end)];
    elseif regexp(S(aSec1Idx+4), '^% *$') == 1
        S(aSec1Idx+4) = "% Copyright " + creationYear;
    end

    % Add current year to the copyright or change it to the current year
    copyrightSYear = regexp(S(aSec1Idx+4), '20\d\d','match');
    copyrightSYear = copyrightSYear{1};
    % Created year and copyright year should be the same
    if ~isequal(copyrightSYear,creationYear{1})
        display(['Creation and copyright year are inconsistent in: ' mFileName])
        continue
    end
    currentYear = char(datetime('today','Format','yyyy'));
    if regexp(S(aSec1Idx+4), '% Copyright 20\d\d( |$)') == 1
        if ~isequal(copyrightSYear, currentYear)
            S(aSec1Idx+4) = regexprep(S(aSec1Idx+4), '% Copyright 20\d\d', ...
                ['% Copyright ' copyrightSYear '-' currentYear]);
        end
    elseif regexp(S(aSec1Idx+4), '% Copyright 20\d\d-20\d\d( |$)') == 1
        copyrightEYear = S{aSec1Idx+4}(18:21);
        if ~isequal(copyrightEYear, currentYear)
            S(aSec1Idx+4) = regexprep(S(aSec1Idx+4), '% Copyright 20\d\d-20\d\d', ...
                ['% Copyright ' copyrightSYear '-' currentYear]);
        end
    else
        continue
    end

    % General checks
    containsHeaderWords = arrayfun(@(x,y) contains(x,y), ...
        S(aSec1Idx:aSec1Idx+4), ["% ------"; "Author"; "E-mail"; "Created"; "Copyright"]);
    if ~all(containsHeaderWords)
        continue
    end
    if any(contains(S(aSec1Idx:aSec1Idx+4), '%  '))
        % Remove leading spaces before the header words
        S(aSec1Idx:aSec1Idx+4) = arrayfun(@(x) regexprep(x,'% +','% '), S(aSec1Idx:aSec1Idx+4));
    end

    aSecLIs = [...
        startsWith(S(aSec1Idx-2), '%') ...
        isempty(S{aSec1Idx-1}) ...
        startsWith(S(aSec1Idx+1), '% Author') ...
        startsWith(S(aSec1Idx+2), '% E-mail: ') ...
        startsWith(S(aSec1Idx+3), '% Created: ')...
        startsWith(S(aSec1Idx+4), '% Copyright ')...
        isempty(S{aSec1Idx+5}) ...
        ];
    if all(aSecLIs)
        mFiles(f).validHeader = true;
    elseif all(aSecLIs == [1 0 1 1 1 1 1])
        S(aSec1Idx-1) = '';
    elseif all(aSecLIs == [0 0 1 1 1 1 1])
        if isempty(S{aSec1Idx-2})
            S(aSec1Idx-2) = '%';
            S(aSec1Idx-1) = '';
        end
    elseif all(aSecLIs == [1 0 1 1 1 1 0]) || all(aSecLIs == [1 1 1 1 1 1 0])
        if regexp(S(aSec1Idx+5), '^% *$') == 1
            S([aSec1Idx-1, aSec1Idx+5]) = '';
        end
    else
        continue
    end

    % No dot after creation date
    if regexp(S(aSec1Idx+3), '% Created: 20\d\d\-\d\d\-\d\d\.','start') == 1
        S(aSec1Idx+3) = regexprep(S(aSec1Idx+3),'\. *$','');
    end
    % Remove multiple spaces after creation date
    S(aSec1Idx+3) = regexprep(S(aSec1Idx+3),' +',' ');
    % If text follows the creation date seperate it by a comma
    if regexp(S(aSec1Idx+3), '% Created: 20\d\d\-\d\d\-\d\d\ .','start') == 1
        S(aSec1Idx+3) = [S{aSec1Idx+3}(1:21) ',' S{aSec1Idx+3}(22:end)];
    end
    % Remove dots at the end of the Created and Copyright line
    if endsWith(S(aSec1Idx+3), '.')
        S(aSec1Idx+3) = S{aSec1Idx+3}(1:end-1);
    end
    if endsWith(S(aSec1Idx+4), '.')
        S(aSec1Idx+4) = S{aSec1Idx+4}(1:end-1);
    end

    if aSec1Idx+6 < length(S)
        while isempty(S{aSec1Idx+6})
            S(aSec1Idx+6) = [];
        end
    end

    %% SEE ALSO section
    if ~isempty(cell2mat(regexp(lower(S), '^% *see also')))
        saSIdx = find(contains(lower(S),'see also'));
        if length(saSIdx) ~=1
            display(['Multiple ''See also'' lines found in: ' mFileName])
            continue
        end
        S(saSIdx) = regexprep(S(saSIdx), '^% *[sS]ee [aA]lso *:* *', '%   See also ');
    end

    %% HISTORY section
    if ~isempty(cell2mat(regexp(lower(S), '^% *history')))
        histSIdx = find(contains(lower(S),'history'));
        if length(histSIdx) ~=1
            display(['Multiple ''History'' sections found in: ' mFileName])
            continue
        end
        histLines = 1;
        while startsWith(S(histSIdx+histLines),'% ')
            histLines = histLines+1;
        end
        S(histSIdx:histSIdx+histLines-1) = [];
    end
    % Check for dates that appear after the header and may be part of a of a "history"
    if ~isempty(cell2mat(regexp(S(aSec1Idx+5:end), '^%.*20\d\d.')))
        display(['Suspicous date found in: ' mFileName])
    end

    %% TODO section
    % There should not be a to-do section. To-dos should be adressed using
    % GitHub Issues.
    if ~isempty(cell2mat(regexp(lower(S), '^%.*(todo|to-do)')))
        display(['''To-do(s)'' found in: ' mFileName])
    end

    %% Save file if changes were performed
    if ~isequal(S, S_BU)
        writelines(S, mFile);
    end
end

% Display remaining invalid files
mFileNames = {mFiles.name}';
validHeaders = {mFiles.validHeader}';
invalidFilesIdx = cellfun(@(x) ~isempty(x) && x==0, validHeaders);
invalidFiles = mFileNames(invalidFilesIdx);
NoIF = length(invalidFiles);
if isempty(invalidFiles)
    invalidFiles='none';
end
disp([num2str(NoIF) ' files with invalid header:'])
disp([num2cell(find(invalidFilesIdx)), invalidFiles])

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [lat, lon, azi, rk] = cassini_inv(lat0, lon0, x, y, ellipsoid)
%CASSINI_INV  Inverse Cassini-Soldner projection
%
%   [LAT, LON] = CASSINI_INV(LAT0, LON0, X, Y)
%   [LAT, LON, AZI, RK] = CASSINI_INV(LAT0, LON0, X, Y, ELLIPSOID)
%
%   performs the inverse Cassini-Soldner projection of points (X,Y) to
%   (LAT,LON) using (LAT0,LON0) as the center of projection.  These input
%   arguments can be scalars or arrays of equal size.  The ELLIPSOID vector
%   is of the form [a, e], where a is the equatorial radius in meters, e is
%   the eccentricity.  If ellipsoid is omitted, the WGS84 ellipsoid (more
%   precisely, the value returned by DEFAULTELLIPSOID) is used.  GEODPROJ
%   defines the projection and gives the restrictions on the allowed ranges
%   of the arguments.  The forward projection is given by CASSINI_FWD.
%
%   AZI and RK give metric properties of the projection at (LAT,LON); AZI
%   is the azimuth of the easting (X) direction and RK is the reciprocal of
%   the northing (Y) scale.  The scale in the easting direction is 1.
%
%   LAT0, LON0, LAT, LON, AZI are in degrees.  The projected coordinates X,
%   Y are in meters (more precisely the units used for the equatorial
%   radius).  RK is dimensionless.
%
%   This routine depends on the MATLAB File Exchange package "Geodesics on
%   an ellipsoid of revolution":
%
%     http://www.mathworks.com/matlabcentral/fileexchange/39108
%
%   See also GEODPROJ, CASSINI_FWD, GEODRECKON, DEFAULTELLIPSOID.

% Copyright (c) Charles Karney (2012) <charles@karney.com>.
%
% This file was distributed with GeographicLib 1.29.

  if nargin < 4, error('Too few input arguments'), end
  if nargin < 5, ellipsoid = defaultellipsoid; end
  try
    [~] = lat0 + lon0 + x + y;
  catch err
    error('lat0, lon0, x, y have incompatible sizes')
  end

  [lat1, lon1, azi0] = geodreckon(lat0, lon0, y, 0, ellipsoid);
  [lat, lon, azi, ~, ~, rk] = ...
      geodreckon(lat1, lon1, x, azi0 + 90, ellipsoid);
end

author: Sable

repo: mcbench-benchmarks


text:
function [irf_record, D_record, gamma_record, struct_irf_record, irf_estimates, D_estimates, gamma_estimates, strshocks_record, strshocks_estimates]...
    =panel2irf(Ymat,Xmat,beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,IRFband,N,n,m,p,k,T,Y,X,signreslabels,FEVDresperiods,data_exo,const,exo,IRFt,strctident,favar,signrestable,signresperiods)

% because there is only one model, there is no need to loop over units
% hence, estimate the IRFs for the model, as with a standard normal-Wishart
% first run the Gibbs sampler to obtain posterior draws
[irf_record]=bear.irf(beta_gibbs,It,Bu,IRFperiods,n,m,p,k);

% If IRFs have been set to an unrestricted VAR (IRFt=1):
if IRFt==1
% run a pseudo Gibbs sampler to obtain records for D and gamma (for the trivial SVAR)
[D_record, gamma_record]=bear.irfunres(n,It,Bu,sigma_gibbs);
struct_irf_record=[];
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(irf_record,n,IRFperiods,IRFband,IRFt,[],[],favar);

% If IRFs have been set to an SVAR with Choleski identification (IRFt=2):
elseif IRFt==2
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record, D_record, gamma_record]=bear.irfchol(sigma_gibbs,irf_record,It,Bu,IRFperiods,n,favar);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,n,IRFperiods,IRFband,IRFt,D_record,gamma_record,favar);

% If IRFs have been set to an SVAR with triangular factorisation (IRFt=3):
elseif IRFt==3
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record, D_record, gamma_record]=bear.irftrig(sigma_gibbs,irf_record,It,Bu,IRFperiods,n,favar);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,n,IRFperiods,IRFband,IRFt,D_record,gamma_record,favar);

% If IRFs have been set to an SVAR with sign restrictions (IRFt=4):
elseif IRFt==4
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record,D_record,gamma_record]=bear.irfres_old(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,m,p,k,signrestable,signresperiods);
%[struct_irf_record,D_record,gamma_record]=bear.irfsignrespanel(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,p,m,k,signrestable,signresperiods);
%[struct_irf_record,D_record,gamma_record,~,~,~,~]=bear.irfres(beta_gibbs,sigma_gibbs,[],[],IRFperiods,n,m,p,k,T,[],[],signreslabels,FEVDresperiods,data_exo,HD,const,exo,strctident,pref,favar,IRFt,It,Bu);

% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,n,IRFperiods,IRFband,IRFt,D_record,gamma_record,favar);

end

% also, if a s structural identification was implemented, compute structural shocks
strshocks_record={};
strshocks_estimates={};
if IRFt~=1
   % because shocks have to be computed for each unit, loop over units
   for ii=1:N
   % run the Gibbs sampler
   strshocks_record(:,:,ii)=bear.strshocks(beta_gibbs,D_record,Ymat(:,:,ii),Xmat(:,:,ii),n,k,It,Bu,favar);
   % obtain point estimates and credibility intervals
   strshocks_estimates(:,:,ii)=bear.strsestimates(strshocks_record(:,:,ii),n,T,IRFband);
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function loc=LMSloc(X)
%Syntax: loc=LMSloc(X)
%_____________________
%
% Calculates the Least Median of Squares (LMS) location parameter
% of the columns of a matrix X. If X is a vector, it returns the LMS
% location parameter of its components. If X is a scalar, it returns
% X.
%
% loc is the LMS estimated vector of locations.
% X is the matrix with the data sets.
%
% Reference:
% Rousseeuw PJ, Leroy AM (1987): Robust regression and outlier detection. Wiley.
%
%
% Alexandros Leontitsis
% Institute of Mathematics and Statistics
% University of Kent at Canterbury
% Canterbury
% Kent, CT2 7NF
% U.K.
%
% University e-mail: al10@ukc.ac.uk (until December 2002)
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% Sep 3, 2001.

if nargin<1 | isempty(X)==1
   error('Not enough input arguments.');
else
   % X must be 2-dimensional
   if ndims(X)>2
      error('Invalid data set.');
   end
   % If X is a row vector make it a column vector
   if size(X,1)==1
      X=X';
   end

   % n is the length of the data set
   n=size(X,1);
end

% For a single data point there is no need to proceed
if n==1
   loc=X;
else
   % Sort the data set to ascending order
   X=sort(X);

   % P. 169: the length of the "half" of the data points
   h=floor(n/2)+1;

   % P. 169: determine the shortest half and compute its midpoint
   j=1:n-h+1;
   range=X(j+h-1,:)-X(j,:);

   % Calculate the location parameter for every column of X
   for i=1:size(X,2);
      f=find(range(:,i)==min(range(:,i)));
      loc(i)=median((X(f+h-1,i)+X(f,i))/2);
   end
end


text:
function [N,Z,M,A,XYZ] = spm_max(X,L)
% Sizes, maxima and locations of local excursion sets
% FORMAT [N Z M A XYZ] = spm_max(X,L)
% X     - values of 3-D field
% L     - locations [x y z]' {in voxels}
%
% N     - size of region {in voxels)
% Z     - Z values of maxima
% M     - location of maxima {in voxels}
% A     - region number
% XYZ   - cell array of voxel locations
%__________________________________________________________________________
%
% spm_max characterizes a point list of voxel values (X) and their
% locations (L) in terms of edge, face and vertex connected subsets,
% returning a maxima- orientated list:  The value of the ith maximum is
% Z(i) and its location is given by M(:,i). A(i) identifies the ith
% maximum with a region. Region A(i) contains N(i) voxels, whose
% coordinates are in a 3-by-N(i) array in XYZ{i}.
%
% See also: spm_bwlabel.m and spm_clusters.m
%__________________________________________________________________________
% Copyright (C) 2003-2011 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: spm_max.m 4384 2011-07-06 17:00:20Z guillaume $

if isempty(L)
    N = []; Z = []; M = []; A = []; XYZ = [];
    return;
end

%-Ensure that L contains exactly integers
%--------------------------------------------------------------------------
L          = round(L);

%-Turn location list to binary 3D volume
%--------------------------------------------------------------------------
dim        = [max(L(1,:)) max(L(2,:)) max(L(3,:))];
vol        = zeros(dim(1),dim(2),dim(3));
index      = sub2ind(dim,L(1,:)',L(2,:)',L(3,:)');
vol(index) = 1;

%-Label each cluster with its own label using an 18 connectivity criterion
% cci = connected components image volume
%--------------------------------------------------------------------------
[cci,num]  = spm_bwlabel(vol,18);

%-Get size (in no. of voxels) for each connected component
% ccs = connected component size
%--------------------------------------------------------------------------
ccs        = histc(cci(:),(0:num) + 0.5);
ccs        = ccs(1:end-1);

%-Get indices into L for voxels that are indeed local maxima (using an 18
% neighbour criterion)
%--------------------------------------------------------------------------
vol(index) = X;
Lindex     = spm_get_lm(vol,L);

M          = L(:,Lindex);
Z          = X(Lindex);
mindex     = sub2ind(dim,L(1,Lindex)',L(2,Lindex)',L(3,Lindex)');
A          = cci(mindex);
N          = ccs(A);

%-Cell array of XYZ locations of voxels in each cluster
%--------------------------------------------------------------------------
if nargout > 4
    xyz(:,index) = sparse(L);
    cci   = sparse(cci(:));
    XYZ = cell(1, max(A));
    for i = 1:max(A)
        XYZ{i} = full(xyz(:,cci == i));
    end
end

author: spm


text:
function fx = p14_fx ( x )

%*****************************************************************************80
%
%% P14_FX evaluates the Camel.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX(*), the value of the function at X.
%
  fx =   1.0 ./ ( ( x - 0.3 ).^2 + 0.01 ) ...
       + 1.0 ./ ( ( x - 0.9 ).^2 + 0.04 ) ...
       + 2.0 * x - 5.2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef SuperEllipseExponentContourPlotter < handle

    properties (Access = private)
        fig
    end

    properties (Access = private)
        fileName
        titleF
        axisAdder
    end

    methods (Access = public)

        function obj = SuperEllipseExponentContourPlotter(cParams)
            obj.init(cParams);
        end

        function plot(obj,x,y,z)
            obj.plotMesh(x,y,z)
            obj.addAxis();
            obj.addTitle();
            obj.print();
        end

        function print(obj)
            fp = contourPrinter(obj.fig);
            filePath = obj.fileName;
            fp.print(filePath);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.fileName = [cParams.fileName,'Contour'];
            obj.titleF   = cParams.title;
            obj.axisAdder = cParams.axisAdder;
        end

        function plotMesh(obj,x,y,z)
            obj.fig = figure();
            connec = obj.obtainConnec(x,y);
            h = trisurf(connec,x,y,z);
            colorbar
            h.FaceColor = 'interp';
            h.EdgeColor = 'black';
            h.EdgeAlpha = 0.3;
            h.FaceAlpha = 1;
            hold on
            plot3(x,y,10*abs(z),'ro','MarkerFaceColor','red','MarkerSize',3)
            view(2)
        end

        function connec = obtainConnec(obj,x,y)
            connec = delaunay(x,y);
            connec = obj.obtainQualityElements(connec,x,y);
        end

        function addAxis(obj)
            obj.axisAdder.add();
        end

        function addTitle(obj)
            tN = ['\textrm{',obj.titleF,'optimal super-ellipse exponent}'];
          %  title(['$',tN,'$'],'interpreter','latex')
        end

    end

    methods (Access = private, Static)

        function newConnec = obtainQualityElements(connec,x,y)
            s.x = x;
            s.y = y;
            s.connec = connec;
            c = ConnecWithQualityComputer(s);
            newConnec = c.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function A = create_3d_linear_interp_matrix(T, varargin)
% A = create_2d_bilinear_interp_matrix(T, [szin])
% T should contain displacements in pixels
% size(T) = [N1, N2, N3, 2]
%
% x_def = A * x(:);

sz = [size(T, 1), size(T, 2), size(T, 3)];

if numel(varargin) == 1
    szin = varargin{1};
else
    szin = sz;
end

npixin = szin(1) * szin(2) * szin(3);
npixout = sz(1) * sz(2) * sz(3);

[m1, m2, m3] = ndgrid(1:sz(1), 1:sz(2), 1:sz(3));
TT = T + cat(4, m1, m2, m3);

r = cell(3, 1);
rk = cell(3, 1);
for i = 1 : 3
    r{i} = floor(TT(:, :, :, i));
    rk{i} = TT(:,:,:, i) - r{i};
end

%
ssh =  [0, 0, 0, 0, 1, 1, 1, 1;...
        0, 1, 0, 1, 0, 1, 1, 0;...
        0, 0, 1, 1, 0, 0, 1, 1;];
A = sparse(npixout, npixin);
for i = 1 : 8
    t = 1 : npixout;
    s = ones(numel(rk{1}), 1);
    rid = cell(3, 1);
    for j = 1 : 3
        rid{j} = r{j} + ssh(j, i);
        if ssh(j, i) == 1
            tmps = rk{j};
        else
            tmps = 1 - rk{j};
        end
        s = s .* tmps(:);
    end
    idx = (rid{1} <= szin(1)) & (rid{1} >= 1) & (rid{2} <= szin(2)) & (rid{2} >= 1) & (rid{3} <= szin(3)) & (rid{3} >= 1);
    for j = 1 : 3
        rid{j} = rid{j}(idx);
    end
    s = s(idx);
    t = t(idx);
    tmp = sparse(t, rid{1} + szin(1) * (rid{2} - 1) + szin(2)*szin(1) * (rid{3} - 1), s, npixout, npixin);
    A = A + tmp;
end

end

repo: pTVreg


text:
% DEMSPGP1D1 Do a simple 1-D regression after Snelson & Ghahramani's example.

% GP

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'spgp1d';
experimentNo = 1;

% load data
[X, y] = mapLoadData(dataSetName);

% Set up model
options = gpOptions('dtc');
options.numActive = 9;

% use the deterministic training conditional.
q = size(X, 2);
d = size(y, 2);

model = gpCreate(q, d, X, y, options);
model.X_u = randn(9, 1)*0.25 - 0.75;
params = gpExtractParam(model);
model = gpExpandParam(model, params);

% Optimise the model.
iters = 1000;
display = 1;

model = gpOptimise(model, display, iters);

% Save the results.
capName = dataSetName;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

demSpgp1dPlot

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef IsosurfaceOptionsDialog < handle
%ISOSURFACEOPTIONSDIALOG Dialog for 3D isosurfaces of intensity images
%
%   Class IsosurfaceOptionsDialog
%
%   Example
%   IsosurfaceOptionsDialog
%
%   See also
%     IsosurfaceOptionsDialog
%

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

%% Properties
properties
    IsosurfaceValue = 0;

    Parent;

    Handles;
end % end properties

%% Constructor
methods
    function obj = IsosurfaceOptionsDialog(parent, varargin)
    % Constructor for IsosurfaceOptionsDialog class

         % call parent constructor to initialize members
        obj = obj@handle();
        obj.Parent = parent;

        % pre-compute initial value of isosurface (from Otsu threshold)
        obj.IsosurfaceValue = imOtsuThreshold(obj.Parent.ImageData);

        % create default figure
        fig = figure(...
            'MenuBar', 'none', ...
            'NumberTitle', 'off', ...
            'HandleVisibility', 'On', ...
            'Name', 'Isosurface Options', ...
            'CloseRequestFcn', @obj.close);
        obj.Handles.Figure = fig;

        setupLayout(fig);

        function setupLayout(hf)
            % Setup layout of figure widgets

            pos = get(hf, 'Position');
            pos(3:4) = [200 180];
            set(hf, 'Position', pos);
            set(hf, 'Units', 'normalized');

            if verLessThan('matlab', 'R2014b')
                mainPanel = uiextras.VBox('Parent', hf);
            else
                mainPanel = uix.VBox('Parent', hf);
            end

            % control panel: use vertical layout
            if verLessThan('matlab', 'R2014b')
                optionsPanel = uiextras.VButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [150 30], ...
                    'VerticalAlignment', 'top', ...
                    'HorizontalAlignment', 'center', ...
                    'Spacing', 5, 'Padding', 5);
            else
                optionsPanel = uix.VButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [150 30], ...
                    'VerticalAlignment', 'top', ...
                    'HorizontalAlignment', 'center', ...
                    'Spacing', 5, 'Padding', 5);
            end

            obj.Handles.ThresholdText = uicontrol(...
                'Style', 'text', 'Parent', optionsPanel, ...
                'String', 'Isosurface Value:', ...
                'HorizontalAlignment', 'Left');

            obj.Handles.IsosurfaceValue = uicontrol(...
                'Style', 'edit', 'Parent', optionsPanel, ...
                'String', num2str(obj.IsosurfaceValue), ...
                'Value', 1, ...
                'HorizontalAlignment', 'Left');

            obj.Handles.ReverseZAxis = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Reverse Z-axis', ...
                'Value', 1, ...
                'HorizontalAlignment', 'Left');

            obj.Handles.RotateOx = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Rotate around X-axis', ...
                'Value', 1, ...
                'HorizontalAlignment', 'Left');

            obj.Handles.Smooth = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Smooth', ...
                'Value', 1, ...
                'HorizontalAlignment', 'Left');

            obj.Handles.ShowAxisLabel = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Show axes label', ...
                'Value', 0, ...
                'HorizontalAlignment', 'Left');

            % add control panel with "Compute" and "Close" buttons
            if verLessThan('matlab', 'R2014b')
                buttonsPanel = uiextras.HButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [70 30], ...
                    'VerticalAlignment', 'top', ...
                    'Spacing', 5, 'Padding', 5, ...
                    'Units', 'normalized', ...
                    'Position', [0 0 1 1]);
            else
                buttonsPanel = uix.HButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [70 30], ...
                    'VerticalAlignment', 'top', ...
                    'Spacing', 5, 'Padding', 5, ...
                    'Units', 'normalized', ...
                    'Position', [0 0 1 1]);
            end

            obj.Handles.ApplyButton = uicontrol('Style', 'PushButton', ...
                'Parent', buttonsPanel, ...
                'String', 'Compute', ...
                'Enable', 'on', ...
                'Callback', @obj.onApplyButtonClicked);
            obj.Handles.CloseButton = uicontrol('Style', 'PushButton', ...
                'Parent', buttonsPanel, ...
                'String', 'Close', ...
                'Callback', @obj.close);

            if verLessThan('matlab', 'R2014b')
                mainPanel.Sizes = [-1 40];
            else
                mainPanel.Heights = [-1 40];
            end

        end
    end

end % end constructors

%% Methods
methods
    function onApplyButtonClicked(obj, hObject, eventdata) %#ok<INUSD>

        % first, disable "compute" button to avoid multiple calls
        set(obj.Handles.ApplyButton, 'enable', 'off');
        hDlg = msgbox({'Computing isosurfaces,', 'Please wait...'}, ...
            'Isosurface computing');

        % extract options from widgets
        str = get(obj.Handles.IsosurfaceValue, 'String');
        [value, flag] = str2num(str);
        if ~flag
            error(['Could not parse isosurface value from: ' str]);
        end
        reverseZAxis = get(obj.Handles.ReverseZAxis, 'Value');
        rotateXAxis = get(obj.Handles.RotateOx, 'Value');
        smooth = get(obj.Handles.Smooth, 'Value');
        showAxesLabel = get(obj.Handles.ShowAxisLabel, 'Value');

        % get image data stored in parent Slicer
        imgData = obj.Parent.ImageData;
        imgSize = obj.Parent.ImageSize;

        % avoid the case of empty image
        if isempty(imgData)
            return;
        end

        % eventually rotate image around X-axis
        if rotateXAxis
            imgData = stackRotate90(imgData, 'x', 1);
            imgSize = imgSize([1 3 2]);
        end

        % compute display settings
        spacing = obj.Parent.VoxelSize;
        origin  = obj.Parent.VoxelOrigin;

        % compute grid positions
        lx = (0:imgSize(1)-1) * spacing(1) + origin(1);
        ly = (0:imgSize(2)-1) * spacing(2) + origin(2);
        lz = (0:imgSize(3)-1) * spacing(3) + origin(3);

        if rotateXAxis
            ly = ly(end:-1:1);
        end

        % determine the color map to use (default is empty)
        cmap = [];
        if ~isColorStack(imgData) && ~isempty(obj.Parent.ColorMap)
            cmap = obj.Parent.ColorMap;
        end
        if isempty(cmap)
            cmap = jet(256);
        end

        % create figure
        hf = figure(); hold on;
        set(hf, 'renderer', 'opengl');

        % compute isosurface
        im = imgData >= value;
        if ~any(im)
            return;
        end

        if smooth
            im = imGaussianFilter(double(im), [5 5 5], 2);
        end

        % display isosurface
        p = patch(isosurface(lx, ly, lz, im, .5));

        % set face color
        set(p, 'FaceColor', cmap(end,:), 'EdgeColor', 'none');

        % compute display extent (add a 0.5 limit around each voxel)
        extent = stackExtent(imgSize, spacing, origin);

        % setup display
        axis equal;
        axis(extent);
        view(3);
        axis('vis3d');

        light;

        if rotateXAxis || reverseZAxis
            set(gca, 'zdir', 'reverse');
        end

        if showAxesLabel
            xlabel('X axis');
            if rotateXAxis
                ylabel('Z Axis');
                zlabel('Y Axis');
            else
                ylabel('Y Axis');
                zlabel('Z Axis');
            end
        end

        % re-enable "apply" button
        set(obj.Handles.ApplyButton, 'enable', 'on');
        if ishandle(hDlg)
            close(hDlg);
        end

    end

end % end methods

%% Figure management
methods
    function close(obj, varargin)
        delete(obj.Handles.Figure);
    end

end

end % end classdef

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function isWrite(I, fName)
	typ = whos('I');
	typ = typ.class;
	f = fopen(fName, 'w');
	getMatlabStr(typ);
	fwrite(f, [1, size(I,2), size(I,1), size(I, 3), getMatlabStr(typ)], 'uint32');
	I = permute(I, [3 2 1]);
	fwrite(f, I(:), typ);
	fclose(f);
end

function i = getMatlabStr(str)
	switch str
		case 'single', i = 0;
		case 'double', i = 1;
		case 'uint8', i = 2;
    	case 'int8', i = 3;
    	case 'uint16', i = 4;
    	case 'int16', i = 5;
    	case 'uint32', i = 6;
    	case 'int32', i = 7;
    	case 'uint64', i = 8;
    	case 'int64', i = 9;
	end
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function e = rmModelSearchFit_twoGaussianDoGfixed(p,Y,Xv,Yv,stim,t,sigmaRatio,betaRatioAlpha)
% rmModelSearchFit_twoGaussianDoGfixed - actual fit function of rmSearchFit
%
% error = rmModelSearchFit(_twoGaussianDoGfixed(p,Y,trends,Xgrid,YGrid,stimulusMatrix);
%
% Basic barebones fit of a single time-series. Error is returned in
% percentage: 100% is RSS of unfitted time-series. This way we can quantify
% the improvement of the fit independend of the variation in the raw
% time-series.
%
% 2006/06 SOD: wrote it.
% 2006/12 SOD: modifications for fmincon, this is litterally called >10000
% times so we cut every corner possible.
% 2009/12 SOD & WZ: adapted for dog fixed model

% make RF (taken from rfGaussian2d)
Xv = Xv - p(1);   % positive x0 moves center right
Yv = Yv - p(2);   % positive y0 moves center up
RF = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(p(3).^2)) );

% make surround
%s2 = p(3).*sigmaRatio;
[tmp eccentricity] = cart2pol(p(1), p(2));
s2 = p(3).*sigmaRatio(1) + eccentricity.*sigmaRatio(2) + sigmaRatio(3);
RF2 = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(s2.^2)) );

% full pRF
betaRatio = (p(3)./s2).^betaRatioAlpha;
RF = RF - betaRatio.*RF2;

% make prediction (taken from rfMakePrediction)
X = [stim*RF t];

% fit - inlining pinv
%b = pinv(X)*Y;
[U,S,V] = svd(X,0);
s = diag(S);
tol = numel(X) * eps(max(s));
r = sum(s > tol);
if (r == 0)
    pinvX = zeros(size(X'));
else
    s = diag(ones(r,1)./s(1:r));
    pinvX = V(:,1:r)*s*U(:,1:r)';
end
b = pinvX*Y;

% compute residual sum of squares (e)
% e = norm(Y - X*abs(b));
if b(1)>0,
    e = norm(Y - X*b);
else
    e = norm(Y).*(1+sum(abs(b(1))));
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function hdrv = hdrvopen(hdrv)
%
%        hdrv = hdrvopen(hdrv)
%
%
%        Input:
%           -hdrv: a HDR video structure.
%
%        Output:
%           -hdrv: a HDR video structure.
%
%        This function opens the video stream for reading frames
%
%     Copyright (C) 2013-17  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(hdrv.streamOpen == 0)
    hdrv.streamOpen = 1;
end

end

author: banterle


text:
function [MatingPool,offspringLoc] = MatingSelection(B,s)
% The mating selection of EAG-MOEA/D

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

    [N,T] = size(B);

    %% Select N subproblems by roulette-wheel selection
    S = sum(s,2) + 1e-6;
    D = S./sum(S) + 0.002;
    D = D./sum(D);
    offspringLoc = RouletteWheelSelection(N,1./D);

    %% Select two parents in each selected subproblem
    MatingPool = zeros(1,2*N);
    for i = 1 : N
        MatingPool([i,i+N]) = B(offspringLoc(i),randperm(T,2));
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function transforms2d(varargin)
%TRANSFORMS2D Description of functions operating on transforms.
%
%   By 'transform' we mean an affine transform. A planar affine transform
%   can be represented by a 3x3 matrix.
%
%   Example
%     % create a translation by the vector [10 20]:
%     T = createTranslation([10 20])
%     T =
%          1     0    10
%          0     1    20
%          0     0     1
%
%     % apply a rotation on a polygon
%     poly = [0 0; 30 0;30 10;10 10;10 20;0 20];
%     trans = createRotation([10 20], pi/6);
%     polyT = transformPoint(poly, trans);
%     % display the original and the rotated polygons
%     figure; hold on; axis equal; axis([-10 40 -10 40]);
%     drawPolygon(poly, 'k');
%     drawPolygon(polyT, 'b');
%
%
%   See also
%   createTranslation, createRotation, createRotation90, createScaling
%   createHomothecy, createLineReflection, createBasisTransform
%   transformPoint, transformVector, transformLine, transformEdge
%   rotateVector, principalAxesTransform, fitAffineTransform2d
%   polynomialTransform2d, fitPolynomialTransform2d

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2008-10-13, using Matlab 7.4.0.287 (R2007a)
% Copyright 2008-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

help('transforms2d');


text:
function MeshMat = in_tess_simnibs(MeshFile, FileFormat)
% IN_TESS_SIMNIBS: Reads a 3D mesh from a gmsh4 file generated with SimNIBS

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Takfarinas Medani, Francois Tadel, 2022-2023

% Read mesh
m = mesh_load_gmsh4(MeshFile);

% Get file name
[fPath, fBase, fExt] = bst_fileparts(MeshFile);

% Convert to bst format
MeshMat = db_template('femmat');
MeshMat.Comment  = fBase;
MeshMat.Vertices = m.nodes(:,1:3);
MeshMat.Elements = double(m.tetrahedra(:,1:4));
MeshMat.Tissue   = double(m.tetrahedron_regions);

% Swap tetrahedrons orientation
MeshMat.Elements = MeshMat.Elements(:, [2 1 3 4]);

switch (FileFormat)
    case 'SIMNIBS3'
        % Replace the eyes with scalp (not used for now)
        % MeshMat.Tissue(MeshMat.Tissue==6) = 5;

        % Default tissue labels
        switch length(unique(MeshMat.Tissue))
            case 3
                MeshMat.TissueLabels = {'brain', 'skull', 'scalp'};
            case 4
                MeshMat.TissueLabels = {'brain', 'csf', 'skull', 'scalp'};
            case 5
                MeshMat.TissueLabels = {'white', 'gray', 'csf', 'skull', 'scalp'};
            otherwise
                uniqueLabels = unique(MeshMat.Tissue);
                for i = 1:length(uniqueLabels)
                     MeshMat.TissueLabels{i} = num2str(uniqueLabels(i));
                end
        end

    case 'SIMNIBS4'
        % Relabel the electrodes and gel
        MeshMat.Tissue(MeshMat.Tissue==100) = 11;
        MeshMat.Tissue(MeshMat.Tissue==500) = 12;
        % Default tissue labels (from file: final_tissues_LUT.txt)
        MeshMat.TissueLabels = {'white', 'gray', 'csf', 'skull', 'scalp', 'eyes', 'compact', 'spongy', 'blood', 'muscle', 'electrode', 'gel'};
end

author: brainstorm-tools


text:
for num =1:50
clearvars -except num;
warning off;
addpath(genpath(cd));

path1 = ['../road/ir/',num2str(num),'.jpg'];
path2 = ['../road/vi/',num2str(num),'.jpg'];
fused_path = ['result/',num2str(num),'.bmp'];

I=double(imread(path1))/255;
V=double(imread(path2))/255;

% image_left = ['../road/ir/',num2str(i),'.jpg'];
% image_right = ['../road/vi/',num2str(i),'.jpg'];
% fused_path = ['result/',num2str(i),'.bmp'];
%
% x{1}=imread(image_left);
% x{2}=imread(image_right);

calc_metric = 1; % Calculate the metrices is time consuming, it is used for quantitative evaluation. Set it to 0 if you do not want to do it.

%%
%The proposed GTF
nmpdef;
pars_irn = irntvInputPars('l1tv');

pars_irn.adapt_epsR   = 1;
pars_irn.epsR_cutoff  = 0.01;   % This is the percentage cutoff
pars_irn.adapt_epsF   = 1;
pars_irn.epsF_cutoff  = 0.05;   % This is the percentage cutoff
pars_irn.pcgtol_ini = 1e-4;
pars_irn.loops      = 5;
pars_irn.U0         = I-V;
pars_irn.variant       = NMP_TV_SUBSTITUTION;
pars_irn.weight_scheme = NMP_WEIGHTS_THRESHOLD;
pars_irn.pcgtol_ini    = 1e-2;
pars_irn.adaptPCGtol   = 1;

tic;
U = irntv(I-V, {}, 4, pars_irn);
t0=toc;

X=U+V;
X=im2gray(X);
imwrite(X,fused_path);
% imwrite(X,['F/GTF/',num2str(num),'.png'],'png');
% if calc_metric, Result = Metric(uint8(abs(I)*255),uint8(abs(V)*255),uint8(abs(X*255))); end

%%
%The laplacian pyramid as a compact image code(1983)
% level=4;
% tic;
% X1 = lp_fuse(I, V, level, 3, 3);       %LP
% t1=toc;
% X1=im2gray(X1);
% % imwrite(X1,['F/1/',num2str(num),'.png'],'png');
% if calc_metric, Result1 = Metric(uint8(abs(I)*255),uint8(abs(V)*255),uint8(abs(X1*255))); end

%%
%Image fusion by a ratio of low pass pyramid(1989)
% tic;
% X2 = rp_fuse(I, V, level, 3, 3);      %RP
% t2=toc;
% X2=im2gray(X2);
% imwrite(X2,['F/2/',num2str(num),'.png'],'png');
% if calc_metric, Result2 = Metric(uint8(abs(I)*255),uint8(abs(V)*255),uint8(abs(X2*255))); end

%%
% Wavelet
% fusion by taking the mean for both approximations and details
% tic;
% X3 = wfusimg(I,V,'db2',5,'mean','mean');
% X3=im2gray(X3);
% t3=toc;
% imwrite(X3,['F/3/',num2str(num),'.png'],'png');
% imwrite(X3,fused_path);
% if calc_metric, Result3 = Metric(uint8(abs(I)*255),uint8(abs(V)*255),uint8(abs(X3*255))); end

%%
%Pixel-and region-based image fusion with complex wavelets(2007)
% [M,N]=size(I);
% I4=imresize(I,[M+mod(M,2) N+mod(N,2)]);
% V4=imresize(V,[M+mod(M,2) N+mod(N,2)]);
% tic;
% X4 = dtcwt_fuse(I4, V4,level);           %DTCWT
% t4=toc;
% X4=im2gray(X4);
% imwrite(X4,['F/4/',num2str(num),'.png'],'png');
% if calc_metric, Result4 = Metric(uint8(abs(I4)*255),uint8(abs(V4)*255),uint8(abs(X4*255))); end
%
% %%
% %Remote sensing image fusion using the curvelet transform(2007)
% tic;
% X5 = curvelet_fuse(I4, V4,level+1);      %CVT
% t5=toc;
% X5=im2gray(X5);
% imwrite(X5,['F/5/',num2str(num),'.png'],'png');
% if calc_metric, Result5 = Metric(uint8(abs(I4)*255),uint8(abs(V4)*255),uint8(abs(X5*255))); end

%%
%Image Fusion technique using Multi-resolution singular Value decomposition(2011)
%apply MSVD
% tic;
% [Y1, U1] = MSVD(I4);
% [Y2, U2] = MSVD(V4);
%
% %fusion starts
% X6.LL = 0.5*(Y1.LL+Y2.LL);
%
% D  = (abs(Y1.LH)-abs(Y2.LH)) >= 0;
% X6.LH = D.*Y1.LH + (~D).*Y2.LH;
% D  = (abs(Y1.HL)-abs(Y2.HL)) >= 0;
% X6.HL = D.*Y1.HL + (~D).*Y2.HL;
% D  = (abs(Y1.HH)-abs(Y2.HH)) >= 0;
% X6.HH = D.*Y1.HH + (~D).*Y2.HH;
%
% %XX = [X.LL, X.LH; X.HL, X.HH];
% U = 0.5*(U1+U2);
%
% %apply IMSVD
% X6 = IMSVD(X6,U);
% t6=toc;
% X6=im2gray(X6);
% imwrite(X6,fused_path);
% if calc_metric, Result6 = Metric(uint8(abs(I4)*255),uint8(abs(V4)*255),uint8(abs(X6*255))); end
%
% %%
% %Image Fusion with Guided Filtering(2013)
% %run('F:\Code\Lichang\16.Image fusion total variation\Image fusion with guided filtering\Demo.m');
% I7=load_images('.\img',1);% the folder of source image
% tic;
% X7 = double(GFF(I7,5,10^-6,5,10^-6));
% %X7=rgb2gray(X7);
% t7=toc;
% imwrite(X7,['F/7/',num2str(num),'.png'],'png');
% %if calc_metric, Result7 = Metric(uint8(abs(I)*255),uint8(abs(V)*255),uint8(X7)); end
%
% %%
% %A general framework for image fusion based on multi-scale transform and sparse representation(2014)
% overlap = 6;
% epsilon=0.1;
% level=4;
% load('D_100000_256_8.mat');
% tic;
% X8= lp_sr_fuse(I,V,level,3,3,D,overlap,epsilon);      %LP-SR
% t8=toc;
% X8=im2gray(X8);
% imwrite(X8,['F/8/',num2str(num),'.png'],'png');
% if calc_metric, Result8 = Metric(uint8(abs(I)*255),uint8(abs(V)*255),uint8(abs(X8*255))); end
end

author: thfylsty


text:
function [x, cost, info, options] = trustregions(problem, x, options)
% Riemannian trust-regions solver for optimization on manifolds.
%
% function [x, cost, info, options] = trustregions(problem)
% function [x, cost, info, options] = trustregions(problem, x0)
% function [x, cost, info, options] = trustregions(problem, x0, options)
% function [x, cost, info, options] = trustregions(problem, [], options)
%
% This is the Riemannian Trust-Region solver (with tCG inner solve), named
% RTR. This solver will attempt to minimize the cost function described in
% the problem structure. It requires the availability of the cost function
% and of its gradient. It will issue calls for the Hessian. If no Hessian
% nor approximate Hessian is provided, a standard approximation of the
% Hessian based on the gradient will be computed. If a preconditioner for
% the Hessian is provided, it will be used.
%
% For a description of the algorithm and theorems offering convergence
% guarantees, see the references below. Documentation for this solver is
% available online at:
%
% http://www.manopt.org/solver_documentation_trustregions.html
%
%
% The initial iterate is x0 if it is provided. Otherwise, a random point on
% the manifold is picked. To specify options whilst not specifying an
% initial iterate, give x0 as [] (the empty matrix).
%
% The two outputs 'x' and 'cost' are the last reached point on the manifold
% and its cost. Notice that x is not necessarily the best reached point,
% because this solver is not forced to be a descent method. In particular,
% very close to convergence, it is sometimes preferable to accept very
% slight increases in the cost value (on the order of the machine epsilon)
% in the process of reaching fine convergence. In practice, this is not a
% limiting factor, as normally one does not need fine enough convergence
% that this becomes an issue.
%
% The output 'info' is a struct-array which contains information about the
% iterations:
%   iter (integer)
%       The (outer) iteration number, or number of steps considered
%       (whether accepted or rejected). The initial guess is 0.
%	cost (double)
%       The corresponding cost value.
%	gradnorm (double)
%       The (Riemannian) norm of the gradient.
%	numinner (integer)
%       The number of inner iterations executed to compute this iterate.
%       Inner iterations are truncated-CG steps. Each one requires a
%       Hessian (or approximate Hessian) evaluation.
%	time (double)
%       The total elapsed time in seconds to reach the corresponding cost.
%	rho (double)
%       The performance ratio for the iterate.
%	rhonum, rhoden (double)
%       Regularized numerator and denominator of the performance ratio:
%       rho = rhonum/rhoden. See options.rho_regularization.
%	accepted (boolean)
%       Whether the proposed iterate was accepted or not.
%	stepsize (double)
%       The (Riemannian) norm of the vector returned by the inner solver
%       tCG and which is retracted to obtain the proposed next iterate. If
%       accepted = true for the corresponding iterate, this is the size of
%       the step from the previous to the new iterate. If accepted is
%       false, the step was not executed and this is the size of the
%       rejected step.
%	Delta (double)
%       The trust-region radius at the outer iteration.
%	cauchy (boolean)
%       Whether the Cauchy point was used or not (if useRand is true).
%   And possibly additional information logged by options.statsfun.
% For example, type [info.gradnorm] to obtain a vector of the successive
% gradient norms reached at each (outer) iteration.
%
% The options structure is used to overwrite the default values. All
% options have a default value and are hence optional. To force an option
% value, pass an options structure with a field options.optionname, where
% optionname is one of the following and the default value is indicated
% between parentheses:
%
%   tolgradnorm (1e-6)
%       The algorithm terminates if the norm of the gradient drops below
%       this. For well-scaled problems, a rule of thumb is that you can
%       expect to reduce the gradient norm by 8 orders of magnitude
%       (sqrt(eps)) compared to the gradient norm at a "typical" point (a
%       rough initial iterate for example). Further decrease is sometimes
%       possible, but inexact floating point arithmetic will eventually
%       limit the final accuracy. If tolgradnorm is set too low, the
%       algorithm may end up iterating forever (or at least until another
%       stopping criterion triggers).
%   maxiter (1000)
%       The algorithm terminates if maxiter (outer) iterations were executed.
%   maxtime (Inf)
%       The algorithm terminates if maxtime seconds elapsed.
%	miniter (3)
%       Minimum number of outer iterations (used only if useRand is true).
%	mininner (1)
%       Minimum number of inner iterations (for tCG).
%	maxinner (problem.M.dim() : the manifold's dimension)
%       Maximum number of inner iterations (for tCG).
%	Delta_bar (problem.M.typicaldist() or sqrt(problem.M.dim()))
%       Maximum trust-region radius. If you specify this parameter but not
%       Delta0, then Delta0 will be set to 1/8 times this parameter.
%   Delta0 (Delta_bar/8)
%       Initial trust-region radius. If you observe a long plateau at the
%       beginning of the convergence plot (gradient norm VS iteration), it
%       may pay off to try to tune this parameter to shorten the plateau.
%       You should not set this parameter without setting Delta_bar.
%	useRand (false)
%       Set to true if the trust-region solve is to be initiated with a
%       random tangent vector. If set to true, no preconditioner will be
%       used. This option is set to true in some scenarios to escape saddle
%       points, but is otherwise seldom activated.
%	kappa (0.1)
%       Inner kappa convergence tolerance.
%	theta (1.0)
%       Inner theta convergence tolerance.
%	rho_prime (0.1)
%       Accept/reject ratio : if rho is at least rho_prime, the outer
%       iteration is accepted. Otherwise, it is rejected. In case it is
%       rejected, the trust-region radius will have been decreased.
%       To ensure this, rho_prime must be strictly smaller than 1/4.
%   rho_regularization (1e3)
%       Close to convergence, evaluating the performance ratio rho is
%       numerically challenging. Meanwhile, close to convergence, the
%       quadratic model should be a good fit and the steps should be
%       accepted. Regularization lets rho go to 1 as the model decrease and
%       the actual decrease go to zero. Set this option to zero to disable
%       regularization (not recommended). See in-code for the specifics.
%   statsfun (none)
%       Function handle to a function that will be called after each
%       iteration to provide the opportunity to log additional statistics.
%       They will be returned in the info struct. See the generic Manopt
%       documentation about solvers for further information. statsfun is
%       called with the point x that was reached last, after the
%       accept/reject decision. See comment below.
%   stopfun (none)
%       Function handle to a function that will be called at each iteration
%       to provide the opportunity to specify additional stopping criteria.
%       See the generic Manopt documentation about solvers for further
%       information.
%   verbosity (2)
%       Integer number used to tune the amount of output the algorithm
%       generates during execution (mostly as text in the command window).
%       The higher, the more output. 0 means silent. 3 and above includes a
%       display of the options structure at the beginning of the execution.
%   debug (false)
%       Set to true to allow the algorithm to perform additional
%       computations for debugging purposes. If a debugging test fails, you
%       will be informed of it, usually via the command window. Be aware
%       that these additional computations appear in the algorithm timings
%       too.
%   storedepth (20)
%       Maximum number of different points x of the manifold for which a
%       store structure will be kept in memory in the storedb. If the
%       caching features of Manopt are not used, this is irrelevant. If
%       memory usage is an issue, you may try to lower this number.
%       Profiling may then help to investigate if a performance hit was
%       incured as a result.
%
% Notice that statsfun is called with the point x that was reached last,
% after the accept/reject decision. Hence: if the step was accepted, we get
% that new x, with a store which only saw the call for the cost and for the
% gradient. If the step was rejected, we get the same x as previously, with
% the store structure containing everything that was computed at that point
% (possibly including previous rejects at that same point). Hence, statsfun
% should not be used in conjunction with the store to count operations for
% example. Instead, you could use a global variable and increment that
% variable directly from the cost related functions. It is however possible
% to use statsfun with the store to compute, for example, alternate merit
% functions on the point x.
%
% See also: steepestdescent conjugategradient manopt/examples

% This file is part of Manopt: www.manopt.org.
% This code is an adaptation to Manopt of the original GenRTR code:
% RTR - Riemannian Trust-Region
% (c) 2004-2007, P.-A. Absil, C. G. Baker, K. A. Gallivan
% Florida State University
% School of Computational Science
% (http://www.math.fsu.edu/~cbaker/GenRTR/?page=download)
% See accompanying license file.
% The adaptation was executed by Nicolas Boumal.
%
% Change log:
%
%   NB April 3, 2013:
%       tCG now returns the Hessian along the returned direction eta, so
%       that we do not compute that Hessian redundantly: some savings at
%       each iteration. Similarly, if the useRand flag is on, we spare an
%       extra Hessian computation at each outer iteration too, owing to
%       some modifications in the Cauchy point section of the code specific
%       to useRand = true.
%
%   NB Aug. 22, 2013:
%       This function is now Octave compatible. The transition called for
%       two changes which would otherwise not be advisable. (1) tic/toc is
%       now used as is, as opposed to the safer way:
%       t = tic(); elapsed = toc(t);
%       And (2), the (formerly inner) function savestats was moved outside
%       the main function to not be nested anymore. This is arguably less
%       elegant, but Octave does not (and likely will not) support nested
%       functions.
%
%   NB Dec. 2, 2013:
%       The in-code documentation was largely revised and expanded.
%
%   NB Dec. 2, 2013:
%       The former heuristic which triggered when rhonum was very small and
%       forced rho = 1 has been replaced by a smoother heuristic which
%       consists in regularizing rhonum and rhoden before computing their
%       ratio. It is tunable via options.rho_regularization. Furthermore,
%       the solver now detects if tCG did not obtain a model decrease
%       (which is theoretically impossible but may happen because of
%       numerical errors and/or because of a nonlinear/nonsymmetric Hessian
%       operator, which is the case for finite difference approximations).
%       When such an anomaly is detected, the step is rejected and the
%       trust region radius is decreased.
%
%   NB Dec. 3, 2013:
%       The stepsize is now registered at each iteration, at a small
%       additional cost. The defaults for Delta_bar and Delta0 are better
%       defined. Setting Delta_bar in the options will automatically set
%       Delta0 accordingly. In Manopt 1.0.4, the defaults for these options
%       were not treated appropriately because of an incorrect use of the
%       isfield() built-in function.

% Verify that the problem description is sufficient for the solver.
if ~canGetCost(problem)
    warning('manopt:getCost', ...
            'No cost provided. The algorithm will likely abort.');
end
if ~canGetGradient(problem)
    warning('manopt:getGradient', ...
            'No gradient provided. The algorithm will likely abort.');
end
if ~canGetHessian(problem)
    warning('manopt:getHessian:approx', ...
            'No Hessian provided. Using an approximation instead.');
end

% Define some strings for display
tcg_stop_reason = {'negative curvature',...
                   'exceeded trust region',...
                   'reached target residual-kappa',...
                   'reached target residual-theta',...
                   'dimension exceeded',...
                   'model increased'};

% Set local defaults here
localdefaults.verbosity = 2;
localdefaults.maxtime = inf;
localdefaults.miniter = 3;
localdefaults.maxiter = 1000;
localdefaults.mininner = 1;
localdefaults.maxinner = problem.M.dim();
localdefaults.tolgradnorm = 1e-6;
localdefaults.kappa = 0.1;
localdefaults.theta = 1.0;
localdefaults.rho_prime = 0.1;
localdefaults.useRand = false;
localdefaults.rho_regularization = 1e3;

% Merge global and local defaults, then merge w/ user options, if any.
localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
if ~exist('options', 'var') || isempty(options)
    options = struct();
end
options = mergeOptions(localdefaults, options);

% Set default Delta_bar and Delta0 separately to deal with additional
% logic: if Delta_bar is provided but not Delta0, let Delta0 automatically
% be some fraction of the provided Delta_bar.
if ~isfield(options, 'Delta_bar')
    if isfield(problem.M, 'typicaldist')
        options.Delta_bar = problem.M.typicaldist();
    else
        options.Delta_bar = sqrt(problem.M.dim());
    end
end
if ~isfield(options,'Delta0')
    options.Delta0 = options.Delta_bar / 8;
end

% Check some option values
assert(options.rho_prime < 1/4, ...
        'options.rho_prime must be strictly smaller than 1/4.');
assert(options.Delta_bar > 0, ...
        'options.Delta_bar must be positive.');
assert(options.Delta0 > 0 && options.Delta0 < options.Delta_bar, ...
        'options.Delta0 must be positive and smaller than Delta_bar.');

% It is sometimes useful to check what the actual option values are.
if options.verbosity >= 3
    disp(options);
end

% Create a store database
storedb = struct();

tic();

% If no initial point x is given by the user, generate one at random.
if ~exist('x', 'var') || isempty(x)
    x = problem.M.rand();
end

%% Initializations

% k counts the outer (TR) iterations. The semantic is that k counts the
% number of iterations fully executed so far.
k = 0;

% initialize solution and companion measures: f(x), fgrad(x)
[fx fgradx storedb] = getCostGrad(problem, x, storedb);
norm_grad = problem.M.norm(x, fgradx);

% initialize trust-region radius
Delta = options.Delta0;

% Save stats in a struct array info, and preallocate
% (see http://people.csail.mit.edu/jskelly/blog/?x=entry:entry091030-033941)
if ~exist('used_cauchy', 'var')
    used_cauchy = [];
end
stats = savestats(problem, x, storedb, options, k, fx, norm_grad, Delta);
info(1) = stats;
info(min(10000, options.maxiter+1)).iter = [];

% ** Display:
if options.verbosity == 2
   fprintf(['%3s %3s      %5s                %5s     ',...
            'f: %e   |grad|: %e\n'],...
           '   ','   ','     ','     ', fx, norm_grad);
elseif options.verbosity > 2
   fprintf('************************************************************************\n');
   fprintf('%3s %3s    k: %5s     num_inner: %5s     %s\n',...
           '','','______','______','');
   fprintf('       f(x) : %e       |grad| : %e\n', fx, norm_grad);
   fprintf('      Delta : %f\n', Delta);
end

% **********************
% ** Start of TR loop **
% **********************
while true

	% Start clock for this outer iteration
    tic();

    % Run standard stopping criterion checks
    [stop reason] = stoppingcriterion(problem, x, options, info, k+1);

    % If the stopping criterion that triggered is the tolerance on the
    % gradient norm but we are using randomization, make sure we make at
    % least miniter iterations to give randomization a chance at escaping
    % saddle points.
    if stop == 2 && options.useRand && k < options.miniter
        stop = 0;
    end

    if stop
        if options.verbosity >= 1
            fprintf([reason '\n']);
        end
        break;
    end

    if options.verbosity > 2 || options.debug > 0
        fprintf('************************************************************************\n');
    end

    % *************************
    % ** Begin TR Subproblem **
    % *************************

    % Determine eta0
    if ~options.useRand
        % Pick the zero vector
        eta = problem.M.zerovec(x);
    else
        % Random vector in T_x M (this has to be very small)
        eta = problem.M.lincomb(x, 1e-6, problem.M.randvec(x));
        % Must be inside trust-region
        while problem.M.norm(x, eta) > Delta
            eta = problem.M.lincomb(x, sqrt(sqrt(eps)), eta);
        end
    end

    % solve TR subproblem
    [eta Heta numit stop_inner storedb] = ...
                     tCG(problem, x, fgradx, eta, Delta, options, storedb);
    srstr = tcg_stop_reason{stop_inner};

    % This is only computed for logging purposes, because it may be useful
    % for some user-defined stopping criteria. If this is not cheap for
    % specific application (compared to evaluating the cost), we should
    % reconsider this.
    norm_eta = problem.M.norm(x, eta);

    if options.debug > 0
        testangle = problem.M.inner(x, eta, fgradx) / (norm_eta*norm_grad);
    end

    % If using randomized approach, compare result with the Cauchy point.
    % Convergence proofs assume that we achieve at least the reduction of
    % the Cauchy point. After this if-block, either all eta-related
    % quantities have been changed consistently, or none of them have
    % changed.
    if options.useRand
        used_cauchy = false;
        % Check the curvature,
        [Hg storedb] = getHessian(problem, x, fgradx, storedb);
        g_Hg = problem.M.inner(x, fgradx, Hg);
        if g_Hg <= 0
            tau_c = 1;
        else
            tau_c = min( norm_grad^3/(Delta*g_Hg) , 1);
        end
        % and generate the Cauchy point.
        eta_c  = problem.M.lincomb(x, -tau_c * Delta / norm_grad, fgradx);
        Heta_c = problem.M.lincomb(x, -tau_c * Delta / norm_grad, Hg);

        % Now that we have computed the Cauchy point in addition to the
        % returned eta, we might as well keep the best of them.
        mdle  = fx + problem.M.inner(x, fgradx, eta) ...
                   + .5*problem.M.inner(x, Heta,   eta);
        mdlec = fx + problem.M.inner(x, fgradx, eta_c) ...
                   + .5*problem.M.inner(x, Heta_c, eta_c);
        if mdle > mdlec
            eta = eta_c;
            Heta = Heta_c; % added April 11, 2012
            used_cauchy = true;
        end
    end

	% Compute the retraction of the proposal
	x_prop  = problem.M.retr(x, eta);

	% Compute the function value of the proposal
	[fx_prop storedb] = getCost(problem, x_prop, storedb);

	% Will we accept the proposed solution or not?
    % Check the performance of the quadratic model against the actual cost.
    rhonum = fx - fx_prop;
    rhoden = -problem.M.inner(x, fgradx, eta) ...
             -.5*problem.M.inner(x, eta, Heta);

    % Heuristic -- added Dec. 2, 2013 (NB) to replace the former heuristic.
    % This heuristic is documented in the book by Conn Gould and Toint on
    % trust-region methods, section 17.4.2.
    % rhonum measures the difference between two numbers. Close to
    % convergence, these two numbers are very close to each other, so
    % that computing their difference is numerically challenging: there may
    % be a significant loss in accuracy. Since the acceptance or rejection
    % of the step is conditioned on the ratio between rhonum and rhoden,
    % large errors in rhonum result in a large error in rho, hence in
    % erratic acceptance / rejection. Meanwhile, close to convergence,
    % steps are usually trustworthy and we should transition to a Newton-
    % like method, with rho=1 consistently. The heuristic thus shifts both
    % rhonum and rhoden by a small amount such that far from convergence,
    % the shift is irrelevant and close to convergence, the ratio rho goes
    % to 1, effectively promoting acceptance of the step.
    % The rationale is that close to convergence, both rhonum and rhoden
    % are quadratic in the distance between x and x_prop. Thus, when this
    % distance is on the order of sqrt(eps), the value of rhonum and rhoden
    % is on the order of eps, which is indistinguishable from the numerical
    % error, resulting in badly estimated rho's.
    % For abs(fx) < 1, this heuristic is invariant under offsets of f but
    % not under scaling of f. For abs(fx) > 1, the opposite holds. This
    % should not alarm us, as this heuristic only triggers at the very last
    % iterations if very fine convergence is demanded.
    rho_reg = max(1, abs(fx)) * eps * options.rho_regularization;
    rhonum = rhonum + rho_reg;
    rhoden = rhoden + rho_reg;

    if options.debug > 0
        fprintf('DBG:     rhonum : %e\n', rhonum);
        fprintf('DBG:     rhoden : %e\n', rhoden);
    end

    % This is always true if a linear, symmetric operator is used for the
    % Hessian (approximation) and if we had infinite numerical precision.
    % In practice, nonlinear approximations of the Hessian such as the
    % built-in finite difference approximation and finite numerical
    % accuracy can cause the model to increase. In such scenarios, we
    % decide to force a rejection of the step and a reduction of the
    % trust-region radius. We test the sign of the regularized rhoden since
    % the regularization is supposed to capture the accuracy to which
    % rhoden is computed: if rhoden were negative before regularization but
    % not after, that should not be (and is not) detected as a failure.
    model_decreased = (rhoden >= 0);

    if ~model_decreased
        srstr = [srstr ', model did not decrease']; %#ok<AGROW>
    end

    rho = rhonum / rhoden;

    if options.debug > 0
        m = @(x, eta) ...
          getCost(problem, x, storedb) + ...
          getDirectionalDerivative(problem, x, eta, storedb) + ...
          .5*problem.M.inner(x, getHessian(problem, x, eta, storedb), eta);
        zerovec = problem.M.zerovec(x);
        actrho = (fx - fx_prop) / (m(x, zerovec) - m(x, eta));
        fprintf('DBG:   new f(x) : %e\n', fx_prop);
        fprintf('DBG: actual rho : %e\n', actrho);
        fprintf('DBG:   used rho : %e\n', rho);
    end

    % Choose the new TR radius based on the model performance
    trstr = '   ';
    % If the actual decrease is smaller than 1/4 of the predicted decrease,
    % then reduce the TR radius.
    if rho < 1/4 || ~model_decreased
        trstr = 'TR-';
        Delta = Delta/4;
    % If the actual decrease is at least 3/4 of the precicted decrease and
    % the tCG (inner solve) hit the TR boundary, increase the TR radius.
    elseif rho > 3/4 && (stop_inner == 1 || stop_inner == 2)
        trstr = 'TR+';
        Delta = min(2*Delta, options.Delta_bar);
    end
    % Otherwise, keep the TR radius constant.

    % Choose to accept or reject the proposed step based on the model
    % performance.
    if model_decreased && rho > options.rho_prime
        accept = true;
        accstr = 'acc';
        x = x_prop;
        fx = fx_prop;
        [fgradx storedb] = getGradient(problem, x, storedb);
        norm_grad = problem.M.norm(x, fgradx);
    else
        accept = false;
        accstr = 'REJ';
    end

    % Make sure we don't use too much memory for the store database
    storedb = purgeStoredb(storedb, options.storedepth);

    % k is the number of iterations we have accomplished.
    k = k + 1;

    % Log statistics for freshly executed iteration.
    % Everything after this in the loop is not accounted for in the timing.
    stats = savestats(problem, x, storedb, options, k, fx, norm_grad, ...
                      Delta, info, rho, rhonum, rhoden, accept, numit, ...
                      norm_eta, used_cauchy);
    info(k+1) = stats; %#ok<AGROW>

    % ** Display:
    if options.verbosity == 2,
        fprintf(['%3s %3s   k: %5d     num_inner: %5d     ', ...
        'f: %e   |grad|: %e   %s\n'], ...
        accstr,trstr,k,numit,fx,norm_grad,srstr);
    elseif options.verbosity > 2,
        if options.useRand && used_cauchy,
            fprintf('USED CAUCHY POINT\n');
        end
		fprintf('%3s %3s    k: %5d     num_inner: %5d     %s\n', ...
				accstr, trstr, k, numit, srstr);
		fprintf('       f(x) : %e     |grad| : %e\n',fx,norm_grad);
		if options.debug > 0
			fprintf('      Delta : %f          |eta| : %e\n',Delta,norm_eta);
		end
		fprintf('        rho : %e\n',rho);
    end
    if options.debug > 0,
        fprintf('DBG: cos ang(eta,gradf): %d\n',testangle);
        if rho == 0
            fprintf('DBG: rho = 0, this will likely hinder further convergence.\n');
        end
    end

end  % of TR loop (counter: k)

% Restrict info struct-array to useful part
info = info(1:k+1);

if (options.verbosity > 2) || (options.debug > 0),
   fprintf('************************************************************************\n');
end
if (options.verbosity > 0) || (options.debug > 0)
    fprintf('Total time is %f [s] (excludes statsfun)\n', info(end).time);
end

% Return the best cost reached
cost = fx;

end

% Routine in charge of collecting the current iteration stats
function stats = savestats(problem, x, storedb, options, k, fx, ...
                           norm_grad, Delta, info, rho, rhonum, ...
                           rhoden, accept, numit, norm_eta, used_cauchy)
    stats.iter = k;
    stats.cost = fx;
    stats.gradnorm = norm_grad;
    stats.Delta = Delta;
    if k == 0
        stats.time = toc();
        stats.rho = inf;
        stats.rhonum = NaN;
        stats.rhoden = NaN;
        stats.accepted = true;
        stats.numinner = NaN;
        stats.stepsize = NaN;
        if options.useRand
            stats.cauchy = false;
        end
    else
        stats.time = info(k).time + toc();
        stats.rho = rho;
        stats.rhonum = rhonum;
        stats.rhoden = rhoden;
        stats.accepted = accept;
        stats.numinner = numit;
        stats.stepsize = norm_eta;
        if options.useRand,
          stats.cauchy = used_cauchy;
        end
    end

    % See comment about statsfun above: the x and store passed to statsfun
    % are that of the most recently accepted point after the iteration
    % fully executed.
    stats = applyStatsfun(problem, x, storedb, options, stats);

end

author: andrewssobral

repo: lrslibrary


text:
classdef ResultsPrinterFactory < handle

    properties (Access = private)
        printer
    end

    methods (Access = public)
        function p = create(obj,resultCase,d)
            obj.createPrinter(resultCase,d);
            p = obj.printer;
        end

    end

    methods (Access = private)

        function createPrinter(obj,resultCase,d)
            switch resultCase
                case 'Elasticity'
                    p = ElasticityResultsPrinter(d);
                case 'ElasticityMicro'
                    p = ElasticityMicroResultsPrinter(d);
                case 'NumericalHomogenizer'
                    p = NumericalHomogenizerPrinter(d);
                case 'TopOptProblem'
                    p = TopOptResultsPrinter(d);
                case 'HomogenizedTensor'
                    p = HomogenizedTensorPrinter(d);
                case 'HomogenizedTensorStressBasis'
                    p = HomogenizedTensorStressBasisPrinter(d);
                case 'ElasticityMicroAndLevelSet'
                    p = ElasticityMicroDefinedByLevelSet(d);
                case {'ShapeFunction'}
                    p = ShapesPrinter(d);
                case {'PerimeterPrinter'}
                    p = PerimeterPrinter(d);
                case {'VectorGauss'}
                    d.fieldName = d.name;
                    p = VectorGaussVariablePrinter(d);
                case {'ScalarGauss'}
                    d.fieldName = d.name;
                    p = ScalarGaussVariablePrinter(d);
                case {'ScalarNodal'}
                    d.fieldName = d.name;
                    p = ScalarNodalVariablePrinter(d);
                case {'DensityGauss'}
                    p = DensityGaussResultsPrinter(d);
                case {'LevelSet'}
                    p = LevelSetResultsPrinter(d);
                otherwise
                    p = [];
            end
            obj.printer = p;
        end

    end

end


text:
function out = iszero( f )
%ISZERO   Check if a SEPARABLEAPPROX is identically zero on its domain.
%   OUT = ISZERO( F ) return 1 if the SEPARABLEAPPROX is exactly the zero function, and
%   0 otherwise.

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

% Get data:
pivots = f.pivotValues;
cols = f.cols;
rows = f.rows;

% Trivial check: If all the pivots are zero, then the SEPARABLEAPPROX is zero:
if ( norm(1./pivots, inf) == 0 )
    out = 1;
    return
end

% Quick check: Evaluate on a meshgrid. If the matrix is nonzero then the
% SEPARABLEAPPROX is nonzero.
dom = f.domain;
x = linspace(dom(1), dom(2), 10);
y = linspace(dom(3), dom(4), 10);
vals = fevalm(f, x, y);
if ( norm( vals, inf ) > 0 )
   out = 0;
   return
end

% Slower check: A pivot may be positive, but the columns or rows may be zero:
rk = length( pivots );
out_cols = zeros( length( pivots ), 1 );
out_rows = zeros( length( pivots ), 1 );
for j = 1:rk
    out_cols( j ) = iszero( cols(:, j) );
    out_rows( j ) = iszero( rows(:, j) );
end
bolslices = ( all(out_cols) || all(out_rows) );
out = bolslices;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun