text:
function f = example(varargin)
% Construct the odf from the dubna data set as example for an SO3FunRBF.

CS = crystalSymmetry('-3m',[4.9 4.9 5.4]);
SS = specimenSymmetry;

% specify file names
fname = {...
  fullfile(mtexDataPath,'PoleFigure','dubna','Q(10-10)_amp.cnv'),...
  fullfile(mtexDataPath,'PoleFigure','dubna','Q(10-11)(01-11)_amp.cnv'),...
  fullfile(mtexDataPath,'PoleFigure','dubna','Q(11-22)_amp.cnv')};

% specify crystal directions
h = {Miller(1,0,-1,0,CS),...
     [Miller(0,1,-1,1,CS),Miller(1,0,-1,1,CS)],... % superposed pole figures
     Miller(1,1,-2,2,CS)};

% specify structure coefficients
c = {1,[0.52 ,1.23],1};

% import data
pf = PoleFigure.load(fname,h,CS,SS,'interface','dubna','superposition',c);

[~,f] = evalc('calcODF(pf)');

end

author: mtex-toolbox

repo: mtex

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


text:
classdef CoherentOrientationSelector < handle

    properties (Access = private)
       mesh
       orientation
    end

    methods (Access = public)

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

        function isCoherent = isOrientationCoherent(obj)
            nnode   = obj.mesh.nnodeElem;
            nElem   = obj.mesh.nelem;
            isCoh   = false(1,nnode,nElem);
            orient  = obj.orientation.fValues;
            bN1 = squeeze(orient(:,1,:));
            for iNode = 1:nnode
                bNi = squeeze(orient(:,iNode,:));
                bN1bNI= dot(bN1,bNi);
                isCoh(1,iNode,:) = (bN1bNI)>0;
            end
            s.fValues  = isCoh;
            s.mesh     = obj.mesh;
            isCoherent = P1DiscontinuousFunction(s);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh        = cParams.mesh;
            obj.orientation = cParams.orientation;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [cl] = mi32cl(mi3)
% Convert volume from cubic miles to centiliters.
% Chad Greene 2012
cl = mi3*416818182540000;

author: Sable

repo: mcbench-benchmarks


text:
classdef HarrisLaplaceFeatureDetector < handle
    %HARRISLAPLACEFEATUREDETECTOR  Class implementing the Harris-Laplace feature detector
    %
    % As described in [Mikolajczyk2004].
    %
    % ## References
    % [Mikolajczyk2004]:
    % > Krystian Mikolajczyk and Cordelia Schmid. "Scale & affine invariant
    % > interest point detectors". International journal of computer vision,
    % > 60(1):63-86, 2004.
    %
    % See also: cv.HarrisLaplaceFeatureDetector.HarrisLaplaceFeatureDetector,
    %  cv.FeatureDetector
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = HarrisLaplaceFeatureDetector(varargin)
            %HARRISLAPLACEFEATUREDETECTOR  The full constructor
            %
            %     obj = cv.HarrisLaplaceFeatureDetector()
            %     obj = cv.HarrisLaplaceFeatureDetector('OptionName',optionValue, ...)
            %
            % ## Options
            % * __NumOctaves__ the number of octaves in the scale-space
            %   pyramid. default 6
            % * __CornThresh__ the threshold for the Harris cornerness
            %   measure. default 0.01
            % * __DOGThresh__ the threshold for the Difference-of-Gaussians
            %   scale selection. default 0.01
            % * __MaxCorners__ the maximum number of corners to consider.
            %   default 5000
            % * __NumLayers__ the number of intermediate scales per octave.
            %   default 4
            %
            % See also: cv.HarrisLaplaceFeatureDetector.detect
            %
            this.id = HarrisLaplaceFeatureDetector_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.HarrisLaplaceFeatureDetector
            %
            if isempty(this.id), return; end
            HarrisLaplaceFeatureDetector_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = HarrisLaplaceFeatureDetector_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.HarrisLaplaceFeatureDetector.empty,
            %  cv.HarrisLaplaceFeatureDetector.load
            %
            HarrisLaplaceFeatureDetector_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.HarrisLaplaceFeatureDetector.clear
            %
            b = HarrisLaplaceFeatureDetector_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.HarrisLaplaceFeatureDetector.load
            %
            HarrisLaplaceFeatureDetector_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.HarrisLaplaceFeatureDetector.save
            %
            HarrisLaplaceFeatureDetector_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.HarrisLaplaceFeatureDetector.save,
            %  cv.HarrisLaplaceFeatureDetector.load
            %
            name = HarrisLaplaceFeatureDetector_(this.id, 'getDefaultName');
        end
    end

    %% Features2D: FeatureDetector
    methods
        function keypoints = detect(this, img, varargin)
            %DETECT  Detects keypoints in an image or image set
            %
            %     keypoints = obj.detect(img)
            %     keypoints = obj.detect(imgs)
            %     [...] = obj.detect(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit grayscale image.
            % * __imgs__ Image set (second variant), cell array of images.
            %
            % ## Output
            % * __keypoints__ The detected keypoints. In the first variant, a
            %   1-by-N structure array. In the second variant of the method,
            %   `keypoints{i}` is a set of keypoints detected in `imgs{i}`.
            %
            % ## Options
            % * __Mask__ A mask specifying where to look for keypoints
            %   (optional). It must be a logical or 8-bit integer matrix with
            %   non-zero values in the region of interest. In the second
            %   variant, it is a cell-array of masks for each input image,
            %   `masks{i}` is a mask for `imgs{i}`. Not set by default.
            %
            % See also: cv.HarrisLaplaceFeatureDetector.HarrisLaplaceFeatureDetector
            %
            keypoints = HarrisLaplaceFeatureDetector_(this.id, 'detect', img, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function y = isempty( x )
y = any( size( x ) == 0 );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function ADEM_lorenz_surprise
% This demo computes the cost-function (negative reward) for a Lorenz
% system; to show cost can be computed easily from value (negative
% surprise or sojourn time). However, value is not a Lyapunov function
% because the flow is not curl-free (i.e., is not irrotational).
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_lorenz_surprise.m 4804 2012-07-26 13:14:18Z karl $

% dynamics and parameters
%-------------------------------------------------------------------------
DEMO     = 0;                          % switch for demo
LOR      = 1;                          % Lorenz vs a linear system

% generative model
%==========================================================================                       % switch for demo
G(1).E.s = 1/4;                        % smoothness
G(1).E.n = 6;                          % smoothness
G(1).E.d = 2;                          % smoothness

% dynamics and parameters
%--------------------------------------------------------------------------
if LOR
    fL  = '[v; 0; 0] + [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]*x/64';
    PL  = [10; -8/3; 32];
    x0  = [1; 1; 24];
    W   = exp(8);
else
    fL  = '[v; 0; 0] + [-P(1) 0 0; 0 -P(2) 0; 0 0 -P(3)]*x/64';
    PL  = [1; 1; 8];
    x0  = [-16; 16; 0];
    W   = diag([32; 32; exp(8)]);
end

% level 1
%--------------------------------------------------------------------------
G(1).x  = x0;
G(1).f  = inline(fL ,'x','v','P');
G(1).g  = inline('x','x','v','P');
G(1).pE = PL;
G(1).V  = exp(8);                           % error precision
G(1).W  = W;                                % error precision

% level 2
%--------------------------------------------------------------------------
G(2).a  = [0;0;0];                          % action variables
G(2).v  = 0;                                % inputs
G(2).V  = exp(16);
G       = spm_DEM_M_set(G);

% space
%--------------------------------------------------------------------------
N        = 66;
if LOR
    x{1} = linspace(-32,32,N);
    x{2} = linspace(-32,32,N);
    x{3} = linspace(  4,64,N);
else
    x{1} = linspace(-8,8,N);
    x{2} = linspace(-8,8,N);
    x{3} = linspace(-1,1,3);
end

% equilibrium density (q0), loss (L) and value (V) functions
%==========================================================================
if DEMO

    % Fokker-Planck operator and equilibrium density
    %----------------------------------------------------------------------
    [M0,q0] = spm_fp(G,x);

    % loss-function and negative surprise (value)
    %----------------------------------------------------------------------
    V    = log(q0);
    L    = spm_unvec(spm_vec(V)'*M0,q0);

    % trim
    %----------------------------------------------------------------------
    q0   = q0(2:end - 1,2:end - 1,2:end - 1);
    L    =  L(2:end - 1,2:end - 1,2:end - 1);
    V    =  V(2:end - 1,2:end - 1,2:end - 1);
    x{1} = x{1}(2:end - 1);
    x{2} = x{2}(2:end - 1);
    x{3} = x{3}(2:end - 1);

    if LOR
        save DEM_lorenz_suprise q0 L V x
    end
else
    load DEM_lorenz_suprise
end

% axes and trajectory
%--------------------------------------------------------------------------
spm_figure('GetWin','DEM');

if LOR, z = 24; else, z = 1; end

i    = 3;
j    = 1:3;
j(i) = [];
T    = 1024;
U.u  = sparse(T,G(1).m);
t    = spm_int_J(G(1).pE,G,U);

% surprise
%--------------------------------------------------------------------------
a    = [x{j(2)}(1) x{j(2)}(end) x{j(1)}(1) x{j(1)}(end)];
subplot(3,2,1)
imagesc(x{j(2)},x{j(1)},V(:,:,z))
hold on, plot(t(:,j(2)),t(:,j(1)),'r'), hold off
axis(a)
axis square xy
title('value','Fontsize',16)

% cost function
%--------------------------------------------------------------------------
subplot(3,2,2)
imagesc(x{j(2)},x{j(1)},L(:,:,z))
hold on, plot(t(:,j(2)),t(:,j(1)),'r'), hold off
axis(a)
axis square xy
title('cost','Fontsize',16)

% equilibrium density
%--------------------------------------------------------------------------
subplot(3,2,3)
imagesc(x{j(2)},x{j(1)},squeeze(max(q0,[],i)))
hold on, plot(t(:,j(2)),t(:,j(1)),'r'), hold off
axis(a)
axis square xy
title('density','Fontsize',16)

% exemplar trajectory
%--------------------------------------------------------------------------
subplot(3,2,4)
plot(t(:,j(2)),t(:,j(1)),'k')
axis(a)
axis square xy
title('trajectory','Fontsize',16)

% evaluate V(t)
%--------------------------------------------------------------------------
for i = 1:T
    [q ii] = min(abs(t(i,1) - x{1}));
    [q ij] = min(abs(t(i,2) - x{2}));
    [q ik] = min(abs(t(i,3) - x{3}));
    v(i)   = V(ii,ij,ik);
end

% and plot
%--------------------------------------------------------------------------
maxV = max(V(:));
subplot(3,1,3)
plot(v,'k'), hold on
plot([1 T],[maxV maxV],'r:'), hold off
title('value','Fontsize',16)
set(gca,'XLim',[1 T]);
box off

author: spm


text:
function [m, state] = str2num (s)
% Uses `try eval(x); catch ...; end instead of `eval(x,y)`, which seems to be buggy
  if (nargin != 1)
    print_usage ();
  elseif (! ischar (s))
    error ("str2num: S must be a string or string array");
  endif
  s(:, end+1) = ";";
  s = sprintf ("m = [%s];", reshape (s', 1, numel (s)));
  state = true;
  try eval(s); catch m = []; state = false; end
  if (ischar (m))
    m = [];
    state = false;
  endif
endfunction

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [ r, s, area ] = node_reference_ql ( )

%*****************************************************************************80
%
%% NODE_REFERENCE_QL returns the basis nodes for a quadratic/linear.
%
%  Reference Element QL:
%
%    |
%    1  4---5---6
%    |  |       |
%    |  |       |
%    S  |       |
%    |  |       |
%    |  |       |
%    0  1---2---3
%    |
%    +--0---R---1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real R(6), S(6), the coordinates of the basis nodes.
%
%    Output, real AREA, the area of the element.
%
  r(1:6) = [ 0.0, 0.5, 1.0, 0.0, 0.5, 1.0 ];
  s(1:6) = [ 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 ];

  area = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [bcUx, bcUy, bcPR, bcdUndt] = KovasznayBC2D(x, y, nx, ny, mapI, mapO, mapW, mapC, time, nu)

% function [bcUx, bcUy, bcPR, bcdUndt] = KovasznayBC2D(x, y, nx, ny, mapI, mapO, mapW, mapC, time, nu)
% Purpose: evaluate boundary conditions for Kovasznay flow

zer = zeros(size(x)); bcUx = zer; bcUy = zer; bcPR = zer; bcdUndt = zer;

lam = (0.5/nu) - sqrt( (0.25/(nu^2)) + 4*pi^2 );

% inflow
xI = x(mapI); yI = y(mapI);
bcUx(mapI)= 1-exp(lam*xI).*cos(2*pi*yI);
bcUy(mapI)= (0.5*lam/pi)*exp(lam*xI).*sin(2*pi*yI);

% outflow
xO = x(mapO); yO = y(mapO);

if(0)
  bcPR(mapO) = .5*(1-exp(2*lam*xO));
  bcUx(mapO)= 1-exp(lam*xO).*cos(2*pi*yO);
  bcUy(mapO)= (0.5*lam/pi)*exp(lam*xO).*sin(2*pi*yO);
else
  bcPR(mapO) = .5*(1-exp(2*lam*xO));

  % Neumann data for each velocity
  bcUx(mapO) = -lam*exp(lam*xO).*cos(2*pi*yO);
  bcUy(mapO) = lam*(0.5*lam/pi)*exp(lam*xO).*sin(2*pi*yO);
end
return

repo: Aero-matlab


text:
function ps = normpdf2( xs, m, C )
% Normal prob. density function (pdf) with arbitrary covariance matrix.
%
% Evaluate the multi-variate density with mean vector m and covariance
% matrix C for the input vector xs.  Assumes that the N datapoints are d
% dimensional.  Then m is dx1 or 1xd, C is dxd, and xs is dxN or NxD where
% N is the number of samples to be evaluated.
%
% USAGE
%  ps = normpdf2( xs, m, C )
%
% INPUTS
%  xs  - points to evaluated (Nxd or dxN)
%  m   - mean vector (dx1 or 1xd)
%  C   - Covariance matrix (dxd)
%
% OUTPUTS
%  ps  - probability density at each x (Nx1)
%
% EXAMPLE
%  ps = normpdf2( randn(10,2), [0 0], eye(2) )
%
% See also NORMPDF
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.30
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get dimensions of data
d=length(m);
if( size(xs,1)~=d ); xs=xs'; end
N=size(xs,2);

if( d==1 ) % fast special case
  ps = 1/sqrt(2*pi*C) * exp(-(xs-m).*(xs-m)/(2*C))';

elseif( rcond(C)<eps ) % if matrix is badly conditioned
  warning('normpdf2: Covariance matrix close to singular.'); %#ok<WNTAG>
  ps = zeros(N,1);

else % get probabilities
  xs = (xs-m(:)*ones(1,N))';
  denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
  mahal = sum( (xs/C).*xs, 2 );
  numer = exp(-0.5*mahal);
  ps = numer/denom;
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function vout=nextchoose(v,w)
%NEXTCHOOSE  find the next choice of k elements from a vector (with possible repeats)
% usage: vout = nextchoose(v,w)
%
% NEXTCHOOSE(V,W), where V is a row or column vector of integers in W,
% returns the next choice of length(V) elements from W by lexical ordering
% (see below) or [] if there is none.  W can also be specified as a
% positive integer N, in which case W=1:N.
%
% NEXTCHOOSE can be used to iterate through all unique k-combinations in
% 1:N but keep in mind that the number is nchoosek(N,length(V))
%
% In the set of all k-combinations of W, with lexical order:
% - V is equivalent to sort(V)
% - pairs are ordered according to the first sorted elements that differ
% - the "first" choice is V0=W(1:k) (if W is sorted)
% - the "last" choice is V1=W(N-k+1:N) (if W is sorted, and N=numel(W))
% - the successor of V is determined according to the rightmost elements
%   such that V(i)=W(j) and j < N-(k-i). Then V(i:k) -> W(j+1:j+k-i+1)
%
% Example 1:
%   nextchoose([1,2,3],4)       % returns [1,2,4]
%   nextchoose([1,2,4],4)       % returns [1,3,4]
%   nextchoose([1,3,4],4)       % returns [2,3,4]
%   nextchoose([2,3,4],4)       % returns []
%
% Example 2: choose 2 elements from a vector with repeats
%   nextchoose([1,2],[1,2,2,3]) % returns [1,3]
%   nextchoose([1,3],[1,2,2,3]) % returns [2,2]
%   nextchoose([2,2],[1,2,2,3]) % returns [2,3]
%   nextchoose([2,3],[1,2,2,3]) % returns []
%
% Example 3: determine all 4-combinations of 1:8
%   n=8;
%   k=4;
%   v=1:k;
%   M=zeros(nchoosek(n,k),k);
%   i=0;
%   while ~isempty(v),
%     i=i+1;
%     M(i,:)=v;
%     v=nextchoose(v,n);
%   end;
%   isequal(M,sortrows(nchoosek(1:n,k))) % =1

% Author: Ben Petschel 25/6/2009
%
% Change history:
%  25/6/2009 - first release (modelled after nextperm.m)
%  29/6/2009 - added support for general choice vectors

vout = sort(v);
k = length(v);

if isscalar(w),
  W = 1:w;
else
  W = sort(w);
end;

% find the rightmost element that can be increased
ind = find(vout<W(end-k+1:end),1,'last');

if isempty(ind),
  % v is the last vector
  vout = [];
else
  % increment vout(ind) and update the tail elements
  j = find(vout(ind)==W,1,'last');
  if isempty(j),
    error('an element of V is not in W');
  end;
  vout(ind:k) = W(j+1:j+k-ind+1);
end;


text:
function test_failed = test_libltfat_fftshift(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

Larr = [1,9,11,110];

for do_complex = 0:1
    complexstring = '';
    if do_complex, complexstring = 'complex'; end
    funname = makelibraryname('fftshift',flags.complexity,do_complex);

    for L = Larr

            if do_complex
                z = cast((1:L)'+1i*(L:-1:1)',flags.complexity);
                zi = complex2interleaved(z);
                zout = randn(size(zi),flags.complexity);

                ziPtr = libpointer(dataPtr,zi);
                zoutPtr = libpointer(dataPtr,zout);
            else
                z = cast((1:L)',flags.complexity);
                zi = z;
                zout = randn(size(zi),flags.complexity);

                ziPtr = libpointer(dataPtr,zi);
                zoutPtr = libpointer(dataPtr,zout);
            end

            trueres = fftshift(z);

            status = calllib('libltfat',funname,ziPtr,L,zoutPtr);

            if do_complex
                res = norm(trueres - interleaved2complex(zoutPtr.Value));
            else
                res = norm(trueres - zoutPtr.Value);
            end

            [test_failed,fail]=ltfatdiditfail(res+status,test_failed,0);
            fprintf(['FFTSHIFT OP L:%3i, %s %s %s %s\n'],L,flags.complexity,complexstring,ltfatstatusstring(status),fail);

            status = calllib('libltfat',funname,ziPtr,L,ziPtr);

            if do_complex
                res = norm(trueres - interleaved2complex(ziPtr.Value));
            else
                res = norm(trueres - ziPtr.Value);
            end

            [test_failed,fail]=ltfatdiditfail(res+status,test_failed,0);
            fprintf(['FFTSHIFT IP L:%3i, %s %s %s %s\n'],L,flags.complexity,complexstring,ltfatstatusstring(status),fail);
    end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% mex CFLAGS="-O3 -march=nocona -Wall" -largeArrayDims mexSOSVMLearn.cpp LaRank.cpp Config.cpp
mex CFLAGS="-O3 -march=nocona -Wall" mexSOSVMLearn.cpp LaRank.cpp Config.cpp

author: flyers

repo: drone-tracking


text:
function model = modelLoadResult(type, dataSet, number, dataLoaderStr)

% MODELLOADRESULT Load a previously saved result.
% FORMAT
% DESC loads a previously saved model result.
% ARG type : the type of model to load.
% ARG dataSet : the name of the data set to load.
% ARG number : the number of the model run to load.
% RETURN model : the saved model.
%
% SEEALSO : modelWriteResult
%
% COPYRIGHT : Neil D. Lawrence, 2009

% MLTOOLS
  if nargin < 4
    dataLoaderStr = 'lvmLoadData';
  end
  dataLoaderHandle = str2func(dataLoaderStr);
  origDataSet = dataSet;
  dataSet(1) = upper(dataSet(1));
  origType = type;

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

  fhandle = [origType 'LoadResult'];
  if exist(fhandle)==2
    % There is load result code, use it.
    fhandle = str2func(fhandle);
    if nargin > 3
      model = fhandle(origDataSet, number, dataLoaderStr);
    else
      model = fhandle(origDataSet, number);
    end
  else
    fhandle = [origType 'Reconstruct'];
    if exist(fhandle)==2
      % There is reconstruct code, use it to reconstruct.
      [Y, lbls] = dataLoaderHandle(origDataSet);
      load(fileName);
      varName = [origType 'Info'];
      eval(['model = ' fhandle '(' varName ', Y);'])
    else
      % No code is provided, just load the file.
      load(fileName);
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef computingOptimal < handle

    methods (Access = public)

        function obj = computingOptimal()
            obj.init()
        end

    end

    methods (Access = private)

        function init(obj)

        end

        function obj = computingOpt()
            fileName = 'OptimalSuperEllipse';
            outputFolder = fullfile(pwd,'Output',fileName);
            gmsFile = [fullfile(outputFolder,fileName),'.msh'];
            createMesh(fileName,outputFolder);
            homog = obj.createNumericalHomogenizer(fileName,gmsFile);
            maxSnorm = obj.computeMaxStressNorm(homog);
        end

        function maxSnorm = computeMaxStressNorm(obj,homog)
            Ch = homog.cellVariables.Ch;
            microProblem = homog.getMicroProblem();
            stress = [cos(pi/4) sin(pi/4) 0]';
            strain = Ch\stress;
            microProblem.element.setVstrain(strain');
            microProblem.computeVariables;
            stresses = microProblem.variables.stress;
            sx = squeeze(stresses(:,1,:));
            sy = squeeze(stresses(:,2,:));
            sxy = squeeze(stresses(:,3,:));
            sNorm = sqrt(sx.*sx + 2*sxy.*sxy + sy.*sy);
            maxSnorm = max(sNorm);
        end

        function createMesh(obj,fileName,outputFolder)
            d = SettingsFreeFemMeshGenerator();
            d.freeFemFileName = 'SmoothRectangle';
            d.hMax  = 0.02;%0.0025;
            d.mxV             = 0.5;
            d.myV             = 0.5;
            d.fileName        = fileName;
            d.printingDir     = outputFolder;
            d.qNorm           = 32;

            fG = FreeFemMeshGenerator(d);
            fG.generate();

        end

        function homog = createNumericalHomogenizer(obj,fileName,gmsFile)
            d.gmsFile = gmsFile;
            d.outFile = fileName;
            d.print   = false;
            d.iter = 0;
            nH = NumericalHomogenizerCreatorFromGmsFile(d);
            homog = nH.getHomogenizer();
        end

    end

end


text:
function [oldvariable,sqrList] = expandforindex(sqrList,candidates,nonconvindex);
%EXPANDFORINDEX Internal function for organizing nonlinear variables

% Expand the incoming variables

if 1
z = sqrList(:,1);

x1_list = find(z==candidates(2));
if candidates(2)==candidates(1)
    x2_list = x1_list;
else
    x2_list = find(z==candidates(1));
end

temp = [];
if isempty(x1_list)
    temp = candidates(2);
else
    temp = sqrList(x1_list,2:end);
end
if isempty(x2_list)
    temp = [temp candidates(1)];
else
    temp = [temp sqrList(x2_list,2:end)];
end
temp = sort(temp(find(any(temp,1))));

oldsqrList = sqrList;
if length(temp)>size(sqrList,2)-1
    % This one cannot already exist
    index = [];

    %sqrList = [sqrList zeros(size(sqrList,1),1+size(sqrList,2)-length(temp));nonconvindex temp(end:-1:1)];
    sqrList = [sqrList zeros(size(sqrList,1),1+length(temp)-size(sqrList,2));nonconvindex temp(end:-1:1)];
else
%    sqrList = [sqrList;nonconvindex temp(end:-1:1) zeros(1,size(sqrList,2)-length(temp)-1)];
%    searchfor = sqrList(end,2:end);

    searchfor = [temp(end:-1:1) zeros(1,size(sqrList,2)-length(temp)-1)];
    searchin = sqrList(:,2:end);

    % Simple Hash
    key = sum(searchfor);
    tbl = sum(searchin,2);

    possible = find(tbl==key);

    index = findrows(searchin(possible,:),searchfor);
    if isempty(index)
         sqrList = [sqrList;nonconvindex temp(end:-1:1) zeros(1,size(sqrList,2)-length(temp)-1)];
    else
         index = possible(index);
    end
end

%searchfor = sqrList(end,2:end);
%index = findrows(sqrList(1:end-1,2:end),searchfor);

if length(index)>1
    index
end

if ~isempty(index)
    oldvariable = sqrList(index,1);
    sqrList = oldsqrList;
else
    oldvariable = [];
end

else

% Expand the list to begin with
if isempty(sqrList)
    sqrList = [nonconvindex candidates];
else
    sqrList = [sqrList; nonconvindex candidates(2) candidates(1) zeros(1,size(sqrList,2)-3)];
end
bottom = size(sqrList,1);

%%FIXz = sqrList(1:bottom,1);
z = sqrList(:,1);

x1_list = find(z==candidates(2));
if candidates(2)==candidates(1)
    x2_list = x1_list;
else
    x2_list = find(z==candidates(1));
end

temp = [];
if isempty(x1_list)
    temp = sqrList(bottom,2);
else
    temp = sqrList(x1_list,2:end);
end
if isempty(x2_list)
    temp = [temp sqrList(bottom,3)];
else
    temp = [temp sqrList(x2_list,2:end)];
end
temp = sort(temp(find(any(temp,1))));
sqrList(bottom,2:2+length(temp)-1)=temp(end:-1:1);%fliplr(sort(temp));

bottom = bottom+1;

searchfor = sqrList(end,2:end);
index = findrows(sqrList(1:end-1,2:end),searchfor);

if ~isempty(index)
    oldvariable = sqrList(index,1);
    sqrList = sqrList(1:end-1,:);
else
    oldvariable = [];
end

end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% pdf_Sobol: Foo function for simulate a Sobol Set
%
% Usage:
%   pdf_Sobol()
%
% Inputs:
%     range            vector [min max] range of the random variable
%
% Output:
%     range            vector [min max] range of the random variable
% ------------------------------------------------------------------------
% See also
%
% Author : Flavio Cannavo'
% e-mail: flavio(dot)cannavo(at)gmail(dot)com
% Release: 1.0
% Date   : 01-02-2011
%
% History:
% 1.0  01-02-2011  First release.
%%

function range = pdf_Sobol(range)

range = range(:)';
% empty function


text:
 Example due to Wang Hee Lin" <engp1622@nus.edu.sg

intra = zeros(2);
intra(1,2) = 1;
inter = zeros(2);
inter(1,1) = 1;

Q = 2; % num hidden states
O = 2; % num observable symbols
ns = [Q O];%number of states
dnodes = 1:2;
%onodes = [1:2]; % only possible with jtree, not hmm
onodes = [2];
bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'observed', onodes);
for i=1:4
  bnet.CPD{i} = tabular_CPD(bnet, i);
end

prior0 = normalise(rand(Q,1));
transmat0 = mk_stochastic(rand(Q,Q));
obsmat0 = mk_stochastic(rand(Q,O));

%engine = smoother_engine(hmm_2TBN_inf_engine(bnet));
engine = smoother_engine(jtree_2TBN_inf_engine(bnet));

ss = 2;%slice size(ss)
ncases = 10;%number of examples
T=10;
max_iter=2;%iterations for EM
cases = cell(1, ncases);
for i=1:ncases
  ev = sample_dbn(bnet, T);
  cases{i} = cell(ss,T);
  cases{i}(onodes,:) = ev(onodes, :);
end
[bnet2, LLtrace] = learn_params_dbn_em(engine, cases, 'max_iter', 4);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [trl, event] = ft_trialfun_hed(cfg)

% FT_TRIALFUN_HED is a trial function that can be used with HED tags. It demonstrates
% some basic functionality for selecting specific events, but mainly serves as an
% example or template for your own conditial trial definitions. For that you would
% copy this function and giuve it your own name, e.g. FT_TRIALFUN_MYEXPERIMENT.
%
% Use this function by calling
%   [cfg] = ft_definetrial(cfg)
% where the configuration structure should contain
%   cfg.dataset           = string with the filename
%   cfg.trialfun          = 'ft_trialfun_hed' % or your own copy
%
% The selection of events and timing of the epochs is specified with
%   cfg.trialdef.regexp     = regular expression that is applied to the HED tags
%   cfg.trialdef.prestim    = number, in seconds
%   cfg.trialdef.poststim   = number, in seconds
%
% See also FT_DEFINETRIAL, FT_TRIALFUN_GENERAL, FT_TRIALFUN_EXAMPLE1,
% FT_TRIALFUN_EXAMPLE2

% read the header information
hdr = ft_read_header(cfg.dataset);

% determine the number of samples before and after the event
prestim  = -round(cfg.trialdef.prestim  * hdr.Fs);
poststim =  round(cfg.trialdef.poststim * hdr.Fs);

% rather than using FT_READ_EVENT which has a limited representation, we assume that
% the data is formatted according to BIDS and includes HED tags

tsvfile  = bids_sidecar(cfg.dataset, 'events', 'tsv');
jsonfile = bids_sidecar(cfg.dataset, 'events', 'json');

tsv = ft_read_tsv(tsvfile);
if ~isempty(jsonfile)
  json = ft_read_json(jsonfile);
else
  json = [];
end

% this will be returned as the complete set of events without selection
event = tsv;

% make a selection based on the HED tags
if any(strcmp(fieldnames(tsv), 'HED'))
  hedtags = tsv.HED;
elseif ~isempty(jsonfile)
  warning('assembling the HED tags on the fly');
  % this requires functions from an external toolbox, see https://hed-examples.readthedocs.io/en/latest/HedMatlabTools.html
  ft_hastoolbox('hedtools', 1);
  tsv = hed_assemble(tsvfile, jsonfile);
  hedtags = tsv.HED;
else
  error('there are no HED tags')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the following code uses a regular expression to select events
% you may want to change this and use another mechanism to define the trials/epochs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% find the HED tags that match the regular expression
sel = find(~cellfun(@isempty, regexp(hedtags, cfg.trialdef.regexp)));

% these are the first three required columnss
trl = table;
trl.begsample = tsv.sample(sel) + prestim;
trl.endsample = tsv.sample(sel) + poststim - 1;
trl.offset = repmat(prestim, size(trl.begsample));

% concatenate all other details to the table
trl = cat(2, trl, tsv(sel,:));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function c = gsp_lanczos_op(G,fi,s,param)
%GSP_LANCZOS_OP Perform the lanczos approximation of the signal s

if nargin < 4, param = struct; end
if ~isfield(param,'verbose'), param.verbose = 1; end;
if ~isfield(param,'order'), param.order = 30; end

Nf = numel(fi);
Nv = size(s,2);
c = zeros(G.N*Nf,Nv);

for jj = 1:Nv

    if sum(abs(s(:,jj)))>eps
        [V,H] = lanczos(G.L, param.order, s(:,jj));

        [Uh, Eh] = eig(H);

        Eh = diag(Eh);
        Eh(Eh<0) = 0;
        fie = gsp_filter_evaluate(fi,Eh);
        V = V*Uh;

        for ii=1:Nf
           c((1:G.N) + G.N*(ii-1),jj) = V * (fie(:, ii) .* (V'*s(:,jj)));
        end
    else
        c(:,jj) = 0;
    end
end

% [V,H] = lanczos(G.L, k, s);
%
% for jj = 1:Nv
%     ind = (1:size(H,1))+size(H,1)*(jj-1);
%     [Uh, Eh] = eig(H(:,ind));
%
%
%     Eh = diag(Eh);
%     fie = gsp_filter_evaluate(fi,Eh);
%     V(:,ind) = V(:,ind)*Uh;
%
%
%     for ii=1:Nf
%        c((1:G.N) + G.N*(ii-1),jj) = V(:,ind) * fie(:, ii) .* (V(:,ind)'*s(:,jj));
%     end
% end

end

function [V,H,orth] = lanczos(A,order,x)

[N,M] = size(x);

% normalization
norm2vec = @(x) (sum(x.^2,1)).^0.5;
q = x./repmat(norm2vec(x),N,1);

% Initialization
hiv =0:order:(order*M-1); % helping indice vector

V = zeros(N,M*order);
V(:,1+hiv) = q;

H = zeros(order+1,order*M);

r = A*q;
H(1,1+hiv) = sum(q .* r, 1 );
r = r - repmat(H(1,1+hiv),N,1).*q;
H(2,1+hiv) = norm2vec(r);

if (nargout > 2)
    orth = zeros(M,1);
    orth(1) = norm(V'*V - M);
end

for k = 2:order

    if (sum(abs(H(k,k-1+hiv))) <= eps)
        H = H(1:k-1,sum_ind(1:k-1,hiv));
        V = V(:,sum_ind(1:k-1,hiv));
        if (nargout > 2)
            orth = orth(1:k-1);
        end
        return;
    end

    H(k-1,hiv+k) = H(k,hiv+k-1);
    v = q;
    q = r./repmat(H(k-1,k+hiv),N,1);
    V(:,k+hiv) = q;

    r = A*q;
    r = r - repmat(H(k-1,k+hiv),N,1).*v;
    H(k,k+hiv) = sum(q .* r, 1 );

    r = r - repmat(H(k,k+hiv),N,1).*q;
    % The next line has to be checked
    r = r - V*(V'*r); % full reorthogonalization
    H(k+1,k+hiv) = norm2vec(r);

    if (nargout > 2)
        orth(k) = [orth, norm(V'*V - M)];
    end
end

H = H(1:order,1:order);

% H = zeros(order+1,order);
%
% q = x/norm(x);
% V(:,1) = q;
%
% r = A*q;
%
% H(1,1) = q'*r;
% r = r - H(1,1)*q ;
% H(2,1) = norm(r);
%
% if (nargout > 2)
%  orth = norm(V'*V - 1);
% end
%
% for k = 2:order
%
%     if (abs(H(k,k-1)) <= 1e-15)
%         H = H(1:k-1,1:k-1);
%         return;
%     end
%
%     H(k-1,k) = H(k,k-1);
%     v = q;
%     q = r/H(k-1,k);
%     V(:,k) = q;
%
%     r = A*q;
%     r = r - H(k-1,k)*v;
%     H(k,k) = q'*r;
%
%     r = r - H(k,k)*q;
%     r = r - V*(V'*r); % full reorthogonalization
%     H(k+1,k) = norm(r);
%
%     if (nargout > 2)
%         orth = [orth, norm(V'*V - eye(k))];
%     end
% end
%
% H = H(1:order,1:order);
%
end

function ind = sum_ind(ind1,ind2)
    ind = repmat(ind1(:),1,numel(ind2)) + repmat((ind2(:))',numel(ind1),1);
    ind = (ind(:))';
end

repo: gspbox


text:
% PCA on Fish data
tic
[coeff,score,latent,tsquared,explained,mu] = pca(M');
toc
code_dir = GetCurrentCodeDir();
save(fullfile(code_dir,'\saved mats\Fish3_full100_PCA.mat'),'coeff','score','latent','tsquared','explained','mu');

cIX_full = cIX;
gIX_full = gIX;

%%
cIX = cIX_full;
gIX = gIX_full;

%%
nPCs = 200;
figure;
im = coeff; % % rows reordered with optimal-leaf-order from dendrogram
% im = im.*repmat(explained',size(coeff,1),1); % PC's weighted by percent-explained
imagesc(im(:,1:nPCs));
% axis equal;
% set(gca,'YTick',1:size(coeff,1),'YTickLabel',ylabels);
xlim([0.5,nPCs+0.5])
xlabel('PC''s')

%% screen for cells that have high PC scores

coeff_crp = coeff(:,1:nPCs); % cropped
cellscore = abs(sum(coeff_crp,2));
[A,IX] = sort(cellscore,'descend');

percCell_PC = 0.1;
numCell = size(M,1);
nCell_PC = round(numCell*percCell_PC);
cIX = cIX(IX(1:nCell_PC));
gIX = ceil((1:length(cIX))'/length(cIX)*min(20,length(cIX)));

author: xiuyechen

repo: FishExplorer


text:
function [xmlStruct, map] = transformXML2Map(fileName)
% Create a MATLAB structure from a given XML file.
% The XML file is first parsed through the "xml2struct" function and
% then transformed into a structure. The content of this structure can be
% found in the [description document](https://github.com/opencobra/cobratoolbox/blob/master/docs/source/notes/simpleMATLABStructure.md)
%
% USAGE:
%
%   [xmlStruct, map] = transformXML2Map(fileName)
%
% INPUT:
%   fileName:       Path to the XML file.
%
% OUTPUTS:
%   xmlStruct:      Structure obtained from the "xml2struct" function.
%                   To be kept for the conversion back to an XML file
%                   of the structure.
%   map:            Matlab structure of the map containing all the
%                   relevant fields usable for checking and correction.
%
%                * sID -  Stoichiometric matrix with `rows = MetabolitesID` and
%                  `columns = ReactionsID` in the same order as in the map
%                  structure. Contains `-1` if the metabolite is a
%                  reactant/substract, `+1` if the metabolite is a product
%                  and `0` if it does not participate in the reaction.
%                * sAlias - Stoichiometric matrix with `rows = MetabolitesAlias` and
%                  `columns = ReactionsID` in the same order as in the map
%                  structure. Contains `-1` if the metabolite is a
%                  reactant/substract, `+1` if the metabolite is a product
%                  and `0` if it does not participate in the reaction.
%                * idAlias - Logical matrix with `rows = MetabolitesID` and
%                  `columns = MetabolitesAlias`. Contains `+1` if the
%                  `MetaboliteID` match with the `MetaboliteAlias` and `0`
%                  if it doesn't.
%
% .. Author: - N.Sompairac - Institut Curie, Paris, 24/07/2017

tic

% Works nicely but has a huge tree structured

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

xmlStruct = xml2struct(fileName);

% Loop over molecules to get the needed information and store it in a
% structure. Molecules refer to each individual node.
for mol = 1:length(xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias)
    map.molAlias{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.Attributes.id;
    map.molID{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.Attributes.species;
    % Check if the info on compartment exists
    if any(strcmp('compartmentAlias', fieldnames(xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.Attributes)))
        map.molCompartAlias{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.Attributes.compartmentAlias;
    else
        map.molCompartAlias{mol, 1} = '';
    end
    map.molXPos{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.x;
    map.molYPos{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.y;
    map.molWidth{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.w;
    map.molHeight{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_bounds.Attributes.h;
    map.molColor{mol, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfSpeciesAliases.celldesigner_colon_speciesAlias{mol}.celldesigner_colon_usualView.celldesigner_colon_paint.Attributes.color;
end
clearvars mol

% Loop over species to get the needed information and store it in a
% structure. Species refer to each individual name.
for spec = 1:length(xmlStruct.sbml.model.listOfSpecies.species)
    map.specID{spec, 1} = xmlStruct.sbml.model.listOfSpecies.species{spec}.Attributes.id;
    map.specMetaID{spec, 1} = xmlStruct.sbml.model.listOfSpecies.species{spec}.Attributes.metaid;
    map.specName{spec, 1} = xmlStruct.sbml.model.listOfSpecies.species{spec}.Attributes.name;
    map.specType{spec, 1} = xmlStruct.sbml.model.listOfSpecies.species{spec}.annotation.celldesigner_colon_extension.celldesigner_colon_speciesIdentity.celldesigner_colon_class.Text;
    % Check if the info on notes exists for the specie
    if any(strcmp('notes', fieldnames(xmlStruct.sbml.model.listOfSpecies.species{spec})))
        map.specNotes{spec, 1} = xmlStruct.sbml.model.listOfSpecies.species{spec}.notes.html.body.Text;
    else
        map.specNotes{spec, 1} = '';
    end
end
clearvars spec

% Loop over reactions to get the need information and store it in a
% structure. Reactions refer to each individual reaction name.
for react = 1:length(xmlStruct.sbml.model.listOfReactions.reaction)
    % Test if the reaction has an ID
    if any(strcmp('id', fieldnames(xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes)))
        map.rxnID{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.id;
    else
        map.rxnID{react, 1} = '';
    end
    map.rxnMetaID{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.metaid;
    % Test if the reaction has a Name
    if any(strcmp('name', fieldnames(xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes)))
        map.rxnName{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.name;
    else
        map.rxnName{react, 1} = '';
    end
    map.rxnType{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_reactionType.Text;
    % Check if info on reversibility exists
    if any(strcmp('reversible', fieldnames(xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes)))
        map.rxnReversibility{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.Attributes.reversible;
    else
        map.rxnReversibility{react, 1} = 'true';
    end
    map.rxnBaseReactantAlias{react, 1} = {};
    map.rxnBaseReactantID{react, 1} = {};
    % Test if there is only 1 base reactant
    if length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant) == 1
        map.rxnBaseReactantAlias{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant.Attributes.alias;
        map.rxnBaseReactantID{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant.Attributes.species;
    else
        % Looping over the multiple base reactants
        for base = 1:length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant)
            map.rxnBaseReactantAlias{react, 1}{base, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant{base}.Attributes.alias;
            map.rxnBaseReactantID{react, 1}{base, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseReactants.celldesigner_colon_baseReactant{base}.Attributes.species;
        end
        clearvars base
    end
    map.rxnBaseProductAlias{react, 1} = {};
    map.rxnBaseProductID{react, 1} = {};
    % Test if there is only 1 base product
    if length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct) == 1
        map.rxnBaseProductAlias{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct.Attributes.alias;
        map.rxnBaseProductID{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct.Attributes.species;
    else
        % Looping over the multiple base products
        for base = 1:length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct)
            map.rxnBaseProductAlias{react, 1}{base, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct{base}.Attributes.alias;
            map.rxnBaseProductID{react, 1}{base, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_baseProducts.celldesigner_colon_baseProduct{base}.Attributes.species;
        end
        clearvars base
    end
    % Test if there are some reactants in the reaction
    if any(strcmp('celldesigner_colon_listOfReactantLinks', fieldnames(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension)))
        map.rxnReactantAlias{react, 1} = {};
        map.rxnReactantID{react, 1} = {};
        % Test if there is only 1 reactant
        if length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink) == 1
            map.rxnReactantAlias{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(1).Attributes.alias;
            map.rxnReactantID{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink(1).Attributes.reactant;
        else
            for x = 1:length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink)
                map.rxnReactantAlias{react, 1}{x, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{x}.Attributes.alias;
                map.rxnReactantID{react, 1}{x, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfReactantLinks.celldesigner_colon_reactantLink{x}.Attributes.reactant;
            end
            clearvars x
        end
    else
        map.rxnReactantAlias{react, 1} = '';
        map.rxnReactantID{react, 1} = '';
    end
    % Test if there are some products in the reaction
    if any(strcmp('celldesigner_colon_listOfProductLinks', fieldnames(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension)))
        map.rxnProductAlias{react, 1} = {};
        map.rxnProductID{react, 1} = {};
        % Test if there is only 1 product
        if length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink) == 1
            map.rxnProductAlias{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(1).Attributes.alias;
            map.rxnProductID{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink(1).Attributes.product;
        else
            for x = 1:length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink)
                map.rxnProductAlias{react, 1}{x, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{x}.Attributes.alias;
                map.rxnProductID{react, 1}{x, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfProductLinks.celldesigner_colon_productLink{x}.Attributes.product;
            end
            clearvars x
        end
    else
        map.rxnProductAlias{react, 1} = '';
        map.rxnProductID{react, 1} = '';
    end
    % Check if there are modifications for the reaction
    if any(strcmp('listOfModifiers', fieldnames(xmlStruct.sbml.model.listOfReactions.reaction{react})))
        % Test if there is only one modification
        if length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification) == 1
            map.rxnModAlias{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).Attributes.aliases;
            map.rxnModID{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).Attributes.modifiers;
            map.rxnModType{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).Attributes.type;
            map.rxnModColor{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).celldesigner_colon_line.Attributes.color;
            map.rxnModWidth{react, 1}{1, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification(1).celldesigner_colon_line.Attributes.width;
        else
            % Loop over the possible modifications
            for mod = 1:length(xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification)
                map.rxnModAlias{react, 1}{mod, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{mod}.Attributes.aliases;
                map.rxnModID{react, 1}{mod, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{mod}.Attributes.modifiers;
                map.rxnModType{react, 1}{mod, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{mod}.Attributes.type;
                map.rxnModColor{react, 1}{mod, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{mod}.celldesigner_colon_line.Attributes.color;
                map.rxnModWidth{react, 1}{mod, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_listOfModification.celldesigner_colon_modification{mod}.celldesigner_colon_line.Attributes.width;
            end
            clearvars mod
        end
    else
        map.rxnModAlias{react, 1} = '';
        map.rxnModID{react, 1} = '';
        map.rxnModType{react, 1} = '';
        map.rxnModColor{react, 1} = '';
        map.rxnModWidth{react, 1} = '';
    end
    map.rxnColor{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_line.Attributes.color;
    map.rxnWidth{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.annotation.celldesigner_colon_extension.celldesigner_colon_line.Attributes.width;

    % Check if the info on notes exists for the reaction
    if any(strcmp('notes', fieldnames(xmlStruct.sbml.model.listOfReactions.reaction{react})))
        map.rxnNotes{react, 1} = xmlStruct.sbml.model.listOfReactions.reaction{react}.notes.html.body.Text;
    else
        map.rxnNotes{react, 1} = '';
    end
end
clearvars react

% Check if there is any information on Compartments
if any(strcmp('celldesigner_colon_compartmentAlias', fieldnames(xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases)))
    % Test is there is only 1 compartment
    if length(xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias) == 1
        map.compartAlias{1, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias(1).Attributes.id;
        map.compartName{1, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias(1).Attributes.compartment;
    else
        % Loop over compartments to get the need information and store it in a
        % structure. Compartments refer to each individual compartment name.
        for compart = 1:length(xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias)
            map.compartAlias{compart, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias{compart}.Attributes.id;
            map.compartName{compart, 1} = xmlStruct.sbml.model.annotation.celldesigner_colon_extension.celldesigner_colon_listOfCompartmentAliases.celldesigner_colon_compartmentAlias{compart}.Attributes.compartment;
        end
        clearvars compart
    end
else
    map.compartAlias = '';
    map.compartName = '';
end
map = getMapMatrices(map);
toc

author: opencobra

repo: cobratoolbox


text:
function [outColor] = nflcolor(teamName, colorType)
%NFLCOLOR Return NFL RGB color values by team or city name
%   COLOR = NFLCOLOR(TEAMNAME,COLORTYPE) returns the 3-element RGB
%   color values of an NFL team specified using the team's city name
%   (Buffalo for example) or by the team name (Bills for example).  The
%   type of color returned is determined by COLORTYPE, COLORTYPE must
%   either be 'team' for the primary team color or 'name for the color used
%   on jersy names.  If COLORTYPE is omitted, the primary team color is
%   output.  Inputs are not case sensitive.  Partial matches are
%   supported.  In cases where more than one match is found, the first
%   match is provided and a warning is generated.
%
%   Examples:
%
%   bengalsBlack = nflcolor('bengals', 'team');
%   bengalsOrange = nflcolor('bengals','name');
%   bengalsOrange = nflcolor('cinci','name');
%   dolphinAqua = nflcolor('miami');
%
%  Color values are from the following web site...
%  http://en.wikipedia.org/wiki/Template:Infobox_NFL_player#Team_colors
%
%  Note:  Giants fans may wish to edit this function and claim New York for
%  themselves. As written, the fuction gives Jets colors when New York is
%  provided as the team name.
%

% Define nfl colors, column 3 is team color, column 4 is name color
nflData = ...
    {'buffalo',        'bills',     '0f4589', 'ffffff';
    'miami',          'dolphins',  '005e6a', 'ffffff';
    'new england',    'patriots',  '0d254c', 'd6d6d6';
    'new york',       'jets',      '313f36', 'ffffff';
    'baltimore',      'ravens',    '2a0365', 'd0b239';
    'cincinnati',     'bengals',   '000000', 'f03a16';
    'cleveland',      'browns',    'de6108', '312821';
    'pittsburgh',     'steelers',  '000000', 'f1c500';
    'houston',        'texans',    '001831', 'ffffff';
    'indianapolis',   'colts',     '003b7b', 'ffffff';
    'jacksonville',   'jaguars',   '00576e', 'd0b239';
    'tennessee',      'titans',    '648fcc', '0d254c';
    'denver',         'broncos',   '10274c', 'df6109';
    'kansas city',    'chiefs',    'b20032', 'f2c800';
    'oakland',        'raiders',   '000000', 'c4c8cb';
    'san diego',      'chargers',  '08214a', 'eec607';
    'dallas',         'cowboys',   'c5ced6', '000080';
    'new york',       'giants',    '192f6b', 'ffffff';
    'philadelphia',   'eagles',    '003b48', 'c4c8cb';
    'washington',     'redskins',  '7d0008', 'ffbe26';
    'chicago',        'bears',     '03182f', 'df6108';
    'detroit',        'lions',     '005da6', 'c4c8cb';
    'green bay',      'packers',   '213d30', 'ffcc00';
    'minnesota',      'vikings',   '3b0160', 'f0bf00';
    'atlanta',        'falcons',   'bd0d18', '000000';
    'carolina',       'panthers',  '000000', '0088d4';
    'new orleans',    'saints',    '000000', 'c9b074';
    'tampa bay',      'buccaneers','b20032', 'ffffff';
    'arizona',        'cardinals', '870619', 'ffffff';
    'st louis',       'rams',      '0d254c', 'c9b074';
    'saint louis',    'rams',      '0d254c', 'c9b074';
    'st. louis',      'rams',      '0d254c', 'c9b074';
    'san francisco',  '49ers',     '840026', 'c9b074';
    'seattle',        'seahawks',  '03182f', '4eae17'};

idx=strmatch(lower(teamName),nflData(:,1:2));

if isempty(idx)
    error('Team name not found');
end

nEntries = size(nflData,1);
if idx(1)>nEntries
    idx(1) = idx(1)-nEntries;
end

if length(idx)>1
    warning('nflcolor:multipleMatches',['Found multiple matches for %s, ' ...
        'using %s'], teamName, nflData{idx(1),2});
end

if nargin<2
    colorType = 'team';
end

switch lower(colorType)
    case 'team'
        colorHex = nflData{idx(1),3};
    case 'name'
        colorHex = nflData{idx(1),4};
    otherwise
        error('Invalid colorType specified');
end

outColor = [hex2dec(colorHex(1:2)),...
    hex2dec(colorHex(3:4)),...
    hex2dec(colorHex(5:6))]/255;


text:
classdef SMMOP2 < PROBLEM
% <multi/many> <real> <large/none> <multimodal> <sparse/none>
% Sparse multi-modal multi-objective optimization problem
% theta --- 0.1 --- Sparsity of the Pareto sets
% np    ---   4 --- Number of the Pareto sets

%------------------------------- Reference --------------------------------
% Y. Tian, R. Liu, X. Zhang, H. Ma, K. C. Tan, and Y. Jin, A
% multipopulation evolutionary algorithm for solving large-scale multimodal
% multiobjective optimization problems, IEEE Transactions on Evolutionary
% Computation, 2021, 25(3): 405-418.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties(Access = private)
        theta = 0.1;    % Sparsity of the Pareto sets
        np    = 4;    	% Number of the Pareto sets
        POS;            % Pareto optimal set for IGDX calculation
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.theta,obj.np] = obj.ParameterSet(0.1,4);
            if isempty(obj.M); obj.M = 2; end
            if isempty(obj.D); obj.D = 100; end
            obj.lower    = [zeros(1,obj.M-1)+0,zeros(1,obj.D-obj.M+1)-1];
            obj.upper    = [zeros(1,obj.M-1)+1,zeros(1,obj.D-obj.M+1)+2];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D] = size(X);
            M     = obj.M;
            S     = ceil(obj.theta*(D-M));
            g     = zeros(N,obj.np);
            for i = 1 : obj.np
            	g(:,i) =  sum(g1(X(:,M+(i-1)*S:M+i*S-1),pi/3),2)+sum(g3(X(:,[M:M+(i-1)*S-1,M+i*S:end]),0),2);
            end
            PopObj = repmat(1+min(g,[],2)/(D-M+1),1,M).*fliplr(cumprod([ones(N,1),X(:,1:M-1)],2)).*[ones(N,1),1-X(:,M-1:-1:1)];
        end
        %% Generate Pareto optimal solutions
        function R = GetOptimum(obj,N)
            % Generate points in Pareto optimal set
            A       = GetPS(obj.D-obj.M+1,obj.np,ceil(obj.theta*(obj.D-obj.M)));
            X       = UniformPoint(N/size(A,1),obj.M-1,'grid');
            obj.POS = [repmat(X,size(A,1),1),A(repmat(1:end,size(X,1),1),:)];
            % Generate points on Pareto front
            R = UniformPoint(N,obj.M);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                a = linspace(0,1,100)';
                R = [a,1-a];
            elseif obj.M == 3
                a = linspace(0,1,10)';
                R = {a*a',a*(1-a'),(1-a)*ones(size(a'))};
            else
                R = [];
            end
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case 'IGDX'
                    score = feval(metName,Population,obj.POS);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopDec = Population.decs;
            A      = GetPS(obj.D-obj.M+1,obj.np,ceil(obj.theta*(obj.D-obj.M)));
            [~,Label]  = min(pdist2(PopDec(:,obj.M:end),A),[],2);
            tempStream = RandStream('mlfg6331_64','Seed',2);
            if obj.M == 2
                for i = 1 : size(A,1)
                    color = rand(tempStream,1,3);
                    Draw(Population(Label==i).objs+(i-1)*0.05,'o','MarkerSize',6,'Marker','o','Markerfacecolor',sqrt(color),'Markeredgecolor',color,{'\it f\rm_1','\it f\rm_2',[]});
                    Draw(obj.PF+(i-1)*0.05,'-','LineWidth',1,'Color',color);
                end
            elseif obj.M == 3
                for i = 1 : size(A,1)
                    color = rand(tempStream,1,3);
                    ax = Draw(Population(Label==i).objs+(i-1)*0.05,'o','MarkerSize',8,'Marker','o','Markerfacecolor',sqrt(color),'Markeredgecolor',color,{'\it f\rm_1','\it f\rm_2','\it f\rm_3'});
                    surf(ax,obj.PF{1}+(i-1)*0.05,obj.PF{2}+(i-1)*0.05,obj.PF{3}+(i-1)*0.05,'EdgeColor',color,'FaceColor','none');
                end
            else
                for i = 1 : size(A,1)
                    Draw(Population(Label==i).objs,'-','Color',rand(tempStream,1,3),'LineWidth',2);
                end
            end
        end
    end
end

function g = g1(x,t)
    g = (x-t).^2;
end

function g = g3(x,t)
    g = 4-(x-t)-4./exp(100*(x-t).^2);
end

function PS = GetPS(D,np,S)
    PS = zeros(np,D);
    for i = 1 : np
        PS(i,(i-1)*S+1:i*S) = 1;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%
% function [Tx,fs] = synsq_cwt_squeeze(Wx, w, t, nv, opt)
%
% Calculates synchrosqueezed transform of f on a logarithmic
% scale.  Used internally by synsq_cwt_fw.
%
% Input:
%   Wx: wavelet transform of x
%   w: estimate of frequency at locations in Wx (see synsq_cwt_fw code)
%   t: time vector
%   nv: number of voices
%   opt: options struct, not currently used
%
% Output:
%   Tx: synchrosqueezed output
%   fs: associated frequencies
%
% Note the multiplicative correction term f in synsq_cwt_squeeze_mex (and in
% the matlab equivalent code commented out), required due to the fact that
% the squeezing integral of Eq. (2.7), in, [1], is taken w.r.t. dlog(a).
% This correction term needs to be included as a factor of Eq. (2.3), which
% we implement here.
%
% A more detailed explanation is available in Sec. III of [2].
% Specifically, this is an implementation of Sec. IIIC, Alg. 1.
% Note the constant multiplier log(2)/nv has been moved to the
% inverse of the normalization constant, as calculated in synsq_adm.m
%
% 1. I. Daubechies, J. Lu, H.T. Wu, "Synchrosqueezed Wavelet Transforms: a
% tool for empirical mode decomposition", 2010.
%
% 2. E. Brevdo, N.S. Fučkar, G. Thakur, and H-T. Wu, "The
% Synchrosqueezing algorithm: a robust analysis tool for signals
% with time-varying spectrum," 2011.
%
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function [Tx,fs] = synsq_cwt_squeeze(Wx, w, t, nv, opt)
    dt = t(2)-t(1);
    dT = t(end)-t(1);

    % Maximum measurable frequency of data
    %fM = 1/(4*dt); % wavelet limit - tested
    fM = 1/(2*dt); % standard
    % Minimum measurable frequency, due to wavelet
    fm = 1/dT; % really
    %fm = 1/(2*dT); % standard

    [na, N] = size(Wx);
    as = 2^(1/nv) .^ [1:1:na]';
    das = [1; diff(as)];
    lfm = log2(fm); lfM = log2(fM);
    fs = 2.^linspace(lfm, lfM, na);
    %dfs = [fs(1) diff(fs)];

    % Harmonics of diff. frequencies but same
    % magnitude have same |Tx|
    dfs = ones(size(fs));

    % The rest of the computation is performed efficiently in MEX
    if norm(Wx,'fro') < eps,
      Tx = zeros(size(Wx));
    else
      Tx = 1/nv * synsq_cwt_squeeze_mex(Wx, w, as, fs, dfs, lfm, lfM);
    end
end

author: HeLiangHIT


text:
% Chapter 10 - Limit Cycles.
% Programs_10b - Phase Portrait (Fig. 10.2).
% Copyright Birkhauser 2013. Stephen Lynch.

% Limit cycle of a van der Pol system.
% IMPORTANT - Programs_10a.m is vectorfield.m.
clear
hold on
% sys=inline('[x(2);-x(1)-5*x(2)*((x(1))^2-1)]','t', 'x');
sys = @(t,x) [x(2);-x(1)-5*x(2)*((x(1))^2-1)];
vectorfield(sys,-3:.3:3,-10:1.3:10);
[t,xs] = ode45(sys,[0 30],[2 1]);
plot(xs(:,1),xs(:,2))
hold off
axis([-3 3 -10 10])
fsize=15;
set(gca,'XTick',-3:1:3,'FontSize',fsize)
set(gca,'YTick',-10:5:10,'FontSize',fsize)
xlabel('x(t)','FontSize',fsize)
ylabel('y(t)','FontSize',fsize)
hold off

% End of Programs_10b.

repo: mcbench-benchmarks


text:
function result = nsga2(opt, varargin)
% Function: result = nsga2(opt, varargin)
% Description: The main flowchart of of NSGA-II. Note:
%   All objectives must be minimization. If a objective is maximization, the
%   objective should be multipled by -1.
%
% Syntax:
%   result = nsga2(opt): 'opt' is generated by function nsgaopt().
%   result = nsga2(opt, param): 'param' can be any data type, it will be
%       pass to the objective function objfun().
%
%   Then ,the result structure can be pass to plotnsga to display the
%   population:  plotnsga(result);
%
% Parameters:
%   opt : A structure generated by funciton nsgaopt().
%   varargin : Additional parameter will be pass to the objective functions.
%       It can be any data type. For example, if you call: nsga2(opt, param),
%       then objfun would be called as objfun(x,param), in which, x is the
%       design variables vector.
% Return:
%   result : A structure contains optimization result.
%
%         LSSSSWC, NWPU
%   Revision: 1.2  Data: 2011-07-26
%*************************************************************************

tStart = tic();
%*************************************************************************
% Verify the optimization model
%*************************************************************************
opt = verifyOpt(opt);

%*************************************************************************
% variables initialization
%*************************************************************************
nVar    = opt.numVar;
nObj    = opt.numObj;
nCons   = opt.numCons;
popsize = opt.popsize;

% pop : current population
% newpop : new population created by genetic algorithm operators
% combinepop = pop + newpop;
pop = repmat( struct(...
    'var', zeros(1,nVar), ...
    'obj', zeros(1,nObj), ...
    'cons', zeros(1,nCons),...
    'rank', 0,...
    'distance', 0,...
    'prefDistance', 0,...       % preference distance used in R-NSGA-II
    'nViol', 0,...
    'violSum', 0),...
    [1,popsize]);

% state: optimization state of one generation
state = struct(...
'currentGen', 1,...         % current generation number
'evaluateCount', 0,...      % number of objective function evaluation
'totalTime', 0,...          % total time from the beginning
'firstFrontCount', 0,...    % individual number of first front
'frontCount', 0,...         % number of front
'avgEvalTime', 0 ...        % average evaluation time of objective function (current generation)
);

result.pops     = repmat(pop, [opt.maxGen, 1]);     % each row is the population of one generation
result.states   = repmat(state, [opt.maxGen, 1]);   % each row is the optimizaiton state of one generation
result.opt      = opt;                              % use for output

% global variables
global STOP_NSGA;   %STOP_NSGA : used in GUI , if STOP_NSGA~=0, then stop the optimizaiton
STOP_NSGA = 0;

%*************************************************************************
% initialize the P0 population
%*************************************************************************
ngen = 1;
pop = opt.initfun{1}(opt, pop, opt.initfun{2:end});
[pop, state] = evaluate(opt, pop, state, varargin{:});
[opt, pop] = ndsort(opt, pop);

% state
state.currentGen = ngen;
state.totalTime = toc(tStart);
state = statpop(pop, state);

result.pops(1, :) = pop;
result.states(1)  = state;

% output
% plotnsga(result, ngen);
opt = callOutputfuns(opt, state, pop);

%*************************************************************************
% NSGA2 iteration
%*************************************************************************
while( ngen < opt.maxGen && STOP_NSGA==0)
    % 0. Display some information
	ngen = ngen+1;
    state.currentGen = ngen;
%     fprintf('\n\n************************************************************\n');
%     fprintf('*      Current generation %d / %d\n', ngen, opt.maxGen);
%     fprintf('************************************************************\n');

    % 1. Create new population
    newpop = selectOp(opt, pop);
    newpop = crossoverOp(opt, newpop, state);
    newpop = mutationOp(opt, newpop, state);
    [newpop, state] = evaluate(opt, newpop, state, varargin{:});

    % 2. Combine the new population and old population : combinepop = pop + newpop
    combinepop = [pop, newpop];

    % 3. Fast non dominated sort
    [opt, combinepop] = ndsort(opt, combinepop);

    % 4. Extract the next population
    pop = extractPop(opt, combinepop);

    % 5. Save current generation results
    state.totalTime = toc(tStart);
    state = statpop(pop, state);

    result.pops(ngen, :) = pop;
    result.states(ngen)  = state;

    % 6. plot current population and output
%     if( mod(ngen, opt.plotInterval)==0 )
%         plotnsga(result, ngen);
%     end
%     if( mod(ngen, opt.outputInterval)==0 )
%         opt = callOutputfuns(opt, state, pop);
%     end

end

% call output function for closing file
opt = callOutputfuns(opt, state, pop, -1);

% close worker processes
if( strcmpi(opt.useParallel, 'yes'))
    matlabpool close
end

toc(tStart);

author: Eric-Bradford

repo: TS-EMO


text:
function f_x = ParFor7(in1)
%PARFOR7
%    F_X = PARFOR7(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.1.
%    16-Jul-2019 11:57:27

u = in1(:,1);
uxx = in1(:,5);
f_x = (u.*9.981631376631412e-1+uxx.*9.965561702440027+u.*uxx.*1.006065455272619)./(u+9.978905414318433);

author: dynamicslab

repo: SINDy-PI


text:
function [H, H_f, H_g] = composeFrames(F,G)

% COMPOSEFRAMES  Compose two 3D frames.
%   H = COMPOSEFRAMES(F,G) composes frames F and G, where frame G is
%   specified in frame F, to get a single frame transform H. Frames are
%   structures with at least the following fields:
%       .t  translation vector
%       .q  orientation quaternion
%
%   The resulting frame H, however, contains the full frame structure.
%
%   [H, H_f, H_g] = COMPOSEFRAMES(...) returns the Jacobians of H.x wrt F.x
%   and G.x.
%
%   See also FRAME, SPLITFRAME, QUATERNION.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout == 1

    H.t = fromFrame(F,G.t);
    H.q = qProd(F.q,G.q);
    H.x = [H.t;H.q];
    H   = updateFrame(H);

else

    [H.t, T_f, T_tg]  = fromFrame(F,G.t);
    [H.q, Q_qf, Q_qg] = qProd(F.q,G.q);

    H.x = [H.t;H.q];
    H   = updateFrame(H);

    H_f(1:3,:)   = T_f;
    H_f(4:7,4:7) = Q_qf;
    H_g(1:3,1:3) = T_tg;
    H_g(4:7,4:7) = Q_qg;

end

return

%%
syms x y z a b c d X Y Z A B C D real
F.x = [x y z a b c d]';
G.x = [X Y Z A B C D]';
F = updateFrame(F);
G = updateFrame(G);
[H, H_f, H_g] = composeFrames(F,G);
simplify(H_f-jacobian(H.x,F.x))
simplify(H_g-jacobian(H.x,G.x))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% patch_mesh.m
%
% Goal: surface mesh rendering
%       need to open figure or subplot before use
%
% Li Shen
% 02/21/2006 - create

function h = patch_mesh(vertices,faces)

h = patch('faces', faces, 'vertices', vertices, ...
	'FaceColor', 'w', 'EdgeColor', 'k', 'FaceAlpha', 1); % transparancy can be adjusted by FaceAlpha
axis image; box on; view(3);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function x=ibm2num(b)
% ibm2num : convert IBM 32 bit floating point format to doubles
%    x=num2ibm(b)
% b is a matrix of uint32
% x is a corresponding matrix of doubles
%
%
% See also num2ibm

%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
%
% (C) Brian Farrelly, 22 October 2001
%  mailto:Brian.Farrelly@nho.hydro.com          Norsk Hydro Research Centre
%  phone +47 55 99 68 74                 (((                  Postboks 7190
%  fax   +47 55 99 69 70                2oooS                 N-5020 Bergen
%  home  +47 55 13 78 49                HYDRO                        Norway
%

x=repmat(NaN,size(b));

sign=bitget(b,32);                            % get sign from first bit
sign=double(sign);

% format hex
exp=bitand(b,uint32(hex2dec('7f000000')));    % get exponent from first byte, last 7 bits
exp=bitshift(exp,-24);
%format long

exp=double(exp)- 64;                          % remove bias from exponent

%format hex
frac=bitand(b,uint32(hex2dec('00ffffff')));   % get mantissa from last 3 bytes
%format long
frac=double(frac);
frac=frac/2^24;

x=(1-2*sign).*16.^exp .* frac;

err = frac==0 & (exp~=-64 | sign~=0);         % bias removal is incorrect for zero
if any(err)
   % TMH 19/06/2003
   disp(['WARNING : ',mfilename,' Invalid zero input --> Sure data are IBM FLOAT formatted ?'])
   return;
   %warning('Invalid zero input in ibm2num for the following:')
   % format hex; disp(b(err)); format
end

err = frac~=0 & (frac<1/16 | frac>=1);
if any(err)
   % TMH 19/06/2003
   disp(['WARNING : ',mfilename,' Invalid mantissa input --> Sure data are IBM FLOAT formatted ?'])
   return;
   % warning('Invalid mantissa input in ibm2num for the following:')
   % format hex; disp(b(err)); format
end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
% CORRCOEF Correlation coefficients.
%    R=CORRCOEF(X) calculates a matrix R of correlation coefficients for
%    an array X, in which each row is an observation and each column is a
%    variable.
%
%    R=CORRCOEF(X,Y), where X and Y are column vectors, is the same as
%    R=CORRCOEF([X Y]).  CORRCOEF converts X and Y to column vectors if they
%    are not, i.e., R=CORRCOEF(X,Y) is equivalent to R=CORRCOEF([X(:) Y(:)])
%    in that case.
%
%    If C is the covariance matrix, C = COV(X), then CORRCOEF(X) is
%    the matrix whose (i,j)'th element is
%
%           C(i,j)/SQRT(C(i,i)*C(j,j)).
%
%    [R,P]=CORRCOEF(...) also returns P, a matrix of p-values for testing
%    the hypothesis of no correlation.  Each p-value is the probability
%    of getting a correlation as large as the observed value by random
%    chance, when the true correlation is zero.  If P(i,j) is small, say
%    less than 0.05, then the correlation R(i,j) is significant.
%
%    [R,P,RLO,RUP]=CORRCOEF(...) also returns matrices RLO and RUP, of
%    the same size as R, containing lower and upper bounds for a 95%
%    confidence interval for each coefficient.
%
%    [...]=CORRCOEF(...,'PARAM1',VAL1,'PARAM2',VAL2,...) specifies additional
%    parameters and their values.  Valid parameters are the following:
%
%        Parameter  Value
%        'Alpha'    A number between 0 and 1 to specify a confidence
%                   level of 100*(1-ALPHA)%.  Default is 0.05 for 95%
%                   confidence intervals.
%        'Rows'     Either 'all' (default) to use all rows, 'complete' to
%                   use rows with no NaN values, or 'pairwise' to compute
%                   R(i,j) using rows with no NaN values in column i or j.
%                   The 'pairwise' option potentially uses different sets
%                   of rows to compute different elements of R, and can
%                   produce a matrix that is indefinite.
%
%    The p-value is computed by transforming the correlation to create a t
%    statistic having N-2 degrees of freedom, where N is the number of rows
%    of X.  The confidence bounds are based on an asymptotic normal
%    distribution of 0.5*log((1+R)/(1-R)), with an approximate variance equal
%    to 1/(N-3).  These bounds are accurate for large samples when X has a
%    multivariate normal distribution.
%
%    Example:  Generate random data having correlation between column 4
%              and the other columns.
%        x = randn(30,4);       % uncorrelated data
%        x(:,4) = sum(x,2);     % introduce correlation
%        [r,p] = corrcoef(x)    % compute sample correlation and p-values
%        [i,j] = find(p<0.05);  % find significant correlations
%        [i,j]                  % display their (row,col) indices
%
%    Class support for inputs X,Y:
%       float: double, single
%
%    See also COV, VAR, STD, PLOTMATRIX.
%
%    Reference page in Doc Center
%       doc corrcoef
%
%    Other functions named corrcoef
%
%       codistributed/corrcoef    gpuArray/corrcoef    ts/corrcoef
%       fints/corrcoef
%

author: jmaih