text:
classdef PREA < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Promising-region based EMO algorithm

%------------------------------- Reference --------------------------------
% J. Yuan, H. Liu, F. Gu, Q. Zhang, and Z. He, Investigating the properties
% of indicators and an evolutionary many-objective algorithm based on a
% promising region, IEEE Transactions on Evolutionary Computation, 2021,
% 25(1): 75-86.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Jiawei Yuan

    methods
        function main(Algorithm,Problem)
            %% Generate the random population
            Population = Problem.Initialization();
            Zmin       = min(Population.objs,[],1);

            %% shift the objective space to R+
            PopObj = Population.objs;
            PopObj = PopObj - repmat(Zmin,Problem.N,1) + 1e-6;

            %% calculate the ratio based indicator matrix
            IMatrix = ones(Problem.N,Problem.N);
            for i=1:1:Problem.N
                Fi             = PopObj(i,:);
                % calculate ratio based indicator value of each individual
                Ir             = PopObj./repmat(Fi,Problem.N,1) - 1;
                InvertIr       = repmat(Fi,Problem.N,1)./PopObj - 1;
                MaxIr          = max(Ir,[],2);
                MinIr          = max(InvertIr,[],2);
                DomInds        = find(MaxIr<=0);
                MaxIr(DomInds) = -MinIr(DomInds);
                IMatrix(i,:)   = MaxIr';
                IMatrix(i,i)   = Inf;
            end

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = MatingStrategy(IMatrix);
                Offspring  = OperatorGAhalf(Problem,Population(MatingPool));
                Zmin       = min([Zmin;Offspring.objs],[],1);
                [Population,IMatrix] = PREA_Update([Population,Offspring],Problem.N,Zmin);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function g =  lfmGradientSigmaH3AV(gamma1_p, gamma1_m, sigma2, t1, ...
    t2, preFactor, mode)

% LFMGRADIENTSIGMAH3AV Gradient of the function h_i(z) with respect \sigma.
% FORMAT
% DESC Computes the gradient of the function h_i(z) with respect to the
% length-scale of the input "force", \sigma.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN g : Gradient of the function with respect to \sigma.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

g = preFactor(1)*lfmavGradientSigmaUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode) ...
    + preFactor(2)*lfmavGradientSigmaUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%kckmeans2 'Performs K-Means clustering with more options '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ckmeans2.pane file
%
% Parameters:
% InputFile: i1 'Input data object', required: 'input data object (to be clustered)'
% InputFile: i2 'Cluster center input object', optional: 'cluster center input object'
% Toggle: map 'Generate output map', default: 0: 'generate output map'
% Toggle: spectrum 'SPECTRUM compatable map segment', default: 0: 'SPECTRUM compatable map segment'
% Integer: n 'Max number of iterations', default: 50000: 'max number of iterations'
% Integer: mergeifcount 'Cluster count is <', default: 1: 'merge with closest if the cluster has less than mergeifcount points'
% Integer: maskifcount 'Mask points in clusters with counts < ', default: 1: 'mask if point is in cluster with less than maskifcount'
% OutputFile: o1 'Cluster number output object', required: 'cluster number output object'
% OutputFile: o2 'Cluster center output object', optional: 'cluster center output object'
% OutputFile: o3 'Cluster variance output object', optional: 'cluster variance output object'
% OutputFile: o4 'Cluster membership count output', optional: 'Cluster membership count output'
% OutputFile: o5 'K-means statistics output (ASCII)', optional: 'K-means statistics output (ASCII)'
% Integer: k 'Initial number of clusters', default: 2: 'initial number of clusters'
%
% Example: [o1, o2, o3, o4, o5] = kckmeans2({i1, i2}, {'i1','';'i2','';'map',0;'spectrum',0;'n',50000;'mergeifcount',1;'maskifcount',1;'o1','';'o2','';'o3','';'o4','';'o5','';'k',2})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ckmeans2 - Performs K-Means clustering with more options
%
%  DESCRIPTION
% "ckmeans2" is the implementation of a clustering algorithm that can be used for unsupervised classification of data. \fIckmeans2\fP has some enhancements over the \fIkkmeans\fP algorithm of the datamanip toolbox. These enhancements provides more variation of parameters for clustering but does not guarantee that better results could be obtained.
% Main differences from "kkmeans" to \fIckmeans2\fP are:
% 1 - "ckmeans2" consider mask data. Masked points will not be taken in account for clustering and/or calculation of means and variances.
% 2 - If the user don't provide an initial cluster center object with the [-i2] parameter, "ckmeans2" will create the initial cluster centers from random sparse values - there is no guarantee that this will improve the convergence (in some cases may even slow down the clustering) but if you use the first N points from the data there is the risk that these points will be similar (e.g. with image data)
% 3 - "ckmeans2" allows merging of clusters based on the number of points in the clusters - this allow clusters with few points  (specified with the parameter [-mergeifcount]) to be joined to the nearest cluster. Clusters can also be joined if their distances is smaller than a threshold ([-mergeifdistance]). Cluster distances are calculated with the Euclidean distance of their centers. Please note that the joining algorithm is not optimum - joining is perfomed after kmeans clustering to avoid conflicts while clustering (slowing convergence), and is not done recursively until a stable state is achieved.
% 4 - "ckmeans2" allows the user to mask points belonging to clusters with few points (specified with the parameter [-maskifcount]). This has the effect of rejecting clusters with few points that can be considered as outliers on the data. A caveat is needed: the kmeans algorithm tries to create similar-sized, spherical-shaped clusters, but few points in a cluster does not necessarily means that this cluster is made of outliers.
% This implementation of "ckmeans2" uses the same basic clustering procedure as the \fIkkmeans\fP implementation in Khoros - meaning that the code was copied partially, specially for creation of that complicated SPECTRUM compatibility thingie. Please refer to the \fIkkmeans\fP help for more information about \fIkkmeans\fP.
% The input file is specified with the parameter [-i1], and optionally initial cluster centers can be passed with the [-i2] parameter. If passed, this data must have dimensions Cx1x1x1xE, where C is the number of classes and E is the number of elements in the [-i1] input data.
% The number of clusters can be specified either with the [-k] parameter or by proportion [-divk] of the number of valid points on the data. The maximum number of iterations is specified with the [-n] parameter.
% The output file [-o1] can have a map segment that will be the center of the cluster for each cluster, if the parameter [-map] is passed. Optionally a Spectrum-compatable map can be created by specifying [-spectrum].
% As with "kkmeans", this routine can create several output files: [-o1] will contain the cluster numbers for each input vector, [-o2] will contain the cluster centers, [-o3] will contain the cluster variances, [-o4] will contain the cluster counts and [-o5]
% will contain information about the clustering (in ASCII). Please refer to the kman page for "kkmeans" for more details about these files.
% This routine is a first step towards a more complete implementation of the Isodata algorithm, which uses the k-means clustering algorithm for clustering and heuristics for post-clustering merging and splitting of clusters.
%
%
%
%  EXAMPLES
% All examples for the Classify toolbox are listed on the Classify Toolbox Manual. For an example of this kroutine, please see the example workspace Classify:workspaces:kmeans
%
%  "SEE ALSO"
% kkmeans (in the Datamanip toolbox), cfuzzycmeans
%
%  RESTRICTIONS
% Complex data types are not supported. If there is mask, expect the mask to be consistent over vectors, i.e. for a vector's elements either all masks are TRUE or FALSE.
% At this implementation the final result will contain the same number of clusters as specified by [-k] or [-divk], even if joining is performed. Future versions should delete the clusters with zero membership.
%
%  REFERENCES
% All references for the Classify toolbox are listed on the Classify Toolbox Manual.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function value = scasum ( n, x, incx )

