text:
function exportModel(model,fileName,exportGeneComplexes,supressWarnings,sortIds)
% exportModel
%   Exports a constraint-based model to an SBML file (L3V1 FBCv2)
%
%   Input:
%   model               a model structure
%   fileName            filename to export the model to. A dialog window
%                       will open if no file name is specified.
%   exportGeneComplexes true if gene complexes (all gene sets linked with
%                       AND relationship) should be recognised and exported
%                       (opt, default false)
%   supressWarnings     true if warnings should be supressed (opt, default
%                       false)
%   sortIds             logical whether metabolites, reactions and genes
%                       should be sorted alphabetically by their
%                       identifiers (opt, default false)
%
%   Usage: exportModel(model,fileName,exportGeneComplexes,supressWarnings,sortIds)
if nargin<2 || isempty(fileName)
    [fileName, pathName] = uiputfile({'*.xml;*.sbml'}, 'Select file for model export',[model.id '.xml']);
    if fileName == 0
        error('You should provide a file location')
    else
        fileName = fullfile(pathName,fileName);
    end
end
fileName=char(fileName);
if nargin<3
    exportGeneComplexes=false;
end
if nargin<4
    supressWarnings=false;
end
if nargin<5
    sortIds=false;
end
if sortIds==true
    model=sortIdentifiers(model);
end

%If no subSystems are defined, then no need to use groups package
if isfield(model,'subSystems')
    modelHasSubsystems=true;
else
    modelHasSubsystems=false;
end

%The default SBML format settings, which are used as input for appropriate
%libSBML functions to generate the blank SBML model structure before using
%exporting in with OutputSBML to xml file
sbmlLevel=3;
sbmlVersion=1;
sbmlPackages={'fbc'};
sbmlPackageVersions=2;
if modelHasSubsystems
    sbmlPackages={sbmlPackages,'groups'};
    sbmlPackageVersions=[sbmlPackageVersions,1];
end

%Check if the "unconstrained" field is still present. This shows if
%exchange metabolites have been removed
if ~isfield(model,'unconstrained')
    model.unconstrained=zeros(numel(model.mets),1);
end

%If model id and name do not exist, make sure that default
%strings are included
if ~isfield(model,'id')
    fprintf('WARNING: The model is missing the "id" field. Uses "blankID". \n');
    model.id='blankID';
end
if ~isfield(model,'name')
    fprintf('WARNING: The model is missing the "name" field. Uses "blankName". \n');
    model.name='blankName';
end

%Check the model structure
if supressWarnings==false
    checkModelStruct(model,false);
end

%Add several blank fields, if they do not exist already. This is to reduce
%the number of conditions below
if ~isfield(model,'compMiriams')
    model.compMiriams=cell(numel(model.comps),1);
end
if ~isfield(model,'inchis')
    model.inchis=cell(numel(model.mets),1);
end
if ~isfield(model,'metFormulas')
    model.metFormulas=cell(numel(model.mets),1);
end
if ~isfield(model,'metMiriams')
    model.metMiriams=cell(numel(model.mets),1);
end
if ~isfield(model,'geneMiriams') && isfield(model,'genes')
    model.geneMiriams=cell(numel(model.genes),1);
end
if ~isfield(model,'geneShortNames') && isfield(model,'genes')
    model.geneShortNames=cell(numel(model.genes),1);
end
if ~isfield(model,'subSystems')
    model.subSystems=cell(numel(model.rxns),1);
end
if ~isfield(model,'eccodes')
    model.eccodes=cell(numel(model.rxns),1);
end
if ~isfield(model,'rxnReferences')
    model.rxnReferences=cell(numel(model.rxns),1);
end
if ~isfield(model,'rxnConfidenceScores')
    model.rxnConfidenceScores=NaN(numel(model.rxns),1);
end
if ~isfield(model,'rxnNotes')
    model.rxnNotes=cell(numel(model.rxns),1);
end
if ~isfield(model,'rxnMiriams')
    model.rxnMiriams=cell(numel(model.rxns),1);
end

if sbmlLevel<3
    %Check if genes have associated compartments
    if ~isfield(model,'geneComps') && isfield(model,'genes')
        if supressWarnings==false
            EM='There are no compartments specified for genes. All genes will be assigned to the first compartment. This is because the SBML structure requires all elements to be assigned to a compartment';
            dispEM(EM,false);
        end
        model.geneComps=ones(numel(model.genes),1);
    end
end

