text:
function transformMap2XML(xmlStruct, map, fileName)
% Creates a new XML file from the information contained in the map
% structure. Uses the function "struct2xml" to transform a matlab
% structure into an XML text format.
%
% USAGE:
%
%   transformMap2XML(xmlStruct, map, fileName)
%
% INPUTS:
%   xmlStruct:      XML structure obtained from the function
%                   "xml2struct". Used by the function "struct2xml"
%                   to obtain the XML file.
%   map:            Matlab structure of the map with the relevant
%                   information. This information is then transfered to
%                   the xmlStruct for the conversion.
%   fileName:       Path and name of the new XML file.
%
% .. Author: - N.Sompairac - Institut Curie, Paris, 24/07/2017

tic

if strcmp('~',fileName(1))
    savedPath=pwd;
    cd('~')
    homePath=pwd;
    fileName = strrep(fileName,'~',homePath);
    cd(savedPath)
end

% Loop over molecules to put the info back to the struct for XML
% conversion
for mol = 1:length(map.molAlias)
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.Attributes.id = map.molAlias{mol};
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.Attributes.species = map.molID{mol};
    % Check if the info on compartment exists
    if ~isempty(map.molCompartAlias{mol})
        xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.Attributes.compartmentAlias = map.molCompartAlias{mol};
    end
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.x = map.molXPos{mol};
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.y = map.molYPos{mol};
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.w = map.molWidth{mol};
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.h = map.molHeight{mol};
    xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_usualView.celldesigner_colon_paint.Attributes.color = map.molColor{mol};
end
clearvars mol

% Loop over species to put the info back to the struct for XML
% conversion
for spec = 1:length(map.specID)
    xmlStruct.sbml.model.listOfSpecies.species{spec}.Attributes.id = map.specID{spec};
    xmlStruct.sbml.model.listOfSpecies.species{spec}.Attributes.metaid = map.specMetaID{spec};
    xmlStruct.sbml.model.listOfSpecies.species{spec}.Attributes.name = map.specName{spec};
    xmlStruct.sbml.model.listOfSpecies.species{spec}.annotation.celldesigner_colon_extension.celldesigner_colon_speciesIdentity.celldesigner_colon_name.Text = map.specName{spec};
    xmlStruct.sbml.model.listOfSpecies.species{spec}.annotation.celldesigner_colon_extension.celldesigner_colon_speciesIdentity.celldesigner_colon_class.Text = map.specType{spec};
    % Check if there are notes for the specie
    if ~isempty(map.specNotes{spec})
        xmlStruct.sbml.model.listOfSpecies.species{spec}.notes.html.body.Text = map.specNotes{spec};
    end
end
clearvars spec

for r = 1:length(map.rxnID)
    lt=length(xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfReactants.speciesReference);
    if lt==1
        if strcmp(xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfReactants.speciesReference.annotation.celldesigner_colon_extension.celldesigner_colon_alias.Text,map.molAlias)
            error('removing only substrate in reaction')
        end
    else
        bool=true(lt,1);
        for s=1:lt
            if ~any(strcmp(xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfReactants.speciesReference{s}.annotation.celldesigner_colon_extension.celldesigner_colon_alias.Text,map.molAlias))
                bool(s)=0;
            end
        end
        xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfReactants.speciesReference =xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfReactants.speciesReference(bool);
    end
%     if ~iscell(xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfReactants.speciesReference)
%         xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfReactants.speciesReference
%     end

    lt=length(xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfProducts.speciesReference);
    if lt==1
        if strcmp(xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfProducts.speciesReference.annotation.celldesigner_colon_extension.celldesigner_colon_alias.Text,map.molAlias)
            error('removing only product in reaction')
        end
    else
        bool=true(lt,1);
        for s=1:lt
            if ~any(strcmp(xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfProducts.speciesReference{s}.annotation.celldesigner_colon_extension.celldesigner_colon_alias.Text,map.molAlias))
                bool(s)=0;
            end
        end
        xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfProducts.speciesReference =xmlStruct.sbml.model.listOfReactions.reaction{r}.listOfProducts.speciesReference(bool);
    end

end

% Loop over reactions to put the info back to the struct for XML
% conversion
for react = 1:length(map.rxnID)
    % Check if there is a reaction ID
    if ~isempty(map.rxnID{react})
        xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.id = map.rxnID{react};
    else
        xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.id = 'No_ID';
    end
    xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.metaid = map.rxnMetaID{react};

    % Check if the reaction has a Name
    if ~isempty(map.rxnName{react})
        xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.name = map.rxnName{react};
        xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_name.Text = map.rxnName{react};
    else
        xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.name = 'No_Name';
        xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_name.Text = 'No_Name';
    end
    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_reactionType.Text = map.rxnType{react};
    if strcmp(map.rxnReversibility{react},'false')
        xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.reversible = map.rxnReversibility{react};
    end
    % Test if there is only 1 base reactant
    if length(map.rxnBaseReactantAlias{react}) == 1
        xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant(1).Attributes.alias = map.rxnBaseReactantAlias{react}{1};
        xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant(1).Attributes.species = map.rxnBaseReactantID{react}{1};
    else
        for x = 1:length(map.rxnBaseReactantAlias{react})
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant{x}.Attributes.alias = map.rxnBaseReactantAlias{react}{x};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant{x}.Attributes.species = map.rxnBaseReactantID{react}{x};
        end
        clearvars x
    end

    % Test if there is only 1 base product
    if length(map.rxnBaseProductAlias{react}) == 1
        xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct(1).Attributes.alias = map.rxnBaseProductAlias{react}{1};
        xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct(1).Attributes.species = map.rxnBaseProductID{react}{1};
    else
        for x = 1:length(map.rxnBaseProductAlias{react})
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct{x}.Attributes.alias = map.rxnBaseProductAlias{react}{x};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct{x}.Attributes.species = map.rxnBaseProductID{react}{x};
        end
    end

    % Test if there are some reactants in the reaction
    if ~isempty(map.rxnReactantAlias{react})
        % Test if there is only 1 reactant
        if length(map.rxnReactantAlias{react}) == 1
            %needed if there is a reduction in number of reactant aliases
            if iscell(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink)
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink=xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{1};
            end
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(1).Attributes.alias = map.rxnReactantAlias{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(1).Attributes.reactant = map.rxnReactantID{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(1).celldesigner_colon_line.Attributes.color = map.rxnColor{react};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(1).celldesigner_colon_line.Attributes.width = map.rxnWidth{react};
        else
            for x = 1:length(map.rxnReactantAlias{react})
                if length(map.rxnReactantAlias{react})== length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink)-1
                    error('inconsistent reactant aliases')
                end
                if iscell(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink)
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{x}.Attributes.alias = map.rxnReactantAlias{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{x}.Attributes.reactant = map.rxnReactantID{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{x}.celldesigner_colon_line.Attributes.color = map.rxnColor{react};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{x}.celldesigner_colon_line.Attributes.width = map.rxnWidth{react};
                else
                    error('should not be cell')
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(x).Attributes.alias = map.rxnReactantAlias{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(x).Attributes.reactant = map.rxnReactantID{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(x).celldesigner_colon_line.Attributes.color = map.rxnColor{react};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(x).celldesigner_colon_line.Attributes.width = map.rxnWidth{react};
                end
            end
            clearvars x
        end
    end

    % Test if there are some products in the reaction
    if ~isempty(map.rxnProductAlias{react})
        % Test if there is only 1 product
        if length(map.rxnProductAlias{react}) == 1
            %needed if there is a reduced number of product aliases
            if iscell(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink)
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink=xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{1};
            end
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(1).Attributes.alias = map.rxnProductAlias{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(1).Attributes.product = map.rxnProductID{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(1).celldesigner_colon_line.Attributes.color = map.rxnColor{react};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(1).celldesigner_colon_line.Attributes.width = map.rxnWidth{react};
        else
            for x = 1:length(map.rxnProductAlias{react})
                if length(map.rxnProductAlias{react})== length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink)-1
                    error('inconsistent product aliases')
                end
                if iscell(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink)
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{x}.Attributes.alias = map.rxnProductAlias{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{x}.Attributes.product = map.rxnProductID{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{x}.celldesigner_colon_line.Attributes.color = map.rxnColor{react};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{x}.celldesigner_colon_line.Attributes.width = map.rxnWidth{react};
                else
                    error('should not be cell')
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(x).Attributes.alias = map.rxnProductAlias{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(x).Attributes.product = map.rxnProductID{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(x).celldesigner_colon_line.Attributes.color = map.rxnColor{react};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(x).celldesigner_colon_line.Attributes.width = map.rxnWidth{react};
                end
            end
            clearvars x
        end
    end

    % Check if there are modifications for the reaction
    if ~isempty(map.rxnModAlias{react})
        % Test if there is only one modification
        if length(map.rxnModAlias{react}) == 1
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).Attributes.aliases = map.rxnModAlias{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).Attributes.modifiers = map.rxnModID{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).Attributes.type = map.rxnModType{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).celldesigner_colon_line.Attributes.color = map.rxnModColor{react}{1};
            xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).celldesigner_colon_line.Attributes.width = map.rxnModWidth{react}{1};
        else
            for m = 1:length(map.rxnModAlias{react})
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{m}.Attributes.aliases = map.rxnModAlias{react}{m};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{m}.Attributes.modifiers = map.rxnModID{react}{m};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{m}.Attributes.type = map.rxnModType{react}{m};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{m}.celldesigner_colon_line.Attributes.color = map.rxnModColor{react}{m};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{m}.celldesigner_colon_line.Attributes.width = map.rxnModWidth{react}{m};
            end
            clearvars m
        end
    end

    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_line.Attributes.color = map.rxnColor{react};
    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_line.Attributes.width = map.rxnWidth{react};

    % Check if there are notes for the reaction
    if ~isempty(map.rxnNotes{react})
        xmlStruct.sbml.model.listOfReactions.reaction{react}.notes.html.body.Text = map.rxnNotes{react};
    end
end
clearvars react

% Check if there is any information on Compartments
if ~isempty(map.compartAlias)
    % Test if there is only 1 compartment
    if length(map.compartAlias) == 1
        xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias(1).Attributes.id = map.compartAlias{1};
        xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias(1).Attributes.compartment = map.compartName{1};
    else
        % Loop over compartments to put the info back to the struct for XML
        % conversion
        for compart = 1:length(map.compartAlias)
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias{compart}.Attributes.id = map.compartAlias{compart};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias{compart}.Attributes.compartment = map.compartName{compart};
        end
        clearvars compart
    end
end

struct2xml(xmlStruct, fileName);
toc


text:
% Lightspeed Toolbox.
% Efficient operations for Matlab programming.
% Version 2.1   10-Oct-2006
% By Tom Minka
% (c) Microsoft Corporation. All rights reserved.
%
% Matrix algebra
%   repmat           - Fast replacement for matlab's repmat.
%   xrepmat          - Matlab's original repmat.
%   row_sum          - Sum for each row.  Faster than 'sum'.
%   scale_rows       - Scale each row of a matrix.
%   scale_cols       - Scale each column of a matrix.
%   solve_triu       - Left division by upper triangular matrix.
%   solve_tril       - Left division by lower triangular matrix.
%   sqdist           - Squared Euclidean and Mahalanobis distance.
%   isposdef         - Check for positive-definiteness.
%   logdet           - log(determinant) for positive definite matrix.
%
% Statistics
%   normpdf          - Multivariate normal density.
%   normpdfln        - Log of multivariate normal density.
%   normcdf          - Normal cumulative distribution.
%   normcdfln        - Log of normal cumulative distribution.
%   normcdflogit     - Logit of normal cumulative distribution.
%   wishpdf          - Wishart probability density function.
%   wishpdfln        - Log of Wishart probability density function.
%   sample           - Sample from categorical distribution.
%   sample_vector    - Sample from multiple categorical distributions.
%   sample_hist      - Sample from multinomial distribution.
%   randbinom        - Sample from binomial distribution.
%   randnorm         - Sample from multivariate normal.
%   randgamma        - Sample from Gamma distribution.
%   randbeta         - Sample from Beta distribution.
%   randwishart      - Sample from Wishart distribution.
%   int_hist         - Histogram of integer values.
%
% Utility
%   logsumexp        - Sum in the log domain.
%   logmulexp        - Matrix multiply in the log domain.
%   ndgridmat        - Matrix of grid points.
%   ind2subv         - Subscript vector from linear index.
%   subv2ind         - Linear index from subscript vector.
%   gammaln          - Fast replacement for matlab's gammaln.
%   digamma          - Derivative of gammaln.
%   trigamma         - Derivative of digamma.
%   ndsum            - Sum over multiple dimensions.
%   ndmax            - Maximum over multiple dimensions.
%   ndlogsumexp      - Sum over multiple dimensions in the log domain.
%   maxdiff          - Maximum difference between structs or arrays.
%   makestruct       - Cell-friendly alternative to STRUCT.
%   setfields        - Set multiple fields of a structure.
%   struct2arglist   - Convert structure to cell array of fields/values.
%   sameobject       - Test if two variables correspond to the same object.
%   find_sameobject  - Find an object in a cell array.
%   toJava           - Convert to Java representation.
%   fromJava         - Convert from Java to Matlab.
%   glob             - Filename expansion via wildcards.
%   globstrings      - String matching via wildcards.
%
% Mutation
%   mutable          - Convert to a mutable object.
%   immutable        - Convert to an ordinary (immutable) object.
%
% Set operations
%   ismember_sorted  - True for member of sorted set.
%   match            - Location of matches in a set.
%   match_sorted     - Location of matches in a sorted set.
%   setdiff_sorted   - Set difference between sorted sets.
%   intersect_sorted - Set intersection between sorted sets.
%   union_sorted     - Set union of sorted sets.
%   union_sorted_rows - Set union of sorted sets of row vectors.
%   duplicated       - Find duplicated rows in a matrix.
%
% Readability
%   rows             - Number of rows.
%   cols             - Number of columns.
%   col_sum          - Sum for each column.
%   setdiag          - Modify the diagonals of a matrix.
%   finddiag         - Index of elements on diagonals.
%   argmax           - Index of maximum element.
%   argmin           - Index of minimum element.
%
% Flop counting
%   flops            - Read/write flop counter.
%   addflops         - Add to flop counter.
%   flops_chol       - Flops for Cholesky decomposition.
%   flops_col_sum    - Flops for column sums.
%   flops_det        - Flops for matrix determinant.
%   flops_digamma    - Flops for gammaln, digamma, and trigamma.
%   flops_div        - Flops for division.
%   flops_exp        - Flops for exponential.
%   flops_inv        - Flops for matrix inversion.
%   flops_mul        - Flops for real matrix multiplication.
%   flops_normpdfln  - Flops for normpdfln.
%   flops_pow        - Flops for raising to real power.
%   flops_randnorm   - Flops for randnorm.
%   flops_row_sum    - Flops for row sums.
%   flops_sample     - Flops for sample(p,n).
%   flops_solve      - Flops for matrix left division.
%   flops_solve_tri  - Flops for triangular left division.
%   flops_spadd      - Flops for sparse matrix addition.
%   flops_spmul      - Flops for sparse matrix multiplication.
%   flops_sqrt       - Flops for square root.
%
% Stand alone programs
%   matfile          - Read/write MAT files.
%   test_flops       - Compare time versus flops for various math operations.
%
% Graphics utilities
%  see graphics/Contents.m
%
% Demos
%   test_repmat,
%   test_solve_tri, ...

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%kmlud 'Compute LU Decomposition of Matrix (Input = L*U)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros mlud.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input matrix'
% OutputFile: o1 'Lower Triangular Matrix', optional: 'Resulting lower triangular output matrix'
% OutputFile: o2 'Upper Triangular Matrix', optional: 'Resulting upper triangular output matrix'
% OutputFile: o3 'Pivot Index Vector', optional: 'Pivot index vector'
%
% Example: [o1, o2, o3] = kmlud(i, {'i','';'o1','';'o2','';'o3',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% mlud - Compute LU Decomposition of Matrix  (Input = L*U)
%
%  DESCRIPTION
% .I mlud
% computes the LU decomposition of a matrix with partial pivoting and
% row interchanges.  If A is the input matrix, "mlud\fR will solve
% the system P*L*U = A for the L and U matrices. The P matrix is
% represented on output by a pivot index vector PI, where row i has been
% interchanged with row PI[i].  U is an upper
% triangular matrix, and L lower triangular with unit diagonal. PI is of length
% min(rows,cols).  The LAPACK routines DGETRF and ZGETRF are used to perform
% the actual decomposition.
%
% After execution of the appropriate LAPACK routine, a diagnostic may be
% printed out indicating that the input matrix is singular and at what point
% this was detected. This will be printed only if KHOROS_NOTIFY is set to
% STANDARD or more verbosity. This message is NOT an error, just information,
% and "mlud\fR will run to completion anyway.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
% LAPACK Users' Guide, E. Anderson et. al., SIAM 1992. ISBN 0-89871-294-7
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kmlud(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kmlud(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o1', '__output';'o2', '__output';'o3', '__output'};
maxval={0,1,1,1};
minval={0,1,1,1};
istoggle=[0,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','OutputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'mlud"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function symb_pvec = sdisplay(pvec,precision)
%SDISPLAY Symbolic display of SDPVAR expression
%
% As variables in YALMIP do not have any actual names internally, the
% command is not guaranteed to recover the variable names you use in the
% work-space correctly. It tries to figure out the names of the variables
% by searching in the work-space for variables involving the same variable
% indicies, but this can fails in many instances.
%
% EXAMPLES
%  sdpvar x y
%  sdisplay(x^2+y^2)
%    ans =
%       'x^2+y^2'
%
% Optionally, we can supply an option to round the data to N decimal digits
%
%  sdisplay(pi*x1^2,4)
%    ans = '3.1416*x^2'

if nargin < 2
    precision = inf;
end

% Support displaying non-sdpvar input
if ~isa(pvec,'sdpvar')
    for r1=1:size(pvec,1)
        for r2=1:size(pvec,2)
            p = full(pvec(r1,r2));
            if isa(p,'double')
                symb_pvec{r1,r2} = num2str2(p,precision);
            else
                symb_pvec{r1,r2} = p;
            end
        end
    end
    if prod(size(symb_pvec))==1 & nargout==0
        display(symb_pvec{1,1});
        clear symb_pvec
    end
    return;
end

% First, some boooring stuff. we need to
% figure out the symbolic names and connect
% these names to YALMIPs variable indicies
W = evalin('caller','whos');

% First, sort the variables available in the work-space based
% on the creation. Early creation means it is more likely the
% relevant variable to display
createTime = [];
for i = 1:size(W,1)
    if strcmp(W(i).class,'sdpvar') || strcmp(W(i).class,'ncvar')
        keep(i) = 1;
        z = evalin('caller',['struct(' W(i).name ').extra.createTime;']);
        createTime = [createTime z];
    else
        keep(i) = 0;
    end
end
W = W(find(keep));
[sorted,index] = sort(createTime);
W = W(index);

global_LinearVariables = depends(pvec);
global_x = recover(global_LinearVariables);
[global_exponent_p,global_ordered_list] = exponents(pvec,global_x);
global_exponent_p = full(global_exponent_p);
global_names = cell(length(global_x),1);

for i = 1:size(W,1)
    if strcmp(W(i).class,'sdpvar') || strcmp(W(i).class,'ncvar')
        % Get the SDPVAR variable
        thevars = evalin('caller',W(i).name);

        % Distinguish 4 cases
        % 1: Sclalar varible x
        % 2: Vector variable x(i)
        % 3: Matrix variable x(i,j)
        % 4: Variable not really defined
        if is(thevars,'scalar') && is(thevars,'linear') && length(getvariables(thevars))==1 & isequal(getbase(thevars),[0 1])
            index_in_p = find(ismember(global_LinearVariables,getvariables(thevars)));

            if ~isempty(index_in_p)
                already = ~isempty(global_names{index_in_p});
                if already
                    already = (isempty(strfind(global_names{index_in_p},'internal')) | isempty(strfind(global_names{index_in_p},'ans')));
                end
            else
                already = 0;
            end

            if ~isempty(index_in_p) & ~already
                % Case 1
                global_names{index_in_p}=W(i).name;
            end
        elseif is(thevars,'lpcone')

            if size(thevars,1)==size(thevars,2)
                % Case 2
                vars = getvariables(thevars);
                indicies = find(ismember(vars,global_LinearVariables));
                for ii = indicies
                    index_in_p = find(ismember(global_LinearVariables,vars(ii)));

                    if ~isempty(index_in_p)
                        already = ~isempty(global_names{index_in_p});
                        if already
                            already = (isempty(strfind(global_names{index_in_p},'internal')) | isempty(strfind(global_names{index_in_p},'ans')));
                        end
                    else
                        already = 0;
                    end

                    if ~isempty(index_in_p) & ~already
                        B = reshape(getbasematrix(thevars,vars(ii)),size(thevars,1),size(thevars,2));
                        [ix,jx,kx] = find(B);
                        ix=ix(1);
                        jx=jx(1);
                        global_names{index_in_p}=[W(i).name '(' num2str(ix) ',' num2str(jx) ')'];
                    end
                end

            else
                % Case 3
                vars = getvariables(thevars);
                indicies = find(ismember(vars,global_LinearVariables));
                for ii = indicies
                    index_in_p = find(ismember(global_LinearVariables,vars(ii)));

                    if ~isempty(index_in_p)
                        already = ~isempty(global_names{index_in_p});
                        if already
                            already = (isempty(strfind(global_names{index_in_p},'internal')) | isempty(strfind(global_names{index_in_p},'ans')));
                        end
                    else
                        already = 0;
                    end

                    if ~isempty(index_in_p) & ~already
                        global_names{index_in_p}=[W(i).name '(' num2str(ii) ')'];
                    end
                end
            end

        elseif is(thevars,'sdpcone')
            % Case 3
            vars = getvariables(thevars);
            indicies = find(ismember(vars,global_LinearVariables));
            for ii = indicies
                index_in_p = find(ismember(global_LinearVariables,vars(ii)));
                if ~isempty(index_in_p)
                    already = ~isempty(global_names{index_in_p});
                    if already
                        already = ~strfind(global_names{index_in_p},'internal');
                    end
                else
                    already = 0;
                end

                if ~isempty(index_in_p) & ~already
                    B = reshape(getbasematrix(thevars,vars(ii)),size(thevars,1),size(thevars,2));
                    [ix,jx,kx] = find(B);
                    ix=ix(1);
                    jx=jx(1);
                    global_names{index_in_p}=[W(i).name '(' num2str(ix) ',' num2str(jx) ')'];
                end
            end

        else
%             % Case 4
%             vars = getvariables(thevars);
%             indicies = find(ismember(vars,global_LinearVariables));
%
%             for i = indicies
%                 index_in_p = find(ismember(global_LinearVariables,vars(i)));
%                 if ~isempty(index_in_p) & isempty(global_names{index_in_p})
%                     global_names{index_in_p}=['internal(' num2str(vars(i)) ')'];
%                 end
%             end

        end
    end
end

for i = 1:length(global_names)
    if length(global_names{i}) == 0
        global_names{i} = ['internal(' num2str(global_LinearVariables(i)) ')'];
    end
end

for pi = 1:size(pvec,1)
    for pj = 1:size(pvec,2)
        p = pvec(pi,pj);
        if isa(p,'double')
            symb_p = num2str2(p,precision);
        else
            symb_p = createSymbolicExpression(p,global_LinearVariables,global_names,precision);
        end
        symb_pvec{pi,pj} = symb_p;
    end
end

if prod(size(symb_pvec))==1 & nargout==0
    display(symb_pvec{1,1});
    clear symb_pvec
end

function symb_p =  createSymbolicExpression(p,global_LinearVariables,global_names,precision)
LinearVariables = depends(p);
x = recover(LinearVariables);
[exponent_p,ordered_list] = exponents(p,x);
exponent_p = full(exponent_p);
[~,map] = ismember(LinearVariables,global_LinearVariables);
names = {global_names{map}};

% Remove 0 constant
symb_p = '';
if size(ordered_list,1)>0
    nummonoms = size(ordered_list,1);
    if full(getbasematrix(p,0)) ~= 0
        symb_p = num2str2(full(getbasematrix(p,0)),precision);
    end
elseif all(exponent_p(1,:)==0)
    symb_p = num2str2(full(getbasematrix(p,0)),precision);
    exponent_p = exponent_p(2:end,:);
    nummonoms = size(exponent_p,1);
else
    nummonoms = size(exponent_p,1);
end

% Loop through all monomial terms
for i = 1:nummonoms
    coeff = full(getbasematrixwithoutcheck(p,i));
    switch coeff
        case 1
            coeff='+';
        case -1
            coeff = '-';
        otherwise
            if isreal(coeff)
                if coeff >0
                    coeff = ['+' num2str2(coeff,precision)];
                else
                    coeff=[num2str2(coeff,precision)];
                end
            else
                coeff = ['+' '(' num2str2(coeff,precision) ')' ];
            end
    end
    if isempty(ordered_list)
        symb_p = [symb_p coeff symbmonom(names,exponent_p(i,:))];
    else
        symb_p = [symb_p coeff symbmonom_noncommuting(names,ordered_list(i,:))];
    end
end
% Clean up some left overs, lazy coding...
symb_p = strrep(symb_p,'+*','+');
symb_p = strrep(symb_p,'-*','-');
if symb_p(1)=='+'
    symb_p = symb_p(2:end);
end
if symb_p(1)=='*'
    symb_p = symb_p(2:end);
end

function s = symbmonom(names,monom)
s = '';
for j = 1:length(monom)
    if abs( monom(j))>0
        if isempty(names{j})
            names{j} = ['internal(' num2str(j) ')'];
        end
        s = [s '*' names{j}];
        if monom(j)~=1
            s = [s '^' num2str(monom(j))];
        end
    end
end

function s = symbmonom_noncommuting(names,monom)
s = '';
j = 1;
while j <= length(monom)
    if abs( monom(j))>0
        if isempty(names{monom(j)})
            names{monom(j)} = ['internal(' num2str(j) ')'];
        end
        s = [s '*' names{monom(j)}];
        power = 1;
        k = j;
        while j<length(monom) & monom(j) == monom(j+1)
            power = power + 1;
            j = j + 1;
            %if j == (length(monom)-1)
            %    j = 5;
            %end
        end
        if power~=1
            s = [s '^' num2str(power)];
        end
    end
    j = j + 1;
end

function s = num2str2(x,precision)
if isinf(precision)
    if isreal(x)
        s = sprintf('%.12g',x);
    else
        s = sprintf('%.12g+%.12gi',real(x),imag(x));
    end
else
    if isreal(x)
        s = sprintf('%.12g',round(x*10^precision)/10^precision);
    else
        s = sprintf('%.12g+%.12gi',round(real(x)*10^precision)/10^precision,round(imag(x)*10^precision)/10^precision);
    end
end
s(s==10)=[];
s(s==32)=[];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Reference performance
accTestRef = 0.9866;
allowedError = 0.001;
method = 'HPOLD';

% Create the algorithm object
algorithmObj = HPOLD();

% Clear parameter struct
clear param;

% Parameter C (Cost)
param.C = 10;

param.k = 10;

% Run the algorithm
info = algorithmObj.fitpredict(train,test,param);

trainCM = confusionmat(info.predictedTrain,train.targets);
testCM = confusionmat(info.predictedTest,test.targets);

accTrain = CCR.calculateMetric(trainCM);
accTest  = CCR.calculateMetric(testCM);

% Report accuracy
fprintf('Performing test for %s\n', method);
fprintf('Accuracy Train %f, Accuracy Test %f\n',accTrain,accTest);

if abs(accTestRef-accTest)<allowedError
    fprintf('Test accuracy matches reference accuracy\n');
else
    warning('Test accuracy does NOT match reference accuracy');
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
% compute bidirectional phase offsets from poor line scanning timing
function BiDiPhase = BiDiPhaseOffsets(data)

[Ly, Lx, nplanes, NT] = size(data);

% lines scanned one direction
yr1 = 2:2:floor(Ly/2)*2;
% lines scanned in other direction
yr2 = 1:2:floor(Ly/2)*2;

% compute phase correlation between lines in x-direction
eps0 = single(1e-6);
Nmax = min(50, NT);
d1 = fft(data(yr1,:,:,1:Nmax),[],2);
d2 = conj(fft(data(yr2,:,:,1:Nmax),[],2));
d1 = d1./(abs(d1) + eps0);
d2 = d2./(abs(d2) + eps0);

cc = ifft(d1 .* d2,[],2);
cc = fftshift(cc, 2);
cc = mean(mean(mean(cc,1),3),4);

% max shift of +/-5 pixels
[cx, ix] = max(cc(floor(Lx/2)+1 + [-5:5]));
ix       = ix - (6);

BiDiPhase = -1 * ix;

author: cortex-lab


text:
function LU = getbounds(F)

K.f = 0;
K.l = 0;
L = [];
U = [];
LU =  yalmip('getbounds',1:yalmip('nvars'));
F = flatten(F);
for i = 1:length(F.clauses)
    if F.clauses{i}.type == 2
        X = F.clauses{i}.data;
        AB = getbase(X);
        K.l = prod(size(X));
        variables = getvariables(X);
        [lb,ub,cand_rows] = find_lp_bounds(AB,K);
        LU(variables,1) = max([lb LU(variables,1)]')';
        LU(variables,2) = min([ub LU(variables,2)]')';
    elseif F.clauses{i}.type == 3
        % FIX : Extract from equalities and binary constraints
    end
end

binary = yalmip('binvariables');
LU(binary,1) = 0;
LU(binary,2) = 1;

% Try to bound some nonlinear terms
% FIX: complete code
[mt,variable_type] = yalmip('monomtable');
quadratic = find(variable_type == 2);
if ~isempty(quadratic)
    M = mt(quadratic,:);
    for i = 1:size(M,1)
        [ii,jj] = find(M(i,:));
        if length(ii) == 1
            LU(quadratic(i),1) = min([0 LU(jj,1)^2]);
            LU(quadratic(i),2) = max([LU(jj,1)^2 LU(jj,2)^2]);
        end
    end
end


text:
function a = i4vec_sort_insert_a ( n, a )

%*****************************************************************************80
%
%% I4VEC_SORT_INSERT_A uses an ascending insertion sort on an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher and Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998, page 11.
%
%  Parameters:
%
%    Input, integer N, the number of items in the vector.
%    N must be positive.
%
%    Input, integer A(N), the array to be sorted.
%
%    Output, integer A(N), the sorted array.
%
  for i = 2 : n

    x = a(i);

    j = i - 1;

    while ( 1 <= j )

      if ( a(j) <= x )
        break
      end

      a(j+1) = a(j);
      j = j - 1;

    end

    a(j+1) = x;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
close all;
clear all;
clc;
rng('default');
png_export = true;
pdf_export = false;
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

mf = spx.graphics.Figures();

% Signal space
N = 256;
% Number of measurements
M = 64;
% Sparsity level
K = 8;
% Number of signals
S = 4;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
% Generate bi-uniform signals
X = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
% Measurement vectors
Y = Phi * X;

solver = spx.pursuit.joint.BasisPursuit(Phi);
result = solver.solve_l2_l1(Y);
% Solution vectors
Z = result.Z;
% Comparison
cs = spx.commons.SparseSignalsComparison(X, Z, K);
cs.summarize();
for s=1:S
    mf.new_figure(sprintf('MMV signal: %d', s));
    subplot(411);
    stem(X(:, s), '.');
    title('Sparse vector');
    subplot(412);
    stem(Z(:, s), '.');
    title('Recovered sparse vector');
    subplot(413);
    stem(abs(X(:, s) - Z(:, s)), '.');
    title('Recovery error');
    subplot(414);
    stem(Y(:, s), '.');
    title('Measurement vector');
end

author: indigits

repo: sparse-plex


text:
function [obj] = ft_convert_units(obj, target, varargin)

% FT_CONVERT_UNITS changes the geometrical dimension to the specified SI unit.
% The units of the input object is determined from the structure field
% object.unit, or is estimated based on the spatial extend of the structure,
% e.g. a volume conduction model of the head should be approximately 20 cm large.
%
% Use as
%   [object] = ft_convert_units(object, target)
%
% The following geometrical objects are supported as inputs
%   electrode or gradiometer array, see FT_DATATYPE_SENS
%   volume conductor, see FT_DATATYPE_HEADMODEL
%   anatomical mri, see FT_DATATYPE_VOLUME
%   segmented mri, see FT_DATATYPE_SEGMENTATION
%   dipole grid definition, see FT_DATATYPE_SOURCE
%
% Possible target units are 'm', 'dm', 'cm ' or 'mm'. If no target units
% are specified, this function will only determine the native geometrical
% units of the object.
%
% See also FT_DETERMINE_UNITS, FT_CONVERT_COORDSYS, FT_DETERMINE_COODSYS

% Copyright (C) 2005-2016, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% This function consists of three parts:
%   1) determine the input units
%   2) determine the requested scaling factor to obtain the output units
%   3) try to apply the scaling to the known geometrical elements in the input object

% the "target" input argument has been made required in Aug 2017
% prior to that it was also possible to use this function to estimate units
% the backward compatibility support can be removed in Aug 2018
if nargin<2
  ft_warning('calling this function only to determine units is deprecated, please use FT_DETERMINE_UNITS instead');
  obj = ft_determine_units(obj);
  return
end

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

% parse input options
feedback = ft_getopt(varargin, 'feedback', false);

if isstruct(obj) && numel(obj)>1
  % deal with a structure array
  for i=1:numel(obj)
    tmp(i) = ft_convert_units(obj(i), target, varargin{:});
  end
  obj = tmp;
  return
elseif iscell(obj) && numel(obj)>1
  % deal with a cell-array
  % this might represent combined EEG, ECoG and/or MEG
  for i=1:numel(obj)
    obj{i} = ft_convert_units(obj{i}, target, varargin{:});
  end
  return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine the units of the input object
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
obj = ft_determine_units(obj);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the scaling factor from the input units to the desired ones
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isequal(obj.unit, target)
  % there is nothing to do
  return
end

if istrue(feedback)
  % give some information about the conversion
  fprintf('converting units from ''%s'' to ''%s''\n', obj.unit, target)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apply the scaling factor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
scale = ft_scalingfactor(obj.unit, target);

% volume conductor model
if isfield(obj, 'r'), obj.r = scale * obj.r; end
if isfield(obj, 'o'), obj.o = scale * obj.o; end
if isfield(obj, 'bnd') && isfield(obj.bnd, 'pnt')
  for i=1:length(obj.bnd)
    obj.bnd(i).pnt = scale * obj.bnd(i).pnt;
  end
end
if isfield(obj, 'bnd') && isfield(obj.bnd, 'pos')
  for i=1:length(obj.bnd)
    obj.bnd(i).pos = scale * obj.bnd(i).pos;
  end
end

% old-fashioned gradiometer array
if isfield(obj, 'pnt1'), obj.pnt1 = scale * obj.pnt1; end
if isfield(obj, 'pnt2'), obj.pnt2 = scale * obj.pnt2; end
if isfield(obj, 'prj'),  obj.prj  = scale * obj.prj;  end

% gradiometer array, electrode array, head shape or dipole grid
if isfield(obj, 'pnt'),        obj.pnt        = scale * obj.pnt;        end
if isfield(obj, 'pos'),        obj.pos        = scale * obj.pos;        end
if isfield(obj, 'chanpos'),    obj.chanpos    = scale * obj.chanpos;    end
if isfield(obj, 'chanposorg'), obj.chanposold = scale * obj.chanposorg; end % pre-2016 version
if isfield(obj, 'chanposold'), obj.chanposold = scale * obj.chanposold; end % 2016 version and later
if isfield(obj, 'coilpos'),    obj.coilpos    = scale * obj.coilpos;    end
if isfield(obj, 'elecpos'),    obj.elecpos    = scale * obj.elecpos;    end

% gradiometer array that combines multiple coils in one channel
if isfield(obj, 'tra') && isfield(obj, 'chanunit')
  % find the gradiometer channels that are expressed as unit of field strength divided by unit of distance, e.g. T/cm
  for i=1:length(obj.chanunit)
    tok = tokenize(obj.chanunit{i}, '/');
    if ~isempty(regexp(obj.chanunit{i}, 'm$', 'once'))
      % assume that it is T/m or so
      obj.tra(i,:)    = obj.tra(i,:) / scale;
      obj.chanunit{i} = [tok{1} '/' target];
    elseif ~isempty(regexp(obj.chanunit{i}, '[T|V]$', 'once'))
      % assume that it is T or V, don't do anything
    elseif strcmp(obj.chanunit{i}, 'unknown')
      % assume that it is T or V, don't do anything
    elseif strcmp(obj.chanunit{i}, 'snr')
      %
    else
      ft_error('unexpected units %s', obj.chanunit{i});
    end
  end % for
end % if

% fiducials
if isfield(obj, 'fid') && isfield(obj.fid, 'pnt'), obj.fid.pnt = scale * obj.fid.pnt; end
if isfield(obj, 'fid') && isfield(obj.fid, 'pos'), obj.fid.pos = scale * obj.fid.pos; end

% dipole grid
if isfield(obj, 'resolution'), obj.resolution = scale * obj.resolution; end

% x,y,zgrid can also be 'auto'
if isfield(obj, 'xgrid') && ~ischar(obj.xgrid), obj.xgrid = scale * obj.xgrid; end
if isfield(obj, 'ygrid') && ~ischar(obj.ygrid), obj.ygrid = scale * obj.ygrid; end
if isfield(obj, 'zgrid') && ~ischar(obj.zgrid), obj.zgrid = scale * obj.zgrid; end

% anatomical MRI or functional volume
if isfield(obj, 'transform')
  H = diag([scale scale scale 1]);
  obj.transform = H * obj.transform;
end

if isfield(obj, 'transformorig')
  H = diag([scale scale scale 1]);
  obj.transformorig = H * obj.transformorig;
end

% remove initial and params structure if they exist
if isfield(obj, 'initial') && ~strcmp(target, 'mm')
  obj = rmfield(obj, 'initial');
  ft_warning('Removing field "initial" because potential transformations of normalised volumes only work if geometrical values are expressed in "mm"');
end

if isfield(obj, 'params') && ~strcmp(target, 'mm')
  ft_warning('Removing field "params" because potential transformations of normalised volumes only work if geometrical values are expressed in "mm"');
  obj = rmfield(obj, 'params');
end

% sourcemodel obtained through mne also has a orig-field with the high
% number of vertices
if isfield(obj, 'orig')
  if isfield(obj.orig, 'pnt')
    obj.orig.pnt = scale * obj.orig.pnt;
  end
  if isfield(obj.orig, 'pos')
    obj.orig.pos = scale * obj.orig.pos;
  end
end

% remember the unit
obj.unit = target;

author: spm


text:
% QQDIAGRAM - Empirical quantile-quantile diagram.
%
% Description:
%               The quantiles (percentiles) of the input distribution Y are plotted (Y-axis)
%               against the corresponding quantiles of the input distribution X.
%               If only X is given, the corresponding quantiles are plotted (Y-axis)
%               against the quantiles of a Gaussian distribution ('Normal plot').
%               Two black dots indicate the lower and upper quartiles.
%               If the data in X and Y belong the same distribution the plot will be linear.
%               In this case,the red and black reference lines (.-.-.-.-) will overlap.
%               This will be true also if the data in X and Y belong to two distributions with
%               the same shape, one distribution being rescaled and shifted with respect to the
%               other.
%               If only X is given, a line is plotted to indicate the mean of X, and a segment
%               is plotted to indicate the standard deviation of X. If the data in X are normally
%               distributed, the red and black reference lines (.-.-.-.-) will overlap.
%
% Usage:
%   >>  ah  =  qqdiagram( x, y, pk );
%
% Inputs:
%   x       - vector of observations
%
% Optional inputs:
%   y       - second vector of observation to compare the first to
%   pk      - the empirical quantiles will be estimated at the values in pk [0..1]
%
% Author: Luca Finelli, CNL / Salk Institute - SCCN, 20 August 2002
%
% Reference: Stahel W., Statistische Datenanalyse, Vieweg, Braunschweig/Wiesbaden, 1995
%
% See also:
%   QUANTILE, SIGNALSTAT, EEGLAB

% Copyright (C) 2002 Luca Finelli, Salk/SCCN, La Jolla, CA
%
% Reference: Stahel, W. Statistische Datenanalyse, Vieweg, Braunschweig/Wiesbaden 1995
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function qqdiagram( x , y, pk )

if nargin < 1
	help qqdiagram;
	return;
end;

if (nargin == 3 && (any(pk > 1) || any(pk < 0)))
    error('qqdiagram(): elements in pk must be between 0 and 1');
end

if nargin==1
	y=x;
	nn=max(1000,10*length(y))+1;
	x=randn(1,nn);
end

if nargin < 3
	nx=sum(~isnan(x));
	ny=sum(~isnan(y));
   	k=min(nx,ny);
    pk=((1:k) - 0.5) ./ k;  % values to estimate the empirical quantiles at
else
    k=length(pk);
end

if nx==k
    xx=sort(x(~isnan(x)));
else
    xx=quantile(x(~isnan(x)),pk);
end

if ny==k
    yy=sort(y(~isnan(y)));
else
    yy=quantile(y(~isnan(y)),pk);
end

% QQ diagram
plot(xx,yy,'+')
hold on

% x-axis range
maxx=max(xx);
minx=min(xx);
rangex=maxx-minx;
xmin=minx-rangex/50;
xmax=maxx+rangex/50;

% Quartiles
xqrt1=quantile(x,0.25); xqrt3=quantile(x,0.75);
yqrt1=quantile(y,0.25); yqrt3=quantile(y,0.75);

plot([xqrt1 xqrt3],[yqrt1 yqrt3],'k-','LineWidth',2); % IQR range

% Drawing the line
sigma=(yqrt3-yqrt1)/(xqrt3-xqrt1);
cy=(yqrt1 + yqrt3)/2;

if nargin ==1
    maxy=max(y);
    miny=min(y);
    rangey=maxy-miny;
	ymin=miny-rangey/50;
	ymax=maxy+rangey/50;

	plot([(miny-cy)/sigma (maxy-cy)/sigma],[miny maxy],'r-.') % the line
    % For normally distributed data, the slope of the plot line
    % is equal to the ratio of the standard deviation of the distributions
	plot([0 (maxy-mean(y))/std(y)],[mean(y) maxy],'k-.') % the ideal line

	xlim=get(gca,'XLim');
	plot([1 1],[ymin  mean(y)+std(y)],'k--')
	plot([1 1],[mean(y)  mean(y)+std(y)],'k-','LineWidth',2)
        % textx = 1.0;
        % texty = mean(y)+3.0*rangey/50.0;
	% text(double(textx), double(texty),' St. Dev.','horizontalalignment','center')
    set(gca,'xtick',get(gca,'xtick'));  % show that vertical line is at 1 sd
	plot([0 0],[ymin  mean(y)],'k--')
	plot(xlim,[mean(y) mean(y)],'k--')
	% text(double(xlim(1)), double(mean(y)+rangey/50),'Mean X')
	plot([xqrt1  xqrt3],[yqrt1 yqrt3],'k.','MarkerSize',10)
	set(gca,'XLim',[xmin xmax],'YLim',[ymin ymax])
	xlabel('Standard Normal Quantiles')
	ylabel('X Quantiles')
else
    cx=(xqrt1 + xqrt3)/2;
    maxy=cy+sigma*(max(x)-cx);
	miny=cy-sigma*(cx-min(x));

	plot([min(x) max(x)],[miny maxy],'r-.'); % the line
    xlabel('X Quantiles');
    ylabel('Y Quantiles');
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [xEst,PEst,logLikes]=singleScanUpdate(xHyp,PHyp,A,algSel1,algSel2,param3)
%%SINGLESCANUPDATE Perform the measurement update step in a single-scan
%                  tracking algorithm that uses Gaussian approximations
%                  to represent the target state before and after the
%                  measurement update given a matrix of likelihoods and a
%                  set of track update hypotheses. The algorithm can be a
%                  global nearest neighbor (GNN) update, a joint
%                  probabilistic data association (JPDA) update, a JPDA*
%                  update, a GNN-JPDA update, or approximate JPDA and
%                  GNN-JPDA updates as well as a naive nearest neighbor
%                  algorithm. In all instances except for the naive nearest
%                  neighbor algorithm, it is assumed that each measurement
%                  can be assigned to at most one target, and each target
%                  can be assigned to at most one measurement.
%
%INPUTS: xHyp An xDimXnumTarXnumHyp set of track states updated with each
%             of the (numHyp-1) measurements for each of the numTar
%             targets, with the last hypothesis being the update for the
%             missed detection hypothesis. The order of the measurement
%             hypotheses is the same as the ordering of the likelihoods in
%             the columns of the A matrix.
%        PHyp An xDimXxDimXnumTarXnumHyp set of covariance matrices for
%             each of the track states for each of the targets updated
%             conditioned on each of the numHyp measurements with the last
%             one being for the missed detection hypothesis.
%           A A matrix of positive likelihoods or likelihood ratios (NOT
%             log-likelihood ratios). A is a numTar X (numMeas+numTar)
%             matrix of all-positive likelihoods or likelihood ratios for
%             assigning the target specified by the row to the measurement
%             specified by the column. Columns > numMeas hold
%             missed-detection likelihoods. Thus, off-diagonal terms for
%             columns > numMeas should be set to 0 and the diagonal terms
%             set to the costs of a missed detection for each given
%             target.
%     algSel1 An optional parameter that selects the assignment algorithm
%             to use. Approximate variants of the algorithms also use
%             algSel2. If this parameter and the next parameter are both
%             omitted, then the algorithm is chosen as described below.
%             Possible values are
%             0) GNN-JPDA
%             1) JPDA
%             2) GNN
%             3) Parallel single-target PDAs
%             4) Naive nearest neighbor
%             5) JPDA*
%             6) Approximate GNN-JPDA
%             7) Approximate JPDA
%             8) Set JPDA
%             9) Naive nearest neighbor JPDA
%             10) Approximate naive nearest neighbor JPDA
%     algSel2 An optional parameter that further specifies the algorithm
%             used when algSel1=6-7. If omitted but algSel1 is specified, a
%             default value of 0 is used. If both algSel1 and algSel2
%             are omitted, the algorithm is chosen as described below.
%             algSel2 chooses the approximation for the beta terms in the
%             function calc2DAssignmentProbsApprox (the approxType input)
%             when using algSel1=6,7; the comments to that function
%             describe the values it can take. If algSel1=8 (Set JPDAF),
%             this is the number of position components in the state (the
%             posDim input in the calcSetJPDAUpdate function).
%      param3 For the case where algSel1=6,7 and algSel2=0, this is the
%             optional delta input to the function
%             calc2DAssignmentProbsApprox. If omitted or an empty matrix
%             is passed, the default value for that function is used.
%
%OUTPUTS: xEst An xDimXnumTar matrix of updated target states.
%         PEst An xDimXxDimXnumtar matrix of updated covariance matrices
%              for the targets.
%      logLike The numTarX1 set of logarithms of the likelihood (ratio)
%              function for the update for each target. If A contains
%              likelihoods, then this is the logarithm of the likelihood;
%              if A constains likelihood ratios, then this is a logarithm
%              of a likelihood ratio, For hard assignments, this is clearly
%              defined (it is just the entry in A assigned to that target).
%              For soft assignments in the JPDA, JPDA*, PDA, and
%              approximate JPDA this is the expected value of the
%              log-likelihood computed as using the beta terms. If the Set
%              JPDAF is chosen, an empty matrix is returned for logLike.
%
%If both algSel1 and algSel2 are omitted, then the algorithm chosen depends
%on the size of A. If A has only one row (one target), then the JPDA is
%used. Otherwise, the approximate GNN-JPDA with algSel2=4 is used (which is
%actually exact if the number of targets is two).
%
%The JPDA is discussed in Chapter 6.2 of [1], and the PDA is discussed in
%various other sections as the JPDA is a generalization of the PDA. The
%JPDA* is discussed in [2]. The concept of the GNN-JPDA is described (but
%not named) in [3]. The GNN-JPDA consists of using the global nearest
%neighbor estimate with a covariance matrix computed as in the JPDA. Thus
%the same approximations that can be used in the JPDA can be used in the
%GNN-JPDA. The Set JPDA is described in [4].
%
%The approximate routines simply change how the target-measurement
%association probabilities are computed. They are described in the comments
%to the function calc2DAssignmentProbsApprox.
%
%Naive nearest neighbor association consists of just choosing the most
%likely assignment of a target to a measurement or missed detection without
%regard to whether one measurmeent is assigned to multiple targets.
%
%One possible way to form the A matrix is to use likelihood ratios going
%into the dimensionless score function from [5]. This is implemented for a
%standard Cartesian/Gaussian model using makeStandardCartOnlyLRMatHyps.
%
%REFERENCES:
%[1] Y. Bar-Shalom, P. K. Willett, and X. Tian, Tracking and Data Fusion.
%    Storrs, CT: YBS Publishing, 2011.
%[2] H. A. Blom and E. A. Bloem, "Probabilistic data association avoiding
%    track coalescence," IEEE Transactions on Automatic Control, vol. 45,
%    no. 2, pp. 247-259, Feb. 2000.
%[3] O. E. Drummond, "Best hypothesis target tracking and sensor fusion,"
%    in Proceedings of SPIE: Signal and Data Processing of Small Targets
%    Conference, vol. 3809, Denver, CO, Oct. 1999, pp. 586-600.
%[4] L. Svensson, D. Svensson, M. Guerriero, and P. Willett, "Set JPDA
%    filter for multitarget tracking," IEEE Transactions on Signal
%    Processing, vol. 59, no. 10, pp. 4677-4691, Oct. 2011.
%[5] Y. Bar-Shalom, S. S. Blackman, and R. J. Fitzgerald, "Dimensionless
%    score function for multiple hypothesis tracking," IEEE Transactions on
%    Aerospace and Electronic Systems, vol. 43, no. 1, pp. 392-400, Jan.
%    2007.
%
%March 2015 David F. Crouse, generalizing the basic JPDAF code of David
%Karnick to many more algorithms, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    numTar=size(A,1);
    numMeas=size(A,2)-numTar;
    numHyp=numMeas+1;

    if(nargin<6)
        param3=[];
    end

    if(nargin<4||isempty(algSel1))
        if(numTar==1)
            algSel1=1;
        else
            algSel1=6;
            algSel2=4;
        end
    end

    if(nargin==4||isempty(algSel2))%If only algSel2 was omitted.
        algSel2=0;
    end
    xDim=size(xHyp,1);

    xEst=zeros(xDim,numTar);
    PEst=zeros(xDim,xDim,numTar);
    switch(algSel1)
        case 0%GNN-JPDA
            %Perform 2D assignment
            ALog=log(A);
            tar2Meas=assign2D(ALog,true);
            logLikes=zeros(numTar,1);
            for curTar=1:numTar
                logLikes(curTar)=ALog(curTar,tar2Meas(curTar));
            end

            %Adjust for the index of the missed detection hypothesis.
            sel=tar2Meas>numMeas;
            tar2Meas(sel)=numMeas+1;

            %Set the target estimate to the ML estimate.
            for curTar=1:numTar
                maxIdx=tar2Meas(curTar);
                xEst(:,curTar)=xHyp(:,curTar,maxIdx);
            end

            %Compute the JPDAF probabilities.
            beta=calc2DAssignmentProbs(A,true);

            %Compute the covariance matrix in the manner of the JPDAF, but
            %about the ML estimate, so the result is a MSE matrix estimate.
            for curTar=1:numTar
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [~,PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P,xEst(:,curTar));
            end
        case 1%JPDA
            %Compute the JPDAF probabilities.
            beta=calc2DAssignmentProbs(A,true);

            %The JPDAF is just the mixture mean and covariance matrix.
            for curTar=1:numTar
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [xEst(:,curTar),PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P);
            end

            %If a likelihood is desired, then find the weighted likelihood.
            if(nargout>2)
                logLikes=logLikefromABeta(A,beta);
            end
        case 2%GNN
            %Perform 2D assignment
            ALog=log(A);
            tar2Meas=assign2D(ALog,true);
            logLikes=zeros(numTar,1);
            for curTar=1:numTar
                logLikes(curTar)=ALog(curTar,tar2Meas(curTar));
            end
            %Adjust for the index of the missed detection hypothesis.
            sel=tar2Meas>numMeas;
            tar2Meas(sel)=numMeas+1;

            for curTar=1:numTar
                maxIdx=tar2Meas(curTar);
                xEst(:,curTar)=xHyp(:,curTar,maxIdx);
                PEst(:,:,curTar)=PHyp(:,:,curTar,maxIdx);
            end
        case 3%Paralle single-target PDAs
            %This is just a bunch of independent PDAFs for each target.
            beta=zeros(numTar,numHyp);
            for curTar=1:numTar
                hypIdx=[1:numMeas,numMeas+curTar];
                beta(curTar,:)=calc2DAssignmentProbs(A(curTar,hypIdx),true);
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [xEst(:,curTar),PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P);
            end

            %If a likelihood is desired, then find the weighted likelihood.
            if(nargout>2)
                logLikes=logLikefromABeta(A,beta);
            end
        case 4%Na�ve nearest neighbor
            %The assignment for each target
            logLikes=zeros(numTar,1);
            for curTar=1:numTar
                [maxVal,maxIdx]=max(A(curTar,:));
                %If the missed detection hypothesis is the most likely.
                if(maxIdx>numMeas)
                    maxIdx=numMeas+1;
                end

                xEst(:,curTar)=xHyp(:,curTar,maxIdx);
                PEst(:,:,curTar)=PHyp(:,:,curTar,maxIdx);

                logLikes(curTar)=log(maxVal);
            end
        case 5%JPDA*
            %This is the same as the JPDA, except the method for computing
            %the betas is different.
            %Compute the JPDAF* probabilities.
            beta=calcStarBetasBF(A);

            %The JPDAF* is just the mixture mean and covariance matrix.
            for curTar=1:numTar
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [xEst(:,curTar),PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P);
            end

            %If a likelihood is desired, then find the weighted likelihood.
            if(nargout>2)
                logLikes=logLikefromABeta(A,beta);
            end
        case 6%Approximate GNN-JPDA
            %This is the same as the GNN-JPDA, except the method for
            %computing the betas is different.

            %Perform 2D assignment
            ALog=log(A);
            tar2Meas=assign2D(ALog,true);
            logLikes=zeros(numTar,1);
            for curTar=1:numTar
                logLikes(curTar)=ALog(curTar,tar2Meas(curTar));
            end

            %Adjust for the index of the missed detection hypothesis.
            sel=tar2Meas>numMeas;
            tar2Meas(sel)=numMeas+1;

            %Set the target estimate to the ML estimate
            for curTar=1:numTar
                maxIdx=tar2Meas(curTar);
                xEst(:,curTar)=xHyp(:,curTar,maxIdx);
            end

            %Compute the approximate JPDAF probabilities.
            beta=calc2DAssignmentProbsApprox(A,algSel2,true,param3);

            %Compute the covariance matrix in the manner of the JPDAF, but
            %about the ML estimate, so the result is a MSE matrix estimate.
            for curTar=1:numTar
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [~,PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P,xEst(:,curTar));
            end
        case 7%Approximate JPDA
            %This is the same as the JPDA, except the method for computing
            %the betas is different.

            %Compute the approximate JPDAF probabilities.
            beta=calc2DAssignmentProbsApprox(A,algSel2,true,param3);

            %The JPDAF* is just the mixture mean and covariance matrix.
            for curTar=1:numTar
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [xEst(:,curTar),PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P);
            end

            %If a likelihood is desired, then find the weighted likelihood.
            if(nargout>2)
                logLikes=logLikefromABeta(A,beta);
            end
        case 8%Set JPDAF
            %This is the same as the JPDA, except the method for computing
            %the betas is different.
            [xEst,PEst]=calcSetJPDAUpdate(xHyp,PHyp,A,AlgSel2);

            if(nargout>2)
                logLikes=[];
            end
        case 9%Na�ve nearest neighbor-JPDA
            %The assignment for each target
            logLikes=zeros(numTar,1);
            for curTar=1:numTar
                [maxVal,maxIdx]=max(A(curTar,:));

                %If the missed detection hypothesis is the most likely.
                if(maxIdx>numMeas)
                    maxIdx=numMeas+1;
                end

                xEst(:,curTar)=xHyp(:,curTar,maxIdx);
                logLikes(curTar)=log(maxVal);
            end

            %Compute the JPDAF probabilities.
            beta=calc2DAssignmentProbs(A,true);

            %Compute the covariance matrix in the manner of the JPDAF, but
            %about the ML estimate, so the result is a MSE matrix estimate.
            for curTar=1:numTar
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [~,PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P,xEst(:,curTar));
            end
        case 10%Approxmate na�ve nearest neighbor-JPDA
            %This is the same as the na�ve nearest neighbor-JPDA, except
            %the method for computing the betas is different.

            %The assignment for each target
            logLikes=zeros(numTar,1);
            for curTar=1:numTar
                [maxVal,maxIdx]=max(A(curTar,:));

                %If the missed detection hypothesis is the most likely.
                if(maxIdx>numMeas)
                    maxIdx=numMeas+1;
                end

                xEst(:,curTar)=xHyp(:,curTar,maxIdx);
                logLikes(curTar)=log(maxVal);
            end

            %Compute the approximate JPDAF probabilities.
            beta=calc2DAssignmentProbsApprox(A,algSel2,true,param3);

            %Compute the covariance matrix in the manner of the JPDAF, but
            %about the ML estimate, so the result is a MSE matrix estimate.
            for curTar=1:numTar
                x=reshape(xHyp(:,curTar,:),[xDim,numHyp,1]);
                P=reshape(PHyp(:,:,curTar,:),[xDim,xDim,numHyp]);
                [~,PEst(:,:,curTar)]=calcMixtureMoments(x,beta(curTar,:),P,xEst(:,curTar));
            end
        otherwise
            error('Unknown algorithm selected.')
    end
end

function logLike=logLikefromABeta(A,beta)
    numTar=size(A,1);
    numHyp=size(beta,2);

    %Create a copy of A with the missed detection likelihoods
    %in the last column and  take the logarithm of it.
    ANew=zeros(numTar,numHyp);
    ANew(:,1:(numHyp-1))=A(:,1:(numHyp-1));
    ANew(:,numHyp)=diag(A(:,numHyp:(numHyp+numTar-1)));
    ANew=log(ANew);

    %Multiply each log-likelihood by the probability of it
    %being correct.
    ANew=ANew.*beta;
    %Get rid of non-finite terms, which should generally only
    %arise in things with zero likelihood.
    sel=~isfinite(ANew);
    ANew(sel)=0;
    %The likelihood is the sum of all of the weighted
    %log-likelihoods.
    logLike=sum(ANew,2);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function isContained=rectContainedInRect(rectMin1,rectMax1,rectMin2,rectMax2)
%RECTCONTAINEDINRECT Determine whether an axis-aligned rectangle (or a
%                    more general axis-aligned hyperrectangle if the number
%                    of dimensions is not two) is completely engulfed by
%                    another rectangle. This does not indicate which
%                    rectangle is contained in the other.
%
%INPUTS: rectMin1 A kX1 or 1Xk vector of the lower bounds of each of the
%                 dimensions of the first k-dimensional hyperrectangle.
%        rectMax1 A kX1 or 1Xk vector of the upper bounds of the first
%                 hyperectangle.
%        rectMin2 A kX1 or 1Xk vector of the lower bounds of the second
%                 hyperrectangle.
%        rectMax2 A kX1 or 1Xk vector of the upper bounds of the second
%                 hyperrectangle.
%
%OUTPUTS: isContained A boolean value that is true if one hyperrectangle
%             is completely contained in the other and is false otherwise.
%
%December 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

k=length(rectMin1);

isContained=true;
for curIdx=1:k
    if(rectMin1(curIdx)<rectMin2(curIdx))
        isContained=false;
        break;
    end

    if(rectMax1(curIdx)>rectMax2(curIdx))
       isContained=false;
       break;
    end
end

if(isContained==true)
    return;
end

isContained=true;
for curIdx=1:k
    if(rectMin2(curIdx)<rectMin1(curIdx))
        isContained=false;
        return;
    end

    if(rectMax2(curIdx)>rectMax1(curIdx))
       isContained=false;
       return;
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
%% A10 WARTHOG FIXED WING AIRCRAFT (A10.m) ///////////////////////////////

% Author: James A. Douthwaite

classdef A10 < fixedWing
    properties
        % Any specific properties of globalhawk
    end
    %  CLASS METHODS
    methods
        % CONSTRUCTOR METHOD
        function obj = A10(varargin)
            % CALL THE SUPERCLASS CONSTRUCTOR
            obj@fixedWing(varargin);  % Create the super class 'fixedWing'

            % //////////////// Check for user overrides ///////////////////
            [obj] = obj.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
    end
end

repo: OpenMAS


text:
function jtrMask3M = jitterMask(mask3M,pctJitter)
% function jtrMask3M = jitterMask(mask3M,pctJitter)
%
% Jitters the input mask3M by translation, rotation and scaling by a factor
% of pctJitter.
%
% Example:
% mask3M = rand(10,10,5)  0.5;
% pctJitter = 5; % percent
% jtrMask3M = jitterMask(mask3M,pctJitter);
%
% APA, 12/2/2018

% global planC
% indexS = planC{end};
%
% structNum = 1;
% scanNum = 1;
% pctJitter = 5;
%
% fullmaks3M = getUniformStr(structNum,planC);
% [rasterSegments, planC, isError] = getRasterSegments(structNum,planC);
% [mask3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);

jtrMask3M = zeros(size(mask3M),'like',mask3M);
for slc = 1:size(mask3M,3)
    fullSiz = size(mask3M(:,:,slc));
    [iV,jV] = find3d(mask3M(:,:,slc));
    maskM = mask3M(min(iV):max(iV),min(jV):max(jV),slc);
    di = max(iV) - min(iV);
    dj = max(jV) - min(jV);
    iCtr = ceil((max(iV) + min(iV))/2);
    jCtr = ceil((max(jV) + min(jV))/2);
    minRow = -di*pctJitter/100;
    maxRow = di*pctJitter/100;
    minCol = -dj*pctJitter/100;
    maxCol = dj*pctJitter/100;
    numRows = minRow + (maxRow-minRow)*rand(1);
    numCols = minCol + (maxCol-minCol)*rand(1);
    minAng = -pctJitter/100*180;
    maxAng = pctJitter/100*180;
    angl = minAng + (maxAng-minAng)*rand(1);
    scl = (100-pctJitter)/100 + 2*pctJitter/100*rand(1);
    maskM = imtranslate(maskM,[numRows,numCols],'nearest','FillValues',0);
    maskM = imrotate(maskM,angl,'nearest');
    maskM = imresize(maskM, scl, 'nearest');
    newSiz = size(maskM);
    iStart = 1;
    jStart = 1;
    iEnd = 0;
    jEnd = 0;
    iMin = iCtr - ceil(newSiz(1)/2);
    if iMin < 0
        iStart = 1-iMin;
        iMin = 1;
    end
    iMax = iCtr + floor(newSiz(1)/2) - 1;
    if iMax > fullSiz(1)
        iEnd = fullSiz(1) - iMax;
        iMax = fullSiz(1);
    end
    jMin = jCtr - ceil(newSiz(2)/2);
    if jMin < 0
        jStart = 1-jMin;
        jMin = 1;
    end
    jMax = jCtr + floor(newSiz(2)/2) - 1;
    if jMax > fullSiz(2)
        jEnd = fullSiz(2) - jMax;
        jMax = fullSiz(2);
    end
    %mask3M(:,:,slc) = 0;
    jtrMask3M(iMin:iMax,jMin:jMax,slc) = maskM(iStart:end-iEnd,jStart:end-jEnd);
end

% fullmaks3M(:,:,uniqueSlices) = mask3M;
% isUniform = 1;
% strname = [planC{indexS.structures}(structNum).structureName, '_Perturbed'];
% planC = maskToCERRStructure(fullmaks3M, isUniform, scanNum, strname, planC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%Method of Rosenbriock using line searches

function [xmin,fmin]=Risenbriock_Method(x0)

epsilon=1e-5;
x_k=x0;
k=0;
n=length(x0);
D=eye(n);
Lambda=zeros(n,1);

while 1
    % orthogonal search
    z_j=x_k;
    for j=1:n
        d_j=D(:,j);
        lambda_j=fminbnd(@(lambda) theta(lambda,z_j,d_j),-10,10);
        Lambda(j)=lambda_j;
        z_j=z_j+lambda_j*d_j;
    end
    x_kplus1=z_j;

    % stop criteria
    if (norm(x_kplus1-x_k)<epsilon) && (objective_fun(x_kplus1)-objective_fun(x_k)<epsilon)
        xmin=x_kplus1;
        fmin=objective_fun(xmin);
        k
        return;
    end

    % Gram-Schmidt procedure
    D=Gram_Schmidt_Procedure(D,Lambda);

    % iteration
    x_k=x_kplus1;
    k=k+1;
end

repo: Optimization-Theory-and-Algorithm


text:
load 'sdae_rbm_mnist_vis.mat';

colors = colormap;

figure;
hold on;
for c = 1:10
    x = H(X_labels == (c-1), 1);
    y = H(X_labels == (c-1), 2);
    rndidx = randperm(length(x));
    x = x(rndidx(1:500));
    y = y(rndidx(1:500));

    plot(x, y, 'x', 'Color', colors(ceil(c/10 * 64), :));
end
hold off;

author: kyunghyuncho

repo: deepmat


text:
function m = TotalMass(j)
global uLINK

if j == 0
   m = 0;
else
   m = uLINK(j).m + TotalMass(uLINK(j).sister) ...
                  + TotalMass(uLINK(j).child);
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function RegMeanSquareSimilarity(handles, command, metric)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS;
indexS = planC{end};

    [originF, spacingF, centerF] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationBaseDataset));

    [originM, spacingM, centerM] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationMovDataset));

    minstep = str2double(get(handles.para.Similarity_minstep, 'string'));
    maxstep = str2double(get(handles.para.Similarity_maxstep, 'string'));
    iternum = str2double(get(handles.para.Similarity_IterNum, 'string'));

    ScaleFactor = str2double(get(handles.para.Similarity_ScaleFactor, 'string'));
    RotationScale = str2double(get(handles.para.Similarity_RS, 'string'));
    TranslationScale = str2double(get(handles.para.Similarity_TS, 'string'));
    InitialAngle = str2double(get(handles.para.Similarity_InitAngle, 'string'));
    RegLevels = str2double(get(handles.para.Similarity_RegLevels, 'string'));

    output = cell(1, 16);

    FImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
    MImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray;

    dimF = size(FImg);
    dimM = size(MImg);
    %generate cliped base dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseTrans'));
    if ~isempty(clipBox)
        FImg = FImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originF(1) = originF(1) + spacingF(1)*double(clipBox(1)-1);
        originF(2) = originF(2) + spacingF(2)*double(uint16(dimF(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseSag'));
    if ~isempty(clipBox)
        FImg = FImg(:, :, clipBox(2):clipBox(4));
        originF(3) = originF(3) + spacingF(3)*double(clipBox(2)-1);
    end

    %generate cliped move dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movTrans'));
    if ~isempty(clipBox)
        MImg = MImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originM(1) = originM(1) + spacingM(1)*double(clipBox(1)-1);
        originM(2) = originM(2) + spacingM(2)*double(uint16(dimM(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movSag'));
    if ~isempty(clipBox)
        MImg = MImg(:, :, clipBox(2):clipBox(4));
        originM(3) = originM(3) + spacingM(3)*double(clipBox(2)-1);
    end

%downsample the input datasets
    dsampled = 0;
    if (get(handles.dsampleCheck,'Value') == get(handles.dsampleCheck,'Max'))

        tic;
        xyScale =  2; zScale = 2;
        spacingF = [spacingF(1)*xyScale spacingF(2)*xyScale spacingF(3)*zScale];
        spacingM = [spacingM(1)*xyScale spacingM(2)*xyScale spacingM(3)*zScale];
        disp('downSampling ...');
%         FImg = imdesample3d(FImg,xyScale,zScale);
%         MImg = imdesample3d(MImg,xyScale,zScale);
%         FImg=GPReduce2D(FImg,1);
%         MImg=GPReduce2(MImg,1);
        FImg = FImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);
        MImg = MImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);

        toc;
        dsampled = 1;
    end

    switch metric
        case 'mean squares'
            metricIndex = 0;
        case 'normalized correlation'
            metricIndex = 1;
    end

    tic;

%flip the source datasets on X for itk coordinate system
    fdim = 1;
    FImg = flipdim(FImg, fdim);
    MImg = flipdim(MImg, fdim);

% transform initializing modes 0:MomentON 1:GeometryOn 2:initail transform On
    initMode = get(handles.InitTrans, 'value');

%prepare the initial transform matrix
    transMB = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).transM;
    transMM = planC{indexS.scan}(stateS.imageRegistrationMovDataset).transM;
    if isempty(transMM), transMM = eye(4); end;
    if isempty(transMB), transMB = eye(4); end;
    transM = inv(transMB) * transMM;

%flip the moving dataset to match the initial transM
    Tf = eye(4);
    flipM = 0;
    if (get(handles.flipX, 'value'))
        MImg = flipdim(MImg, 2);
        Tf = [-1 0 0 2*centerM(1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
        flipM = 2;
    end
    if (get(handles.flipY, 'value'))
        MImg = flipdim(MImg, 1);
        Tf = [1 0 0 0; 0 -1 0 2*centerM(2); 0 0 1 0; 0 0 0 1];
        flipM = 1;
    end
    if (get(handles.flipZ, 'value'))
        MImg = flipdim(MImg, 3);
        Tf = [1 0 0 0; 0 1 0 0; 0 0 -1 2*centerM(3); 0 0 0 1];
        flipM = 3;
    end
    transM = transM * inv(Tf);
    rotM = transM(1:3, 1:3); transV = transM(1:3, 4);
    rotM = rotM';
    transV = -transV;

%run the registration method
    if strcmpi(command, 'mr')
        try
            [im, Rotation, Offset, TrackData] = MultiRes_Similarity3D(int16(FImg), originF, spacingF, ...
                                                int16(MImg), originM, spacingM, ...
                                                minstep, maxstep, iternum, ...
                                                ScaleFactor, RotationScale, TranslationScale, InitialAngle, ...
                                                RegLevels, metricIndex, rotM, transV, initMode);
        catch
            [im, Rotation, Offset, TrackData] = MultiRes_Similarity3D_64(int16(FImg), originF, spacingF, ...
                                                int16(MImg), originM, spacingM, ...
                                                minstep, maxstep, iternum, ...
                                                ScaleFactor, RotationScale, TranslationScale, InitialAngle, ...
                                                RegLevels, metricIndex, rotM, transV, initMode);
        end
    else
        try
            [im, Rotation, Offset, TrackData] = Similarity3D(int16(FImg), originF, spacingF, ...
                                                int16(MImg), originM, spacingM, ...
                                                minstep, maxstep, iternum, ...
                                                ScaleFactor, RotationScale, TranslationScale, InitialAngle, ...
                                                rotM, transV, initMode);
        catch
            [im, Rotation, Offset, TrackData] = Similarity3D_64(int16(FImg), originF, spacingF, ...
                                                int16(MImg), originM, spacingM, ...
                                                minstep, maxstep, iternum, ...
                                                ScaleFactor, RotationScale, TranslationScale, InitialAngle, ...
                                                rotM, transV, initMode);
        end
    end

    im = flipdim(im, fdim);

    toc;
    RegTime = num2str(ceil(toc/60));

    output{1} = ['Rotation Angle Z = ' num2str(asin(Rotation(1, 2))*45.0/atan(1.0))];
    output{2} = ['Translation X = ' num2str(Offset(1))];
    output{3} = ['Translation Y = ' num2str(Offset(2))];
    output{4} = ['Translation Z = ' num2str(Offset(3))];
    output{5} = ['Number of Iterations = ' num2str(Offset(7))];
    output{6} = ['Best Value = ' num2str(Offset(8))];

    output{7} = ['Rotation Center X = ' num2str(Offset(9))];
    output{8} = ['Rotation Center Y = ' num2str(Offset(10))];
    output{9} = ['Rotation Center Z = ' num2str(Offset(11))];
    output{10} = ['Scale = ' num2str(Offset(12))];
    output{11} = ['Offset X = ' num2str(Offset(13))];
    output{12} = ['Offset Y = ' num2str(Offset(14))];
    output{13} = ['Offset Z = ' num2str(Offset(15))];

    set(handles.OutputList, 'string', output);

%update the transM;
    Tv = [Offset(1) Offset(2) Offset(3)];
    Cv = [Offset(7) Offset(8) Offset(9)];
    Cv(3) = -Cv(3); %CERR z axis is opposite to the ITK.

    rot = reshape(Rotation, 3,3);
    offset = [Offset(13) Offset(14) Offset(15)]';

%     rot(2,1) = -rot(2,1); rot(1,2) = -rot(1,2);
%     rot(1,3) = -rot(1,3); rot(3,1) = -rot(3,1);
%     rot(2,3) = -rot(2,3); rot(3,2) = -rot(3,2);
%
%     offset(1) = -offset(1);
%     offset(2) = -offset(2);
%     offset(3) = -offset(3);

    TM = eye(4);
    TM(:,4) = [offset; 1];

    RM = eye(4);
    RM(1:3, 1:3) = rot;

    newTransform = inv(TM*RM);
    newTransform =  transMB * newTransform * Tf;

    scanSetM = stateS.imageRegistrationMovDataset;
    scanSetF = stateS.imageRegistrationBaseDataset;
    planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = newTransform;

% save the resampled dataset
    if (get(handles.saveCheck, 'value'))
        if (~dsampled)
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = im;
        else
            %imReg = resample(scanSetF, scanSetM); % need to up-sample im ????
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = imReg;
        end
        controlFrame('refresh');
    end

    sliceCallBack('refresh');

end