%*****************************************************************************80
%
%% SCASUM takes the sum of the absolute values of a complex vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for FORTRAN usage,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, pages 308-323, 1979.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, complex X(*), the vector.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Output, real VALUE, the sum of the absolute values.
%
  value = sum ( abs ( real ( x(1:incx:1+(n-1)*incx) ) ) ...
              + abs ( imag ( x(1:incx:1+(n-1)*incx) ) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function[x,y]=ellcurves(varargin)
%ELLCURVES  Returns curves corresponding to specified ellipse properties.
%
%   For direct plotting of ellipses, see ELLPLOT, which calls ELLCURVES.
%
%   ZC=ELLCURVES(KAPPA,LAMBDA,THETA,ZRES) returns complex-valued curves ZC
%   tracing out the periphery of ellipses with amplitude KAPPA, linearity
%   LAMBDA, and orientation THETA, located at complex positions ZRES.
%
%   All input arguments are arrays of the same length, say N.  By default,
%   ZC is calculated at 32 locations around the periphery.  Thus ZC will be
%   a complex-valued matrix with 32 rows and N columns.
%
%   [XC,YC]=ELLCURVES(...) alternately returns the real and imaginary parts
%   of ZC, corresponding to the X- and Y- components of the ellipses.
%
%   ELLCURVES(KAPPA,LAMBDA,THETA,PHI,ZRES) with five input arguments begins
%   each ellipse at phase PHI.  For most applications, the starting phase
%   does not matter.  The default value is to begin with a phase of zero.
%
%   ELLCURVES(... ,'npoints',M) uses M points around the ellipse periphery.
%
%   ELLCURVES(... ,'aspect',AR) with AR=[XAR YAR] multiplies the X-signal
%   and Y-signal by XAR and YAR, respectively, for plotting purposes.
%
%   ELLCURVES is called by ELLPLOT, and is also useful in applications,
%   e.g. analysis of model fields within elliptical contours.
%
%   The curves can be plotted with PLOT(ZC).
%
%   See also ELLIPSEPLOT, ELLSIG, INELLIPSE.
%
%   Usage: zc=ellcurves(kappa,lambda,theta,zres);
%          zc=ellcurves(kappa,lambda,theta,phi,zres);
%          zc=ellcurves(kappa,lambda,theta,phi,zres,'npoints',64);
%          [xc,yc]=ellcurves(kappa,lambda,theta,phi,zres,'npoints',64);
%          zc=ellcurves(kappa,lambda,theta,phi,zres,'aspect',ar,'npoints',64);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014--2018 J.M. Lilly --- type 'help jlab_license' for details

na=length(varargin);

k=varargin{1};
l=varargin{2};
theta=varargin{3};

ar=[1 1];
npoints=32;

if na>4
    for i=1:3
        if ischar(varargin{end-1})
            if strcmpi(varargin{end-1}(1:3),'npo')
                npoints=varargin{end};
            elseif strcmpi(varargin{end-1}(1:3),'asp')
                ar=varargin{end};
            end
            na=na-2;
            varargin=varargin(1:end-2);
        end
    end
end

if length(varargin)==5
    phi=varargin{4};
    x=varargin{5};
else
    x=varargin{4};
    if iscell(k)
        for i=1:length(k)
             phi{i}=zeros(size(k{i}));
        end
    else
        phi=zeros(size(k));
    end
end

if ~iscell(k)
    z=ellcurves_one(k,l,theta,phi,x,ar,npoints);
else
    for i=1:length(k)
        z{i,1}=ellcurves_one(k{i},l{i},theta{i},phi{i},x{i},ar,npoints);
    end
end

if nargout==2
    if iscell(z)
        x=cellreal(z);
        y=cellimag(z);
    else
        x=real(z);
        y=imag(z);
    end
else
    x=z;
end

function[z]=ellcurves_one(kappa,lambda,theta,phi,x,ar,npoints)
phio=linspace(0,2*pi,npoints);
%x=conj(x');

z=zeros(length(phio),length(kappa));
for i=1:length(phio)  %Loop over phio, very fast
    [xc,yc]=ellsig(kappa,lambda,theta,phi+phio(i),'real');
    %vsize(xc,yc,x)
    z(i,:)=ar(1)*xc+sqrt(-1)*ar(2)*yc+x;
end
%vsize(xc,yc,x)

%if ~isempty(strfind(str,'not'))
%    z=permute(z,[2 1]);  %This is like 6 times faster than conjugate transpose
%end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [] = generateMetaboReport(modelProperties,reportDir,orgNames)
% This function generates
F = fieldnames(modelProperties);

for i = 1 : length(F)
    clear html
    fid=fopen('reportTemplate.html');
    tline = fgetl(fid);
    cnt =1;
    while ischar(tline)
        disp(tline)
        tline = fgetl(fid);
        disp(tline)
        % replace space holders in template html file
        if ~isnumeric(tline)
            if ~exist('orgNames','var')
                ModelName = regexprep(F{i},'_', ' ');
            else
                if strcmp(F{i},orgNames{i})
                    ModelName = regexprep(F{i},'_', ' ');
                else
                    ModelName = regexprep(orgNames{i},'_', ' ');
                end
            end
            tline= regexprep(tline,'###ModelName',ModelName);
            % Basic information
            tline= regexprep(tline,'###OverallScore',num2str(round(modelProperties.(F{i}).modelProp2.Scores.Overall,2)));
            tline= regexprep(tline,'###NumMets',num2str(modelProperties.(F{i}).modelProp2.m));
            tag = '###ListNumMets'; list = modelProperties.(F{i}).modelProp2.Details.metabolites; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRxns',num2str(modelProperties.(F{i}).modelProp2.n));
            tag = '###ListNumRxns'; list = modelProperties.(F{i}).modelProp2.Details.reactions; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);
            tline= regexprep(tline,'###NumGenes',num2str(modelProperties.(F{i}).modelProp2.genes));
            tline= regexprep(tline,'###NumComp',num2str(modelProperties.(F{i}).modelProp2.compartments));
            tline= regexprep(tline,'###MetCov',num2str(round(modelProperties.(F{i}).modelProp2.metCov,2)));
            % Metabolic information
            tline= regexprep(tline,'###NumMUnique',num2str(modelProperties.(F{i}).modelProp2.metUnique));
            tag = '###ListMetUnique'; list = modelProperties.(F{i}).modelProp2.Details.metabolites_unique; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetNoCharge',num2str(modelProperties.(F{i}).modelProp2.MissingMetCharge));

            tline= regexprep(tline,'###NumMetNoForm',num2str(modelProperties.(F{i}).modelProp2.MissingMetFormulae));
            tag = '###ListNumMetNoForm'; list = modelProperties.(F{i}).modelProp2.Details.MissingMetFormulae; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###MedComp',num2str(modelProperties.(F{i}).modelProp2.MediumMets));
            tag = '###ListMedComp'; list = modelProperties.(F{i}).modelProp2.Details.MediumMets; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            % Reaction information
            tline= regexprep(tline,'###NumMetRxns',num2str(modelProperties.(F{i}).modelProp2.MetabolicRxns));
            tag = '###ListNumMetRxns'; list = modelProperties.(F{i}).modelProp2.Details.MetabolicRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumTransRxns',num2str(modelProperties.(F{i}).modelProp2.TransportRxns));
            tag = '###ListNumTransRxns'; list = modelProperties.(F{i}).modelProp2.Details.TransportRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumExRxns',num2str(modelProperties.(F{i}).modelProp2.ExchangeRxns));
            tag = '###ListNumExRxns'; list = modelProperties.(F{i}).modelProp2.Details.ExchangeRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumDmRxns',num2str(modelProperties.(F{i}).modelProp2.DemandRxns));
            tag = '###ListNumDmRxns'; list = modelProperties.(F{i}).modelProp2.Details.DemandRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumSinkRxns',num2str(modelProperties.(F{i}).modelProp2.SinkRxns));
            tag = '###ListNumSinkRxns'; list = modelProperties.(F{i}).modelProp2.Details.SinkRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumBmRxns',num2str(modelProperties.(F{i}).modelProp2.BiomassRxns));
            tag = '###ListNumBmRxns'; list = modelProperties.(F{i}).modelProp2.Details.BiomassRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            % GPR information
            tline= regexprep(tline,'###NumRxnWoGpr',num2str(round(modelProperties.(F{i}).modelProp2.RxnsWithoutGpr,2)));
            tag = '###ListNumRxnWoGpr'; list = modelProperties.(F{i}).modelProp2.Details.RxnsWithoutGpr; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumTRxnWoGpr',num2str(round(modelProperties.(F{i}).modelProp2.TRxnsWithoutGpr,2)));
            tag = '###ListNumTRxnWoGpr'; list = modelProperties.(F{i}).modelProp2.Details.TRxnsWithoutGpr; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            % Network topology
            tline= regexprep(tline,'###NumBlockedR',num2str(round(modelProperties.(F{i}).modelProp2.BlockedRxns,2)));
            tag = '###ListNumBlockedR'; list = modelProperties.(F{i}).modelProp2.Details.BlockedRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumDeadM',num2str(round(modelProperties.(F{i}).modelProp2.DeadendsMets,2)));
            tag = '###ListNumDeadM'; list = modelProperties.(F{i}).modelProp2.Details.DeadendsMets; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumStoichCycleRxns',num2str(round(modelProperties.(F{i}).modelProp2.StoichCycleRxns,2)));
            tag = '###ListNumStoichCycleRxns'; list = modelProperties.(F{i}).modelProp2.Details.StoichCycleRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            % Matrix Conditioning
            tline= regexprep(tline,'###NumMinMaxCoeff',num2str(round(modelProperties.(F{i}).modelProp2.maxminCoeff,2)));
            tag = '###ListNumMinMaxCoeff'; list = modelProperties.(F{i}).modelProp2.Details.maxminCoeff; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRank',num2str(round(modelProperties.(F{i}).modelProp2.Rank,2)));

            % Consistency
            tline= regexprep(tline,'###NumConsRxns',num2str(round(modelProperties.(F{i}).modelProp2.ConsRxns,2)));
            tag = '###ListNumConsRxns'; list = modelProperties.(F{i}).modelProp2.Details.InconsRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumBalRxns',num2str(round(modelProperties.(F{i}).modelProp2.BalancedMassRxns,2)));
            tag = '###ListNumBalRxns'; list = modelProperties.(F{i}).modelProp2.Details.UnBalancedMassRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumChargeRxns',num2str(round(modelProperties.(F{i}).modelProp2.BalancedChargeRxns,2)));
            tag = '###ListNumChargeRxns'; list = modelProperties.(F{i}).modelProp2.Details.UnBalancedChargeRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetConnec',num2str(round(modelProperties.(F{i}).modelProp2.MetConn,2)));
            tag = '###ListNumMetConnec'; list = modelProperties.(F{i}).modelProp2.Details.MetConn; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumUnboundR',num2str(round(modelProperties.(F{i}).modelProp2.UnboundedFlux,2)));
            tag = '###ListNumUnboundR'; list = modelProperties.(F{i}).modelProp2.Details.UnboundedFluxRxns; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###ScoreConsistency',num2str(round(modelProperties.(F{i}).modelProp2.Scores.Consistency,2)));

            % Annotation - Metabolites
            tline= regexprep(tline,'###NumMAnnoAny',num2str(round(modelProperties.(F{i}).modelProp2.MetWAnno,2)));
            tag = '###ListNumMAnnoAny'; list = modelProperties.(F{i}).modelProp2.Details.metWOAnno; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetPubCh',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetPubChemID,2)));
            tag = '###ListNumMetPubCh'; list = modelProperties.(F{i}).modelProp2.Details.missingmetPubChemID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetKegg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetKEGGID,2)));
            tag = '###ListNumMetKegg'; list = modelProperties.(F{i}).modelProp2.Details.missingmetKEGGID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetSeed',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetSEEDID,2)));
            tag = '###ListNumMetSeed'; list = modelProperties.(F{i}).modelProp2.Details.missingmetSEEDID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetIKey',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetInchiKey,2)));
            tag = '###ListNumMetIKey'; list = modelProperties.(F{i}).modelProp2.Details.missingmetInchiKey; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetIString',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetInchiString,2)));
            tag = '###ListNumMetIString'; list = modelProperties.(F{i}).modelProp2.Details.missingmetInchiString; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetChebi',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetChEBIID,2)));
            tag = '###ListNumMetChebi'; list = modelProperties.(F{i}).modelProp2.Details.missingmetChEBIID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetHmdb',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetHMDBID,2)));
            tag = '###ListNumMetHmdb'; list = modelProperties.(F{i}).modelProp2.Details.missingmetHMDBID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetReac',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetReactomeID,2)));
            tag = '###ListNumMetReac'; list = modelProperties.(F{i}).modelProp2.Details.missingmetReactomeID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetMetNetX',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetMetaNetXID,2)));
            tag = '###ListNumMetMetNetX'; list = modelProperties.(F{i}).modelProp2.Details.missingmetMetaNetXID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetBigg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetBiGGID,2)));
            tag = '###ListNumMetBigg'; list = modelProperties.(F{i}).modelProp2.Details.missingmetBiGGID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetBioc',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetmetBioCycID,2)));
            tag = '###ListNumMetBioc'; list = modelProperties.(F{i}).modelProp2.Details.missingmetBioCycID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            % Annotation - Metabolites - Conformity
            tline= regexprep(tline,'###NumMetCPubCh',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetPubChemID,2)));
            tag = '###ListNumMetCPubCh'; list = modelProperties.(F{i}).modelProp2.Details.missingmetPubChemID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCKegg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetKEGGID,2)));
            tag = '###ListNumMetCKegg'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetKEGGID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCSeed',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetSEEDID,2)));
            tag = '###ListNumMetCSeed'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetSEEDID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCIKey',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetInchiKey,2)));
            tag = '###ListNumMetCIKey'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetInchiKey; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCIString',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetInchiString,2)));
            tag = '###ListNumMetCIString'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetInchiString; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCChebi',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetChEBIID,2)));
            tag = '###ListNumMetCChebi'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetChEBIID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCHmdb',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetHMDBID,2)));
            tag = '###ListNumMetCHmdb'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetHMDBID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCReac',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetReactomeID,2)));
            tag = '###ListNumMetCReac'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetReactomeID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCMetNetX',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetMetaNetXID,2)));
            tag = '###ListNumMetCMetNetX'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetMetaNetXID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCBigg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetBiGGID,2)));
            tag = '###ListNumMetCBigg'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetBiGGID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetCBioc',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetConfmetBioCycID,2)));
            tag = '###ListNumMetCBioc'; list = modelProperties.(F{i}).modelProp2.Details.AnnoMetNonConfmetBioCycID; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###ScoreAnnotationMetabolites',num2str(round(modelProperties.(F{i}).modelProp2.Scores.AnnotationMetabolites,2)));

            % Reaction Annotation
            tline= regexprep(tline,'###NumRAnnoAny',num2str(round(modelProperties.(F{i}).modelProp2.rxnWAnno,2)));
            tag = '###ListNumRAnnoAny'; list = modelProperties.(F{i}).modelProp2.Details.rxnWOAnno; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRRhea',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnRheaID,2)));
            tag = '###ListNumRRhea'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnRheaID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRRhea',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnRheaID,2)));
            tag = '###ListNumCRRhea'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnRheaID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRKegg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnKEGGID,2)));
            tag = '###ListNumRKegg'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnKEGGID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRKegg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnKEGGID,2)));
            tag = '###ListNumCRKegg'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnKEGGID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRSeed',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnSEEDID,2)));
            tag = '###ListNumRSeed'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnSEEDID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRSeed',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnSEEDID,2)));
            tag = '###ListNumCRSeed'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnSEEDID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRMetX',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnMetaNetXID,2)));
            tag = '###ListNumRMetX'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnMetaNetXID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRMetX',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnMetaNetXID,2)));
            tag = '###ListNumCRMetX'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnMetaNetXID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRBigg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnBiGGID,2)));
            tag = '###ListNumRBigg'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnBiGGID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRBigg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnBiGGID,2)));
            tag = '###ListNumCRBigg'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnBiGGID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRReac',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnReactomeID,2)));
            tag = '###ListNumRReac'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnReactomeID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRReac',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnReactomeID,2)));
            tag = '###ListNumCRReac'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnReactomeID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRBren',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnBRENDAID,2)));
            tag = '###ListNumRBren'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnBRENDAID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRBren',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnBRENDAID,2)));
            tag = '###ListNumCRBren'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnBRENDAID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRBioc',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnBioCycID,2)));
            tag = '###ListNumRBioc'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnBioCycID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCRBioc',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnBioCycID,2)));
            tag = '###ListNumCRBioc'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnBioCycID; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumREcn',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnrxnECNumbers,2)));
            tag = '###ListNumREcn'; list = modelProperties.(F{i}).modelProp2.Details.missingrxnECNumbers; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCREcn',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnConfrxnECNumbers,2)));
            tag = '###ListNumCREcn'; list = modelProperties.(F{i}).modelProp2.Details.AnnoRxnNonConfrxnECNumbers; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###ScoreAnnotationReactions',num2str(round(modelProperties.(F{i}).modelProp2.Scores.AnnotationReactions,2)));

            % Gene Annotation

            tline= regexprep(tline,'###NumGAnnoAny',num2str(round(modelProperties.(F{i}).modelProp2.geneWAnno,2)));
            tag = '###ListNumGAnnoAny'; list = modelProperties.(F{i}).modelProp2.Details.geneWOAnno; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGRefSeq',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneRefSeqID,2)));
            tag = '###ListNumGRefSeq'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneRefSeqID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGUniProt',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneUniprotID,2)));
            tag = '###ListNumGUniProt'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneUniprotID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGEcoGene',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneEcoGeneID,2)));
            tag = '###ListNumGEcoGene'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneEcoGeneID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGKegg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneKEGGID,2)));
            tag = '###ListNumGKegg'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneKEGGID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGNCBIG',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneEntrezID,2)));
            tag = '###ListNumGNCBIG'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneEntrezID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGNCBIP',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneNCBIProteinID,2)));
            tag = '###ListNumGNCBIP'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneNCBIProteinID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGCCDS',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneCCDSID,2)));
            tag = '###ListNumGCCDS'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneCCDSID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGHPRD',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneHPRDID,2)));
            tag = '###ListNumGHPRD'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneHPRDID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGASAP',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGenegeneASAPID,2)));
            tag = '###ListNumGASAP'; list = modelProperties.(F{i}).modelProp2.Details.missinggeneASAPID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGRefSeq',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneRefSeqID,2)));
            tag = '###ListNumCGRefSeq'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneRefSeqID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGUniProt',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneUniprotID,2)));
            tag = '###ListNumCGUniProt'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneUniprotID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGEcoGene',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneEcoGeneID,2)));
            tag = '###ListNumCGEcoGene'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneEcoGeneID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGKegg',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneKEGGID,2)));
            tag = '###ListNumCGKegg'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneKEGGID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGNCBIG',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneEntrezID,2)));
            tag = '###ListNumCGNCBIG'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneEntrezID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGNCBIP',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneNCBIProteinID,2)));
            tag = '###ListNumCGNCBIP'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneNCBIProteinID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGCCDS',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneCCDSID,2)));
            tag = '###ListNumCGCCDS'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneCCDSID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGHPRD',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneHPRDID,2)));
            tag = '###ListNumCGHPRD'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneHPRDID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumCGASAP',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneConfgeneASAPID,2)));
            tag = '###ListNumCGASAP'; list = modelProperties.(F{i}).modelProp2.Details.AnnoGeneNonConfgeneASAPID; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###ScoreAnnotationGenes',num2str(round(modelProperties.(F{i}).modelProp2.Scores.AnnotationGenes,2)));

            % SBO Annotation
            tline= regexprep(tline,'###NumMSBOAny',num2str(round(modelProperties.(F{i}).modelProp2.metWSBO,2)));
            tag = '###ListNumMSBOAny'; list = modelProperties.(F{i}).modelProp2.Details.metWOSBO; type = 'met';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumMetSBO247',num2str(round(modelProperties.(F{i}).modelProp2.AnnoMetSBO0000247,2)));

            tline= regexprep(tline,'###NumRSBOAny',num2str(round(modelProperties.(F{i}).modelProp2.rxnWSBO,2)));
            tag = '###ListNumRSBOAny'; list = modelProperties.(F{i}).modelProp2.Details.rxnWOSBO; type = 'rxn';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumRSBO176',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnSBO0000176,2)));
            tline= regexprep(tline,'###NumRSBO176',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnSBO0000176,2)));
            tline= regexprep(tline,'###NumRSBO185',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnSBO0000185,2)));
            tline= regexprep(tline,'###NumRSBO627',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnSBO0000627,2)));
            tline= regexprep(tline,'###NumRSBO628',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnSBO0000628,2)));
            tline= regexprep(tline,'###NumRSBO632',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnSBO0000632,2)));
            tline= regexprep(tline,'###NumRSBO629',num2str(round(modelProperties.(F{i}).modelProp2.AnnoRxnSBO0000629,2)));

            tline= regexprep(tline,'###NumGSBOAny',num2str(round(modelProperties.(F{i}).modelProp2.geneWSBO,2)));
            tag = '###ListNumGSBOAny'; list = modelProperties.(F{i}).modelProp2.Details.geneWOSBO; type = 'none';
            tline = writeListItem(tline,tag,list,type);

            tline= regexprep(tline,'###NumGSBO243',num2str(round(modelProperties.(F{i}).modelProp2.AnnoGeneSBO0000243,2)));

            tline= regexprep(tline,'###ScoreAnnotationSBO',num2str(round(modelProperties.(F{i}).modelProp2.Scores.AnnotationSBO,2)));

        end
        html{cnt} = tline;
        cnt = cnt + 1;
    end
    fclose(fid);
    html = html';
    ModelName = regexprep(ModelName,' ','_');
    fid =fopen([reportDir filesep 'modelreport_' ModelName '.html'], 'w');

    for j = 1 : length(html)-1
        fprintf(fid,strcat(html{j},'\n'));
    end

    fclose(fid);
end

function tlineOut = writeListItem(tline,tag,list,type)
%type = 'rxn';
%type = 'met';
if  regexp(tline,tag)
    %  tline= regexprep(tline,'###ListMedComp',num2str(modelProperties.(F{i}).modelProp2.MediumMets+2));
    tlineOut = '';
    for k = 1 : length(list)
        m = list{k};
        if strcmp(type,'rxn')
            tlineOut = strcat(tlineOut,'<li class="list-group-item"> <a href = https://vmh.life/#reaction/',m,' target = "_blank">',m,'</a>','</li>','\n');
        elseif strcmp(type,'met')
            if contains(m,'[')
                mpart = split(m,'[');
                mpart = mpart{1};
            else
                mpart = m;
            end
            tlineOut = strcat(tlineOut,'<li class="list-group-item"> <a href = https://vmh.life/#metabolite/',mpart,' target = "_blank">',m,'</a>','</li>','\n');
        elseif strcmp(type,'none')
            tlineOut = strcat(tlineOut,'<li class="list-group-item">',m,'</li>','\n');
        end
    end
else
    tlineOut = tline;
end

author: opencobra

repo: cobratoolbox


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [c,fc] = FInvSABR4_1(x, f)
% This function computes the inverse of a SABR cumulative distribution
% Since it is based on FSABR2 it is very fast and works for vectors!
f2 = @(t) (x - FSABR2_1(t,f));

% Use simple bisection
a = 0.00001 * ones(1,length(x)); % left starting points
b = 1*ones(1,length(x));        % right starting points

eps = 1e-3;                     % accuracy level
k = 0;                          % init iteration counter

c=(b+a)/2;
fc = f2(c);
% Ind = -eps <= fc & fc <= eps;
while(k <= 200)
    k = k+1;                            % increase iteration counter
%     if sum(Ind) == length(x)
%         break;
%     else
%         b(fc(Ind)<0) = 0.5*(b(fc(Ind)<0)+a(fc(Ind)<0));               % update right boundary
%         a(fc(Ind)>=0) = 0.5*(b(fc(Ind)>=0)+a(fc(Ind)>=0));            % update left boundary
%     end
%     c(Ind) = 0.5*(b(Ind)+a(Ind));
%     fc = f2(c);
%     Ind = -eps <= fc & fc <= eps;

    if( (-eps <= min(fc)) & (max(fc) <= eps))
        break;                          % end if all values have been found
    else
        b(fc<0) = 0.5*(b(fc<0)+a(fc<0));               % update right boundary
        a(fc>=0) = 0.5*(b(fc>=0)+a(fc>=0));            % update left boundary
    end
    c = 0.5*(b+a);
    fc = f2(c);                         % evaluate function at new c
end

end

author: Sable

repo: mcbench-benchmarks


text:
function [Mz_z,Mz_xy,F,ref_eff,Mx_xy,My_xy]=simRf(rf,rephase_factor,prephase_factor)
%simRf Simulate an RF pulse with the given pulse shape.
%   [Mz_z,Mz_xy,F,ref_eff,Mx_xy,My_xy]=simRf(pulse,prephase_factor,rephase_factor)
%   Performs a rapid RF pulse simulation based on the rotation formalism.
%   The algorithm is optimized by using quaternions to represent rotations.
%   The compulsory parameter 'rf' is the Pulseq RF pulse. Optional
%   parameter 'rephase_factor' is needed in several cases e.g. to correclty
%   visualize the phase of the magnetization for slice-selective
%   excitation. Another optional parameter 'prephase_factor' is an
%   experimental parameter useful for simulating refocusing pulses or
%   spoiling needed.
%   Return values:
%     Mz_z,Mz_xy:  z and xy comnponents of the magnetisation after the pulse
%                  assuming the unit magnetization was aligned with z before
%                  the pulse. Useful for assessing excitation RF pulses.
%     F:           frequency axis in Hz
%     ref_eff:     Refocusing efficiency of the pulse as a complex value.
%                  Magnitude of ref_eff seems to closely follow Mz_z. Phase
%                  of ref_eff is related to the effective phase of the RF
%                  pulse, e.g. the axis of the planar flip.
%     Mx_xy,My_xy: xy magnetizations after the RF pulse assuming the unit
%                  magnetization was aligned with x or y axis prior to the
%                  pulse, respectively. Useful for detailed analyses of
%                  refocusing pulses.
%
%   The implementation was inspired by the example by Dr. Tony Stoecker
%   (https://github.com/stoeckert/mr-simu-example-ismrm19)
%   The algorithm was rewritten to quaternions and vectorized for
%   performance by MZ
%

bw_mul=4;    % simulation bandwidth (multiplier of the pulse bandwidth)
df=1;        % spectral resolution [Hz]
dt=10e-6;    % (re-)sampling interval

if nargin < 2
    if isfield(rf,'use') && strcmp(rf.use,'refocusing')
        rephase_factor = 0;
    else
        rephase_factor = -(rf.shape_dur-mr.calcRfCenter(rf))/rf.shape_dur;
    end
end

if nargin < 3
    prephase_factor = 0;
end

[bw,f0,spectrum,FF,rfs,tt]=mr.calcRfBandwidth(rf,0.5,df*10,dt);

T     = (1:round(rf.shape_dur/dt))*dt-0.5*dt;                           % timesteps axis [s]
F     = 2*pi*linspace(f0-bw_mul*bw/2,f0+bw_mul*bw/2,bw/df)';               % offset frequencies [rad/s]

shapea = interp1(rf.t, 2*pi*rf.signal.*exp(1i*(rf.phaseOffset+2*pi*rf.freqOffset*rf.t)),T,'linear',0);

% intialize result vectors
M_ROT=zeros(size(F));
Z_ROT=zeros(size(F));
sf=size(F);
q=zeros(sf(1),4);
q(:,1)=1; % init rotation quaternions

% prephaser / left spoiler
W = -F*dt*length(T)*prephase_factor;     % effective field rotation angle
Q = [cos(W/2) zeros(sf) zeros(sf) sin(W/2)];
q=quat_multiply(q,Q);

% RF pulse simulation
for j=1:length(T)
    W = -dt*sqrt(abs(shapea(j))^2+F.^2); % effective field rotation angles
    n = dt * [real(shapea(j))*ones(sf) imag(shapea(j))*ones(sf) F]./abs(W); % effective field rotation axes
    Q = [cos(W/2) sin(W/2).*n];
    q=quat_multiply(q,Q);
end

% rephaser / right spoiler / refocusing pulse
W = -F*dt*length(T)*rephase_factor;     % effective field rotation angle
Q = [cos(W/2) zeros(sf) zeros(sf) sin(W/2)];
q=quat_multiply(q,Q);

% export results
F=F/(2*pi);
m=zeros(sf(1),4);

% excitation: start with M0=M_z
m(:,4)=1;
m0rf=quat_multiply(quat_conj(q),quat_multiply(m,q));
Mz_z=m0rf(:,4);
Mz_xy=m0rf(:,2)+1i*m0rf(:,3);

% refocusing: start both with M0=M_x and them M0=M_y
m=zeros(sf(1),4);
m(:,2)=1;
Mx_xy=quat_multiply(quat_conj(q),quat_multiply(m,q));
Mx_xy=Mx_xy(:,2)+1i*Mx_xy(:,3);
m=zeros(sf(1),4);
m(:,3)=1;
My_xy=quat_multiply(quat_conj(q),quat_multiply(m,q));
My_xy=My_xy(:,2)+1i*My_xy(:,3);
ref_eff=(Mx_xy+My_xy*1i)/2;
end

function qout = quat_multiply( q, r )
%  quat_multiply: Calculate the product of two quaternions.

% Calculate vector portion of quaternion product
% vec = s1*v2 + s2*v1 + cross(v1,v2)
vec = [q(:,1).*r(:,2) q(:,1).*r(:,3) q(:,1).*r(:,4)] + ...
         [r(:,1).*q(:,2) r(:,1).*q(:,3) r(:,1).*q(:,4)]+...
         [ q(:,3).*r(:,4)-q(:,4).*r(:,3) ...
           q(:,4).*r(:,2)-q(:,2).*r(:,4) ...
           q(:,2).*r(:,3)-q(:,3).*r(:,2)];

% Calculate scalar portion of quaternion product
% scalar = s1*s2 - dot(v1,v2)
scalar = q(:,1).*r(:,1) - q(:,2).*r(:,2) - ...
             q(:,3).*r(:,3) - q(:,4).*r(:,4);

qout = [scalar  vec];
end

function q = quat_conj( q )
%  quat_conj Calculate the conjugate of a quaternion.
q(:,2:4) = -q(:,2:4);
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function circle_trajectory %(p1, p2, p3)

close all

p1=rand(1,3);p2=rand(1,3);p3=rand(1,3);
[pc,r,v1,v2] = circlefit3d(p1,p2,p3);
plot3(p1(:,1),p1(:,2),p1(:,3),'ro');hold on;plot3(p2(:,1),p2(:,2),p2(:,3),'go');plot3(p3(:,1),p3(:,2),p3(:,3),'bo');

x1 = (p1-pc);
x1 = x1/norm(x1);
phi1 = dot(v1,x1);
x2 = (p2-pc);
x2 = x2/norm(x2);
phi2 = dot(v1,x2);

angle = [phi1 phi2];
%angle = 0:0.1:2*pi;
%angle = pi/4:0.1:pi/2;
for i=1:length(angle)
        a = angle(i);
        x = pc(:,1)+cos(a)*r.*v1(:,1)+sin(a)*r.*v2(:,1);
        y = pc(:,2)+cos(a)*r.*v1(:,2)+sin(a)*r.*v2(:,2);
        z = pc(:,3)+cos(a)*r.*v1(:,3)+sin(a)*r.*v2(:,3);
        plot3(x,y,z,'r+');
end
axis equal;grid on;rotate3d on;

% fit a circle with three points

%[pm, r, v1, v2] = circlefit3d(p1', p2', p3');

repo: ARTE


text:
%  Figure 3.16      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.16
%  fig3_14.m
%  Example 3.23
clf;
num=[2 1];
den=[1 3 2];
t=0:0.1:6;
y=impulse(num,den,t);
plot(t,y,'-')
grid
title('Fig. 3.16  Example 3.23 system impulse response.')
xlabel('Time (sec)')
ylabel('h(t)')
% grid
nicegrid


text:
%%%
%> @brief  computes the pose vector v from an homogeneous transformation A
%> param   A homogeneous transformation
%> return  v pose vector
%> @author Giorgio Grisetti
%%%
function v = t2v(A)
% T2V homogeneous transformation to vector
v(1:2,1) = A(1:2,3);
v(3,1) = atan2(A(2,1), A(1,1));
end

repo: graphslam


text:
function [labels rts] = graph_connected_components(C)
% C - connection matrix
% labels =[1 1 1 2 2 3 3 ...]  lenght(labels)=L, label for each vertex
% labels(i) is order number of connected component, i is vertex number
% rts - roots, numbers of started vertex in each component

L=size(C,1); % number of vertex

% Breadth-first search:
labels=zeros(1,L); % all vertex unexplored at the begining
rts=[];
ccc=0; % connected components counter
while true
    ind=find(labels==0);
    if ~isempty(ind)
        fue=ind(1); % first unexplored vertex
        rts=[rts fue];
        list=[fue];
        ccc=ccc+1;
        labels(fue)=ccc;
        while true
            list_new=[];
            for lc=1:length(list)
                p=list(lc); % point
                cp=find(C(p,:)); % points connected to p
                cp1=cp(labels(cp)==0); % get only unexplored vertecies
                labels(cp1)=ccc;
                list_new=[list_new cp1];
            end
            list=list_new;
            if isempty(list)
                break;
            end
        end
    else
        break;
    end
end


text:
function f = solharm(l, m, varargin)
% SOLHARM Complex-valued, solid harmonic of degree L and order M.
%
%   F = SOLHARM(L, M) returns the degree L, order M real-valued
%   solid harmonic on the ball.  F is normalized so that its two-norm
%   over the sphere is 1.
%
%   F = SOLHARM(L, M, 'complex') returns the degree L, order M complex-valued
%   solid harmonic on the ball.  F is normalized so that its two-norm
%   over the sphere is 1.
%
% See also spherefun.sphharm.

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

if nargin == 2
    % Complex solid harmonics
    F = complex_solharm(l,abs(m));
    % Determine if the cos or sin term should be added:
    pos = abs(max(0, sign(m+1)));
    % Compute the real spherical harmonic:
    % Real part of pos = 1, imaginary part if pos = 0
    F(:,1,:) = (-1)^(1-pos)*F(:,2*abs(m)+1,:);
    F = (-1i)^(1-pos)*F/(1+(m~=0));
    f = ballfun(F,'coeffs');

elseif nargin > 2 && strcmp(varargin{1},'complex')
    % Complex solid harmonics
    F = complex_solharm(l,m)/sqrt(1+(abs(m)>0));
    f = ballfun(F,'coeffs');

else
    % Don't know what to do.
    error('CHEBFUN:BALLFUN:solharm:input', ...
                'Unrecognized inputs.')
end
end

function F = complex_solharm(l,m)
% SOLHARM Complex-valued, solid harmonic of degree L and order M.
%
%   F = SOLHARM(L, M) returns the degree L, order M complex-valued
%   solid harmonic on the ball.  F is normalized so that its two-norm
%   over the sphere is 1.

if ( l < abs(m) )
    error('CHEBFUN:BALLFUN:solHarm', 'Degree must be >= order for solid harmonic ');
end

% abs(m) <= l
% Compute P^abs(m)_l
S = [l+1,2*abs(m)+1,2*l+1];
Plm = normalized_legendre(l,abs(m));
if m < 0
   Plm = (-1)^m*Plm;
end
% Normalize the solid harmonic so that its two-norm over the ball is 1
% This is the normalization constant so that the integral of r^l is 1
Plm = Plm*sqrt(2*l+3);

% Compute the chebyshev coefficients of r^l
r = chebpts(l+1, [-1 1]);
monomial = r.^l;
monomial = chebtech2.vals2coeffs(monomial);

% Return the Spherical Harmonic Y^m_l
F = zeros(S);
F(:,abs(m)+m+1,:) = reshape(monomial*Plm.',l+1,1,2*l+1);
end

function Pold = normalized_legendre(l_max,m_max)
%NORMALIZED_LEGENDRE   Generate the fully normalized associated Legendre
% polynomials P^m_max_l_max
%
% Normalization for the spherical harmonics is \tilde{P}^m_l =
% P^m_l/(2*sqrt(pi)).

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

%% Below is the implementation the Modified Forward Column (MFC) method
% described in the Holmes and Featherstones paper (2002)
% It computes P^m_n/u^m by a stable recurrence to avoid numerical errors
% near the poles

% Discretization number of the polynomial
p = 2*l_max+1;

% Interpolation points
th = pi*trigpts(p);

% Precompute vector cos(th)
CosTh = cos(th);

%% Compute P_m_max^m_max / u^m_max

% Initialize P^0_0 / u^0
Pold = ones(p, 1);

% Compute P^m_m/u^m
for m = 1:m_max
    % Compute Pm^m with Pm-1^m-1
    Pold = sqrt((2*m+1)/(2*m-(m==1)))*Pold;
end

% Initialize the recurrence (Pm^m-1 does not exist)
Poldold = zeros(p, 1);

%% Compute P^m_l / u^m with the recurrence formula, m_max+1 <= l <= l_max
for l = m_max+1:l_max
    anm = sqrt((4*l^2-1)/((l-m_max)*(l+m_max)));
    bnm = sqrt((2*l+1)*(l+m_max-1)*(l-m_max-1)/((l-m_max)*(l+m_max)*(2*l-3)));
    % Compute the normalized associated legendre polynomial P^m_l/u^m
    Pl = anm*CosTh.*Pold - bnm*Poldold;

    % Update the polynomials for the recurrence
    Poldold = Pold;
    Pold = Pl;
end

% Normalize the polynomial and recover associated Legendre polynomials
% Divide by sqrt(2) if m_max > 0
Pold = (-1)^m_max*sin(th).^m_max.*Pold/sqrt(4*pi);

% FFT to recover the coefficients
Pold = trigtech.vals2coeffs(Pold);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function compute_print_results(masks, seg_obj, gt_fldr)
% You can use this function to compute accuracy scores for the segment
% proposals generated against some GT. This function computes the overlap
% for each segment to all GTs. Then prints the results and stores them. The
% results are stored at the path defined by scores_dirpath defined in
% internal_params.m. The file stores some collated scores including
% overlap, and covering (collated_scores). It also stores the number of
% segments generated by parametric min-cut, and the number remaining after
% filtration (init_num_segs and final_num_segs). Also overlaps of all
% segments against the GT is stored (Q). The function finally displays the
% results in a pretty table :)
%
% For instance, if you are working with PASCAL VOC (e.g. stored at
%  ~/pascalvoc12), you can do the following:
%       [masks, seg_obj, total_time] = ...
%           rigor_obj_segments('~/pascalvoc12/JPEGImages/2007_000033.jpg');
%       compute_print_results(masks, seg_obj, ...
%           '~/pascalvoc12/SegmentationObject');
%
% You can also run compute_print_results.m by loading the data saved to
% disk by rigor_obj_segments.m, and supplying masks, and seg_obj to the
% function.
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    img_name = seg_obj.input_info.img_name;

    % compute the quality of each segment given the GT
    [Q, collated_scores] = SvmSegm_segment_quality(img_name, gt_fldr, ...
                                                   masks, 'overlap');

    % collect scores for all segments
    overlap_scores = [Q.q];
    tp = cell2mat(arrayfun(@(q) q.extra_info.tp, Q, 'UniformOutput',false));
    fp = cell2mat(arrayfun(@(q) q.extra_info.fp, Q, 'UniformOutput',false));
    fn = cell2mat(arrayfun(@(q) q.extra_info.fn, Q, 'UniformOutput',false));

    num_gt = numel(Q);
    num_sgmntrs = numel(seg_obj.segm_params.graph_methods);

    % # of segments returned by each segmenter and the remaining number
    % after the filteration step
    init_num_segs = seg_obj.num_segs.init_segs;
    final_num_segs = seg_obj.num_segs.after_clustering_FINAL;

    % get the name of each segmenter
    sgmntr_names = seg_obj.segm_params.graph_methods;
    sgmntr_names = cellfun(@(s) s(1:strfind(s, 'Graph')-1), ...
                           sgmntr_names, 'UniformOutput',false);

    % compute the scores over different segmenters individually
    start_idx = 1;
    avg_best_overlap_segmenters = zeros(1, num_sgmntrs);
    collective_overlap_segmenters = zeros(1, num_sgmntrs);
    best_seg_overlap_segmenters = zeros(num_gt, num_sgmntrs);
    best_seg_idx_segmenters = zeros(num_gt, num_sgmntrs);
    % iterate over each segmenter
    for sgmntr_idx = 1:num_sgmntrs
        % get the scores for just the current segmenter
        segmenter_overlap = ...
            overlap_scores(start_idx:start_idx + ...
                                     final_num_segs(sgmntr_idx)-1,:);
        % compute the best scores for each GT for the current segmenter
        [best_seg_overlap, best_seg_idx] = max(segmenter_overlap);
        avg_best_overlap_segmenters(sgmntr_idx) = mean(best_seg_overlap);
        best_seg_overlap_segmenters(:,sgmntr_idx) = best_seg_overlap;
        best_seg_idx_segmenters(:,sgmntr_idx) = ...
                                            start_idx + best_seg_idx - 1;

        % compute the collective overlap for all GT's for the current
        % segmenter
        best_ind = sub2ind(size(tp), ...
                           best_seg_idx_segmenters(:,sgmntr_idx)', ...
                           1:num_gt);
        total_tp = tp(best_ind);
        total_union = fp(best_ind) + fn(best_ind) + tp(best_ind);
        collective_overlap_segmenters(sgmntr_idx) = ...
            sum(total_tp) / sum(total_union);

        start_idx = start_idx + final_num_segs(sgmntr_idx);
    end

    % store the scores computed
    collated_scores.collective_overlap_segmenters = ...
                                            collective_overlap_segmenters;
    collated_scores.avg_best_overlap_segmenters = ...
                                            avg_best_overlap_segmenters;
    collated_scores.best_seg_overlap_segmenters = ...
                                            best_seg_overlap_segmenters;
    collated_scores.best_seg_idx_segmenters = best_seg_idx_segmenters;

    % print the table of scores
    print_scores_table([Q.gt_seg_szs], sgmntr_names, init_num_segs, ...
                       final_num_segs, collated_scores);

    % save the scores
    scores_filepath = seg_obj.return_scores_filepath();

    dir_name = fileparts(scores_filepath);
    if(~exist(dir_name, 'dir'))
        mkdir(dir_name);
    end
    save(scores_filepath, 'Q', 'sgmntr_names', 'init_num_segs', ...
                          'final_num_segs', 'collated_scores');

    % print a summary line, useful for noting down results/performance
    print_summary(seg_obj, collated_scores);
end

function print_scores_table(gt_seg_szs, sgmntr_names, init_num_segs, ...
                            final_num_segs, cscores)
    % get the row headers
    row_txts = arrayfun(@(i,gt_sz) sprintf('GT %2d (%5d)', i, gt_sz), ...
                        1:length(gt_seg_szs), gt_seg_szs, ...
                        'UniformOutput', false);
    row_txts = [{'', 'Init. # Segs', 'Final # Segs', 'Avg. Best Overlap', ...
                'Collective Overlap'}, row_txts];

    % get the column headers
    col_txts = [{''}, sgmntr_names, {'Overall'}];

    % compute the width based on the header lengths
    first_col_width = max(cellfun(@length, row_txts)) + 1;
    other_col_width = max(cellfun(@length, col_txts)) + 1;

    % create a divider line
    total_width = 2 + first_col_width + ...
                  (length(col_txts) - 1)*(other_col_width + 1);
    divider_line = [repmat('-', 1, total_width), '\n'];

    fprintf(divider_line);

    % printf first column
    print_line(row_txts{1}, col_txts(2:end), 's', ...
               first_col_width, other_col_width, length(sgmntr_names)+1);

    fprintf(divider_line);

    % print initial number of segments
    print_line(row_txts{2}, [init_num_segs, sum(init_num_segs)], 'd', ...
               first_col_width, other_col_width);

    % print final number of segments
    print_line(row_txts{3}, [final_num_segs, sum(final_num_segs)], 'd', ...
               first_col_width, other_col_width, length(sgmntr_names)+1);

    fprintf(divider_line);

    % print average best overlap
    print_line(row_txts{4}, [cscores.avg_best_overlap_segmenters, ...
                             cscores.avg_best_overlap], '.4f', ...
               first_col_width, other_col_width, length(sgmntr_names)+1);

    % print collective best overlap
    print_line(row_txts{5}, [cscores.collective_overlap_segmenters, ...
                             cscores.collective_overlap], '.4f', ...
               first_col_width, other_col_width);

    fprintf(divider_line);

    % print stats for each GT
    for idx = 1:size(cscores.best_seg_overlap_segmenters,1)
        [best_val, best_idx] = ...
            max(cscores.best_seg_overlap_segmenters(idx,:));
        print_line(row_txts{5+idx}, ...
                   [cscores.best_seg_overlap_segmenters(idx,:), ...
                    best_val], '.4f', ...
                   first_col_width, other_col_width, best_idx);
    end

    fprintf(divider_line);
end

function print_line(header, vals, spec, first_col_width, ...
                    other_col_width, strong_idx)
    % controls whether a certain value is going to be bold
    strong = false(1, length(vals));
    if exist('strong_idx','var')
        strong(strong_idx) = true;
    end
    strong_o = repmat({'<strong>'}, 1, length(vals));
    strong_e = repmat({'</strong>'}, 1, length(vals));
    strong_o(~strong) = {''};
    strong_e(~strong) = {''};

    % print header
    eval(['fprintf(''|%-', num2str(first_col_width), 's'', header)']);

    if spec == 's'
        % in case of string
        temp = [strong_o; vals; strong_e];
    else
        temp = [strong_o; num2cell(vals); strong_e];
    end
    eval(['fprintf(''|%s%-', num2str(other_col_width), spec '%s'', temp{:})']);

    fprintf('|\n');
end

function print_summary(seg_obj, cscores)
    for idx = 1:numel(seg_obj.segm_params.graph_methods)
        fprintf(1, '% 4d (% 4d)', ...
                seg_obj.num_segs.init_segs(idx), ...
                seg_obj.num_segs.after_clustering_FINAL(idx));
        if idx < numel(seg_obj.segm_params.graph_methods)
            fprintf(1, ', ');
        else
            fprintf(1, ' - ');
        end
    end
    best_gt_scores = max(cscores.best_seg_overlap_segmenters,[],2);
    for idx = 1:length(best_gt_scores)
        fprintf(1, '%.4f', best_gt_scores(idx));
        if idx < length(best_gt_scores)
            fprintf(1, ' + ');
        else
            fprintf(1, ' => ');
        end
    end
    for idx = 1:length(cscores.avg_best_overlap_segmenters)
        fprintf(1, '%.4f', cscores.avg_best_overlap_segmenters(idx));
        if idx < length(cscores.avg_best_overlap_segmenters)
            fprintf(1, ' + ');
        else
            fprintf(1, ' => ');
        end
    end
    fprintf('%.4f\n', cscores.avg_best_overlap);
end

author: Cloud-CV

repo: object-proposals


text:
function [ha] = km22ha(km2)
% Convert area from square kilometers to hectares.
% Chad A. Greene 2012
ha = km2*100;


text:
function [ eigenvalues ] = ml_kernel_grid_search(X,options,kpars)
%ML_KERNEL_GRID_SEARCH
%
%   input -----------------------------------------------------------------
%
%       o X           : (N x D),            original data N samples
%                                           dimension D.
%
%       o options     : struct,             same as for ml_projection.
%
%
%       o kpars       : (1 x T) or (2 x T), depending on if you are use the
%                                           gaussina or polynomial kernel.
%
%   output ----------------------------------------------------------------
%
%       o eigenvalues : (M x T),             T eigenvalues of dimension M.
%
%
%
%

[m,T] = size(kpars);
disp(['num parameters: ' num2str(T)]);
disp(' ');

eigenvalues = [];

for t=1:T

    options.kpar = [];
    if m==1
        options.kpar(1)     = kpars(1,t);
    else
        options.kpar(1)     = kpars(1,t);
        options.kpar(2)     = kpars(2,t);
    end

    disp([num2str(t) '/' num2str(T)]);

    [~,mapping]             = ml_projection(X,options);
    eigenvalues = [eigenvalues;mapping.L(:)'];

end
disp(' ');

end

author: epfl-lasa


text:
function out = spm_run_reorient(varargin)
% SPM job execution function
% takes a harvested job data structure and call SPM functions to perform
% computations on the data.
% Input:
% job    - harvested job data structure (see matlabbatch help)
% Output:
% out    - computation results, usually a struct variable.
%__________________________________________________________________________
% Copyright (C) 2006-2014 Wellcome Trust Centre for Neuroimaging

% Volkmar Glauche
% $Id: spm_run_reorient.m 6078 2014-06-30 18:10:33Z guillaume $

job = varargin{1};
if isfield(job.transform,'transprm')
    job.transform.transM = spm_matrix(job.transform.transprm);
elseif isfield(job.transform,'transF')
    load(char(job.transform.transF), 'M');
    job.transform.transM = M;
end

K = numel(job.srcfiles);
spm_progress_bar('Init', K, 'Reorient', 'Images completed');
if isempty(job.prefix)

    % read and write separately, so duplicates get harmlessly overwritten
    MM = zeros(4, 4, K);
    for k = 1:K
        MM(:, :, k) = spm_get_space(job.srcfiles{k});
    end
    for k = 1:K
        spm_get_space(job.srcfiles{k}, job.transform.transM * MM(:, :, k));
        spm_progress_bar('Set',k);
    end
    out.files = job.srcfiles;

else

    out.files = cell(size(job.srcfiles));
    for k = 1:K
        V       = spm_vol(job.srcfiles{k});
        X       = spm_read_vols(V);
        V.mat   = job.transform.transM * V.mat;
        V.fname = spm_file(V.fname, 'prefix',job.prefix);
        spm_write_vol(V,X);
        out.files{k} = [V.fname ',' num2str(V.n(1))];
        spm_progress_bar('Set',k);
    end

end
spm_progress_bar('Clear');

author: spm


text:
function [ o, x, w ] = en_r2_07_2 ( n )

%*****************************************************************************80
%
%% EN_R2_07_2 implements the Stroud rule 7.2 for region EN_R2.
%
%  Discussion:
%
%    The rule has order O = 2^(N+1) + 4 * N^2.
%
%    The rule has precision P = 7.
%
%    EN_R2 is the entire N-dimensional space with weight function
%
%      w(x) = exp ( - x1^2 - x2^2 ... - xn^2 )
%
%    The rule requires 3 <= N.
%
%    The reference has a typographical error in the description of this rule.
%    The formula:
%
%      (t,t,t,...,t)FS
%
%    should read
%
%      (t,t,0,...,0)FS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%    3 <= N.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EN_R2_07_2 - Fatal error!\n' );
    fprintf ( 1, '  3 <= N is required.\n' );
    error ( 'EN_R2_07_2 - Fatal error!' );
  end

  o = 2^(n+1) + 4 * n^2;
  volume = sqrt ( pi^n );

  rho1 = sqrt ( ( n + 2 - sqrt ( 2 * ( n + 2 ) ) ) / 2 );
  rho2 = sqrt ( ( n + 2 + sqrt ( 2 * ( n + 2 ) ) ) / 2 );
  a1 = ( n + 2 + sqrt ( 2 * ( n + 2 ) ) ) / 2 / ( n + 2 );
  a2 = ( n + 2 - sqrt ( 2 * ( n + 2 ) ) ) / 2 / ( n + 2 );

  r = 1.0;
  s = sqrt ( 1 / n );
  t = sqrt ( 1 / 2 );
  b = ( 8 - n ) * volume / n / ( n + 2 ) / ( n + 4 );
  c = n^3 * volume / 2^n / n / ( n + 2 ) / ( n + 4 );
  d = 4 * volume / n / ( n + 2 ) / ( n + 4 );

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  k = 0;
%
%  2 * 2 * N points.
%
  for i = 1 : n
    k = k + 1;
    x(i,k) = - rho1 * r;
    w(k) = a1 * b;
    k = k + 1;
    x(i,k) = - rho2 * r;
    w(k) = a2 * b;
    k = k + 1;
    x(i,k) = + rho1 * r;
    w(k) = a1 * b;
    k = k + 1;
    x(i,k) = + rho2 * r;
    w(k) = a2 * b;
  end
%
%  2 * 2^N points.
%
  k = k + 1;
  x(1:n,k) = - rho1 * s;
  w(k) = a1 * c;
  k = k + 1;
  x(1:n,k) = - rho2 * s;
  w(k) = a2 * c;
  more = 1;
  while ( more )
    more = 0;
    for i = n : -1 : 1
      if ( x(i,k) < 0.0 )
        k = k + 1;
        x(1:n,k) =     x(1:n,k-2);
        x(i,k)     =   abs ( x(i,k) );
        x(i+1:n,k) = - abs ( x(i+1:n,k) );
        w(k) = a1 * c;
        k = k + 1;
        x(1:n,k) =     x(1:n,k-2);
        x(i,k)     =   abs ( x(i,k) );
        x(i+1:n,k) = - abs ( x(i+1:n,k) );
        w(k) = a2 * c;
        more = 1;
        break;
      end
    end
  end
%
%  2 * 4 * ( N * ( N - 1 ) / 2 ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - rho1 * t;
      x(j,k) = - rho1 * t;
      w(k) = a1 * d;
      k = k + 1;
      x(i,k) = - rho1 * t;
      x(j,k) = + rho1 * t;
      w(k) = a1 * d;
      k = k + 1;
      x(i,k) = + rho1 * t;
      x(j,k) = - rho1 * t;
      w(k) = a1 * d;
      k = k + 1;
      x(i,k) = + rho1 * t;
      x(j,k) = + rho1 * t;
      w(k) = a1 * d;
      k = k + 1;
      x(i,k) = - rho2 * t;
      x(j,k) = - rho2 * t;
      w(k) = a2 * d;
      k = k + 1;
      x(i,k) = - rho2 * t;
      x(j,k) = + rho2 * t;
      w(k) = a2 * d;
      k = k + 1;
      x(i,k) = + rho2 * t;
      x(j,k) = - rho2 * t;
      w(k) = a2 * d;
      k = k + 1;
      x(i,k) = + rho2 * t;
      x(j,k) = + rho2 * t;
      w(k) = a2 * d;
    end
  end

  return
end


text:
% CircularSparseMat
%   Multidimensional circular sparse matrix. Builds on SparseMat adding circular
%   wrapping of indices
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/18/2010
classdef CircularSparseMat < edu.stanford.covert.util.SparseMat
    properties (SetAccess=protected)
        circularDims
    end

    %constructor
    methods
        %five ways to instatiate:
        %(1) x = CircularSparseMat()
        %    Creates an empty sparse matrix of size 0x2 with no dimensions
        %    wrapping
        %
        %(2) x = CircularSparseMat(mat)
        %    Casts matrix mat to sparse matrix with no dimensions wrapping
        %
        %(3) x = CircularSparseMat(mat, circularDims)
        %    Casts matrix mat to sparse matrix with selected dimensions wrapping
        %
        %(4) x = CircularSparseMat(subs, vals, sz)
        %    Creates sparse matrix with size sz and non-zeros with subscripts
        %    subs and values vals with no dimensions wrapping
        %
        %(5) x = CircularSparseMat(subs, vals, sz, circularDims)
        %    Creates sparse matrix with size sz and non-zeros with subscripts
        %    subs and values vals with selected dimensions wrapping
        function this = CircularSparseMat(varargin)
            import edu.stanford.covert.util.CircularSparseMat;

            switch nargin
                case 4, supervarargin = varargin(1:3);
                    circDims = varargin{end};

                    subs = varargin{1};
                    siz  = varargin{3};
                    if ~isempty(subs)
                        for i = 1:numel(circDims)
                            subs(:, circDims(i)) = mod(subs(:, circDims(i)) - 1, siz(circDims(i))) + 1;
                        end
                        supervarargin{1} = subs;
                    end
                case 0,
                case 1, supervarargin = varargin(1);
                    circDims = varargin{1}.circularDims;
                case 2, supervarargin = varargin(1);
                    circDims = varargin{end};
                case 3, supervarargin = varargin(1:3);
                    circDims = zeros(1,0);
                otherwise
                    throw(MException('CircularSparseMat:error','no constructor matches the calling signature'));
            end

            this = this@edu.stanford.covert.util.SparseMat(supervarargin{:});

            validateattributes(circDims, {'numeric'}, {'positive', 'real', 'integer', '<=', this.ndims});
            if numel(circDims) > 1 && ~all(diff(sort(circDims)))
                throw(MException('CircularSparseMat:invalidDimensions', 'circularDims must be positive real intergers in the range 1..ndims(this)'));
            end

            if numel(circDims) == 1
                this.circularDims = circDims;
            else
                this.circularDims = sort(circDims);
            end
        end

        function this = normalize(this, reSize, reFind, reSort, oldDims)
            this = this.normalize@edu.stanford.covert.util.SparseMat(reSize, reFind, reSort, oldDims);

            %inlined is ismember optimized for small number of circular dimensions
            if isscalar(this.circularDims)
                this.circularDims = find(oldDims == this.circularDims);
            else
                tfs = false(size(oldDims));
                for i = 1:numel(oldDims)
                    tfs(i) = any(this.circularDims == oldDims(i));
                end
                this.circularDims = find(tfs);
            end
        end

        function value = isDimCircular(this, dim)
            value = any(this.circularDims == dim);
        end

        function tf = isequal(A, B)
            tf = isequal@edu.stanford.covert.util.SparseMat(A, B) && ...
                isequal(A.circularDims, B.circularDims);
        end

        function C = tprod(A, B, Aid, Bid)
            Aod = setdiff(1:A.ndims, Aid);
            Bod = setdiff(1:B.ndims, Bid);

            circDims = [];
            if isa('CircularSparseMat',A)
                [tfs, idxs] = ismember(A.circularDims, Aod);
                circDims = idxs(tfs);
            end
            if isa('CircularSparseMat',A)
                [tfs, idxs] = ismember(B.circularDims, Bod);
                circDims = [circDims idxs(tfs)+numel(Aod)];
            end

            C = tprod@edu.stanford.covert.util.SparseMat(A, B);
            C.circularDims = circDims;
        end

        function this = subsasgn(this, s, rhs)
            %Case 1: Dot reference to properties
            if strcmp(s.type,'.') || strcmp(s.type, '{}')
                this = this.subsasgn@edu.stanford.covert.util.SparseMat(s, rhs);
                return;
            end

            if numel(s.subs)==1
                %Case 2: Subscripts of elements
                subs = s.subs{1};

                if size(subs,2) ~= this.ndims
                    throw(MException('SparseMat:invalidDimensions','Subscripts must have number of columns equal to dimensions of matrix.'));
                end
                for i=1:numel(this.circularDims)
                    subs(:,this.circularDims(i)) = mod(subs(:,this.circularDims(i))-1, this.siz(this.circularDims(i)))+1;
                end

                s.subs{1}=subs;
            else
                %Case 3: Subscripts ranges
                for i=1:numel(this.circularDims)
                    if this.circularDims(i) > numel(s.subs)
                        break;
                    end

                    subs = s.subs{this.circularDims(i)};

                    if isnumeric(subs)
                        subs = mod(subs-1, this.siz(this.circularDims(i)))+1;
                    end

                    s.subs{this.circularDims(i)} = subs;
                end
            end

            this = this.subsasgn@edu.stanford.covert.util.SparseMat(s, rhs);
        end

        function value = subsref(this, s)
            %Case 1: Dot reference to properties
            if strcmp(s.type,'.') || strcmp(s.type, '{}')
                value = this.subsref@edu.stanford.covert.util.SparseMat(s);
                return;
            end

            if numel(s.subs)==1
                %Case 2: Subscripts of elements
                subs = s.subs{1};

                if size(subs,2) ~= this.ndims
                    throw(MException('SparseMat:invalidDimensions','Subscripts must have number of columns equal to dimensions of matrix.'));
                end
                for i=1:numel(this.circularDims)
                    subs(:,this.circularDims(i)) = mod(subs(:,this.circularDims(i))-1, this.siz(this.circularDims(i)))+1;
                end

                s.subs{1}=subs;
            else
                %Case 3: Subscripts ranges
                for i=1:numel(this.circularDims)
                    if this.circularDims(i) > numel(s.subs)
                        break;
                    end

                    subs = s.subs{this.circularDims(i)};

                    if isnumeric(subs)
                        subs = mod(subs-1, this.siz(this.circularDims(i)))+1;
                    end

                    s.subs{this.circularDims(i)} = subs;
                end
            end

            value = this.subsref@edu.stanford.covert.util.SparseMat(s);
        end

        function display(this)
            if isempty(this.circularDims)
                circDims = 'no circular dimensions';
            else
                circDims = sprintf('circular dimensions {%d', this.circularDims(1));
                if numel(this.circularDims)>1
                    circDims = [circDims sprintf(',%d', this.circularDims(2:end))];
                end
                circDims = [circDims '}'];
            end

            fprintf('%d',this.siz(1));
            fprintf('x%d',this.siz(2:end));
            fprintf(' %s of type %s with %s and %d non-zeros\n', ...
                class(this), valueClass(this), circDims, this.nnz);
            for i=1:this.nnz
                fprintf('(%d',this.subs(i,1));
                fprintf(',%d',this.subs(i,2:end));
                fprintf(')\t\t%.4f\n',this.vals(i));
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
clear all;
randn('seed',0);
format compact;
param.num_threads=-1; % all cores (-1 by default)
param.verbose=false;   % verbosity, false by default
param.lambda=0.1; % regularization parameter
param.it0=1;      % frequency for duality gap computations
param.max_it=100; % maximum number of iterations
param.L0=0.1;
param.tol=1e-5;
param.intercept=false;
param.pos=false;

graph.eta_g=[1 1 1 1 1];
graph.groups=sparse([0 0 0 1 0;
                     0 0 0 0 0;
                     0 0 0 0 0;
                     0 0 0 0 0;
                     0 0 1 0 0]);   % g5 is included in g3, and g2 is included in g4
graph.groups_var=sparse([1 0 0 0 0;
                         1 0 0 0 0;
                         1 0 0 0 0 ;
                         1 1 0 0 0;
                         0 1 0 1 0;
                         0 1 0 1 0;
                         0 1 0 0 1;
                         0 0 0 0 1;
                         0 0 0 0 1;
                         0 0 1 0 0]); % represents direct inclusion relations

X=randn(100,10);
param.verbose=true;
%X=eye(10);
X=X-repmat(mean(X),[size(X,1) 1]);
X=mexNormalize(X);
Y=randn(100,1);
Y=Y-repmat(mean(Y),[size(Y,1) 1]);
Y=mexNormalize(Y);
W0=zeros(size(X,2),size(Y,2));
% Regression experiments
% 100 regression problems with the same design matrix X.
fprintf('\nVarious regression experiments\n');
param.compute_gram=true;
fprintf('\nFISTA + Regression graph\n');
param.loss='square';
param.regul='graph';
tic
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
t=toc;
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));

fprintf('\nADMM + Regression graph\n');
param.admm=true;
param.lin_admm=true;
param.c=1;
param.delta=1;
tic
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
t=toc;
fprintf('mean loss: %f, stopping criterion: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));

param.admm=false;
param.max_it=5;
param.it0=1;
tic
[W optim_info]=mexFistaGraph(Y,X,W,graph,param);
t=toc;
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));

% works also with non graph-structured regularization. graph is ignored
fprintf('\nFISTA + Regression Fused-Lasso\n');
param.regul='fused-lasso';
param.lambda2=0.01;
param.lambda3=0.01; %
tic
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
t=toc;
fprintf('mean loss: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),t,mean(optim_info(4,:)));

fprintf('\nFISTA + Regression graph with intercept \n');
param.intercept=true;
param.regul='graph';
tic
[W optim_info]=mexFistaGraph(Y,[X ones(size(X,1),1)],[W0; zeros(1,size(W0,2))],graph,param);
t=toc;
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));
param.intercept=false;

% Classification
fprintf('\nOne classification experiment\n');
Y=2*double(randn(100,size(Y,2)) > 0)-1;
fprintf('\nFISTA + Logistic + graph-linf\n');
param.regul='graph';
param.loss='logistic';
param.lambda=0.01;
tic
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
t=toc;
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));
% can be used of course with other regularization functions, intercept,...

% Multi-Class classification
Y=double(ceil(5*rand(100,size(Y,2)))-1);
param.loss='multi-logistic';
param.regul='graph';
fprintf('\nFISTA + Multi-Class Logistic + graph \n');
tic
nclasses=max(Y(:))+1;
W0=zeros(size(X,2),nclasses*size(Y,2));
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
t=toc;
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));
% can be used of course with other regularization functions, intercept,...

% Multi-Task regression
Y=randn(100,size(Y,2));
Y=Y-repmat(mean(Y),[size(Y,1) 1]);
Y=mexNormalize(Y);
param.compute_gram=false;
param.verbose=true;   % verbosity, false by default
W0=zeros(size(X,2),size(Y,2));
param.loss='square';
fprintf('\nFISTA + Regression multi-task-graph \n');
param.regul='multi-task-graph';
param.lambda2=0.01;
tic
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
toc
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));

% Multi-Task Classification
fprintf('\nFISTA + Logistic + multi-task-graph \n');
param.regul='multi-task-graph';
param.lambda2=0.01;
param.loss='logistic';
Y=2*double(randn(100,size(Y,2)) > 0)-1;
tic
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
toc
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));
% Multi-Class + Multi-Task Regularization

param.verbose=false;
fprintf('\nFISTA + Multi-Class Logistic +multi-task-graph \n');
Y=double(ceil(5*rand(100,size(Y,2)))-1);
param.loss='multi-logistic';
param.regul='multi-task-graph';
tic
nclasses=max(Y(:))+1;
W0=zeros(size(X,2),nclasses*size(Y,2));
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);
t=toc;
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));
% can be used of course with other regularization functions, intercept,...

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab