text:
function transformFullMap2XML(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:
%
%   transformFullMap2XML(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

    % 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 species
        if ~isempty(map.specNotes{spec})
            xmlStruct.sbml.model.listOfSpecies.species{spec}.notes.html.body.Text = map.specNotes{spec};
        end
    end
    clearvars spec

    % Check if there is any information on Compartments
    if ~isempty(map.complexAlias)
        % Loop over complexes to put the info back to the struct for XML
        % conversion
        for k = 1:length(map.complexAlias)
            map.complexAlias{k, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.Attributes.id;
            map.complexID{k, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.Attributes.species;
            % Test if info on compartment exists
            if ~isempty(map.complexCompartAlias{k})
                xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.Attributes.compartmentAlias = map.complexCompartAlias{k};
            end
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.celldesigner_colon_bounds.Attributes.x = map.complexXPos{k};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.celldesigner_colon_bounds.Attributes.y = map.complexYPos{k};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.celldesigner_colon_bounds.Attributes.w = map.complexWidth{k};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.celldesigner_colon_bounds.Attributes.h = map.complexHeight{k};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfComplexSpeciesAliases.celldesigner_colon_complexSpeciesAlias{k}.celldesigner_colon_usualView.celldesigner_colon_paint.Attributes.color = map.complexColor{k};
        end
    end

    % Test if included species exist in the map
    if ~isempty(map.specIncID)
        % Loop over included species to get the need information and
        % store it in a structure. Included species refer to each
        % individual species included inside a k.
        for included = 1:length(map.specIncID)
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfIncludedSpecies.celldesigner_colon_species{included}.Attributes.id = map.specIncID{included};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfIncludedSpecies.celldesigner_colon_species{included}.Attributes.name = map.specIncName{included};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfIncludedSpecies.celldesigner_colon_species{included}.celldesigner_colon_annotation.celldesigner_colon_complexSpecies.Text = map.specIncCplxID{included};
            xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfIncludedSpecies.celldesigner_colon_species{included}.celldesigner_colon_annotation.celldesigner_colon_speciesIdentity.celldesigner_colon_class.Text = map.specIncType{included};
            % Test if there are notes for the included species
            if ~isempty(map.specIncNotes)
                xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfIncludedSpecies.celldesigner_colon_species{included}.celldesigner_colon_notes.html.body.Text = map.specIncNotes{included};
            end
        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};
        xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.reversible = map.rxnReversibility{react};
        % 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
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink.Attributes.alias = map.rxnReactantAlias{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink.Attributes.reactant = map.rxnReactantID{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink.celldesigner_colon_line.Attributes.type = map.rxnReactantLineType{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink.celldesigner_colon_line.Attributes.color = map.rxnReactantLineColor{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink.celldesigner_colon_line.Attributes.width = map.rxnReactantLineWidth{react}{1};
            else
                for x = 1:length(map.rxnReactantAlias{react})
                    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.type = map.rxnReactantLineType{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.rxnReactantLineColor{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{x}.celldesigner_colon_line.Attributes.width = map.rxnReactantLineWidth{react}{x};
                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
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink.Attributes.alias = map.rxnProductAlias{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink.Attributes.product = map.rxnProductID{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink.celldesigner_colon_line.Attributes.type = map.rxnProductLineType{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink.celldesigner_colon_line.Attributes.color = map.rxnProductLineColor{react}{1};
                xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink.celldesigner_colon_line.Attributes.width = map.rxnProductLineWidth{react}{1};
            else
                for x = 1:length(map.rxnProductAlias{react})
                    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.type = map.rxnProductLineType{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.rxnProductLineColor{react}{x};
                    xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{x}.celldesigner_colon_line.Attributes.width = map.rxnProductLineWidth{react}{x};
                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

end

author: opencobra

repo: cobratoolbox


text:
function cvt_write ( ndim, n, batch, seed_init, seed, init_string, it_max, ...
  it_fixed, it_num, it_diff, energy, sample_string, sample_num, r, ...
  file_out_name )

%*****************************************************************************80
%
%% CVT_WRITE writes a CVT dataset to a file.
%
%  Discussion:
%
%    The initial lines of the file are comments, which begin with a
%    '#' character.
%
%    Thereafter, each line of the file contains the NDIM-dimensional
%    components of the next entry in the set.
%
%  Modified:
%
%    17 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NDIM, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer BATCH, sets the maximum number of sample points
%    generated at one time.  It is inefficient to generate the sample
%    points 1 at a time, but memory intensive to generate them all
%    at once.  You might set BATCH to min ( SAMPLE_NUM, 10000 ), for instance.
%
%    Input, integer SEED_INIT, the initial random number seed.
%
%    Input, integer SEED, the current random number seed.
%
%    Input, character INIT_STRING(*), specifies how the initial
%    generators are chosen:
%    filename, by reading data from a file;
%    'GRID', picking points from a grid;
%    'HALTON', from a Halton sequence;
%    'RAND', using MATLAB RAND function;
%    'UNIFORM', using a simple uniform RNG;
%
%    Input, integer IT_MAX, the maximum number of iterations allowed.
%
%    Input, integer IT_FIXED, the maximum number of iterations to take
%    with a fixed set of sample points.
%
%    Input, integer IT_NUM, the actual number of iterations taken.
%
%    Input, real IT_DIFF, the L2 norm of the change
%    in the CVT coordinates on the last iteration.
%
%    Input, real ENERGY, the discrete "energy", divided
%    by the number of sample points.
%
%    Input, string SAMPLE_STRING, specifies how the region is sampled:
%    'GRID', picking points from a grid;
%    'HALTON', from a Halton sequence;
%    'RAND', using MATLAB RAND function;
%    'UNIFORM', using a simple uniform RNG;
%
%    Input, integer SAMPLE_NUM, the number of sampling points used on
%    each iteration.
%
%    Input, real R(NDIM,N), the points.
%
%    Input, string FILE_OUT_NAME, the name of
%    the output file.
%
  file_out_unit = fopen ( file_out_name, 'w' );

  if ( file_out_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_WRITE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file:\n' );
    fprintf ( 1, '  "%s".\n', file_out_name );
    error ( 'CVT_WRITE - Fatal error!' );
  end

  fprintf ( file_out_unit, '#  %s\n', file_out_name );
  fprintf ( file_out_unit, '#  created by CVT_WRITE.M\n' );
  fprintf ( file_out_unit, '#\n' );
  fprintf ( file_out_unit, '#  Spatial dimension NDIM =   %12d\n', ndim );
  fprintf ( file_out_unit, '#  Number of points N =       %12d\n', n );
  fprintf ( file_out_unit, '#  Initial SEED_INIT =        %12d\n', seed_init );
  fprintf ( file_out_unit, '#  Current SEED =             %12d\n', seed );
  fprintf ( file_out_unit, '#  INIT =                    "%s"\n', init_string );
  fprintf ( file_out_unit, '#  Max iterations IT_MAX =    %12d\n', it_max );
  fprintf ( file_out_unit, '#  IT_FIXED (fixed samples) = %12d\n', it_fixed );
  fprintf ( file_out_unit, '#  Iterations IT_NUM =        %12d\n', it_num );
  fprintf ( file_out_unit, '#  Difference IT_DIFF =       %12f\n', it_diff );
  fprintf ( file_out_unit, '#  CVT ENERGY =               %12f\n', energy );
  fprintf ( file_out_unit, '#  SAMPLE =                  "%s"\n', sample_string );
  fprintf ( file_out_unit, '#  SAMPLE_NUM =               %12d\n', sample_num );
  fprintf ( file_out_unit, '#  Sample BATCH size =        %12d\n', batch );
  fprintf ( file_out_unit, '#  EPSILON (unit roundoff ) = %e\n', eps );
  fprintf ( file_out_unit, '#\n' );

  for j = 1 : n
    for i = 1 : ndim
      fprintf ( file_out_unit, '  %10f', r(i,j) );
    end
    fprintf ( file_out_unit, '\n' );
  end

  fclose ( file_out_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fileName = gpWriteResult(model, dataSet, number)

% GPWRITERESULT Write a GP result.
% FORMAT
% DESC writes a GP result.
% ARG model : the model to write.
% ARG dataSet : the name of the data set to write.
% ARG number : the number of the GP run to write.
% RETURN fileName : the file name used to write.
%
% SEEALSO : gpLoadResult
%
% COPYRIGHT : Neil D. Lawrence, 2009

% GP

dataSet(1) = upper(dataSet(1));
type = model.type;
type(1) = upper(type(1));
fileName = ['dem' dataSet type num2str(number)];

[kern, noise, gpInfo] = gpDeconstruct(model);

save(fileName, 'kern', 'noise', 'gpInfo');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ npart, a, rank ] = partn_successor ( n, nmax, npart, a, rank )

%*****************************************************************************80
%
%% PARTN_SUCCESSOR computes partitions whose largest part is NMAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Input, integer NMAX, the maximum size of any part of the
%    partition.  1 <= NMAX <= N.
%
%    Input/output, integer NPART, the number of parts of the
%    partition.  1 <= NPART <= N.
%
%    Input/output, integer A(N), contains the partition.
%    A(1) through A(NPART) contain the nonzero integers which
%    sum to N.
%
%    Input/output, integer RANK, the rank.
%    If RANK = -1 on input, then the routine understands that this is
%    the first call, and that the user wishes the routine to supply
%    the first element in the ordering, which has RANK = 0.
%    In general, the input value of RANK is increased by 1 for output,
%    unless the very last element of the ordering was input, in which
%    case the output value of RANK is 0.
%

%
%  Return the first element.
%
  if ( rank == -1 )
    a(1) = nmax;
    npart = n + 1 - nmax;
    a(2:npart) = 1;
    rank = 0;
    return
  end
%
%  Check.
%
  ierror = partn_sf_check ( n, nmax, npart, a );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARTN_SUCCESSOR - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  IERROR = %d\n', ierror );
    error ( 'PARTN_SUCCESSOR - Fatal error!' );
  end
%
%  If there are at least two parts, and the next to last is not NMAX,
%  then rob the last part and pay the next to the last part.
%  Then, if the next to last part is too big, swap it leftwards.
%
  if ( 1 < npart )

    if ( a(npart-1) < nmax )

      a(npart) = a(npart) - 1;
      a(npart-1) = a(npart-1) + 1;
      index = npart - 1;

      while ( 1 )

        if ( index <= 1 )
          break
        end

        if ( a(index) <= a(index-1) )
          break
        end

        temp       = a(index-1);
        a(index-1) = a(index);
        a(index)   = temp;

        index = index - 1;

      end
%
%  Sum the tail.
%
      temp = sum ( a(index+1:npart) );
%
%  Spread the sum as 1's.
%
      npart = index + temp;
      a(index+1:npart) = 1;
      rank = rank + 1;
      return

    end
%
%  Otherwise, we've reached the last item.
%  Return the first one.
%
  else

    npart = n + 1 - nmax;
    a(1) = nmax;
    a(2:npart) = 1;
    rank = 0;
    return

  end

  return
end


text:
function E = edges4connected(height,width)

% EDGE4CONNECTED Creates edges where each node
%   is connected to its four adjacent neighbors on a
%   height x width grid.
%   E - a vector in which each row i represents an edge
%   E(i,1) --> E(i,2). The edges are listed is in the following
%   neighbor order: down,up,right,left, where nodes
%   indices are taken column-major.
%
%   (c) 2008 Michael Rubinstein, WDI R&D and IDC
%   $Revision$
%   $Date$
%

N = height*width;
I = []; J = [];
% connect vertically (down, then up)
is = [1:N]'; is([height:height:N])=[];
js = is+1;
I = [I;is;js];
J = [J;js;is];
% connect horizontally (right, then left)
is = [1:N-height]';
js = is+height;
I = [I;is;js];
J = [J;js;is];

E = [I,J];

end

author: Sable

repo: mcbench-benchmarks


text:
function matches_refined = refine_matches(...
    part, model, part_ind, model_ind, C_init, options)

F = part_ind;
G = model_ind;

k = options.k;

W = zeros(k);
for i=1:k
    for j=1:k
        slope = 1;
        direction = [1 slope];
        direction = direction./norm(direction);
        W(i,j) = exp(-0.03*sqrt(i.^2 + j.^2))*norm(cross([direction 0], [i,j, 0]-[1 1 0]));
    end
end
d=ones(1,k);
D = repmat(d,k,1);

mu1 = 1e-2; % diagonal mask
mu2 = 1e1; % orthogonality

for iter=1:options.refine_iters

    A = part.evecs'*part.S*F;
    B = model.evecs'*model.S*G;

    manifold = euclideanfactory(k,k);
    problem = {};

    problem.M = manifold;

    problem.cost = @(C) (...
        sum(sum((C*A-B).^2).^0.5) + ...
        mu1 * norm(C.*W,'fro')^2 + ...
        mu2 * (norm(C'*C,'fro')^2 - sum(diag(C'*C).^2) + sum((diag(C'*C) - d').^2) ));

    problem.egrad = @(C) (...
        norm_21_gradient(C,A,B) + ...
        mu1 * 2 * C.*W.*W + ...
        mu2 * 4*(C*C'*C - C.*D ));

    options.verbosity = 2;
%     options.maxiter = 5e3;
%     C_refined = conjugategradient(problem, C_init, options);
    options.maxiter = 3e2;
    C_refined = trustregions(problem, C_init, options);

%     figure,colormap(bluewhitered)
%     subplot(121),imagesc(C_init),colorbar,axis image
%     subplot(122),imagesc(C_refined),colorbar,axis image

    [matches_refined, ~] = flann_search(...
        model.evecs', ...
        C_refined*part.evecs', ...
        1, struct());

%     [matches_init, ~] = flann_search(...
%         model.evecs', ...
%         C_init*part.evecs', ...
%         1, struct());

%     colors = create_colormap(model,model);
%     figure
%     subplot(231), colormap(colors), plot_scalar_map(model, 1:model.n), axis off, view([0 90]), freeze_colors
%     subplot(232), colormap(colors(matches_init,:)), plot_scalar_map(part, 1:part.n), axis off, view([0 90]), freeze_colors
%     subplot(233), colormap(colors(matches_refined,:)), plot_scalar_map(part, 1:part.n), axis off, view([0 90])
%     subplot(234), colormap(colors), plot_scalar_map(model, 1:model.n), axis off, view([-180 -90]), freeze_colors
%     subplot(235), colormap(colors(matches_init,:)), plot_scalar_map(part, 1:part.n), axis off, view([-180 -90]), freeze_colors
%     subplot(236), colormap(colors(matches_refined,:)), plot_scalar_map(part, 1:part.n), axis off, view([-180 -90])

    C_init = C_refined;

    fps = fps_euclidean(part.VERT, 1e3, randi(part.n));
    F = sparse(fps, 1:length(fps), 1, part.n, length(fps));
    G = sparse(matches_refined(fps), 1:length(fps), 1, model.n, length(fps));

end

end

author: OshriHalimi


text:
function accuracy = src(TrainSet, TestSet, train_num, test_num, class_num, lambda, options)
% Sparse representation classification (SRC) algorithm
%
% Inputs:
%       TrainSet            train sets of size dxn, where d is dimension and n is number of sets
%       TestSet             test sets of size dxn, where d is dimension and n is number of sets
%       test_num            numner of test sets
%       class_num           numner of classes
%       lambda              regularization paramter
% Output:
%       accuracy            classification accurary
%
% References:
%       J. Wright, A. Yang, A. Ganesh, S. Sastry, and Y. Ma,
%       "Robust face recognition via sparse representation,"
%       IEEE Transaction on Pattern Analysis and Machine Intelligence, vol.31, no.2, pp.210-227, 2009.
%
%
% Created by H.Kasai on July 06, 2017

    % extract options
    if ~isfield(options, 'verbose')
        verbose = false;
    else
        verbose = options.verbose;
    end

    if ~isfield(options, 'eigenface')
        eigenface = true;
    else
        eigenface = options.eigenface;
    end

    if ~isfield(options, 'eigenface_dim')
        eigenface_dim = train_num;
    else
        eigenface_dim = options.eigenface_dim;
    end

    % generate eigenface
    if eigenface
        [disc_set, ~, ~] = Eigenface_f(TrainSet.X, eigenface_dim);

        % project on subspace
        TrainSet.X  =  disc_set' * TrainSet.X;
        TestSet.X   =  disc_set' * TestSet.X;
    end

    % normalize data to l2-norm
    [TrainSet.X, ~] = data_normalization(TrainSet.X, TrainSet.y, 'std');
    [TestSet.X, ~] = data_normalization(TestSet.X, TestSet.y, 'std');

    % prepare class array
    classes = unique(TrainSet.y);

    % prepare predicted label array
    identity = zeros(1, test_num);

    for i = 1 : test_num

        y = TestSet.X(:, i);

        if 0
            % calculate sparse code
            %tau = max(1e-4*max(abs(TrainSet.X'*y)),sigma*sqrt(log(train_num)));
            tau = lambda;

            in = [];
            in.tau = tau;
            delx_mode = 'mil'; % mil or qr
            in.delx_mode = delx_mode;
            in.debias = 0;
            in.verbose = 0;
            in.plots = 0;

            out = l1homotopy(TrainSet.X, y, in);
            xp = out.x_out;
        elseif 1

            P = inv(TrainSet.X'*TrainSet.X+0.001*eye(size(TrainSet.X,2)))*TrainSet.X';
            x0 = P*y;

            maxIteration = 5000;
            isNonnegative = false;
            lambda = 1e-2; %5e-3;
            tolerance = 0.05;
            STOPPING_GROUND_TRUTH = -1;
            STOPPING_DUALITY_GAP = 1;
            STOPPING_SPARSE_SUPPORT = 2;
            STOPPING_OBJECTIVE_VALUE = 3;
            STOPPING_SUBGRADIENT = 4;
            stoppingCriterion = STOPPING_GROUND_TRUTH;
            [xp, iterationCount] = SolveHomotopy(TrainSet.X, y, ...
                            'maxIteration', maxIteration,...
                            'isNonnegative', isNonnegative, ...
                            'stoppingCriterion', stoppingCriterion, ...
                            'groundtruth', x0, ...
                            'lambda', lambda, ...
                            'tolerance', tolerance);

        elseif 0
            xp = l1_ls(TrainSet.X, y, 1e-3);
        else

            param.lambda = lambda;
            param.lambda2 =  0;
            param.mode = 0;
            xp = full(mexLasso(y, TrainSet.X, param));
        end

        % prepare residual array
        residuals = zeros(1, class_num);

        % calculate residual for each class
        for j = 1 : class_num
            idx = find(TrainSet.y == classes(j));
            %residuals(j) = norm(y-TrainSet.X(:,idx)*xp(idx))/sum(xp(idx).*xp(idx));
            residuals(j) = norm(y - TrainSet.X(:,idx)*xp(idx));
        end

        % calculate the predicted label with minimum residual
        [~, label] = min(residuals);
        identity(i) = label;

        if verbose
            correct = (label == TestSet.y(1, i));
            fprintf('# SRC: test:%04d, predict class: %03d --> ground truth :%03d (%d)\n', i, label, TestSet.y(1, i), correct);
        end

        %identity(i) = src_based_classifier(xp, TrainSet.X, TrainSet.y, TestSet.X(:, i), TestSet.y(1, i), classes, i, class_num, 'SRC', verbose);

    end

    % calculate accuracy
    correct_num = sum(identity == TestSet.y);
    accuracy = correct_num/test_num;
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [res, runTime] = run_tracker(seqLen, detections)

%% setting parameters for tracking
c_en      = 12;     %% 10 birth cost, has no influence on the result
c_ex      = 9;     %% 10 death cost, has no influence on the result
c_ij      = 0;     %% 0 transition cost
betta     = 0.2;   %% 0.2 betta, increase will have less tracks, for every single detection
max_it    = inf;    %% inf max number of iterations (max number of tracks)
thr_cost  = 15;     %% 18 max acceptable cost for a track (increase it to have more tracks.), for every tracklet 19.8

time_start = tic;
%% Run object/human detector on all frames.
frameNums = 1:seqLen;
dres = greedy_detect_generator(detections, frameNums);

%% Running tracking algorithms
dres_dp_nms   = tracking_dp(dres, c_en, c_ex, c_ij, betta, thr_cost, max_it, 1);
dres_dp_nms.r = -dres_dp_nms.id;
runTime = toc(time_start);

%% save the tracking result
bboxes_tracked = dres2bboxes(dres_dp_nms, numel(frameNums));  %% we are visualizing the "DP with NMS in the lop" results. Can be changed to show the results of DP or push relabel algorithm.
res = saveResTxt(bboxes_tracked);

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function filter = helperInitIMMFilter(detection)
% This is a helper function and may be removed in a future release.
% This function initializes an IMM filter for tracking using lidar data.

% Copyright 2019 The MathWorks, Inc.

% Create individual filters for constant turn rate and constant velocity.
trackingfilters = cell(2,1);
trackingfilters{1} = initCTCuboidFilter(detection);
trackingfilters{2} = initCVCuboidFilter(detection);

% Construct IMM and inform it about switching between models.
filter = trackingIMM(trackingfilters,'ModelConversionFcn',@switchimmCuboid);
filter.MeasurementNoise = detection.MeasurementNoise;
filter.TransitionProbabilities = 0.99;

end

function x = switchimmCuboid(modelType1,x1,modelType2,x2)
% This function simply removes helper and Cuboid from the model names and
% uses the switchimm function for switching kinematic state. It appends the
% non-kinematic part after the switch.
modelType1 = strrep(modelType1,'Cuboid','');
modelType1 = strrep(modelType1,'helper','');
modelType2 = strrep(modelType2,'Cuboid','');
modelType2 = strrep(modelType2,'helper','');

isCov = ~isvector(x1);
n1 = size(x1,1);
n2 = size(x2,1);
if ~isCov
    kinematicPartX1 = x1(1:(n1-4));
    shapePartX1 = x1((n1-3):n1);
    kinematicPartX2 = x2(1:(n2-4));
else
    kinematicPartX1 = x1(1:(n1-4),1:(n1-4));
    shapePartX1 = x1((n1-3):n1,(n1-3):n1);
    kinematicPartX2 = x1(1:(n2-4),1:(n2-4));
end

kinematicPartX = switchimm(modelType1,kinematicPartX1,modelType2,kinematicPartX2);

if isCov
    x = blkdiag(kinematicPartX,shapePartX1);
else
    x = [kinematicPartX;shapePartX1];
end

end

function filter = initCVCuboidFilter(detection)
    % This function initializes a constant-velocity cuboid filter from a
    % detection report.
    % detection contains measurement with the following convention:
    % [x;y;z;l;w;h];

    posIndex = [1 3 5];
    dimIndex = [8 9 10];
    yawIndex = 7;

    meas = detection.Measurement;
    measCov = detection.MeasurementNoise;
    dataType = class(meas);
    [pos,posCov,dim,dimCov,yaw,yawCov] = helperInverseLidarModel(meas,measCov);

    % Assemble state and state covariances
    state = zeros(10,1,dataType);
    state(posIndex) = pos;
    state(dimIndex) = dim;
    state(yawIndex) = yaw;
    cov = 100*eye(10,dataType);
    cov(posIndex,posIndex) = posCov;
    cov(dimIndex,dimIndex) = dimCov;
    cov(yawIndex,yawIndex) = yawCov;

    % processNoise. Acceleration and omega
    Q = eye(4);
    Q(4,4) = 2.5;

    % Use a UKF for capture non-linearity.
    filter = trackingUKF(@helperConstvelCuboid,@helperCvmeasCuboid,state,...
         'StateCovariance',cov,...
        'HasAdditiveProcessNoise',false,...
        'ProcessNoise',Q,'MeasurementNoise',detection.MeasurementNoise,...
        'Alpha',0.01);

    setMeasurementSizes(filter,6,6);
end

function filter = initCTCuboidFilter(detection)
    % This function initializes a constant-turn rate cuboid filter from a
    % detection report.
    % detection contains measurement with the following convention:
    % [x;y;z;l;w;h];

    posIndex = [1 3 6];
    dimIndex = [9 10 11];
    yawIndex = 8;

    meas = detection.Measurement;
    measCov = detection.MeasurementNoise;
    dataType = class(meas);
    [pos,posCov,dim,dimCov,yaw,yawCov] = helperInverseLidarModel(meas,measCov);

    % Assemble state and state covariances
    state = zeros(11,1,dataType);
    state(posIndex) = pos;
    state(dimIndex) = dim;
    state(yawIndex) = yaw;
    cov = 100*eye(11,dataType);
    cov(posIndex,posIndex) = posCov;
    cov(dimIndex,dimIndex) = dimCov;
    cov(yawIndex,yawIndex) = yawCov;

    % processNoise
    Q = eye(4);
    Q(4,4) = 25;

    % Use a UKF for capture non-linearity.
    filter = trackingUKF(@helperConstturnCuboid,@helperCtmeasCuboid,state,...
        'StateCovariance',cov,...
        'HasAdditiveProcessNoise',false,...
        'ProcessNoise',Q,...
        'MeasurementNoise',detection.MeasurementNoise,...
        'Alpha',0.01);

    setMeasurementSizes(filter,6,6);
end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
function solution = solveCobraEP(EPproblem, varargin)
% Solves the following optimisation problem:
%
% minimize   osense*(c.*d)'x + d.*x'(log(x) -1) + (1/2)*x'*Q*x
%   x
%
% subject to    A*x    <=> b         : y
%               lb <= x  <=  ub      : z
%
% or
% subject to    blc <= A*x <= buc   : y
%               lb <= x  <=  ub     : z
%
%
% However, when EPproblem.P is present, the following optimisation problem
% is solved:
%
% minimize   osense*(c.*d)'*x  + (d.*x)'*log(x./q) + (1/2)*x'*Q*x = osense*(c.*d)'*x  + (d.*x)'*log(x) - (d.*x)'*log(q) + (1/2)*x'*Q*x
%   x,q
%
% subject to    A*x    <=> b         : y
%               P*x - q =  0         : r
%               lb <= x  <=  ub      : z
%
% or
% subject to    blc <= A*x <= buc   : y
%               P*x  -q =  0        : r
%               lb <= x  <=  ub     : z
%
%
%
% USAGE:
%    solution = solveCobraEP(EPproblem, varargin)
%
% INPUT:
%    EPproblem:     Structure containing the following fields describing the EP problem to be solved
%                     * .A  - m x n Linear constraint matrix
%                     * .c  - n x 1 Linear objective coeff vector
%                     * .lb - n x 1 Lower bound vector
%                     * .ub - n x 1 Upper bound vector
%                     * .d  - n x 1 Non-negative vector indicating the non-negative
%                                   variables whose entropy is maximised. If d(i)==0
%                                   then there is only a linear objective on x(i).
%                     * .osense - Linear objective sense (-1 means maximise, 1 means minimise)
%
%                   With either the following fields
%                     * .b - m x 1 right hand side vector i.e. A <=> b
%                     * .csense - m x 1 string containting the constraint sense for
%                                 each row in A ('E', equality, 'G' greater than, 'L' less than).
%
%                   Or with the following fields
%                     * .blc - m x 1  left hand side vector i.e.  blc <= A*x
%                     * .buc - m x 1 right hand side vector i.e.  A*x <= buc
%
%
% OPTIONAL INPUTS:
%    EPproblem:     Structure containing the following fields describing the EP problem to be solved
%                     * .P - p x n matrix with entries {0,1}, such that q_i := P_{i,:}*x is the sums
%                            of the x corresponding to nonzero columns of the ith row of P, i.e. P_{i,:}.
%                            Used for normalised entropy maximisation.
%
%                     * .Q - positive semidefinite matrix for quadratic part of objective (see above)
%
%    varargin:      Additional parameters either as parameter struct, or as
%                   parameter/value pairs. A combination is possible, if
%                   the parameter struct is either at the beginning or the
%                   end of the optional input.
%                   All fields of the struct which are not COBRA parameters
%                   (see `getCobraSolverParamsOptionsForType`) for this
%                   problem type will be passed on to the solver in a
%                   solver specific manner. Some optional parameters which
%                   can be passed to the function as parameter value pairs,
%                   or as part of the options struct are listed below:
%                  'verify'
%                  'printLevel'
%                  'debug'
%                  'feasTol'
%                  'optTol'
%                  'solver'
%
%
%    printLevel:    Printing level
%
%                     * 0 - Silent (Default)
%                     * 1 - Warnings and Errors
%                     * 2 - Summary information
%                     * 3 - More detailed information
%                     * > 10 - Pause statements, and maximal printing (debug mode)
%
%   solver:         Optimisation solver used, {('mosek'),'pdco'}
%   feasTol:        Feasibility tolerance
%   optTol:         Optimality tolerance
%
% OUTPUT:
%    solution:      Structure containing the following fields describing a LP solution:
%                     * .obj:          Objective value
%                     *.objLinear      osense*c'*x;
%                     *.objEntropy     d.*x'*(log(x) -1);
%                     *.objQuadratic   (1/2)*x'*Q*x;
%                     * .full:         Primal solution vector
%                     * .slack:        bl = A*x + s = bu
%                     * .rcost:        Reduced costs, dual solution to :math:`lb <= x <= ub`
%                     * .dual:         dual solution to constraints :math: `A*x ('E' | 'G' | 'L') b`
%
%                     * .solver:       Solver used to solve EP problem
%                     * .stat:         Solver status in standardized form
%                       * 0 - Infeasible problem
%                       * 1 - Optimal solution
%                       * 2 - Unbounded solution
%                       * 3 - Almost optimal solution
%                       * -1 - Some other problem (timelimit, numerical problem etc)
%                     * .origStat:         Original status returned by the specific solver
%                     * .origStatText:     Original status text returned by the specific solver
%                     * .time:         Solve time in seconds
%
% OPTIONAL OUTPUT (from conic optimisation with mosek):
%  solution.auxPrimal:  auxiliary primal variable
%  solution.auxRcost:   dual to auxiliary primal variable
%  solution.coneF:      affine constraint matrix
%  solution.coneDual:   dual to affine constraints
%  solution.dualNorm:   dual to the probability normalisation constraint
%
% OPTIONAL OUTPUT (from optimisation with pdco):
%  solution.d1:  primal regularisation parameter, see pdco.m
%  solution.d2:  dual regularisation parameter, see pdco.m
%
% EXAMPLE:
%
% NOTE: This code is a draft version released for the ELIXIR Fluxomic course and is not yet published and not to be redistributed without express permission of the author.
%
% Author(s): Ronan M.T. Fleming, 2021

[problemTypeParams, solverParams] = parseSolverParameters('EP', varargin{:});

if ~isfield(problemTypeParams,'debug')
    problemTypeParams.debug = 1;
end

% Remove outer function specific parameters to avoid crashing solver interfaces
% Default EP parameters are removed within solveCobraEP, so are not removed here
solverParams = mosekParamStrip(solverParams);

if any(EPproblem.lb>EPproblem.ub)
    error('EPproblem.lb>EPproblem.ub');
end

% assume constraint A*v = b if csense not provided
if isfield(EPproblem, 'csense')
    bool = EPproblem.csense == 'E' | EPproblem.csense == 'G' | EPproblem.csense == 'L';
    if any(~bool)
        error('Incorrect formulation of EPproblem.csense \n%s','EPproblem.csense must be an m x 1 character array containing the constraint sense {''E'',''L'',''G''} corresponding to each row of EPproblem.A')
    end
end

if isequal(problemTypeParams.solver,'mosek')
    if ~(isfield(EPproblem,'blc') || isfield(EPproblem,'blc'))
        % blc <= A*x <= buc
        EPproblem.blc = EPproblem.b;
        EPproblem.blc(EPproblem.csense == 'L',1) = -inf;
        EPproblem.buc = EPproblem.b;
        EPproblem.buc(EPproblem.csense == 'G',1) = inf;
        %remove other specification of constraints to avoid conflict
        EPproblem = rmfield(EPproblem,'csense');
        EPproblem = rmfield(EPproblem,'b');
    end
end

%% if in debug mode, test to see if the LP part of the problem is feasible
if problemTypeParams.debug
    switch problemTypeParams.solver
        case 'pdco'
            solutionLP2 = solveCobraLP(EPproblem);
            if problemTypeParams.printLevel>2
                disp(solutionLP2)
            end

        case 'mosek'
            %https://docs.mosek.com/8.1/toolbox/solving-linear.html
            if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_PFEAS')
                solverParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol;
            end
            if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.')
                solverParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol;
            end
            %If the feasibility tolerance is changed by the solverParams
            %struct, this needs to be forwarded to the cobra Params for the
            %final consistency test!
            if isfield(problemTypeParams,'MSK_DPAR_INTPNT_TOL_PFEAS')
                solverParams.feasTol = solverParams.MSK_DPAR_INTPNT_TOL_PFEAS;
            end
            [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,solverParams,'minimize');

            %parse mosek result structure
            [solutionLP2.stat,solutionLP2.origStat,x,y,w] = parseMskResult(res);
%             if stat ==1
%                 f=c'*x;
%                 % slack for blc <= A*x <= buc
%                 s = b - A * x; % output the slack variables
%             else
%                 f = NaN;
%                 s = NaN*ones(size(A,1),1);
%             end

            switch solutionLP2.stat
                case 0
                    solution = solutionLP2;
                    message = ['solveCobraEP: LP part of EPproblem is infeasible according to solveCobraLP with ' problemTypeParams.solver '.'];
                    warning(message)

                    return
                case 2
                    solution = solutionLP2;
                    message = ['solveCobraEP: LP part of EPproblem is unbounded according to solveCobraLP with ' problemTypeParams.solver '.'];
                    warning(message)

                    return
                case 1
                    message =['solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with ' problemTypeParams.solver '.'];
                    fprintf('%s\n',message)
                otherwise
                    error('inconclusive solveCobraLP')
            end
            messages = cellstr(message);
    end
end

if ~isfield(EPproblem, 'b')
    EPproblem.b = zeros(size(EPproblem.A, 1), 1);
end

if ~isfield(EPproblem,'d')
    EPproblem.d = zeros(size(EPproblem.A,2),1);
end
if ~isfield(EPproblem,'sumFluxes')
    EPproblem.sumFluxes = [];
end
if ~isfield(EPproblem,'sumConc')
    EPproblem.sumConc = [];
end
if ~isfield(EPproblem,'sumConc0')
    EPproblem.sumConc0 = [];
end
if ~isfield(EPproblem,'Q')
    if isfield(EPproblem,'F')
        % solveCobraQP uses F for a positive semidefinite matrix
        % we use Q instead, because F is used to denote the matrix of
        % affine constraints arising from conic reformulation
        EPproblem.Q = EPproblem.F;
        EPproblem = rmfield(EPproblem,'F');
    end
end

[A,lb,ub,c,osense,d] = ...
    deal(EPproblem.A,EPproblem.lb,EPproblem.ub,EPproblem.c,EPproblem.osense,EPproblem.d);

if isfield(EPproblem,'b')
    b = EPproblem.b;
else
    b = zeros(size(A,1),1);
end

if isfield(EPproblem,'csense')
    csense = EPproblem.csense;
    if isfield(EPproblem,'blc')
        error('Ambiguous specification of EP problem to have EPproblem.blc and EPproblem.csense')
    end
end

if isfield(EPproblem,'blc')
    blc = EPproblem.blc;
    buc = EPproblem.buc;
    if any(blc> buc)
        error('EPproblem.blc must be less than or equal to EPproblem.buc')
    end
    if isfield(EPproblem,'csense')
        error('Ambiguous specification of EP problem to have EPproblem.blc and EPproblem.csense')
    end
end

[mlt,nlt]=size(A);

switch problemTypeParams.solver
    case 'pdco'
        % solves optimization problems of the form
        %
        %    minimize    phi(x) + 1/2 norm(D1*x)^2 + 1/2 norm(r)^2
        %      x,r
        %    subject to  A*x + D2*r = b,   bl <= x <= bu,   r unconstrained,

        %pdco only works with equality constraints and box constraints so
        %any other linear constraints need to be reformulated in terms of
        %slack variables
        %indl = find(csense == 'L'); %  A*x + s =   b
        %indg = find(csense == 'G'); % -A*x + s = - b

        if ~any(csense == 'L' | csense == 'G')
            Aeq  =  A;
            beq  =  b;
            lbeq = lb;
            ubeq = ub;
            ceq  =  osense*c;
            deq  = d;
            if isfield(EPproblem,'Q')
                Q = EPproblem.Q;
            end
        else
            Aeq = A;
            Aeq(csense == 'G',:) = -1*Aeq(csense == 'G',:);
            beq = b;
            beq(csense == 'G',:) = -1*beq(csense == 'G',:);
            K = speye(mlt);
            K = K(:,csense == 'L' | csense == 'G');
            Aeq = [Aeq K];
            nSlacks = nnz(csense == 'L' | csense == 'G');
            lbeq = [lb ; zeros(nSlacks,1)];
            ubeq = [ub ; inf*ones(nSlacks,1)];
            ceq  = [osense*c  ; zeros(nSlacks,1)];
            deq  = [d  ; zeros(nSlacks,1)];

            if isfield(EPproblem,'Q')
                %extend Q matrix to account for slack variables
                Q = sparse(size(Aeq,2),size(Aeq,2));
                Q(1:nlt,1:nlt) = EPproblem.Q;

            end
        end

        %add regularisation in case its not positive definite
        if isfield(EPproblem,'Q')
            try
                R = chol(Q);
                clear R;
            catch ME
                fprintf('%s\n',ME.message)
                Q = Q + diag(sparse(ones(size(Q,1),1)*1e-4));
            end
        end

        if isfield(solverParams,'d1')
            d1 = solverParams.d1;
        else
            d1 = 1e-4;
        end
        if isfield(solverParams,'d2')
            d2 = solverParams.d2;
        else
            d2 = 1e-4;
        end
        if isfield(solverParams,'x0')
            x0 = solverParams.x0;
        else
            x0 = ones(size(Aeq,2),1);
        end
        if isfield(solverParams,'y0')
            y0 = solverParams.y0;
        else
            y0 = ones(size(Aeq,1),1);
        end
        if isfield(solverParams,'z0')
            z0 = solverParams.z0;
        else
            z0 = ones(size(Aeq,2),1);
        end
        if isfield(solverParams,'xsize')
            xsize = solverParams.xsize;
        else
            xsize = 1;
        end
        if isfield(solverParams,'zsize')
            zsize = solverParams.zsize;
        else
            zsize = 1;
        end

        %TODO - still have no idea what the best parameters for pdco are
        options = pdcoSet;
        %options.mu0       = 1; %very small only for entropy function
        options.mu0       = 0; %pdco chooses its own
        options.FeaTol    = problemTypeParams.feasTol;
        options.OptTol    = problemTypeParams.optTol;
        %   If getting linesearch failures, slacken tolerances
        %   i.e. Linesearch failed (nf too big)
        %options.FeaTol    = 1e-6; %%Ecoli core working at 1e-7
        %options.OptTol    = 1e-6;
        %        options.StepSame  = 0; %(allow different primal and dual steps)

        %%%%%%
        %Additional parameter specifications by Ronan
        %increasing to 0.99 reduced the number of iterations required
        %options.StepTol   = 0.9;
        % needed more than 30 iterations when xsize & zsize not tuned set
        options.MaxIter   = 200;
        options.Method = 1;

        if 0
            %options from Michael's pdcotestENTROPY
            xsize = 5/nlt;               % A few elements of x are much bigger than 1/n.
            xsize = min(xsize,1);      % Safeguard for tiny problems.
            zsize = 1;                 % This makes y (sic) about the right size.
            % 10 makes ||y|| even closer to 1,
            % but for some reason doesn't help.

            x0min = xsize;             % Applies to scaled x1, x2
            z0min = zsize;             % Applies to scaled z1, z2

            en    = ones(size(Aeq,2),1);
            x0    = en*xsize;          %
            y0    = zeros(size(Aeq,1),1);
            z0    = en*z0min;          % z is nominally zero (but converges to mu/x)

            d1    = 0;                 % gamma. 1e-3 is normal.  0 seems fine for entropy
            d2    = 1e-3;              % delta

            options = pdcoSet;
            options.MaxIter      =    50;
            options.FeaTol       =  1e-6;
            options.OptTol       =  1e-6;
            options.x0min        = x0min;  % This applies to scaled x1, x2.
            options.z0min        = z0min;  % This applies to scaled z1, z2.
            options.mu0          =  1e-0;  % 09 Dec 2005: BEWARE: mu0 = 1e-5 happens
            %    to be ok for the entropy problem,
            %    but mu0 = 1e-0 is SAFER IN GENERAL.

            options.Method       =     1;  % 1=Chol  2=QR  3=LSQR
            options.LSMRatol1    =  1e-3;
            options.LSMRatol2    =  1e-6;
            options.wait         =     0;
        end

        if 0
            %Michael Saunders suggestion for badly scaled/almost infeasible
            %problem
            options.FeaTol = 1e-8;
            options.OptTol = 1e-8;
            options.Method = 2;
            options.d1    = 1e-4;
            options.d2    = 1e-4;
            options.xsize = 1000;
            options.zsize = 1e+8;
        end

        %set the objective
        if isfield(EPproblem,'Q')
            objHandle = @(x) entropyQPObj(x,ceq,deq,Q);
        else
            objHandle = @(x) entropyObj(x,ceq,deq);
        end
        options.Print = problemTypeParams.printLevel-1;

        saveAndDebug=0;
        if saveAndDebug
            clearvars -except entropyhandle Aeq beq lbeq ubeq d1 d2 options x0 y0 z0 xsize zsize ceq deq saveAndDebug
        end
        tic;
        [x,y,z,inform,~,~,~] = ...
            pdco(objHandle,Aeq,beq,lbeq,ubeq,d1,d2,options,x0,y0,z0,xsize,zsize);

        logx = zeros(length(x),1);
        logx(deq~=0) = reallog(x(deq~=0));     % error if negative
        grad = ceq + deq.*logx;

        if problemTypeParams.printLevel > 2 || problemTypeParams.debug
            % determine the residuals
            fprintf('\n%s\n','KKT with pdco signs:')
            fprintf('%8.2g %s\n',norm(Aeq*x - beq,inf), '|| Aeq*x - beq ||_inf');
            fprintf('%8.2g %s\n',norm(Aeq*x - beq + (d2^2)*y,inf), '|| Aeq*x - beq + (d2^2)*y ||_inf');

            %gradient may differ depending on the solver
            % z includes (d1.^2).*x from the primal regularization.
            fprintf('%8.2g %s\n',norm(grad  - Aeq' * y - z,inf), '|| grad - Aeq''*y - z ||_inf');
            fprintf('%8.2g %s\n',norm(grad  - Aeq' * y - z - (d1^2)*x,inf), '|| grad - Aeq''*y - z - (d1^2)*x ||_inf');
        end
        if saveAndDebug
            save([datestr(now,29) '_pdco_problem_debug.mat'])
            return
        end

        solution.time = toc;

        % inform = 0 if a solution is found;
        %        = 1 if too many iterations were required;
        %        = 2 if the linesearch failed too often;
        %        = 3 if the step lengths became too small;
        %        = 4 if Cholesky said ADDA was not positive definite.
        if (inform == 0)
            solution.stat = 1;

            if ~any(csense == 'L' | csense == 'G')
                slack = zeros(mlt,1);
            else
                slack = zeros(mlt,1);
                slack(csense == 'L' | csense == 'G') = x(nlt+1:end);
                slack(csense == 'G') = -slack(csense == 'G');
                %important to flip the signs of the dual variables for any
                %greater than constraint
                y(csense == 'G') = -y(csense == 'G');
            end

            solution.slack = slack;
            solution.full = x(1:nlt,1);
            solution.dual = -y;
            solution.rcost = -z(1:nlt,1);
            solution.origStat = inform;

            %objective
            logx = zeros(size(A,2),1);
            logx(d~=0) = reallog(solution.full(d~=0));     % error if negative
            if isfield(EPproblem,'Q')
                solution.obj  = c'*solution.full + (d.*solution.full)'*logx + (1/2)*solution.full'*EPproblem.Q*solution.full;
                grad = c + d.*logx + EPproblem.Q*solution.full;
            else
                solution.obj  = c'*solution.full + (d.*solution.full)'*logx;
                grad = c + d.*logx;
            end
            Aty = -A'*y;

            if problemTypeParams.printLevel > 2 || problemTypeParams.debug
                fprintf('\n%s\n','KKT with Rockafellar signs:')
                fprintf('%8.2g %s\n',norm(A*solution.full + solution.slack - b,inf), '|| A*x + s - b ||_inf');
                fprintf('%8.2g %s\n',norm(A*solution.full + solution.slack - b - (d2^2)*solution.dual,inf), '|| A*x + s - b - (d2^2)*y ||_inf');
                res2 = grad  + A'*solution.dual + solution.rcost;
                fprintf('%8.2g %s\n',norm(res2,inf), '|| grad + A''*y + z ||_inf');

                fprintf('%8.2g %s\n',norm(grad  + A'*solution.dual + solution.rcost - (d1^2)*solution.full,inf), '|| grad + A''*y + z - (d1^2)*x ||_inf');
                if problemTypeParams.debug
                    %res2 = grad  + Aty + solution.rcost;
                    res2 = grad  + A'*solution.dual + solution.rcost;
                    solution.T = table(res2,c,d.*logx, -A'*y, solution.rcost,...
                        'VariableNames',{'total','c','dlogx','Aty','z'});
                end
                if any(~isfinite(res2))
                    warning('Infinite variables in dual optimality condition')
                    solution.Tinf=solution.T(~isfinite(solution.T.total),:);
                    ind = find(~isfinite(solution.T.total));
                    solution.Tinf.ind = ind;
                end
            end

            solution.d1=d1;
            solution.d2=d2;

        elseif (inform == 1 || inform == 2 || inform == 3)
            solution.stat = 0;
            solution.obj = NaN;
        else
            solution.stat = -1;
            solution.obj = NaN;
        end
        solution.origStat = inform;

        %update parameters for testing optimality criterion
        problemTypeParams.feasTol = options.FeaTol;
        problemTypeParams.optTol = options.OptTol;
        %                     * .full:         Full LP solution vector
        %                     * .obj:          Objective value
        %                     * .rcost:        Reduced costs, dual solution to :math:`lb <= v <= ub`
        %                     * .dual:         dual solution to `A*v ('E' | 'G' | 'L') b`
        %                     * .solver:       Solver used to solve LP problem
        %                     * .algorithm:    Algorithm used by solver to solve LP problem
        %                     * .stat:         Solver status in standardized form

    case 'mosek'
        %%
        %         https://docs.mosek.com/modeling-cookbook/expo.html
        %         https://docs.mosek.com/modeling-cookbook/qcqo.html#conic-reformulation
        %         min  (d.*x)'*(log(x./y) + c)  + (1/2)*x'*Q*x
        %         s.t. l <= A[x;y] <= u
        %
        %         Assuming Q is positive semidefinite, there exists an F such that Q = F'*F
        %
        %         min  (d.*x)'*(log(x./y) + c)  + (1/2)*x'*(R'*R)*x
        %         s.t. l <= A[x;y] <= u
        %
        %         where d,c,A,l,u,Q are data and x,y are variables, is equivalent to
        %
        %         min   d*e + d*c*x  + q
        %         s.t.   x*log(x/y)    <= e
        %                (1/2)*x'*Q'*x <= s
        %         l <= [A, 0, 0]*[x;y;e;s] <= u
        %
        %         which is equivalent to:
        %
        %         min   d*e + d*c*x + q
        %         s.t.   (y, x, -e) \in K_{exp}     Exponential cone % MSK_CT_PEXP
        %                (1, s, Fx) \in Q^{k+2}_{r} Quadratic cone % MSK_CT_QUAD
        %         l <= A[x;y] <= u
        %
        %         Such a problem could be formulated using the Affine conic constraints, as shown in the following code:

        % subject to    blc <=  A*x      <= buc   : y
        %               0   <= -P*x +  q <= 0     : r
        %               lb  <=    x      <=  ub   : z
        %             -inf  <=    q      <=  inf  : z

        nExpCone  = nnz(d);
        expCone1 = (nExpCone>0)+0;

        if isfield(EPproblem,'Q')
            quadRows  = any(EPproblem.Q,2);
            quadCols  = any(EPproblem.Q,1)';
            quadBool  = quadRows | quadCols;
            nQuadCone = nnz(quadBool);
            R = sparse(nQuadCone,size(A,2));
            %cholesky factorisation of Q
            R(:,quadBool) = chol(EPproblem.Q(quadBool,quadBool));
        else
            R=[];
            nQuadCone = 0;
        end

        if isfield(EPproblem,'P')
            p = size(EPproblem.P,1);
            P = EPproblem.P;
            pBool=(sum(P,1)~=0)'; %identify normalised variables
        else
            p = 0;
            P = [];
            pBool = false(size(d));
        end
        varNotNorm= any(d & ~pBool)+0;
        quadCone1 = (nQuadCone>0)+0;

        Om1 = sparse(size(A,1),varNotNorm);
        Omp = sparse(size(A,1),p);
        Omd = sparse(size(A,1),nExpCone);
        Omq = sparse(size(A,1),nQuadCone);
        Oz1 = sparse(size(A,1),quadCone1);

        Ox1 = sparse(p,varNotNorm);
        Ip  = spdiag(ones(p,1));
        Opd = sparse(p,nExpCone);
        Opq = sparse(p,nQuadCone);
        Op1 = sparse(p,quadCone1);

        prob.a = [...
           %x,   1,   p,   e,   1,    q;
            A, Om1, Omp, Omd, Oz1, Omq;
            P, Ox1, -Ip, Opd, Op1, Opq];

        %This should not be here!!!!!!!!!!!!!1
%         if isfield(EPproblem,'csense')
%             blc(csense == 'E',1) = b(csense == 'E');
%             buc(csense == 'E',1) = b(csense == 'E');
%             blc(csense == 'G',1) = b(csense == 'G');
%             buc(csense == 'G',1) = inf;
%             blc(csense == 'L',1) = -inf;
%             buc(csense == 'L',1) = b(csense == 'L');
%         else
%             blc = b;
%             buc = b;
%         end

        prob.blc = [blc;zeros(p,1)];
        prob.buc = [buc;zeros(p,1)];

        %remove the normalisation constant from the optimality conditions
        %ed = double(d~=0);

        %maximise the variable corresponding to the exponential cone
        %minimise the variable corresponding to the quadratic cone
        %                x,                  1,         p,       e,          1,                       q;
        prob.c = [osense*c;zeros(varNotNorm,1);zeros(p,1);-d(d~=0);zeros(quadCone1,1);ones(nQuadCone,quadCone1)];

        %fix the non-normalised variables corresponding to the exponential cone to one y = 1
        %            x,                  1,                                                          p,                    e, 1,                     q
        prob.blx = [lb; ones(varNotNorm,1);                                                 zeros(p,1);-inf*ones(nExpCone,1); ones(quadCone1,1);-inf*ones(nQuadCone,1)];
        prob.bux = [ub; ones(varNotNorm,1); EPproblem.sumFluxes; EPproblem.sumConc; EPproblem.sumConc0; inf*ones(nExpCone,1); ones(quadCone1,1); inf*ones(nQuadCone,1)];

        % Specify conic part of the problem
        % https://docs.mosek.com/9.2/toolbox/data-types.html#cones
        if problemTypeParams.printLevel>1 || problemTypeParams.debug
            [~, res] = mosekopt('symbcon');
        else
            [~, res] = mosekopt('symbcon echo(0)');
        end
%         For affine conic constraints Fx+g∈, where =1×⋯×s, cones is a list consisting of s
%         concatenated cone descriptions. If a cone requires no additional parameters (quadratic, rotated quadratic, exponential, zero) then its description is
%         [type,len]
%         where type is the type (conetype) and len is the length (dimension). The length must be present.
        %cone type
        prob.cones(1:2:2*nExpCone) = res.symbcon.MSK_CT_PEXP;
        nCone = nExpCone+nQuadCone;
        if nQuadCone>0
            prob.cones(2*nExpCone+1:2:2*nCone) = res.symbcon.MSK_CT_QUAD;
        end
        %dimensions of cone
        prob.cones(2:2:2*nCone) = 3;

        %Conic problem with affine conic constraints
        %https://docs.mosek.com/9.2/toolbox/data-types.html#equation-doc-notation-conic
        %f (double[][]) – The matrix of affine conic constraints. It must be a sparse matrix.

        % Primal exponential cone.
        %  max e <= x * log(x/1),      x >= 0  <=> [ 1; x; e] \in K_{exp}
        %  min e >= x * log(1/x),      x >= 0  <=> [ 1; x;-e] \in K_{exp}
        % x3 <= x2 * log(x1 /x2), x1, x2 >= 0  <=> [x1;x2;x3] \in K_{exp}
        % x1 = 1 (no normalisation)
        % x2 = x
        % x3 = e

        % Primal exponential cone.
        % -e <=  x * log(y /x),    y,  x >= 0  <=> [y;  x;-e] \in K_{exp}
        %  e >=  x * log(x /y),    y,  x >= 0  <=> [y;  x; e] \in K_{exp}
        % x3 <= x2 * log(x1 /x2), x1, x2 >= 0  <=> [x1;x2;x3] \in K_{exp}
        % x1 = 1 or y (normalisation)
        % x2 = x
        % x3 = e
        Id=speye(nExpCone);
        Idn=speye(size(A,2));
        Idn = Idn(d~=0,:);
        Od1 = sparse(nExpCone,varNotNorm);
        Odp = sparse(nExpCone,p);
        Od  = sparse(nExpCone,nExpCone);
        Odn = sparse(nExpCone,size(A,2));

        if varNotNorm
            %entropy maximisation without normalisation
            Id1 = ones(size(A,2),varNotNorm);
            Id1(pBool)=0;
            Id1=Id1(d~=0); %zero out normalised variables from y=1;
        else
            %relative entropy maximsation with normalisation
            Id1 = sparse(nExpCone,varNotNorm);
        end
        Ox1 = sparse(nQuadCone,varNotNorm);
        Oz1 = sparse(nExpCone,quadCone1);

        if isfield(EPproblem,'P')
            Idp = P(:,d~=0)';
            Oqp = sparse(nQuadCone,size(Idp,2));
        else
            Idp = [];
            Oqp = [];
        end

        % Quadratic cone
        % https://docs.mosek.com/modeling-cookbook/cqo.html#convex-quadratic-sets
        % (1/2)*x'*(R'*R)*x <= q            <=>  [ q; 1; R*x] \in Q^{k+2}_{r} Quadratic cone
        % (1/2)*x3'*(R'*R)*x3 <= x1, x2 = 1 <=>  [x1;x2;R*x3] \in Q^{k+2}_{r} Quadratic cone
        % x1 =   q
        % x2 =   1
        % x3 = R*x
        Iq  =  speye(nQuadCone);
        Iq1 =   ones(nQuadCone,quadCone1);
        Oq1 = sparse(nQuadCone,quadCone1);
        Oq  = sparse(nQuadCone,nQuadCone);
        Oqn = sparse(nQuadCone,size(A,2));

        %two cones
        Odq = sparse(nExpCone,nQuadCone);
        Oqd = sparse(nQuadCone,nExpCone);

        F = [...
            %  x,   1,  p,   e,  1,    q;
            Odn, Id1, Idp,  Od, Oz1, Odq;  % exp cone    x1  = 1 or y (if normalisation)
            Oqn, Ox1, Oqp, Oqd, Oq1,  Iq;  % quad cone   x1  = q
            Idn, Od1, Odp,  Od, Oz1, Odq;  % exp cone    x2  = x
            Oqn, Ox1, Oqp, Oqd, Iq1,  Oq;  % quad cone   x2  = 1
            Odn, Od1, Odp,  Id, Oz1, Odq;  % exp cone    x3  = e
              R, Ox1, Oqp, Oqd, Oq1,  Oq]; % quad cone R*x3  = F3*x

        %permute the rows of F to form (x1, x2, x3) triples for each cone
        prob.f = sparse(size(F,1),size(F,2));
        prob.f(1:3:(3*nCone),:) = F((1:nCone)',:);
        prob.f(2:3:(3*nCone),:) = F((nCone+1:2*nCone)',:);
        prob.f(3:3:(3*nCone),:) = F((2*nCone+1:3*nCone)',:);

        %g (double[]) – The constant term of affine conic constraints. If not present or g==[] it is assumed g=0
        prob.g = zeros(size(prob.f,1),1);

        if nlt==1 && 0 %TODO implement for general small problem
            %names on all of the variables, used for debugging small
            %problems
            prob.names.var = cell(size(prob.a,2)+size(prob.f,1),1);
            if 1
                prob.names.var{1,1} = 'f';
                prob.names.var{2,1} = 'r';
                prob.names.var{3,1} = 'vA';
                prob.names.var{4,1} = 'vB';
                prob.names.var{5,1} = 's';
                prob.names.var{6,1} = 'tf';
                prob.names.var{7,1} = 'tr';
                prob.names.var{8,1} = 'xf';
                prob.names.var{9,1} = 'xr';
                prob.names.var{10,1} = 'yf';
                prob.names.var{11,1} = 'yr';
                prob.names.var{12,1} = 'mtf';
                prob.names.var{13,1} = 'mtr';
                %rearrange the names of the variables according to ind
                prob.names.var(8:13)=prob.names.var(7+ind);
            else
                for i=1:nlt
                    prob.names.var{i,1} = ['f' int2str(i)];
                    prob.names.var{nlt+i,1} = ['r' int2str(i)];
                    prob.names.var{2*nlt+1+i,1} = ['tf' int2str(i)];
                    prob.names.var{3*nlt+1+i,1} = ['tr' int2str(i)];
                    prob.names.var{4*nlt+1+i,1} = ['xf' int2str(i)];
                    prob.names.var{5*nlt+1+i,1} = ['xr' int2str(i)];
                    prob.names.var{6*nlt+1+i,1} = ['yf' int2str(i)];
                    prob.names.var{7*nlt+1+i,1} = ['yr' int2str(i)];
                    prob.names.var{8*nlt+1+i,1} = ['mtf' int2str(i)];
                    prob.names.var{9*nlt+1+i,1} = ['mtr' int2str(i)];
                end
            end

            %print out the problem to diagnose problems manually
            mosekopt('write(problem.opf)',prob)
        end

        %set default mosek parameters for this type of problem
        paramMosek=mosekParamSetEFBA;

        if ~isfield(solverParams,'MSK_DPAR_INTPNT_CO_TOL_PFEAS')
            if isfield(solverParams,'MSK_DPAR_INTPNT_CO_TOL_PFEAS')
                paramMosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS = solverParams.feasTol;
            else
                paramMosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS = problemTypeParams.feasTol;
            end
        end
        if ~isfield(solverParams,'MSK_DPAR_INTPNT_CO_TOL_DFEAS')
            if isfield(solverParams,'MSK_DPAR_INTPNT_CO_TOL_DFEAS')
                paramMosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS = solverParams.optTol;
            else
                paramMosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS = problemTypeParams.optTol;
            end
        end

        % only set the print level if not already set via solverParams structure
        if ~isfield(solverParams, 'MSK_IPAR_LOG')
            switch problemTypeParams.printLevel
                case 0
                    echolev = 0;
                case 1
                    echolev = 3;
                case 2
                    paramMosek.MSK_IPAR_LOG_INTPNT = 1;
                    paramMosek.MSK_IPAR_LOG_SIM = 1;
                    echolev = 3;
                otherwise
                    echolev = 0;
            end
            if echolev == 0
                paramMosek.MSK_IPAR_LOG = 0;
                cmd = ['minimize echo(' int2str(echolev) ')'];
            else
                cmd = 'minimize';
            end

        end
        %overised if in debug mode
        if problemTypeParams.debug
            cmd = 'minimize';
        end

        if problemTypeParams.debug && 0
            probBeforeMosekopt = prob;
            save('probBeforeMosekopt','probBeforeMosekopt');
        end

        %param = updateStructData(param,solverParams);

        %call mosek exponential cone solver
        tic;
        if 0
            %default
            [~,res]=mosekopt('minimize',prob);
        else
            [~,res]=mosekopt(cmd,prob,paramMosek);
        end
        solution.time = toc;

        %parse mosek result structure
        [stat,origStat,x,y,z,s,doty] = parseMskResult(res,prob.a,prob.blc,prob.buc,problemTypeParams.printLevel,paramMosek);

        solution.stat = stat;
        solution.origStat = origStat;
        switch stat
            case 1
                %check for zeros in variables within entropy functions
                zeroxBool = x(1:length(d))==0 & d~=0;
                if any(zeroxBool)
                    warning([num2str(nnz(zeroxBool)) ' optimal values that equal zero within entropy objective(s)'])
                    ind = find(zeroxBool);
                    fprintf('%8s %8s %8s\n','xl','x','xu')
                    for i=1:length(ind)
                        fprintf('%8.4g %8.4g %8.4g\n',prob.blx(ind(i)),x(ind(i)),prob.bux(ind(i)));
                    end
                end

                if problemTypeParams.printLevel > 1
                    % Problem definition here: https://docs.mosek.com/9.2/toolbox/prob-def-affine-conic.html
                    fprintf('%s\n','Optimality conditions (numerical)')
                    % Guide to interpreting the solution summary: https://docs.mosek.com/9.2/toolbox/debugging-log.html#continuous-problem
                    fprintf('%8.2g %s\n',norm(prob.a(prob.blc==prob.buc,:)*x - prob.blc(prob.blc==prob.buc),inf), '|| A*x - b ||_inf');
                    val = norm(prob.c - prob.a'*y - z - prob.f'*doty,inf);
                    fprintf('%8.2g %s\n',val, '|| c - A''*y - z - F''*doty ||_inf');
                    if val>1e-6 || problemTypeParams.debug
                        solution.T0 = table(prob.c - prob.a'*y - z - prob.f'*doty,prob.c, prob.a'*y, z,prob.f'*doty,'VariableNames',{'tot','c','Aty','z','Ftdoty'});
                    end
                    %fprintf('%8.2g %s\n',norm(prob.c - prob.f'*s,inf), '|| c - F''s ||_inf');
                    fprintf('%8.2g %s\n',norm(-y + res.sol.itr.slc - res.sol.itr.suc,inf), '|| -y + res.sol.itr.slc - res.sol.itr.suc ||_inf');
                    %fprintf('%8.2g %s\n',prob.c'*x - prob.b'*y, ' c''*x -b''*y');

                    fprintf('%8.2g %s\n',(prob.f*x + prob.g)'*doty, '(F*x + g)''*s >= 0');
                end

                %%% Reorder
                % Dual variables to affine conic constraints, based on original order of rows in F matrix
                y_K = zeros(length(doty),1);
                y_K(1:nCone,1) = doty(1:3:3*nCone);
                y_K(nCone+1:2*nCone,1) = doty(2:3:3*nCone);
                y_K(2*nCone+1:3*nCone,1) = doty(3:3:3*nCone);

                %check with the original order of the affine cone constraints
                val = norm(prob.c - prob.a'*y - z - F'*y_K,inf);
                if problemTypeParams.printLevel > 1
                    fprintf('%8.2g %s\n',val, '|| c - A''*y - z - F''*y_K ||_inf');
                end
                if val>1e-6 || problemTypeParams.debug
                    solution.T = table(prob.c - prob.a'*y - z - F'*y_K,prob.c, prob.a'*y, z,prob.f'*doty,F'*y_K,'VariableNames',{'tot','c','Aty','z','Ftdoty','Fty_K'});
                end

                if problemTypeParams.printLevel > 1
                    x1 = F(1:nCone,:)*x;
                    x2 = F(nCone+1:2*nCone,:)*x;
                    x3 = F(2*nCone+1:3*nCone,:)*x;

                    if nExpCone>0
                        fprintf('\n%s\n','Primal exponential cone:')
                        fprintf('%8.2g %s\n',min(x1(1:nExpCone) - x2(1:nExpCone).*exp(x3(1:nExpCone)./x2(1:nExpCone))), 'min(x1 - x2*exp(x3/x2)) >= 0');
                        %https://docs.mosek.com/modeling-cookbook/expo.html#entropy
                        fprintf('%8.2g %s\n',max(x3(1:nExpCone) - x2(1:nExpCone).*log(x1(1:nExpCone)./x2(1:nExpCone))), 'max(x3 - x2*log(x1/x2)) <= 0');
                        fprintf('%8.2g %s\n',norm(x3(1:nExpCone) - x2(1:nExpCone).*log(x1(1:nExpCone)./x2(1:nExpCone)),inf), '|| x3 - x2*log(x1/x2) ||_inf for exp cones');

                        if nExpCone<=5 && isfield(prob,'names') && 0
                            %TODO complete for general input
                            for i=1:nExpCone %exp cones first
                                fprintf('%7.2g\t%s\n',norm(x3(i) - x2(i).*log(x1(i)/x2(i)),inf),...
                                    ['|| ' prob.names.var{2*n+1+i} ' - ' prob.names.var{i} '.* log(' prob.names.var{i} ' / ' prob.names.var{2*n+1} ') ||_inf']);
                            end
                        end
                    end

                    if nQuadCone>0
                        fprintf('\n%s\n','Primal quadratic cone:')
                        % (1/2)*x'*(R'*R)*x <= q
                        % x1 =   q
                        % x2 =   1
                        % x3 = R*x
                        fprintf('%8.2g %s\n',min(x1 - (1/2)*(x3'*x3)), 'min(x1 - (1/2)*x3''*x3)) >= 0');
                    end

                    y1_K = y_K(1:nCone);          % 1 should  be non-negative
                    y2_K = y_K(nCone+1:2*nCone);  % x
                    y3_K = y_K(2*nCone+1:3*nCone);% e should  be non-positive

                    fprintf('\n%s\n','Dual exponential cone:')
                    % https://docs.mosek.com/9.2/toolbox/prob-def-affine-conic.html
                    % This is moseks convention to the dual exponential cone
                    fprintf('%7.2g\t%s\n',min(y1_K(1:nExpCone) + y3_K(1:nExpCone).*exp(y2_K(1:nExpCone)./y3_K(1:nExpCone))/exp(1)), 'min(y1_k + y3_k.*exp(y2_K./y3_K)/exp(1))  >= 0');
                end

                solution.full = x(1:size(A,2));
                %switch to Rockafellar signs
                solution.dual = -y(1:size(A,1));
                solution.dualNorm = -y(size(A,1)+1:size(A,1)+p);
                solution.rcost = -z(1:size(A,2)+p);
                solution.slack = s;

                %need to zero out the NaN due to log(0) for some variables
                logSolutionFull = real(log(solution.full));
                logSolutionFull(~isfinite(logSolutionFull))=0;
                if isfield(EPproblem,'Q')
                    solution.obj = EPproblem.c'*solution.full + (EPproblem.d.*solution.full)'*(logSolutionFull -1) + (1/2)*solution.full'*EPproblem.Q*solution.full;
                    solution.objLinear = EPproblem.c'*solution.full;
                    solution.objEntropy = -(EPproblem.d.*solution.full)'*(logSolutionFull -1);
                    solution.objQuadratic = (1/2)*solution.full'*EPproblem.Q*solution.full;
                else
                    solution.obj = EPproblem.c'*solution.full + (EPproblem.d.*solution.full)'*(logSolutionFull -1);
                    solution.objLinear = EPproblem.c'*solution.full;
                    solution.objEntropy = -(EPproblem.d.*solution.full)'*(logSolutionFull -1);
                    solution.objQuadratic = 0;

                end

                posRcost = solution.rcost>0;
                negRcost = solution.rcost<0;
                blx = prob.blx(1:size(A,2));
                bux = prob.bux(1:size(A,2));
                solution.lagRcost = sum(solution.rcost(negRcost)'*blx(negRcost) + solution.rcost(posRcost)'*bux(posRcost));

                %pass back the F matrix to check biochemical optimality criteria
                solution.coneF = F;
                solution.auxPrimal = x(size(A,2)+p+1:end);
                solution.auxRcost = -z(size(A,2)+p+1:end);
                solution.coneDual = -y_K;

                % variable to determine the residual 1
                b = prob.blc;
                % variable to determine the residual 2
                grad = prob.c - F'*y_K;
                grad = grad(1:size(A,2)+p,1);
                Aty = -prob.a'*y;
                Aty = Aty(1:size(A,2)+p,1);%strictly this is [A;P]'*y

            otherwise

                doty = NaN*ones(size(prob.f,1),1);
        end
    otherwise
        error([problemTypeParams.solver ' is an unrecognised solver'])
end

switch solution.stat
    case 0
        switch problemTypeParams.solver
            case 'pdco'
                %infeasible, debug the situtation
                disp(solution.origStat)
                %solution.origStat: 'PRIMAL_INFEASIBLE_CER'
                solutionLP = solveCobraLP(EPproblem);
            case 'mosek'
                %https://docs.mosek.com/8.1/toolbox/solving-linear.html
                if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_PFEAS')
                    solverParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol;
                end
                if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.')
                    solverParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol;
                end
                %If the feasibility tolerance is changed by the solverParams
                %struct, this needs to be forwarded to the cobra Params for the
                %final consistency test!
                if isfield(problemTypeParams,'MSK_DPAR_INTPNT_TOL_PFEAS')
                    solverParams.feasTol = solverParams.MSK_DPAR_INTPNT_TOL_PFEAS;
                end

                % only set the print level if not already set via solverParams structure
                if ~isfield(solverParams, 'MSK_IPAR_LOG')
                    switch problemTypeParams.printLevel
                        case 0
                            echolev = 0;
                        case 1
                            echolev = 3;
                        case 2
                            solverParams.MSK_IPAR_LOG_INTPNT = 1;
                            solverParams.MSK_IPAR_LOG_SIM = 1;
                            echolev = 3;
                        otherwise
                            echolev = 0;
                    end
                end
                if echolev == 0
                    solverParams.MSK_IPAR_LOG = 0;
                    cmd = ['minimize echo(' int2str(echolev) ')'];
                else
                    cmd = 'minimize';
                end
                [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,solverParams,cmd);

                [statLP,origStat,x,y,z,s,doty] = parseMskResult(res,EPproblem.A,EPproblem.blc,EPproblem.buc,problemTypeParams.printLevel);

        end

        switch statLP
            case 1
                message =['solveCobraEP: EPproblem with ' problemTypeParams.solver ' is infeasible, but corresponding LPproblem is feasible according to solveCobraLP with ' problemTypeParams.solver];
                warning(message)
            otherwise
                message = ['solveCobraEP: EPproblem with ' problemTypeParams.solver ' is infeasible, because corresponding LPproblem is infeasible according to solveCobraLP with ' problemTypeParams.solver];
                warning(message)
        end
        if exist('messages','var')
            if isfield(solution,'messages')
                solution.messages = [messages;solution.messages;message];
            else
                solution.messages = [messages;message];
            end
        else
            solution.messages = cellstr(message);
        end

    case 1
        % check the optimality conditions for various solvers
        if ~isempty(solution.slack) && ~isempty(solution.full)
            % determine the residual 1
            switch problemTypeParams.solver
                case 'pdco'
                    feasTol = 1e-3;
                    res1 = A*solution.full + solution.slack - b;
                    res1(~isfinite(res1))=0;
                case 'mosek'
                    feasTol = problemTypeParams.feasTol * 1e2;
                    res1 = A(blc==buc,:)*solution.full - blc(blc==buc);
            end
            tmp1 = norm(res1, inf);

            % evaluate the optimality condition 1
            if tmp1 > feasTol
                if strcmp(problemTypeParams.solver,'pdco')
                    res1b = norm(A*solution.full + solution.slack - b + (d2^2)*y,inf);
                    tmp1b = norm(res1b, inf);
                    if tmp1b > feasTol
                        displayError = 1;
                    else
                        displayError = 0;
                        warning(['[' problemTypeParams.solver '] Primal optimality condition in solveCobraEP only approximately satisfied, residual = ' num2str(tmp1) ', regularised residual = ' num2str(tmp1b) ', while problem feasTol = ' num2str(feasTol) '.  origStat = ' solution.origStat])
                    end
                else
                    %TODO - debug why solver reporting optimal but unscaled seems less so.
                    displayError = 0;
                end
                if displayError
                    %disp(solution.origStat)
                    fprintf('%s\n',['[' problemTypeParams.solver '] Primal optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp1) ', while problem feasTol = ' num2str(feasTol) '.  origStat = ' solution.origStat])
                end
            else
                if problemTypeParams.printLevel > 0
                    fprintf(['\n > [' problemTypeParams.solver '] Primal optimality condition in solveCobraEP satisfied.']);
                end
            end
        end

        %gradient may differ depending on the solver
        res2 = grad  + Aty + solution.rcost;
        tmp2 = norm(res2, inf);

        if 0
            optTol = problemTypeParams.optTol * 1e2;
        else
            optTol = 5e-5;
        end
        % evaluate the optimality condition 2
        if tmp2 > optTol
            disp(solution.origStat)
            if ~(length(A)==1 && strcmp(problemTypeParams.solver,'pdco')) %todo, why does pdco choke on small A?
                warning(['[' problemTypeParams.solver '] Dual   optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp2) ', while problem optTol = ' num2str(optTol)])
            end
        else
            if problemTypeParams.printLevel > 0
                fprintf(['\n > [' problemTypeParams.solver '] Dual   optimality condition in solveCobraEP satisfied.\n']);
            end
        end
end

solution.solver=problemTypeParams.solver;

end

function [obj,grad,hess] = entropyObj(x,c,d)

n=length(x);
logx = zeros(n,1);
e = double(d~=0);
logx(d~=0) = reallog(x(d~=0));     % error if negative
obj  = c'*x + (d.*x)'*(logx - e);
grad = c + d.*logx;
hess = d./x;
hess(d==0)=0;
%hess = diag(hess); % not necessary as pdco knows to treat vector as a
%diagonal of a hessian.
end

function [obj,grad,hess] = entropyQPObj(x,c,d,Q)

n=length(x);
logx = zeros(n,1);
e = double(d~=0);
logx(d~=0) = reallog(x(d~=0));     % error if negative
obj  = c'*x + (d.*x)'*(logx - e) + 1/2*x'*Q*x;
grad = c + d.*logx + Q*x;
hess = d./x;
hess(d==0)=0;
hess = spdiags(hess,0,n,n) + Q;
end


text:
% dim_red_type: dimension reduction type: one of the strings: '', 'pca', 'ica'
% new_dim: relevant only if dim_red_type is not ''
function lmm(numClusters, is_sampled, dim_red_type, new_dim, numOfEMIter, numOfCpus)

    fv_init;

    if ~exist('numOfEMIter', 'var')
        numOfEMIter = 10;
    end

    if ~exist('numOfCpus', 'var')
        numOfCpus = 4;
    end

    word2vec_sqrt = false;

    mult_factor = LMM_HGLMM_MULT_FACTOR;

    output(1, 'LMM_HGLMM_MULT_FACTOR = %d\n', LMM_HGLMM_MULT_FACTOR);

    vectors_file_name = add_data_dir_base('GoogleNews_vectors_norm.mat');
    output_file_name = add_data_dir_base(sprintf('GoogleNews_norm_lmm_%d.mat', numClusters));

    if word2vec_sqrt
        vectors_file_name = fname_concat(vectors_file_name, '_sqrt');
        output_file_name = fname_concat(output_file_name, '_sqrt');
    end

    if is_sampled
        vectors_file_name = fname_concat(vectors_file_name, '_sampled');
        output_file_name = fname_concat(output_file_name, '_sampled');
    end

    if strcmp(dim_red_type,'pca') || strcmp(dim_red_type,'ica')
        post_fix = sprintf('_%s_%d', dim_red_type, new_dim);
        vectors_file_name = fname_concat(vectors_file_name, post_fix);
        output_file_name = fname_concat(output_file_name, post_fix);
    end

    output(1, 'reading vectors file %s\n', vectors_file_name);
    load(vectors_file_name);
    output(1, 'done\n');

    output(1, 'Calculating LMM, %d clusters\n', numClusters);

    % multiply by LMM_HGLMM_MULT_FACTOR=10 due to numerical issue
    % transpose because the rows should be the samples
    vectors_ = mult_factor * vectors_';

    % create a random initialization
    initPriors = (1.0 / numClusters) * ones(numClusters,1);
    muInit=vectors_(randsample(size(vectors_,1),numClusters),:);
    init_b=repmat(sum(abs(vectors_ - repmat(median(vectors_),size(vectors_,1),1))) / size(vectors_,1),numClusters,1);

    % calculate LMM

    output(1, 'numOfCpus = %d\n', numOfCpus);
    output(1, 'numOfEMIter = %d\n', numOfEMIter);

    [n_means,n_covariances,n_priors,samplesWeightsOut]=LMMMex(double(vectors_), double(muInit), double(init_b), double(initPriors), numOfEMIter, 0.01, numOfCpus);

    save(output_file_name, 'n_means', 'n_covariances', 'n_priors', 'mult_factor', '-v7.3');

    output(1, 'saved to file: %s\n', output_file_name);
end

author: layumi

repo: Image-Text-Embedding


text:
function options = reg_fine_rigid(rgb1, I1, wl, options)
%REG_FINE_RIGID Summary of this function goes here
%   Detailed explanation goes here

t_start = tic;
% optimize the objective function by block coordinate descent
options = reg_fine_rigid_internal(rgb1, I1, wl, options);

disp(['Elapsed time for fine-scale rigid registration is ',num2str(toc(t_start))]);

function options = reg_fine_rigid_internal(rgb1, I1, wl, options)
I = double(rgb1) / 255;

last_pos = NaN(1,6);

for iter = 1:20
    try
        t_start_iter = tic;

        % optimize w.r.t. degree, t, s
        [degree, t, s] = optimize_wrt_T_s(I, I1, wl, options);
        options.degree = degree;
        options.t = t;
        options.s = s;

        % optimize w.r.t. lambda
        sigma = optimize_wrt_sigma(I, I1, options);
        options.sigma = sigma;

        % check convergence
        if all([degree,t,s,sigma] == last_pos)
            disp(['Iteration ends at ',num2str(iter)]);
            break;
        else
            last_pos = [degree,t,s,sigma];
        end
        disp(['Iteration ',num2str(iter),' lasts ',num2str(toc(t_start_iter))]);
    catch err
        if strcmp(err.identifier, 'Registration:DerivativeBoundary')
            disp(['Stop optimization since boundary touch when ', ...
                'calculating derivatives']);
            break;
        else
            throw(err);
        end
    end
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function f2 = p04_f2 ( x )

%*****************************************************************************80
%
%% P04_F2 evaluates the second derivative for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    LE Scales,
%    Introduction to Non-Linear Optimization,
%    Springer, 1985.
%
%  Parameters:
%
%    Input, real X, the values of the variables.
%
%    Output, real F2, the second derivative.
%
  f2 = exp ( x ) + 0.02 / x / x / x;

  return
end