%Convert ids to SBML-convenient format. This is to avoid the data loss when
%unsupported characters are included in ids. Here we are using part from
%convertSBMLID, originating from the COBRA Toolbox
model.rxns=regexprep(model.rxns,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
model.mets=regexprep(model.mets,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
model.comps=regexprep(model.comps,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
if isfield(model,'genes')
    problemGenes=find(~cellfun('isempty',regexp(model.genes,'([^0-9_a-zA-Z])')));
    originalGenes=model.genes(problemGenes);
    replacedGenes=regexprep(model.genes(problemGenes),'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
    model.genes(problemGenes)=replacedGenes;
    for i=1:numel(problemGenes)
        model.grRules = regexprep(model.grRules, ['(^|\s|\()' originalGenes{i} '($|\s|\))'], ['$1' replacedGenes{i} '$2']);
    end
end

%Generate an empty SBML structure
modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
modelSBML.metaid=model.id;
modelSBML.id=regexprep(model.id,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
modelSBML.name=model.name;

if isfield(model,'annotation')
    if isfield(model.annotation,'note')
        modelSBML.notes=['<notes>
',regexprep(model.annotation.note,'
|
',''),'
</notes>'];
    end
else
    modelSBML.notes='<notes>
This file was generated using the exportModel function in RAVEN Toolbox 2 and OutputSBML in libSBML
</notes>';
end

modelSBML.annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.id '">'];
if isfield(model,'annotation')
    nameString='';
    if isfield(model.annotation,'familyName')
        if ~isempty(model.annotation.familyName)
            nameString=['<vCard:Family>' model.annotation.familyName '</vCard:Family>'];
        end
    end
    if isfield(model.annotation,'givenName')
        if ~isempty(model.annotation.givenName)
            nameString=[nameString '<vCard:Given>' model.annotation.givenName '</vCard:Given>'];
        end
    end
    email='';
    if isfield(model.annotation,'email')
        if ~isempty(model.annotation.email)
            email=['<vCard:EMAIL>' model.annotation.email '</vCard:EMAIL>'];
        end
    end
    org='';
    if isfield(model.annotation,'organization')
        if ~isempty(model.annotation.organization)
            org=['<vCard:ORG rdf:parseType="Resource"><vCard:Orgname>' model.annotation.organization '</vCard:Orgname></vCard:ORG>'];
        end
    end
    if ~isempty(nameString) || ~isempty(email) || ~isempty(org)
        modelSBML.annotation=[modelSBML.annotation '<dc:creator><rdf:Bag><rdf:li rdf:parseType="Resource">'];
        if ~isempty(nameString)
            modelSBML.annotation=[modelSBML.annotation '<vCard:N rdf:parseType="Resource">' nameString '</vCard:N>'];
        end
        modelSBML.annotation=[modelSBML.annotation email org '</rdf:li></rdf:Bag></dc:creator>'];
    end
end
modelSBML.annotation=[modelSBML.annotation '<dcterms:created rdf:parseType="Resource">'...
    '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF>'...
    '</dcterms:created>'...
    '<dcterms:modified rdf:parseType="Resource">'...
    '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF>'...
    '</dcterms:modified>'];

if isfield(model,'annotation')
    if isfield(model.annotation,'taxonomy')
        modelSBML.annotation=[modelSBML.annotation '<bqbiol:is><rdf:Bag><rdf:li rdf:resource="https://identifiers.org/taxonomy/' regexprep(model.annotation.taxonomy,'taxonomy/','') '"/></rdf:Bag></bqbiol:is>'];
    end
end
modelSBML.annotation=[modelSBML.annotation '</rdf:Description></rdf:RDF></annotation>'];

%Prepare compartments
for i=1:numel(model.comps)
    %Add the default values, as these will be the same in all entries
    if i==1
        if isfield(modelSBML.compartment, 'sboTerm')
            modelSBML.compartment(i).sboTerm=290;
        end
        if isfield(modelSBML.compartment, 'spatialDimensions')
            modelSBML.compartment(i).spatialDimensions=3;
        end
        if isfield(modelSBML.compartment, 'size')
            modelSBML.compartment(i).size=1;
        end
        if isfield(modelSBML.compartment, 'constant')
            modelSBML.compartment(i).constant=1;
        end
        if isfield(modelSBML.compartment, 'isSetSize')
            modelSBML.compartment(i).isSetSize=1;
        end
        if isfield(modelSBML.compartment, 'isSetSpatialDimensions')
            modelSBML.compartment(i).isSetSpatialDimensions=1;
        end
    end
    %Copy the default values to the next entry as long as it is not the
    %last one
    if i<numel(model.comps)
        modelSBML.compartment(i+1)=modelSBML.compartment(i);
    end

    if isfield(modelSBML.compartment,'metaid')
        if regexp(model.comps{i},'^[^a-zA-Z_]')
            EM='The compartment IDs are in numeric format. For the compliance with SBML specifications, compartment IDs will be preceded with "c_" string';
            dispEM(EM,false);
            model.comps(i)=strcat('c_',model.comps(i));
        end
        modelSBML.compartment(i).metaid=model.comps{i};
    end
    %Prepare Miriam strings
    if ~isempty(model.compMiriams{i})
        [~,sbo_ind] = ismember('sbo',model.compMiriams{i}.name);
        if sbo_ind > 0
            modelSBML.compartment(i).sboTerm=str2double(regexprep(model.compMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
            % remove the SBO term from compMiriams so the information is
            % not duplicated in the "annotation" field later on
            model.compMiriams{i}.name(sbo_ind) = [];
            model.compMiriams{i}.value(sbo_ind) = [];
        end
    end
    if ~isempty(model.compMiriams{i}) && isfield(modelSBML.compartment(i),'annotation')
        modelSBML.compartment(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.comps{i} '">'];
        modelSBML.compartment(i).annotation=[modelSBML.compartment(i).annotation '<bqbiol:is><rdf:Bag>'];
        modelSBML.compartment(i).annotation=[modelSBML.compartment(i).annotation getMiriam(model.compMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
    end
    if isfield(modelSBML.compartment, 'name')
        modelSBML.compartment(i).name=model.compNames{i};
    end
    if isfield(modelSBML.compartment, 'id')
        modelSBML.compartment(i).id=model.comps{i};
    end

end

%Begin writing species
for i=1:numel(model.mets)
    %Add the default values, as these will be the same in all entries
    if i==1
        if isfield(modelSBML.species, 'sboTerm')
            modelSBML.species(i).sboTerm=247;
        end
        if isfield(modelSBML.species, 'initialAmount')
            modelSBML.species(i).initialAmount=1;
        end
        if isfield(modelSBML.species, 'initialConcentration')
            modelSBML.species(i).initialConcentration=0;
        end
        if isfield(modelSBML.species, 'isSetInitialAmount')
            modelSBML.species(i).isSetInitialAmount=1;
        end
        if isfield(modelSBML.species, 'isSetInitialConcentration')
            modelSBML.species(i).isSetInitialConcentration=1;
        end
    end
    %Copy the default values to the next entry as long as it is not the
    %last one
    if i<numel(model.mets)
        modelSBML.species(i+1)=modelSBML.species(i);
    end

    if isfield(modelSBML.species,'metaid')
        modelSBML.species(i).metaid=['M_' model.mets{i}];
    end
    if isfield(modelSBML.species, 'name')
        modelSBML.species(i).name=model.metNames{i};
    end
    if isfield(modelSBML.species, 'id')
        modelSBML.species(i).id=['M_' model.mets{i}];
    end
    if isfield(modelSBML.species, 'compartment')
        modelSBML.species(i).compartment=model.comps{model.metComps(i)};
    end
    if isfield(model,'unconstrained')
        if model.unconstrained(i)
            modelSBML.species(i).boundaryCondition=1;
        end
    end
    if isfield(modelSBML.species, 'fbc_charge') && isfield(model,'metCharges')
        if ~isnan(model.metCharges(i))
            modelSBML.species(i).fbc_charge=model.metCharges(i);
            modelSBML.species(i).isSetfbc_charge=1;
        else
            modelSBML.species(i).isSetfbc_charge=0;
        end
    end
    if ~isempty(model.metMiriams{i})
        [~,sbo_ind] = ismember('sbo',model.metMiriams{i}.name);
        if sbo_ind > 0
            modelSBML.species(i).sboTerm=str2double(regexprep(model.metMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
            % remove the SBO term from metMiriams so the information is
            % not duplicated in the "annotation" field later on
            model.metMiriams{i}.name(sbo_ind) = [];
            model.metMiriams{i}.value(sbo_ind) = [];
        end
    end
    if isfield(modelSBML.species,'annotation')
        if ~isempty(model.metMiriams{i}) || ~isempty(model.metFormulas{i})
            hasInchi=false;
            if ~isempty(model.metFormulas{i})
                %Only export formula if there is no InChI. This is because
                %the metFormulas field is populated by InChIs if available
                if ~isempty(model.inchis{i})
                    hasInchi=true;
                end
                if hasInchi==false
                    modelSBML.species(i).fbc_chemicalFormula=model.metFormulas{i};
                end
            end
            if ~isempty(model.metMiriams{i}) || hasInchi==true
                modelSBML.species(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_M_' model.mets{i} '">'];
                modelSBML.species(i).annotation=[modelSBML.species(i).annotation '<bqbiol:is><rdf:Bag>'];
                if ~isempty(model.metMiriams{i})
                    modelSBML.species(i).annotation=[modelSBML.species(i).annotation getMiriam(model.metMiriams{i})];
                end
                if hasInchi==true
                    modelSBML.species(i).annotation=[modelSBML.species(i).annotation '<rdf:li rdf:resource="https://identifiers.org/inchi/InChI=' regexprep(model.inchis{i},'^InChI=','') '"/>'];
                    modelSBML.species(i).fbc_chemicalFormula=char(regexp(model.inchis{i}, '/(\w+)/', 'tokens', 'once'));
                end
                modelSBML.species(i).annotation=[modelSBML.species(i).annotation '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
            end
        end
    end
end

if isfield(model,'genes')
    for i=1:numel(model.genes)
        %Add the default values, as these will be the same in all entries
        if i==1
            if isfield(modelSBML.fbc_geneProduct, 'sboTerm')
                modelSBML.fbc_geneProduct(i).sboTerm=243;
            end
        end
        %Copy the default values to the next index as long as it is not the
        %last one
        if i<numel(model.genes)
            modelSBML.fbc_geneProduct(i+1)=modelSBML.fbc_geneProduct(i);
        end

        if isfield(modelSBML.fbc_geneProduct,'metaid')
            modelSBML.fbc_geneProduct(i).metaid=model.genes{i};
        end
        if ~isempty(model.geneMiriams{i})
            [~,sbo_ind] = ismember('sbo',model.geneMiriams{i}.name);
            if sbo_ind > 0
                modelSBML.fbc_geneProduct(i).sboTerm=str2double(regexprep(model.geneMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
                % remove the SBO term from compMiriams so the information is
                % not duplicated in the "annotation" field later on
                model.geneMiriams{i}.name(sbo_ind) = [];
                model.geneMiriams{i}.value(sbo_ind) = [];
            end
        end
        if ~isempty(model.geneMiriams{i}) && isfield(modelSBML.fbc_geneProduct(i),'annotation')
            modelSBML.fbc_geneProduct(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.genes{i} '">'];
            modelSBML.fbc_geneProduct(i).annotation=[modelSBML.fbc_geneProduct(i).annotation '<bqbiol:is><rdf:Bag>'];
            modelSBML.fbc_geneProduct(i).annotation=[modelSBML.fbc_geneProduct(i).annotation getMiriam(model.geneMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
        end
        if isfield(modelSBML.fbc_geneProduct, 'fbc_id')
            modelSBML.fbc_geneProduct(i).fbc_id=model.genes{i};
        end
        if isfield(modelSBML.fbc_geneProduct, 'fbc_label') && isfield(model,'geneShortNames')
            if isempty(model.geneShortNames{i})
                modelSBML.fbc_geneProduct(i).fbc_label=model.genes{i};
            else
                modelSBML.fbc_geneProduct(i).fbc_label=model.geneShortNames{i};
            end
        end
    end
    if exportGeneComplexes==true
        %Also add the complexes as genes. This is done by splitting grRules
        %on "or" and adding the ones which contain several genes
        geneComplexes={};
        if isfield(model,'grRules')
            %Only grRules which contain " and " can be complexes
            uniqueRules=unique(model.grRules);
            I=cellfun(@any,strfind(uniqueRules,' and '));
            uniqueRules(~I)=[];
            uniqueRules=strrep(uniqueRules,'(','');
            uniqueRules=strrep(uniqueRules,')','');
            uniqueRules=strrep(uniqueRules,' and ',':');
            for i=1:numel(uniqueRules)
                genes=regexp(uniqueRules(i),' or ','split');
                genes=genes{1}(:);
                %Check which ones are complexes
                I=cellfun(@any,strfind(genes,':'));
                geneComplexes=[geneComplexes;genes(I)];
            end
        end
        geneComplexes=unique(geneComplexes);
        if ~isempty(geneComplexes)
            %Then add them as genes. There is a possiblity that a complex
            %A&B is added as separate from B&A. This is not really an issue
            %so this is not dealt with
            for i=1:numel(geneComplexes)
                modelSBML.fbc_geneProduct(numel(model.genes)+i)=modelSBML.fbc_geneProduct(1);
                if isfield(modelSBML.fbc_geneProduct,'metaid')
                    modelSBML.fbc_geneProduct(numel(model.genes)+i).metaid=geneComplexes{i};
                end
                if isfield(modelSBML.fbc_geneProduct,'fbc_id')
                    modelSBML.fbc_geneProduct(numel(model.genes)+i).fbc_id=geneComplexes{i};
                else
                    modelSBML.fbc_geneProduct(i).fbc_label=modelSBML.fbc_geneProduct(i).fbc_id;
                end
            end
        end
    end
end

%Generate a list of unique fbc_bound names
totalValues=[model.lb; model.ub];
totalNames=cell(size(totalValues,1),1);

listUniqueValues=unique(totalValues);

for i=1:length(listUniqueValues)
    listUniqueNames{i,1}=['FB',num2str(i),'N',num2str(abs(round(listUniqueValues(i))))]; % create unique flux bound IDs.
    ind=find(ismember(totalValues,listUniqueValues(i)));
    totalNames(ind)=listUniqueNames(i,1);
end

for i=1:length(listUniqueNames)
    %Add the default values, as these will be the same in all entries
    if i==1
        if isfield(modelSBML.parameter, 'constant')
            modelSBML.parameter(i).constant=1;
        end
        if isfield(modelSBML.parameter, 'isSetValue')
            modelSBML.parameter(i).isSetValue=1;
        end
    end
    %Copy the default values to the next index as long as it is not the
    %last one
    if i<numel(listUniqueNames)
        modelSBML.parameter(i+1)=modelSBML.parameter(i);
    end
    modelSBML.parameter(i).id=listUniqueNames{i};
    modelSBML.parameter(i).value=listUniqueValues(i);
end

for i=1:numel(model.rxns)
    %Add the default values, as these will be the same in all entries
    if i==1
        if isfield(modelSBML.reaction, 'sboTerm')
            modelSBML.reaction(i).sboTerm=176;
        end
        if isfield(modelSBML.reaction, 'isSetFast')
            modelSBML.reaction(i).isSetFast=1;
        end
    end
    %Copy the default values to the next index as long as it is not the
    %last one
    if i<numel(model.rxns)
        modelSBML.reaction(i+1)=modelSBML.reaction(i);
    end

    if isfield(modelSBML.reaction,'metaid')
        modelSBML.reaction(i).metaid=['R_' model.rxns{i}];
    end

    %Export notes information
    if (~isnan(model.rxnConfidenceScores(i)) || ~isempty(model.rxnReferences{i}) || ~isempty(model.rxnNotes{i}))
        modelSBML.reaction(i).notes='<notes>';
        if ~isnan(model.rxnConfidenceScores(i))
            modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '
Confidence Level: ' num2str(model.rxnConfidenceScores(i)) '
'];
        end
        if ~isempty(model.rxnReferences{i})
            modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '
AUTHORS: ' model.rxnReferences{i} '
'];
        end
        if ~isempty(model.rxnNotes{i})
            modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '
NOTES: ' model.rxnNotes{i} '
'];
        end
        modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '</notes>'];
    end

    % Export SBO terms from rxnMiriams
    if ~isempty(model.rxnMiriams{i})
        [~,sbo_ind] = ismember('sbo',model.rxnMiriams{i}.name);
        if sbo_ind > 0
            modelSBML.reaction(i).sboTerm=str2double(regexprep(model.rxnMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
            % remove the SBO term from rxnMiriams so the information is not
            % duplicated in the "annotation" field later on
            model.rxnMiriams{i}.name(sbo_ind) = [];
            model.rxnMiriams{i}.value(sbo_ind) = [];
        end
    end

    %Export annotation information from rxnMiriams
    if (~isempty(model.rxnMiriams{i}) && isfield(modelSBML.reaction(i),'annotation')) || ~isempty(model.eccodes{i})
        modelSBML.reaction(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_R_' model.rxns{i} '">'];
        modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation '<bqbiol:is><rdf:Bag>'];
        if ~isempty(model.eccodes{i})
            eccodes=regexp(model.eccodes{i},';','split');
            for j=1:numel(eccodes)
                modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation  '<rdf:li rdf:resource="https://identifiers.org/ec-code/' regexprep(eccodes{j},'ec-code/|EC','') '"/>'];
            end
        end
        modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation getMiriam(model.rxnMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
    end

    if isfield(modelSBML.reaction, 'name')
        modelSBML.reaction(i).name=model.rxnNames{i};
    end
    if isfield(modelSBML.reaction, 'id')
        modelSBML.reaction(i).id=['R_' model.rxns{i}];
    end

    %Add the information about reactants and products
    involvedMets=addReactantsProducts(model,modelSBML,i);
    for j=1:numel(involvedMets.reactant)
        if j<numel(involvedMets.reactant)
            modelSBML.reaction(i).reactant(j+1)=modelSBML.reaction(i).reactant(j);
        end
        modelSBML.reaction(i).reactant(j).species=involvedMets.reactant(j).species;
        modelSBML.reaction(i).reactant(j).stoichiometry=involvedMets.reactant(j).stoichiometry;
        modelSBML.reaction(i).reactant(j).isSetStoichiometry=involvedMets.reactant(j).isSetStoichiometry;
        modelSBML.reaction(i).reactant(j).constant=involvedMets.reactant(j).constant;
    end
    if numel(involvedMets.reactant)==0
        modelSBML.reaction(i).reactant='';
    end
    for j=1:numel(involvedMets.product)
        if j<numel(involvedMets.product)
            modelSBML.reaction(i).product(j+1)=modelSBML.reaction(i).product(j);
        end
        modelSBML.reaction(i).product(j).species=involvedMets.product(j).species;
        modelSBML.reaction(i).product(j).stoichiometry=involvedMets.product(j).stoichiometry;
        modelSBML.reaction(i).product(j).isSetStoichiometry=involvedMets.product(j).isSetStoichiometry;
        modelSBML.reaction(i).product(j).constant=involvedMets.product(j).constant;
    end
    if numel(involvedMets.product)==0
        modelSBML.reaction(i).product='';
    end
    %Export reversibility information. Reactions are irreversible by
    %default
    if model.rev(i)==1
        modelSBML.reaction(i).reversible=1;
    end
    if isfield(model, 'rxnComps')
        modelSBML.reaction(i).compartment=model.comps{model.rxnComps(i)};
    end
    if isfield(model, 'grRules')
        modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association=model.grRules{i};
    end
    modelSBML.reaction(i).fbc_lowerFluxBound=totalNames{i};
    modelSBML.reaction(i).fbc_upperFluxBound=totalNames{length(model.lb)+i};
end

%Prepare subSystems Code taken from COBRA functions getModelSubSystems,
%writeSBML, findRxnsFromSubSystem under GNU General Public License v3.0,
%license file in readme/GPL.MD. Code modified for RAVEN
if modelHasSubsystems
    modelSBML.groups_group.groups_kind = 'partonomy';
    modelSBML.groups_group.sboTerm = 633;
    tmpStruct=modelSBML.groups_group;

    rxns=strcat('R_',model.rxns);
    if ~any(cellfun(@iscell,model.subSystems))
        if ~any(~cellfun(@isempty,model.subSystems))
            subSystems = {};
        else
            subSystems = setdiff(model.subSystems,'');
        end
    else
        orderedSubs = cellfun(@(x) columnVector(x),model.subSystems,'UniformOUtput',false);
        subSystems = setdiff(vertcat(orderedSubs{:}),'');
    end
    if isempty(subSystems)
        subSystems = {};
    end
    if ~isempty(subSystems)
        %Build the groups for the group package
        groupIDs = strcat('group',cellfun(@num2str, num2cell(1:length(subSystems)),'UniformOutput',false));
        for i = 1:length(subSystems)
            cgroup = tmpStruct;
            if ~any(cellfun(@iscell,model.subSystems))
                present = ismember(model.subSystems,subSystems{i});
            else
                present = cellfun(@(x) any(ismember(x,subSystems{i})),model.subSystems);
            end
            groupMembers = rxns(present);
            for j = 1:numel(groupMembers)
                cMember = tmpStruct.groups_member;
                cMember.groups_idRef = groupMembers{j};
                if j == 1
                    cgroup.groups_member = cMember;
                else
                    cgroup.groups_member(j) = cMember;
                end
            end
            cgroup.groups_id = groupIDs{i};
            cgroup.groups_name = subSystems{i};
            if i == 1
                modelSBML.groups_group = cgroup;
            else
                modelSBML.groups_group(i) = cgroup;
            end
        end
    end
end

%Prepare fbc_objective subfield

modelSBML.fbc_objective.fbc_type='maximize';
modelSBML.fbc_objective.fbc_id='obj';

ind=find(model.c);

if isempty(ind)
    modelSBML.fbc_objective.fbc_fluxObjective.fbc_coefficient=0;
else
    for i=1:length(ind)
        %Copy the default values to the next index as long as it is not the
        %last one
        if i<numel(ind)
            modelSBML.reaction(i+1)=modelSBML.reaction(i);
        end
        values=model.c(model.c~=0);
        modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction=modelSBML.reaction(ind(i)).id;
        modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient=values(i);
        modelSBML.fbc_objective(i).fbc_fluxObjective.isSetfbc_coefficient=1;
    end
end

modelSBML.fbc_activeObjective=modelSBML.fbc_objective.fbc_id;

fbcStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/fbc/version',num2str(sbmlPackageVersions(1))];
if modelHasSubsystems
    groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))];
    modelSBML.namespaces=struct('prefix',{'','fbc','groups'},...
    'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
    fbcStr,groupStr});
else
    modelSBML.namespaces=struct('prefix',{'','fbc'},...
    'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
    fbcStr});
end

if sbmlPackageVersions(1) == 2
    modelSBML.fbc_strict=1;
end

modelSBML.rule=[];
modelSBML.constraint=[];

[ravenDir,prevDir]=findRAVENroot();
fileName=checkFileExistence(fileName,1,true,false);
cd(fullfile(ravenDir,'software','libSBML'));
OutputSBML(modelSBML,fileName,1,0,[1,0]);
cd(prevDir);
end

function modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions)
%Returns the blank SBML model structure by using appropriate libSBML
%functions. This creates structure by considering three levels

sbmlFieldNames=getStructureFieldnames('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
sbmlDefaultValues=getDefaultValues('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);

for i=1:numel(sbmlFieldNames)
    modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
    sbmlSubfieldNames=getStructureFieldnames(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
    sbmlSubfieldValues=getDefaultValues(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
    if ~strcmp(sbmlFieldNames{1,i},'event') && ~strcmp(sbmlFieldNames{1,i},'functionDefinition') && ~strcmp(sbmlFieldNames{1,i},'initialAssignment')
        for j=1:numel(sbmlSubfieldNames)
            modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j})=sbmlSubfieldValues{1,j};
            sbmlSubsubfieldNames=getStructureFieldnames(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
            sbmlSubsubfieldValues=getDefaultValues(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
            if ~strcmp(sbmlSubfieldNames{1,j},'modifier') && ~strcmp(sbmlSubfieldNames{1,j},'kineticLaw')
                for k=1:numel(sbmlSubsubfieldNames)
                    %'compartment' and 'species' fields are not supposed to
                    %have their standalone structures if they are subfields
                    %or subsubfields
                    if ~strcmp(sbmlSubfieldNames{1,j},'compartment') && ~strcmp(sbmlSubfieldNames{1,j},'species')
                        modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k})=sbmlSubsubfieldValues{1,k};
                    end
                    %If it is fbc_association in the third level, we need
                    %to establish the fourth level, since libSBML requires
                    %it
                    if strcmp(sbmlSubsubfieldNames{1,k},'fbc_association')
                        fbc_associationFieldNames=getStructureFieldnames('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
                        fbc_associationFieldValues=getDefaultValues('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
                        for l=1:numel(fbc_associationFieldNames)
                            modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k}).(fbc_associationFieldNames{1,l})=fbc_associationFieldValues{1,l};
                        end
                    end
                end
            end
        end
    end
    if ~isstruct(modelSBML.(sbmlFieldNames{1,i}))
        modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
    end
end

modelSBML.unitDefinition.id='mmol_per_gDW_per_hr';

unitFieldNames=getStructureFieldnames('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
unitDefaultValues=getDefaultValues('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);

kinds={'mole','gram','second'};
exponents=[1 -1 -1];
scales=[-3 0 0];
multipliers=[1 1 1*60*60];

for i=1:numel(unitFieldNames)
    modelSBML.unitDefinition.unit(1).(unitFieldNames{1,i})=unitDefaultValues{1,i};
    for j=1:3
        modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=unitDefaultValues{1,i};
        if strcmp(unitFieldNames{1,i},'kind')
            modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=kinds{j};
        elseif strcmp(unitFieldNames{1,i},'exponent')
            modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=exponents(j);
        elseif strcmp(unitFieldNames{1,i},'scale')
            modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=scales(j);
        elseif strcmp(unitFieldNames{1,i},'multiplier')
            modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=multipliers(j);
        end
    end
end
end

function miriamString=getMiriam(miriamStruct)
%Returns a string with list elements for a miriam structure ('<rdf:li
%rdf:resource="https://identifiers.org/go/GO:0005739"/>' for example). This
%is just to speed up things since this is done many times during the
%exporting

miriamString='';
if isfield(miriamStruct,'name')
    for i=1:numel(miriamStruct.name)
        miriamString=[miriamString '<rdf:li rdf:resource="https://identifiers.org/' miriamStruct.name{i} '/' miriamStruct.value{i} '"/>'];
    end
end
end

function [tmp_Rxn]=addReactantsProducts(model,sbmlModel,i)
%This function provides reactants and products for particular reaction. The
%function was 'borrowed' from writeSBML in COBRA toolbox, lines 663-679

met_idx = find(model.S(:, i));
tmp_Rxn.product=[];
tmp_Rxn.reactant=[];
for j_met=1:size(met_idx,1)
    tmp_idx = met_idx(j_met,1);
    sbml_tmp_species_ref.species = sbmlModel.species(tmp_idx).id;
    met_stoich = model.S(tmp_idx, i);
    sbml_tmp_species_ref.stoichiometry = abs(met_stoich);
    sbml_tmp_species_ref.isSetStoichiometry=1;
    sbml_tmp_species_ref.constant=1;
    if (met_stoich > 0)
        tmp_Rxn.product = [ tmp_Rxn.product, sbml_tmp_species_ref ];
    else
        tmp_Rxn.reactant = [ tmp_Rxn.reactant, sbml_tmp_species_ref];
    end
end
end

function vecT = columnVector(vec)
% Code below taken from COBRA Toolbox under GNU General Public License v3.0
% license file in readme/GPL.MD.
%
% Converts a vector to a column vector
%
% USAGE:
%
%   vecT = columnVector(vec)
%
% INPUT:
%   vec:     a vector
%
% OUTPUT:
%   vecT:    a column vector

[n, m] = size(vec);

if n < m
    vecT = vec';
else
    vecT = vec;
end
end

author: SysBioChalmers

repo: RAVEN


text:
function options = get_nmf_default_options()
% Function to set default options
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai
%

    options.not_store_infos     = false;
    options.metric_type         = 'euc';
    options.metric_param        = [];
    options.max_epoch           = 100;
    options.f_opt               = -Inf;
    options.tol_optgap          = 1.0e-12;
    options.verbose             = 0;
    options.store_sol           = false;
    options.batch_size          = 1;
    options.permute_on          = false;
    options.gnd                 = [];
    options.x_init_robust       = false;
    options.calc_symmetry       = false;
    options.calc_clustering_acc = false;
    options.clustering_gnd      = [];
    options.clustering_classnum = 0;
    options.clustering_eval_num = 5;
    options.disp_freq           = 0;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [ spect, frqs ] = LTASS( speech_folder_OR_vec, nfft, fs )
% Computes the Long-Term Average Speech Spectrum from a folder of speech files or vector of speech samples
%
% Syntax:	[ spect, frqs ] = LTASS( speech_folder_OR_vec, nfft )
%
% Inputs:
% 	speech_folder_OR_vec - The path to the folder containing the speech
%                          files OR a vector of concatenated speech signals
% 	nfft - The number of FFT points used to compute the LTASS
% 	fs - The sampling frequency to use (if not provided then the sampling
%        frequency of the file is used)
%
% Outputs:
% 	spect - The LTASS spectrum
% 	frqs - The frequency vector for the spectrum

% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Copyright: Jacob Donley 2017
% Date: 17 June 2016
% Revision: 0.4 (30 March 2017)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isa(speech_folder_OR_vec,'char') % if a character array (string)
    % Just incase this function tries to call getAllFiles within a class folder we
    % should create a function handle for getAllFiles regardless
    inf = dbstack('-completenames');
    funcName = 'getAllFiles';
    funcPath = inf.file;
    classDirs = getClassDirs(funcPath);
    getAllFiles_ = str2func([classDirs funcName]);

    %% Start LTASS
    files = getAllFiles_(speech_folder_OR_vec);
    speech=[];
    F = length(files);
    for file = 1:F
        try
            [audioSig,fs_] = audioread(files{file});
            if nargin < 3, fs = fs_; end
            audioSig = audioSig ./ rms(audioSig(:));
        catch err
            if strcmp(err.identifier, 'MATLAB:audiovideo:audioread:FileTypeNotSupported')
                continue; % Skip unsupported files
            end
        end
        speech = [speech; audioSig];
    end
    if nargin < 2
        nfft = numel(speech);
    end
    if logical(mod(nfft,2)) % if isodd( nfft )
        nfft = nfft-1; % Force nfft to be even so that pwelch returns normalised frequencies [0,...,1]
    end

else
    speech = speech_folder_OR_vec;
end

%%
win_=rectwin(nfft);
ovlap = 0;

[spect,frqs]=pwelch(speech,win_,nfft*ovlap,nfft,fs,'power'); % Power spectrum
spect = sqrt(spect); % Magnitude spectrum

end

function classDirs = getClassDirs(FullPath)
classDirs = '';
classes = strfind(FullPath,'+');
for c = 1:length(classes)
    clas = FullPath(classes(c):end);
    stp = strfind(clas,filesep);
    classDirs = [classDirs  clas(2:stp(1)-1) '.'];
end
end

author: jtkim-kaist

repo: Speech-enhancement


text:
function e = horzcat(varargin)
% overloads [e1,e2,e3..]

e = cat(2,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
load sstest2_results.mat % TM T1 TU f
whos

index = UFget ;

tmax = max (max (TM))
tmin = min (TM (find (TM > 0)))

nmat = length (f) ;
k = nmat ;

    for kind = 1:4

	subplot (2,4,kind) ;
	r = TM (1:k,kind) ./ T1 (1:k,kind) ;
	rmin = min (r) ;
	rmax = max (r) ;
	loglog (TM (1:k,kind), r, 'o', ...
	    [tmin tmax], [1 1], 'r-', ...
	    [tmin tmax], [1.1 1.1], 'r-', ...
	    [tmin tmax], [1/1.1 1/1.1], 'r-', ...
	    [tmin tmax], [2 2], 'g-', ...
	    [tmin tmax], [1.5 1.5], 'g-', ...
	    [tmin tmax], [1/1.5 1/1.5], 'g-', ...
	    [tmin tmax], [.5 .5], 'g-' );
	if (k > 2)
	    axis ([tmin tmax rmin rmax]) ;
	    set (gca, 'XTick', [1e-5 1e-4  1e-3 1e-2 1e-1 1 10]) ;
	    set (gca, 'YTick', [.5 1/1.5 1/1.1 1 1.1 1.5 2]) ;
	end
	xlabel ('MATLAB time') ;
	ylabel ('MATLAB/SM time') ;
	if (kind == 1)
	    title ('real*real') ;
	elseif (kind == 2)
	    title ('complex*real') ;
	elseif (kind == 3)
	    title ('real*complex') ;
	elseif (kind == 4)
	    title ('complex*complex') ;
	end

    end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%  Copyright (C) 2011-2012 by Jun He, Laura Balzano, and Arthur Szlam
%
%  This file is part of the GRASTA library.
%  It is provided without any warranty of fitness
%  for any purpose. You can redistribute this file
%  and/or modify it under the terms of the GNU
%  Lesser General Public License (LGPL) as published
%  by the Free Software Foundation, either version 3
%  of the License or (at your option) any later version.
%  (see http://www.opensource.org/licenses for more info)

function [ Unew, STATUSnew, OPTSnew ] = grasta_stream( y_Omega, idx, U0, STATUS,OPTIONS, OPTS)
%  grasta_stream is the streaming version of GRASTA which can be used for
%  online data processing
%
%   [1] Jun He, Laura Balzano, and John C.S. Lui. Online Robust Subspace Tracking from Partial Information
%       http://arxiv.org/abs/1109.3827
%   [2] Jun He, Laura Balzano, and Arthur Szlam. Incremental gradient on the grassmannian for online foreground
%       and background separation in subsampled video. In IEEE Conference on Computer Vision and Pattern Recognition
%       (CVPR), June 2012.
%
%  Usage:
%       [ Unew, STATUSnew, OPTSnew ] = grasta_stream( y_Omega, idx, U0, STATUS, OPTIONS, OPTS)
%   Inputs:
%       y_Omega: current partial observation of a full data vector
%       idx:     the observed indices of y_Omega
%       U0:      the current estimated subspace
%       STATUS:  the current status of GRASTA
%           last_mu: \mu_{t-1} for adaptive step rule
%           step_scale: the estimated step_scale constant for adaptive
%       step-size rule
%           lelve: the current level of the "Multi-Level" adaptive rule
%           last_gamma and last_w: previous gradient = last_gamma*last_w'
%
%       OPTIONS: the options of GRASTA
%           CONSTANT_STEP: default 0, will use multi-level step-size rule; > 0 will
%           use constant step-size rule
%           DIM_M : the ambient dimension
%           RANK  : the estimated low-rank of the underlying system
%           rho : ADMM constant step
%           ITER_MAX: the max_iter for ADMM solver
%           TOL : the stopping tolerance of admm_srp / mex_srp
%
%       OPTS:    the current options for the subproblem of GRASTA
%           refer to mex_srp.cpp or admm_srp.m
%   Outputs:
%       Unew:    the updated subspace
%       STATUSnew: the updated running status
%       OPTSnew: the updated options for the subproblem of GRASTA
%
%  Author: Jun He, Laura Balzano
%  Email:  hejun.zz@gmail.com, sunbeam@ece.wisc.edu
%  Date:   Sept. 01, 2012
%
%

LEVEL_FACTOR        = 2;

if isfield(OPTIONS,'MIN_MU'),
    MIN_MU = OPTIONS.MIN_MU;
else
    MIN_MU = 1;
end

if isfield(OPTIONS,'MAX_MU'),
    MAX_MU = OPTIONS.MAX_MU;
else
    MAX_MU = 15;
end

if isfield(OPTIONS,'USE_MEX'),
    USE_MEX = OPTIONS.USE_MEX;
else
    USE_MEX = 1;
end

if isfield(OPTIONS,'CONSTANT_STEP'),
    CONSTANT_STEP = OPTIONS.CONSTANT_STEP;
else
    CONSTANT_STEP = 0;
end

if isfield(OPTIONS,'DIM_M'),
    DIM_M = OPTIONS.DIM_M;
else
    error('Should specify OPTIONS.DIM_M data ambient dimension!!!\n');
end

if isfield(OPTIONS,'ITER_MIN'),
    MIN_ITER = OPTIONS.ITER_MIN;
else
    MIN_ITER = 5;
end

if isfield(OPTIONS,'ITER_MAX'),
    ITER_MAX = OPTIONS.ITER_MAX;
else
    ITER_MAX = 60;
end

if isfield(OPTIONS,'TOL'),
    TOL = OPTIONS.TOL;
else
    TOL = 1e-6;
end

if isfield(OPTIONS,'MAX_LEVEL'),
    MAX_LEVEL = OPTIONS.MAX_LEVEL;
else
    MAX_LEVEL = 20;
end

if isfield(OPTIONS,'RANK'),
    RANK = OPTIONS.RANK;
else
    error('Should specify OPTIONS.RANK!!!\n');
end

if isfield(OPTIONS,'QUIET'),
    QUIET = OPTIONS.QUIET;
else
    QUIET = 0;
end

DEFAULT_MU_HIGH     = (MAX_MU-1)/2;
DEFAULT_MU_LOW      = MIN_MU + 2;

% If we first call GRASTA_stream then do the following initilization
if STATUS.init == 0,
    STATUS.init         = 1; % Do not enter this initial part any longer
    STATUS.curr_iter    = 0; % For debug

    STATUS.last_mu      = MIN_MU;
    STATUS.level        = 0;
    STATUS.step_scale   = 0;
    STATUS.last_w       = zeros(RANK,1);
    STATUS.last_gamma   = zeros(DIM_M,1);

    OPTS.TOL        = TOL;
    OPTS.MAX_ITER   = MIN_ITER;  % the max iteration of ADMM at level=0
    OPTS.QUIET      = 1;

    if isfield(OPTIONS,'rho'),
        OPTS.RHO = OPTIONS.rho;
    else
        OPTS.RHO = 1.8;
    end

    U0 = orth(randn(DIM_M,RANK));
end

%%%%%%%%%%%%%%%
% main framework of GRASTA
U_Omega = U0(idx,:);

if USE_MEX,
    [s_t, w, ldual] = mex_srp(U_Omega, y_Omega, OPTS);
else
    [s_t, w, ldual] = admm_srp(U_Omega, y_Omega, OPTS);
end

gamma_1 = ldual;% + OPTS.RHO*(U_Omega*w + s_t - y_Omega);
UtDual_omega = U_Omega' * gamma_1;
gamma_2 = U0 * UtDual_omega;
gamma = zeros(DIM_M,1);
gamma(idx) = gamma_1;
gamma = gamma - gamma_2;

gamma_norm = norm(gamma);
w_norm     = norm(w);
sG = gamma_norm * w_norm;

% Here we use the adaptive step-size rule for SGD. The adaptive can work
% well for both dynamic and static subspace tracking tasks

% 1. determine the step scale from the first observation
if ~STATUS.step_scale,
    STATUS.step_scale = 0.5*pi*(1+MIN_MU)/sG;

    if ~QUIET,
        fprintf('Level 0: %.2e\n',STATUS.step_scale);
    end
end

% 2. inner product of previous grad and current grad
grad_ip = trace(STATUS.last_w * (STATUS.last_gamma' * gamma) * w');

%%% avoid inner product too large
normalization = norm(STATUS.last_gamma * STATUS.last_w','fro') * norm(gamma * w','fro');
if normalization == 0,
    grad_ip_normalization = 0;
else
    grad_ip_normalization = grad_ip/normalization;
end
%%%

% 3. if the two consecutive grad in the same direction, we take a larger
% step along the gradient direction, otherwise take a small step along the
% gradient direction
STATUS.last_mu = max(STATUS.last_mu + sigmoid(-grad_ip_normalization) , MIN_MU);

if CONSTANT_STEP > 0,
    t = CONSTANT_STEP;
else
    % should not take a step larger than pi/2
    t = STATUS.step_scale * LEVEL_FACTOR^(-STATUS.level) * sG / (1+STATUS.last_mu);
    if t>=pi/3,
        t = pi/3;
    end

    % Adjust the level
    bShrUpd = 0;
    if STATUS.last_mu <= MIN_MU,
        if STATUS.level > 1,
            bShrUpd = 1;
            STATUS.level = STATUS.level - 1;

            if ~QUIET,
                fprintf('multi-level adaption - decreasing, t:%.2e, vectors: %d, level: %d\n',...
                    t,STATUS.curr_iter, STATUS.level);
            end
            STATUS.curr_iter  = 0;
        end

        STATUS.last_mu    = DEFAULT_MU_LOW;
    elseif STATUS.last_mu > MAX_MU,
        if STATUS.level < MAX_LEVEL,
            bShrUpd = 1;
            STATUS.level = STATUS.level + 1;

            if ~QUIET,
                fprintf('multi-level adaption - increasing, t:%.2e, vectors: %d, level: %d\n',...
                    t, STATUS.curr_iter,  STATUS.level);
            end
            STATUS.curr_iter  = 0;
            STATUS.last_mu    = DEFAULT_MU_HIGH;
        else
            STATUS.last_mu    = MAX_MU;
        end
    end

    if bShrUpd,
        if STATUS.level>=0 && STATUS.level <4,          % [0,4)
            OPTS.MAX_ITER = MIN_ITER;
        elseif STATUS.level>=4 && STATUS.level <7,      % [4,7)
            OPTS.MAX_ITER = min(MIN_ITER*2, ITER_MAX);
        elseif STATUS.level>=7 && STATUS.level <10,     % [7,10)
            OPTS.MAX_ITER = min(MIN_ITER*4, ITER_MAX);
        elseif STATUS.level>=10 && STATUS.level <14,    % [10,14)
            OPTS.MAX_ITER = min(MIN_ITER*8, ITER_MAX);
        else
            OPTS.MAX_ITER = ITER_MAX;               % [14,...)
        end

        if ~QUIET,
            fprintf('Will use %d ADMM iterations in level %d\n',OPTS.MAX_ITER, STATUS.level);
        end
    end

end

% 4. update the gradient for further step size update
STATUS.last_gamma  = gamma;
STATUS.last_w      = w;

STATUS.grad_ip = grad_ip_normalization; % just for debugging

% Take the gradient step along Grassmannian geodesic.
alpha = w/w_norm;
beta  = gamma/gamma_norm;
step = (cos(t)-1)*U0*(alpha*alpha')  - sin(t)*beta*alpha';

U0 = U0 + step;

%%

STATUS.s_t   = s_t;
STATUS.w     = w;
STATUS.ldual = ldual;
STATUS.SCALE = 1;
STATUS.curr_iter = STATUS.curr_iter + 1;

STATUS.grasta_t = t;

%
%%%%%%%%%%%%%%%%%%%%%%

Unew = U0;
STATUSnew = STATUS;
OPTSnew = OPTS;
end

%% Function of Sigmoid

function fval = sigmoid(x)
FMIN = -1; FMAX = 1;
omega = 0.1;

fval = FMIN + (FMAX - FMIN)/(1 - (FMAX/FMIN)*exp(-x/omega));
end

author: andrewssobral

repo: lrslibrary


text:
function Population = EnvironmentalSelection(Population,N,Z,Zmin)
% The environmental selection of NSGA-III

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

    if isempty(Zmin)
        Zmin = ones(1,size(Z,2));
    end

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);
    Next = FrontNo < MaxFNo;

    %% Select the solutions in the last front
    Last   = find(FrontNo==MaxFNo);
    Choose = LastSelection(Population(Next).objs,Population(Last).objs,N-sum(Next),Z,Zmin);
    Next(Last(Choose)) = true;
    % Population for next generation
    Population = Population(Next);
end

function Choose = LastSelection(PopObj1,PopObj2,K,Z,Zmin)
% Select part of the solutions in the last front

    PopObj = [PopObj1;PopObj2] - repmat(Zmin,size(PopObj1,1)+size(PopObj2,1),1);
    [N,M]  = size(PopObj);
    N1     = size(PopObj1,1);
    N2     = size(PopObj2,1);
    NZ     = size(Z,1);

    %% Normalization
    % Detect the extreme points
    Extreme = zeros(1,M);
    w       = zeros(M)+1e-6+eye(M);
    for i = 1 : M
        [~,Extreme(i)] = min(max(PopObj./repmat(w(i,:),N,1),[],2));
    end
    % Calculate the intercepts of the hyperplane constructed by the extreme
    % points and the axes
    Hyperplane = PopObj(Extreme,:)\ones(M,1);
    a = 1./Hyperplane;
    if any(isnan(a))
        a = max(PopObj,[],1)';
    end
    % Normalization
    PopObj = PopObj./repmat(a',N,1);

    %% Associate each solution with one reference point
    % Calculate the distance of each solution to each reference vector
    Cosine   = 1 - pdist2(PopObj,Z,'cosine');
    Distance = repmat(sqrt(sum(PopObj.^2,2)),1,NZ).*sqrt(1-Cosine.^2);
    % Associate each solution with its nearest reference point
    [d,pi] = min(Distance',[],1);

    %% Calculate the number of associated solutions except for the last front of each reference point
    rho = hist(pi(1:N1),1:NZ);

    %% Environmental selection
    Choose  = false(1,N2);
    Zchoose = true(1,NZ);
    % Select K solutions one by one
    while sum(Choose) < K
        % Select the least crowded reference point
        Temp = find(Zchoose);
        Jmin = find(rho(Temp)==min(rho(Temp)));
        j    = Temp(Jmin(randi(length(Jmin))));
        I    = find(Choose==0 & pi(N1+1:end)==j);
        % Then select one solution associated with this reference point
        if ~isempty(I)
            if rho(j) == 0
                [~,s] = min(d(N1+I));
            else
                s = randi(length(I));
            end
            Choose(I(s)) = true;
            rho(j) = rho(j) + 1;
        else
            Zchoose(j) = false;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Example
% An example script, it is used to show how to use the arcs
% functions of the KUKA iiwa matlab toolbox
% First start the server on the KUKA iiwa controller
% Then run the following script in Matlab

% Note you have 60 seconds to connect to server after starting the
% application (MatlabToolboxServer) from the smart pad on the robot controller.

% Mohammad SAFEEA, 3rd of May 2017

close all,clear all;clc;
warning('off');
ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
else

      %% Go to initial configuration
      relVel=0.25; % over ride relative joint velocities

      pos={0, -pi / 180 * 10, 0, -pi / 180 * 100, pi / 180 * 90,pi / 180 * 90, 0};   % initial cofiguration

      movePTPJointSpace( t_Kuka , pos, relVel); % go to home position
      %% Move in an arc, the orientation of EEF changes while performing the motion,
      % The function utilized (movePTPCirc1OrintationInter)
      % f2 is the final frame, at which the arc motion ends
      % f1 is an intermidiate frame, through wich the robot passes while
      % performing the motion.

      f1=getEEFPos( t_Kuka );
      f2=f1;
      r=75;
      f1{2}=f1{2}+r;
      f1{3}=f1{3}-r;
      f1{6}=f1{6}+pi/8;

      f2{3}=f2{3}-2*r;
      f2{6}=f2{6}+pi/2;

      vel=150; % linear velocity of end-effector mm/sec
      movePTPCirc1OrintationInter( t_Kuka , f1,f2, vel)

          %% Move robot in joint space to some initial configuration
pinit={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % joint angles of initial confuguration
relVel=0.15; % the relative velocity
movePTPJointSpace( t_Kuka , pinit, relVel); % point to point motion in joint space

%% Move EEF -100 mm in Z direction
deltaX=0.0;deltaY=0;deltaZ=-100.; % relative displacemnets of end-effector
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
movePTPLineEefRelBase( t_Kuka , Pos, vel);

%% Store the current position in the memory
Cen=getEEFPos(t_Kuka); % Concider the current position as the center of the arcs

%% Move EEF 50mm in X direction
deltaX=100;deltaY=0;deltaZ=0.;
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
movePTPLineEefRelBase( t_Kuka , Pos, vel);
%% Store the current position in the memory
circle_Starting_Point=getEEFPos( t_Kuka ); % Consider the current point as circle starting point

%% Move in an arc, the arc is drawn on an incliend plane
% using the function ((movePTPArc_AC))
theta=pi/2; % the angle subtended by the arc at the center ((c))
k=[1;1;1]; % normal vector of the plane, on which the circle is drawn
c=[Cen{1};Cen{2};Cen{3}]; % the center of the arc
vel=100; % the motion velocity mm/sec
movePTPArc_AC(t_Kuka,theta,c,k,vel)

%% Go back to ((circle_Starting_Point)) coordinates
vel=150;
movePTPLineEEF( t_Kuka , circle_Starting_Point, vel);
%% Move in an arc, the arc is drawn in XY plane
% using the function ((movePTPArcXY_AC))
theta=1.98*pi; % the angle subtended by the arc at the center ((c))
c=[Cen{1};Cen{2}]; % the XY coordinate of the center of the arc
vel=150; % the motion velocity mm/sec
movePTPArcXY_AC(t_Kuka,theta,c,vel)

      %% turn off the server
       net_turnOffServer( t_Kuka );

       fclose(t_Kuka);
end

warning('on');

repo: KST-Kuka-Sunrise-Toolbox


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = (X' * X) / m;  %'
[U, S, V] = svd(Sigma);

% =========================================================================

end

author: gopaczewski

repo: coursera-ml


text:
%% Optimization of a simple (Rosenbrock) function, with no constraints
% The unconstrained solution is at [1,1]
rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;

% With no constraints, operation simply passes through
% directly to fminsearch. The solution should be [1 1]
xsol = fminsearchbnd(rosen,[3 3])

%% Full lower and upper bound constraints which will all be inactive
xsol = fminsearchbnd(rosen,[3 3],[-1 -1],[4 4])

%% Only lower bound constraints
xsol = fminsearchbnd(rosen,[3 3],[2 2])

%% Only upper bound constraints
xsol = fminsearchbnd(rosen,[-5 -5],[],[0 0])

%% Dual constraints
xsol = fminsearchbnd(rosen,[2.5 2.5],[2 2],[3 3])

%% Dual constraints, with an infeasible starting guess
xsol = fminsearchbnd(rosen,[0 0],[2 2],[3 3])

%% Mixed constraints
xsol = fminsearchbnd(rosen,[0 0],[2 -inf],[inf 3])

%% Provide your own fminsearch options
opts = optimset('fminsearch');
opts.Display = 'iter';
opts.TolX = 1.e-12;

n = [10,5];
H = randn(n);
H=H'*H;
Quadraticfun = @(x) x*H*x';

% Global minimizer is at [0 0 0 0 0].
% Set all lower bound constraints, all of which will
% be active in this test.
LB = [.5 .5 .5 .5 .5];
xsol = fminsearchbnd(Quadraticfun,[1 2 3 4 5],LB,[],opts)

%% Exactly fix one variable, constrain some others, and set a tolerance
opts = optimset('fminsearch');
opts.TolFun = 1.e-12;

LB = [-inf 2 1 -10];
UB = [ inf  inf 1  inf];
xsol = fminsearchbnd(@(x) norm(x),[1 3 1 1],LB,UB,opts)

%% All the standard outputs from fminsearch are still returned
[xsol,fval,exitflag,output] = fminsearchbnd(@(x) norm(x),[1 3 1 1],LB,UB)

author: Sable

repo: mcbench-benchmarks


text:
function [version_num, version_name, version_date] = check_matlab_version()
% [version_num, version_name, version_date] = check_matlab_version()
%
% version_date will be in datenum format
% Compare to a date like this: version_date > datenum('01-Jan-2016')
%
% [version_num, version_name, version_date] = check_matlab_version;
% if version_date > datenum('01-Jan-2016'), disp('You are running matlab 2016 or later!'), else, disp('old version!'), end

which_ver = ver('matlab');

version_num = str2double(which_ver.Version);

version_name = which_ver.Release;

version_date = datenum(which_ver.Date);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function out = CO_TranslateShape(y,shape,d,howToMove)
% CO_TranslateShape  Statistics on datapoints inside geometric shapes across the time series.
%
% Inputs specify a shape and its size, and a method for moving this shape
% through the time domain.
%
% This is usually more informative in an embedding space (CO_Embed2_...), but
% here we do it just in the temporal domain (_t_).
%
% In the future, could perform a similar analysis with a soft boundary, some
% decaying force function V(r), or perhaps truncated...?
%
%---INPUTS:
% y, the input time series
% shape, the shape to move about the time-domain (e.g., 'circle')
% d, a parameter specifying the size of the shape (e.g., d = 2)
% howToMove, a method specifying how to move the shape about, e.g., 'pts'
%               places the shape on each point in the time series.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check inputs:
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(shape)
    shape = 'circle';
end
if nargin < 3 || isempty(d)
    d = 2; % a default distance d = 2
end
if nargin < 4 || isempty(howToMove)
    howToMove = 'pts'; % by default, places shapes on each timepoint
end

% ------------------------------------------------------------------------------
% Preliminaries:
% ------------------------------------------------------------------------------
N = length(y); % the length of the time series

% y must be a column vector, transpose it if it's a row vector
if size(y,2) > size(y,1)
    y = y';
end

% Add a time index
ty = [(1:N)', y]; % has increasing integers as time in the first column

%-------------------------------------------------------------------------------
% Generate the statistics on the number of points inside the shape as it is
% translated across the time series
%-------------------------------------------------------------------------------
switch howToMove
    case 'pts' % Place shapes on each timepoint (excluding a range at start and end)
        switch shape
            case 'circle' % uses a circle of radius 'd'
                r = d;
                w = floor(r); % only consider a window radius w (these are the
                            %    only points that could possibly be inside)
                rnge = 1+w:N-w;
                NN = length(rnge); % number of admissible points
                np = zeros(NN,1); % number of points
                for i = 1:NN
                    win = ty(rnge(i)-w:rnge(i)+w,:); % create window
                    difwin = win - ones(2*w+1,1)*ty(rnge(i),:);
                    np(i) = sum(sum(difwin.^2,2) <= r^2); % number of points enclosed in shape
                end
            case 'rectangle'
                % uses a rectangle of half-width d (integer), height as double the current time
                % series value (on either side of origin)
                w = d;
                rnge = 1+d:N-d;
                NN = length(rnge); % number of admissible points
                np = zeros(NN,1); % number of points
                for i = 1:NN
                    np(i) = sum(abs(y(rnge(i)-d:rnge(i)+d)) <= abs(y(i)));
                end
        otherwise
            error('Unknown shape ''%s''',shape)
        end
    otherwise
        error('Unknown setting for ''howToMove'' input: ''%s''',howToMove)
end

%-------------------------------------------------------------------------------
%% Output statistics
%-------------------------------------------------------------------------------

% Statistics on number of hits inside the shape:
out.max = max(np); % max hits
out.std = std(np); % std of hits
out.mean = mean(np); % mean number of hits

% Count the hits:
histnp = arrayfun(@(x)sum(np==x),unique(np));

% Compute mode of the histogram:
[out.npatmode,out.mode] = max(histnp);
out.npatmode = out.npatmode/NN;

% Output all stats:
if 2*w + 1 >= 1; out.ones = mean(np==1); end
if 2*w + 1 >= 2; out.twos = mean(np==2); end
if 2*w + 1 >= 3; out.threes = mean(np==3); end
if 2*w + 1 >= 4; out.fours = mean(np==4); end
if 2*w + 1 >= 5; out.fives = mean(np==5); end
if 2*w + 1 >= 6; out.sixes = mean(np==6); end
if 2*w + 1 >= 7; out.sevens = mean(np==7); end
if 2*w + 1 >= 8; out.eights = mean(np==8); end
if 2*w + 1 >= 9; out.nines = mean(np==9); end
if 2*w + 1 >= 10; out.tens = mean(np==10); end
if 2*w + 1 >= 11; out.elevens = mean(np==11); end

% -----
% Stationarity of the statistics in 2,3, & 4 segments of the time series:

out.statav2_m = SY_SlidingWindow(np,'mean','std',2,1);
out.statav2_s = SY_SlidingWindow(np,'std','std',2,1);
out.statav3_m = SY_SlidingWindow(np,'mean','std',3,1);
out.statav3_s = SY_SlidingWindow(np,'std','std',3,1);
out.statav4_m = SY_SlidingWindow(np,'mean','std',4,1);
out.statav4_s = SY_SlidingWindow(np,'std','std',4,1);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function obj = addVariable(obj, label, nodes, varargin)
    % Adds NLP decision variables
    %
    % Parameters:
    % label: the label name (column) of the variable @type char
    % nodes: the node list of the variable @type rowvec
    % varargin: variable input arguments for property values non-empty
    % NlpVariables @copydoc NlpVariable::updateProp
    %
    % @see NlpVariable, removeVariable

    varnames = obj.OptVarTable.Properties.VariableNames;

    %     if ismember(label,varnames)
    %         warning('The NLP variable (%s) already exists.\n Overwriting the existing NLP variable.', label);
    %     end

    if ischar(nodes)
        switch nodes
            case 'first'
                node_list = 1;
            case 'last'
                node_list = obj.NumNode;
            case 'except-first'
                node_list = 2:obj.NumNode;
            case 'except-last'
                node_list = 1:obj.NumNode-1;
            case 'except-terminal'
                node_list = 2:obj.NumNode-1;
            case 'all'
                node_list = 1:obj.NumNode;
            case 'cardinal'
                node_list = 1:2:obj.NumNode;
            case 'interior'
                node_list = 2:2:obj.NumNode-1;
            otherwise
                error('Unknown node type.');
        end
    else
        if ~isnumeric(nodes)
            error(['The node must be specified as a list of node indices or following supported characters:\n',...
                '%s'],implode({'first','last','all','except-first','except-last','except-terminal', 'cardinal', 'interior'},','));
        else
            node_list = nodes;
        end
    end

    % create empty NlpVariable array
    var = repmat(NlpVariable(),obj.NumNode,1);
    for i=node_list
        var(i) = NlpVariable(varargin{:});
    end
    % add to the decision variable table
    obj.OptVarTable.(label) = var;

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function mObject3 = timeFrequencySTRAIGHTmorphingExt(mObject1,mObject2,mRate,mixMethod);
%   Morphing based on STRAIGHT parameters
%   mObject3 = timeFrequencySTRAIGHTmorphing(mObject1,mObject2,mRate,mixMethod);

%   Designed and coded by Hideki Kawahara
%   28/Feb./2005
%   Copyright(c) 2005, Hideki Kawahara
%   14/March/2005 bug fix on sampling frequency
%   01/Oct./2005 bug fix on similarity check
%   04/Oct./2005 partial morphing extension
%   18/Oct./2005 direct differential manipulation and API cange
%   29/Jan./2006 bug fix on boundary conditions
%   24/Oct./2006 modificaton of definition

switch nargin
    case 0
        mObject3.morphingObject = createMobject;
        mixRate.F0 = 0;
        mixRate.spectrum = 0;
        mixRate.aperiodicity = 0;
        mixRate.timeCoordinate = 0;
        mixRate.freqCoordinate = 0;
        mObject3.mixRate = mixRate;
        mObject3.mixMethods = {'linear','log','differentialLogarithm'};
        return
end;
if ~isfield(mObject1,'vuv')
    mObject1.vuv = (mObject1.F0>0);
elseif length(mObject1.vuv) == 0
    mObject1.vuv = (mObject1.F0>0);
end;
if ~isfield(mObject2,'vuv')
    mObject2.vuv = (mObject2.F0>0);
elseif length(mObject2.vuv) == 0
    mObject2.vuv = (mObject2.F0>0);
end;
mObject3 = checkForSimilarity(mObject1,mObject2);
mixRate = checkForMorphingConditions(mRate);
mObject1 = checkForIntegrity(mObject1);
mObject2 = checkForIntegrity(mObject2);
fs = mObject1.samplingFrequency;
if length(mObject3) ==0;return;end;
dtFrame = mObject1.frameUpdateInterval;
endLocation1 = (length(mObject1.F0)-1)*dtFrame; % in ms
endLocation2 = (length(mObject2.F0)-1)*dtFrame; % in ms
timeAnchor1 = [0;mObject1.anchorTimeLocation;endLocation1];
timeAnchor2 = [0;mObject2.anchorTimeLocation;endLocation2];
locationOn1 = (0:length(mObject1.F0)-1)*dtFrame;
locationOn2 = (0:length(mObject2.F0)-1)*dtFrame;
mapFrom1to2 = interp1(timeAnchor1,timeAnchor2,locationOn1);
[nr1,nc1] = size(mObject1.spectrogram);
[nr2,nc2] = size(mObject2.spectrogram);

%---- initialize frequency mapping function
fmapFrom1to2OnTime1 = generateFrequencyMap(mObject1,mObject2);

%---- mixing on mObject1's time axis
nAxis1 = length(locationOn1);
nAxis2 = length(locationOn2);
morphedF0 = zeros(nAxis1,1);
morphedvuv = zeros(nAxis1,1);
morphedAp = zeros(nr1,nAxis1);
morphedSgram = zeros(nr1,nAxis1);
weightSumF0 = zeros(nAxis1,1);
for ii=1:nAxis1
    mappedIndexOn2 = mapFrom1to2(ii)/dtFrame+1;
    iFloor = floor(mappedIndexOn2);
    iFraction = mappedIndexOn2-iFloor;
    fIndex = floor(fmapFrom1to2OnTime1(:,ii)/fs*2*(nr1-1))+1;
    dAp = iFraction*(mObject2.aperiodicityIndex(:,min(iFloor+1,nAxis2))-mObject2.aperiodicityIndex(:,min(iFloor,nAxis2)));
    ap2on2faxis = mObject2.aperiodicityIndex(:,min(iFloor,nAxis2))+dAp;
    ap2on1faxis = ap2on2faxis(fIndex);
    morphedAp(:,ii) = (1-mixRate.aperiodicity)*mObject1.aperiodicityIndex(:,ii)+mixRate.aperiodicity*ap2on1faxis; %04/Oct/2005 HK
    switch mixMethod
        case 'linear'
            dSgram = iFraction*(mObject2.spectrogram(:,min(iFloor+1,nAxis2))-mObject2.spectrogram(:,min(iFloor,nAxis2)));
            sgram2on2faxis = mObject2.spectrogram(:,min(iFloor,nAxis2))+dSgram;
            sgram2on1faxis = sgram2on2faxis(fIndex);
            morphedSgram(:,ii) = (1-mixRate.spectrum)*mObject1.spectrogram(:,ii)+mixRate.spectrum*sgram2on1faxis;
        case 'log'
            dSgram = iFraction*(log(mObject2.spectrogram(:,min(iFloor+1,nAxis2)))-log(mObject2.spectrogram(:,min(iFloor,nAxis2))));
            sgram2on2faxis = log(mObject2.spectrogram(:,min(iFloor,nAxis2)))+dSgram;
            sgram2on1faxis = sgram2on2faxis(fIndex);
            tmp = (1-mixRate.spectrum)*log(mObject1.spectrogram(:,ii))+mixRate.spectrum*sgram2on1faxis;
            morphedSgram(:,ii) = exp(tmp);
        case 'differentialLogarithm'
            dSgram = iFraction*(mObject2.spectrogram(:,min(iFloor+1,nAxis2))-mObject2.spectrogram(:,min(iFloor,nAxis2)));
            sgram2on2faxis = mObject2.spectrogram(:,min(iFloor,nAxis2))+dSgram;
            sgram2on1faxis = sgram2on2faxis(fIndex);
            tmp = (1-mixRate.spectrum)*log(mObject1.spectrogram(:,ii))+mixRate.spectrum*sgram2on1faxis;
            morphedSgram(:,ii) = exp(tmp);
    end;
    if mObject1.F0(ii)>0
        morphedF0(ii) = (1-mixRate.F0)*log(mObject1.F0(ii));
        weightSumF0(ii) = (1-mixRate.F0);
    end;
    if (mObject2.F0(iFloor)>0) & (mObject2.F0(min(iFloor+1,nAxis2))>0)
        dF0 = iFraction*(log(mObject2.F0(min(iFloor+1,nAxis2)))-log(mObject2.F0(min(iFloor,nAxis2))));
        morphedF0(ii) = mixRate.F0*(log(mObject2.F0(min(iFloor,nAxis2)))+dF0)+morphedF0(ii);
        weightSumF0(ii) = weightSumF0(ii)+mixRate.F0;
    end;
    morphedvuv(ii) = ((mObject1.vuv(ii)*abs(1-mixRate.F0)+abs(mixRate.F0)*mObject2.vuv(min(iFloor,nAxis2)))>0);
end;
morphedF0(weightSumF0>0) = exp(morphedF0(weightSumF0>0)./weightSumF0(weightSumF0>0));

%----- mapping back onto morphed time axis
timeAnchorMorph = (1-mixRate.timeCoordinate)*timeAnchor1 + mixRate.timeCoordinate*timeAnchor2;
locationOnMorph = (0:(timeAnchorMorph(end)/dtFrame))*dtFrame;
mapFormMorphTo1 = interp1(timeAnchorMorph,timeAnchor1,locationOnMorph);
nAxisMorph = length(locationOnMorph);
morphedApOnMorph = zeros(nr1,nAxisMorph);
morphedSgramOnMorph = zeros(nr1,nAxisMorph);
morphedF0onMorph = zeros(nAxisMorph,1);
morphedVUVonMorph = zeros(nAxisMorph,1);
%----- set place holders
mObject3.samplingFrequency = fs;
mObject3.F0 = morphedF0onMorph;
mObject3.vuv = morphedVUVonMorph;
mObject3.aperiodicityIndex = morphedApOnMorph;
mObject3.spectrogram = morphedSgramOnMorph;
mObject3.anchorTimeLocation = timeAnchorMorph(2:end-1);
mObject3.anchorFrequency = (1-mixRate.freqCoordinate)*mObject1.anchorFrequency+mixRate.freqCoordinate*mObject2.anchorFrequency;
%------ nitialize frequency mapping function
fmapFromMorphto1OnTimeMorph = generateFrequencyMap(mObject3,mObject1);
for ii=1:nAxisMorph
    mappedIndexOn1 = mapFormMorphTo1(ii)/dtFrame+1;
    iFloor = floor(mappedIndexOn1);
    iFraction = mappedIndexOn1-iFloor;
    fIndex = floor(fmapFromMorphto1OnTimeMorph(:,ii)/fs*2*(nr1-1))+1;
    morphedApOnMorph(:,ii) = morphedAp(fIndex,iFloor) ...
        +iFraction*(morphedAp(fIndex,min(iFloor+1,nAxis1))-morphedAp(fIndex,iFloor));
    morphedSgramOnMorph(:,ii) = morphedSgram(fIndex,iFloor) ...
        +iFraction*(morphedSgram(fIndex,min(iFloor+1,nAxis1))-morphedSgram(fIndex,iFloor));
    if (morphedF0(iFloor)>0) & (morphedF0(min(iFloor+1,nAxis1))>0)
        dF0 = iFraction*(morphedF0(min(iFloor+1,nAxis1))-morphedF0(iFloor));
        morphedF0onMorph(ii) = morphedF0(iFloor)+dF0;
    end;
    morphedVUVonMorph(ii) = morphedvuv(iFloor);
end;
mObject3.F0 = morphedF0onMorph;
mObject3.vuv = morphedVUVonMorph;
mObject3.aperiodicityIndex = morphedApOnMorph;
mObject3.spectrogram = morphedSgramOnMorph;
mObject3.anchorTimeLocation = timeAnchorMorph(2:end-1);
%mObject3.anchorFrequency = (1-mRate)*mObject1.anchorFrequency+mRate*mObject2.anchorFrequency;
%mObject3 = fmapFromMorphto1OnTimeMorph; % This line is a dummy.
return;

%%% ------ Internal function to check for object's similarity
function mObject3 = checkForSimilarity(mObject1,mObject2)
mObject3 = [];
if mObject1.samplingFrequency ~= mObject2.samplingFrequency;mObject3 = [];return;end;
if mObject1.frameUpdateInterval ~= mObject2.frameUpdateInterval;mObject3 = [];return;end;
if length(mObject1.anchorTimeLocation) ~= length(mObject2.anchorTimeLocation);mObject3 = [];return;end;
nAnchor = length(mObject1.anchorTimeLocation);
for ii=1:nAnchor % check for similarity of anchor structure
    frequencyAnchor1 = mObject1.anchorFrequency(ii,:)';% 01/Oct./2005 by HK
    frequencyAnchor2 = mObject2.anchorFrequency(ii,:)';% 01/Oct./2005 by HK
    if (sum(frequencyAnchor1>0) ~= sum(frequencyAnchor2>0)) | ...
            (sum(frequencyAnchor1<0) ~= sum(frequencyAnchor2<0))
        display('Warning!! Object structures are inconsistent!'); % 01/Oct./2005 by HK
        return;
    end;
end;
mObject3 = createMobject;
m0bject3.samplingFrequency = mObject1.samplingFrequency;
m0bject3.frameUpdateInterval = mObject1.frameUpdateInterval;
return;

%%%--------
function mixRate = checkForMorphingConditions(mRate);
%   04/Oct./2005 added by HK

if ~isstruct(mRate)
    mixRate.F0 = mRate;
    mixRate.spectrum = mRate;
    mixRate.aperiodicity = mRate;
    mixRate.timeCoordinate = mRate;
    mixRate.freqCoordinate = mRate;
    return;
end;
mixRate.F0 = mRate.F0;
mixRate.spectrum = mRate.spectrum;
mixRate.aperiodicity = mRate.aperiodicity;
mixRate.timeCoordinate = mRate.timeCoordinate;
mixRate.freqCoordinate = mRate.freqCoordinate;
return;

%%%--------
function fmapFrom1to2OnTime1 = generateFrequencyMap(mObject1,mObject2);

dtFrame = mObject1.frameUpdateInterval;
endLocation1 = (length(mObject1.F0)-1)*dtFrame; % in ms
timeAnchor1 = [0;mObject1.anchorTimeLocation;endLocation1];
locationOn1 = (0:length(mObject1.F0)-1)*dtFrame;
fs = mObject1.samplingFrequency;
[nr1,nc1] = size(mObject1.spectrogram);
nAnchor = length(mObject1.anchorTimeLocation);
fmapFrom1to2 = zeros(nr1,nAnchor);
frequencyAxis = (0:nr1-1)'/(nr1-1)*fs/2;
numberOfFrequencyAnchors = zeros(nAnchor,1);
for ii=1:nAnchor
    frequencyAnchor1 = mObject1.anchorFrequency(ii,:)';
    frequencyAnchor1 = [0;frequencyAnchor1(frequencyAnchor1>0);fs/2];
    numberOfFrequencyAnchors(ii) = length(frequencyAnchor1(frequencyAnchor1>0));
    frequencyAnchor2 = mObject2.anchorFrequency(ii,:)';
    frequencyAnchor2 = [0;frequencyAnchor2(frequencyAnchor2>0);fs/2];
    fmapFrom1to2(:,ii) = interp1(frequencyAnchor1,frequencyAnchor2,frequencyAxis);
end;
for ii=1:nAnchor
    if numberOfFrequencyAnchors(ii) == 1
        if numberOfFrequencyAnchors(min(ii+1,nAnchor)) > 1
            fmapFrom1to2(:,ii) = fmapFrom1to2(:,min(ii+1,nAnchor));
        elseif numberOfFrequencyAnchors(max(ii-1,1)) > 1
            fmapFrom1to2(:,ii) = fmapFrom1to2(:,max(ii-1,1));
        end;
    end;
end;
fmapFrom1to2 = [fmapFrom1to2(:,1) fmapFrom1to2 fmapFrom1to2(:,nAnchor)];
fmapFrom1to2OnTime1 = interp1(timeAnchor1,fmapFrom1to2',locationOn1)';
return;

%%%-------
function cleanedUpObject = checkForIntegrity(inputObject);

maximumIndex = max([length(inputObject.F0), ...
    size(inputObject.spectrogram,2) ...
    size(inputObject.aperiodicityIndex,2)]);
if length(inputObject.F0) < maximumIndex
    inputObject.F0 = [inputObject.F0(:);inputObject.F0(end)*ones(maximumIndex - length(inputObject.F0),1)];
end;
if size(inputObject.spectrogram,2) < maximumIndex
    numberOfFillIn = maximumIndex-size(inputObject.spectrogram,2);
    inputObject.spectrogram = [inputObject.spectrogram inputObject.spectrogram(:,end)*ones(1,numberOfFillIn)];
end;
if size(inputObject.aperiodicityIndex,2) < maximumIndex
    numberOfFillIn = maximumIndex-size(inputObject.aperiodicityIndex,2);
    inputObject.aperiodicityIndex = [inputObject.aperiodicityIndex inputObject.aperiodicityIndex(:,end)*ones(1,numberOfFillIn)];
end;
cleanedUpObject = inputObject;

author: HidekiKawahara


text:
function [ all_bboxes_gt ] = get_grount_truth_bboxes_from_voc( voc_path, image_set, voc_year, with_hard_samples, cache_dir )
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

if ~exist('cache_dir','var')
    cache_dir = '.';
end

VOCopts            = initVOCOpts( voc_path, voc_year );
fullimage_set_name = ['voc_', voc_year, '_' image_set];
addpath([voc_path, filesep, 'VOCcode']);
cache_file = [cache_dir, filesep, 'gt_bbox_', fullimage_set_name];
if with_hard_samples, cache_file = [cache_file, '_with_hard_samples']; end
cache_file = [cache_file, '.mat'];

if exist(cache_file,'file')
    all_bboxes_gt = loadGroundTruthBBoxes(cache_file);
else
    class_to_id     = containers.Map(VOCopts.classes, 1:length(VOCopts.classes));
    image_ids       = textread(sprintf(VOCopts.imgsetpath, image_set), '%s');
    num_imgs        = length(image_ids);
    all_bboxes_gt   = cell(num_imgs, 1);
    for img_idx = 1:num_imgs
        if mod(img_idx,500) == 0, fprintf('Load gt bboxes::%s %d/%d\n', fullimage_set_name, img_idx, num_imgs); end
        all_bboxes_gt{img_idx} = getGroundTruthBBoxes(sprintf(VOCopts.annopath, image_ids{img_idx}), class_to_id, with_hard_samples);
    end
    saveGroundTruthBBoxes(cache_file, all_bboxes_gt);
end

end

function gt_bboxes = getGroundTruthBBoxes(filename, class_to_id, with_hard_samples)
try
    voc_rec       = PASreadrecord(filename);
    valid_objects = 1:length(voc_rec.objects(:));
    is_difficult  = cat(1, voc_rec.objects(:).difficult);
   if ~with_hard_samples
       valid_objects = valid_objects(~is_difficult);
       is_difficult  = is_difficult(~is_difficult);
   end

   gt_bboxes = single(cat(1, voc_rec.objects(valid_objects).bbox));
   gt_class_idx = class_to_id.values({voc_rec.objects(valid_objects).class});
   gt_class_idx = single(cat(1, gt_class_idx{:}));
   gt_bboxes    = single([gt_bboxes, gt_class_idx, is_difficult(:)]);
catch
    gt_bboxes = zeros(0, 6, 'single');
end

end

function all_bboxes_gt = loadGroundTruthBBoxes(filename)
load(filename, 'all_bboxes_gt');
end

function saveGroundTruthBBoxes(filename, all_bboxes_gt)
save(filename, 'all_bboxes_gt', '-v7.3');
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function [transformationMatrix] = ...
    optical_flow_linear_registration3d(moving, fixed, varargin)
% OPTICAL_FLOW_LINEAR_REGISTRATION3D Estimates a transformation matrix using optical flow
%
% INPUT ARGUMENTS
% moving                - Moving iamge
% fixed                 - Fixed image
%
% OUTPUT ARGUMENTS
% transformationMatrix  - Estimated displacement field
%
% OPTIONAL INPUT ARGUMENTS
% 'transformationModel'     - Transformation model for estimating the
%                             displacement field
%                             'translation', 'affine', 'non-rigid' (default)
%
% 'multiModal'              - Set whether to perform multi-modal or
%                             uni-modal image registration
%                             false (default), true
%
% 'numberOfChannels'        - Number of channels to use in when computing
%                             the entropy (based on channel coding). This
%                             is only relevant if multiModal is set to
%                             true.
%                             Default value is 8

% Copyright (c) 2012
% danne.forsberg@outlook.com
%
% 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/>.

%% Setup default parameters
% translation, affine
transformationModel = 'affine';

% multi-modal
multiModal = false;

% number of channels, only valid for multi-modal registration
numberOfChannels = 8;

% Overwrites default parameter
for k=1:2:length(varargin)
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

% Initialize transformation matrix
transformationMatrix = eye(4);

[b_moving(:,:,:,1) b_moving(:,:,:,2) b_moving(:,:,:,3)] = gradient(moving);
[b_fixed(:,:,:,1) b_fixed(:,:,:,2) b_fixed(:,:,:,3)] = gradient(fixed);

if multiModal
    b = b_fixed(:,:,:,2);
    b(:,:,:,2) = b_fixed(:,:,:,1);
    b(:,:,:,3) = b_fixed(:,:,:,3);

    [delta_c mask] = estimate_delta_c(fixed,moving,numberOfChannels);
    delta_c(mask ~= 1) = 0;
    mask = repmat(mask,[1 1 3]);
    b(mask ~= 1) = 0;
else
    b = (b_moving(:,:,:,2) + b_fixed(:,:,:,2))/2;
    b(:,:,:,2) = (b_moving(:,:,:,1) + b_fixed(:,:,:,1))/2;
    b(:,:,:,3) = (b_moving(:,:,:,3) + b_fixed(:,:,:,3))/2;

    delta_c = fixed - moving;
end

[G, h] = build_G_h_linear3d(b, delta_c, transformationModel);

% Solve the equation system
switch transformationModel
    case 'translation'
        d = G \ h;
        transformationMatrix(1:3,4) = d;
    case {'rigid','affine'}
        p = G \ h;
        transformationMatrix(1:3,1:4) = [1+p(1) p(2) p(3) p(10);...
                                         p(4) 1+p(5) p(6) p(11);...
                                         p(7) p(8) 1+p(9) p(12)];
end

author: fordanic

repo: image-registration


text:
function chebyshev_polynomial_test03 ( )

%*****************************************************************************80
%
%% CHEBYSHEV_POLYNOMIAL_TEST03 tests T_POLYNOMIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_POLYNOMIAL_TEST03:\n' );
  fprintf ( 1, '  T_POLYNOMIAL_VALUES stores values of\n' );
  fprintf ( 1, '  the Chebyshev polynomials.\n' );
  fprintf ( 1, '  T_POLYNOMIAL evaluates the polynomial.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                        Tabulated                 Computed\n' );
  fprintf ( 1, '     N        X           T(n,x)                    T(n,x)                     Error\n' );

  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx1 ] = t_polynomial_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2_vec = t_polynomial ( 1, n, x );
    fx2 = fx2_vec(1,n+1);
    e = fx1 - fx2;

    fprintf ( 1, '  %4d  %12f  %24.16e  %24.16e  %8.2g\n', n, x, fx1, fx2, e );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% EBSD
%
%%


text:
function gui_LoadResult
% function pf_load
%   Opens a open dialog to load result data
%   and invokes the display dialog
%
% Danijel Schorlemmer
% August 22, 2002

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Invoke the open dialog
[newfile, newpath] = uigetfile('*.mat', 'Load calculated data');
% Cancel pressed?
if newfile == 0
  return;
end
% Everything ok?
newfile = [newpath newfile];
if length(newfile) > 1
  % Load the data
  load(newfile, 'vResults');
  % Open dialog
  hResultFig = gui_result2(vResults);
end

author: CelsoReyes


text:
function r = mldivide(p,q)
%MLDIVIDE     Implements  p \ q  for univariate polynomials (same as q/p)
%

% written  11/02/05     S.M. Rump
%

  r = q / p;

repo: OpenMAS


text:
function gencost = modcost(gencost, alpha, modtype)
%MODCOST  Modifies generator costs by shifting or scaling (F or X).
%   NEWGENCOST = MODCOST(GENCOST, ALPHA)
%   NEWGENCOST = MODCOST(GENCOST, ALPHA, MODTYPE)
%
%   For each generator cost F(X) (for real or reactive power) in
%   GENCOST, this function modifies the cost by scaling or shifting
%   the function by ALPHA, depending on the value of MODTYPE, and
%   and returns the modified GENCOST. Rows of GENCOST can be a mix
%   of polynomial or piecewise linear costs. ALPHA can be a scalar,
%   applied to each row of GENCOST, or an NG x 1 vector, where each
%   element is applied to the corresponding row of GENCOST.
%
%   MODTYPE takes one of the 4 possible values (let F_alpha(X) denote the
%   the modified function):
%       'SCALE_F' (default) : F_alpha(X)         == F(X) * ALPHA
%       'SCALE_X'           : F_alpha(X * ALPHA) == F(X)
%       'SHIFT_F'           : F_alpha(X)         == F(X) + ALPHA
%       'SHIFT_X'           : F_alpha(X + ALPHA) == F(X)

%   MATPOWER
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define named indices into data matrices
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

if nargin < 3
    modtype = 'SCALE_F';
end

[ng, m] = size(gencost);

if ng ~= 0
    if length(alpha) ~= ng
        if length(alpha) == 1 && ng > 1     %% scalar, make it a col vector
            alpha = alpha * ones(ng, 1);
        else
            error('modcost: ALPHA must be a scalar or col vector with NG rows');
        end
    elseif size(alpha, 2) ~= 1
        alpha = alpha';                     %% convert row vector to col vector
    end

    ipwl = find(gencost(:, MODEL) == PW_LINEAR);
    ipol = find(gencost(:, MODEL) == POLYNOMIAL);
    npwl = length(ipwl);
    npol = length(ipol);
    c = gencost(ipol, COST:m);

    switch modtype
        case 'SCALE_F',
            if npol
                gencost(ipol, COST:m) = spdiags(alpha(ipol), 0, npol, npol) * c;
            end
            if npwl
                gencost(ipwl, COST+1:2:m) = spdiags(alpha(ipwl), 0, npwl, npwl) * ...
                    gencost(ipwl, COST+1:2:m);
            end
        case 'SCALE_X',
            for k = 1:length(ipol)
                n = gencost(ipol(k), NCOST);
                for i = 1:n
                    gencost(ipol(k), COST+i-1) = c(k, i) / alpha(ipol(k))^(n-i);
                end
            end
            if npwl
                gencost(ipwl, COST:2:m-1) = spdiags(alpha(ipwl), 0, npwl, npwl) * ...
                    gencost(ipwl, COST:2:m-1);
            end
        case 'SHIFT_F',
            for k = 1:length(ipol)
                n = gencost(ipol(k), NCOST);
                gencost(ipol(k), COST+n-1) = alpha(ipol(k)) + c(k, n);
            end
            if npwl
                gencost(ipwl, COST+1:2:m) = spdiags(alpha(ipwl), 0, npwl, npwl) * ...
                    ones(length(ipwl), (m+1-COST)/2) + gencost(ipwl, COST+1:2:m);
            end
        case 'SHIFT_X',
            for k = 1:length(ipol)
                n = gencost(ipol(k), NCOST);
                gencost(ipol(k), COST:COST+n-1) = polyshift(c(k, 1:n)', alpha(ipol(k)))';
            end
            if npwl
                gencost(ipwl, COST:2:m-1) = spdiags(alpha(ipwl), 0, npwl, npwl) * ...
                    ones(length(ipwl), (m+1-COST)/2) + gencost(ipwl, COST:2:m-1);
            end
        otherwise
            error('modcost: ''%s'' is not a valid modtype\n', modtype);
    end
end

%%-----  POLYSHIFT  -----
function d = polyshift(c, a)
%POLYSHIFT  Returns the coefficients of a horizontally shifted polynomial.
%
%   D = POLYSHIFT(C, A) shifts to the right by A, the polynomial whose
%   coefficients are given in the column vector C.
%
%   Example: For any polynomial with n coefficients in c, and any values
%   for x and shift a, the f - f0 should be zero.
%       x = rand;
%       a = rand;
%       c = rand(n, 1);
%       f0 = polyval(c, x)
%       f  = polyval(polyshift(c, a), x+a)

n = length(c);
d = zeros(size(c));
A = (-a * ones(n, 1)) .^ ((0:n-1)');
b = ones(n, 1);
for k = 1:n
    d(n-k+1) = b' * ( c(n-k+1:-1:1) .* A(1:n-k+1) );
    b = cumsum(b(1:n-k));
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower