text:
function h = showsolutionpoly(node,elem,u,varargin)
%% SHOWSOLUTIONPOLY plots the solution u on a polygonal mesh in 2-D.
%
%    showmeshpoly(node,elem, u) displays u on a polygonal mesh. elem{j} is a cell
%    array with j-th cell representing the j-th element.
%
%    showmeshpoly(node,elem,u,viewangle) changes the display angle. The
%    deault view angle on planar meshes is view(2) and view(3) for surface
%    meshes.
%
%   Example:
%
%     load polygonunstruct
%     u = @(p) sin(2*pi*p(:,1)).*cos(2*pi*p(:,2));
%     uI = u(node);
%     showsolutionpoly(node,elem,uI);
%
%   See also showmesh, showboundary3, showmeshpoly.
%
% Added by Shuhao Cao.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

V = [node, u];
elemVertexNumber = cellfun('length',elem); % the number of vertices per element
for Nv = min(elemVertexNumber):max(elemVertexNumber)
    idx = (elemVertexNumber == Nv); % index of elements with Nv vertices
    elemNv = cell2mat(elem(idx));
    if size(elemNv,2) == 1
        elemNv = reshape(elemNv,[Nv,nnz(idx)])';
    end
    h = patch('faces',elemNv,'vertices',V);
    hold on;
    set(h,    'FaceVertexCData',u,...
              'Facecolor','interp',...
              'Edgecolor','k',...
              'Edgealpha',1/log(size(node,1)),...
              'FaceAlpha',0.65);
   %%
   if (nargin>3) && ~isempty(varargin) % set display property
       if isnumeric(varargin{1})
           view(varargin{1});
           if nargin>4
               set(h,varargin{2:end});
           end
       else
           set(h,varargin{1:end});
       end
   end
end
view(3); axis tight; grid on;
camproj('perspective');
hold off
%%

repo: ifem


text:
function r4_uniform_01_test ( )

%*****************************************************************************80
%
%% R4_UNIFORM_01_TEST tests R4_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  R4_UNIFORM_01 produces a sequence of random values.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  SEED = %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ r, seed ] = r4_uniform_01 ( seed );
    fprintf ( 1, '  %2d  %14g\n', i, r );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef svd

methods(Static)

    function [rank, max_gap] = mahdi_rank(singular_values)
        % rank accordingly to mahdi 2013 heuristic
        [ min_val , ind_min ] = min( diff( singular_values(1:end-1) ) ) ;
        rank = ind_min;
        max_gap = -min_val;
    end

    function rank = vidal_rank(singular_values,kappa)
        if nargin < 2
            kappa = .1;
        end
        % Rank used in Vidal papers
        n = length(singular_values);
        % we will try rank values between 1 to n-1.
        % an array to store criterion value
        criterion = zeros(1, n-1);
        for(r=1:n-1)
            num = singular_values(r+1)^2;
            den = sum(singular_values(1:r).^2);
            criterion(r)=num/den + kappa*r;
        end
        %criterion
        % find the rank with minimum value of the criterion
        [min_value, index]=min(criterion);
        rank = index;
    end

    function result = low_rank_approx(X, r)
        % return low rank approximation of X
        [U S V] = svd(X);
        % keep the first r left singular vectors
        U = U(:, 1:r);
        % keep the first r singular values
        S = S(1:r, 1:r);
        % keep the first r right singular vectors
        V = V(:, 1:r);
        % Return the approximation
        result = U * S * V';
    end

    function [result, basis] = low_rank_projection(X, r)
        % Projects X to a low dimensional space
        % X is NxS
        % result is RxS
        % basis is [NxR] orthonormal basis

        % Compute the SVD
        [U, ~, ~] = svd(X, 0);
        % Choose the low rank basis
        basis = U(:, 1:r);
        % Compute coefficients in this basis
        result = basis' * X;
    end

    function result = low_rank_basis(X, r)
        % Returns the ON basis for low rank approximation
        [U S V] = svd(X, 'econ');
        result = U(:, 1:r);
    end

    function result = low_rank_bases(X, counts, r)
        % low rank bases for individual subspaces
        K = length(counts);
        % bases cell array
        result = cell(1, K);
        [start_indices, end_indices] = spx.cluster.start_end_indices(counts);
        for k=1:K
            ss = start_indices(k);
            ee = end_indices(k);
            XX = X(:, ss:ee);
            basis = spx.la.svd.low_rank_basis(XX, r);
            result{k} = basis;
        end
    end

    function [result, r] = mahdi_rank_basis(X)
        [U S V] = svd(X, 'econ');
        sv = diag(S);
        r = spx.la.svd.mahdi_rank(sv);
        % r = r + 1;
        result = U(:, 1:r);
    end

    function [result, ranks] = mahdi_rank_bases(X, counts)
        % low rank bases for individual subspaces
        K = length(counts);
        % bases cell array
        result = cell(1, K);
        [start_indices, end_indices] = spx.cluster.start_end_indices(counts);
        ranks = zeros(1, K);
        for k=1:K
            ss = start_indices(k);
            ee = end_indices(k);
            XX = X(:, ss:ee);
            [basis,r] = spx.la.svd.mahdi_rank_basis(XX);
            result{k} = basis;
            ranks(k) = r;
        end
    end

    function [result, r] = vidal_rank_basis(X, kappa)
        if nargin < 3
            kappa = 0.1;
        end
        [U S V] = svd(X, 'econ');
        sv = diag(S)';
        r = spx.la.svd.vidal_rank(sv, kappa);
        % r = r + 1;
        result = U(:, 1:r);
    end

    function [result, ranks] = vidal_rank_bases(X, counts, kappa)
        if nargin < 3
            kappa = 0.1;
        end
        % low rank bases for individual subspaces
        K = length(counts);
        % bases cell array
        result = cell(1, K);
        [start_indices, end_indices] = spx.cluster.start_end_indices(counts);
        ranks = zeros(1, K);
        for k=1:K
            ss = start_indices(k);
            ee = end_indices(k);
            XX = X(:, ss:ee);
            [basis,r] = spx.la.svd.vidal_rank_basis(XX, kappa);
            result{k} = basis;
            ranks(k) = r;
        end
    end

end

end

author: indigits

repo: sparse-plex


text:
function [x1, x2, fval1, fval2, beta1, beta2] = computeLaunchWindows(bodyInfo, launchLong, launchLat, targetInc, targetRAAN, windowSearchStartUT)
%computeLaunchWindows Summary of this function goes here
%   Detailed explanation goes here
    %d = "dummy"

    dSMA = bodyInfo.radius + 1;
    dEcc = 0.0;
    dInc = targetInc;
    dRAAN = targetRAAN;
    dArg = 0.0;
    dTru = 0.0;
    gmu = bodyInfo.gm;

    [dRVect,dVVect]=getStatefromKepler(dSMA, dEcc, dInc, dRAAN, dArg, dTru, gmu);
    dVVectMag = norm(dVVect);
    dHVect = cross(dRVect,dVVect);
    dPoint = [0;0;0];

    distFunc = @(ut) getDistBetweenPlaneAndLaunchSite(ut, launchLat, launchLong, bodyInfo.radius, bodyInfo, dHVect, dPoint);
    rotPeriod = bodyInfo.rotperiod;

    options = optimset('TolX',1E-8);

    lb1 = windowSearchStartUT;
    ub1 = windowSearchStartUT + rotPeriod;
    [x1,fval1] = fminbnd(distFunc,lb1,ub1, options);

    lb2 = lb1;
    ub2 = x1-1;
    [x21,fval21] = fminbnd(distFunc,lb2,ub2, options);

    lb3 = x1+1;
    ub3 = ub1;
    if(ub3<lb3)
        ub3 = abs(ub3-lb3)+lb3;
    end
    [x22,fval22] = fminbnd(distFunc,lb3,ub3, options);

    if(fval21 < fval22)
        x2 = x21;
        fval2 = fval21;
    else
        x2 = x22;
        fval2 = fval22;
    end

    sinBeta = cos(targetInc)/cos(launchLat);
    beta1 = AngleZero2Pi(getLaunchAz(x1, sinBeta, dHVect, distFunc, bodyInfo, dVVectMag));
    beta2 = AngleZero2Pi(getLaunchAz(x2, sinBeta, dHVect, distFunc, bodyInfo, dVVectMag));
end

function [d, rVectECI] = getDistBetweenPlaneAndLaunchSite(ut, lat, long, alt, bodyInfo, planeNormalVect, planePoint)
    rVectECI = getInertialVectFromLatLongAlt(ut, lat, long, alt, bodyInfo, [NaN;NaN;NaN]);

    v = rVectECI - planePoint;
    n = normVector(planeNormalVect); %unit vector
    d = abs(dot(v,n)); %distance between plane and point
end

function beta = getLaunchAz(x, sinBeta, hVect, distFunc, bodyInfo, dVVectMag)
    [~,rVect] = distFunc(x);
    vVect = dVVectMag*normVector(cross(normVector(hVect),normVector(rVect)));
    beta = asin(sinBeta);
    if(vVect(3) < 0)
        beta = pi - beta;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function test_suite = test_survival_coxph

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_SURVIVAL_COXPH
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testPredictionsCoxph(testCase)
  verifyVarsEqual(testCase, getName(), {'Ef','Varf'}, ...
    'AbsTol', 0.1)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function y = truerand(varargin)
% TRUERAND Returns truly random integers using random.org's Random Integer
%   Generator. According to random.org, the numbers are generated based on
%   atmospheric noise and skew-corrected to generate uniform numbers. The
%   generated numbers have been shown to pass the NIST tests for RNGs.
%
%   The range of numbers is -1e9 to 1e9 and the maximum number of values
%   that can be generated is 10,000
%
% USAGE:
%  truerand(rows,cols,min,max) returns a matrix of size rows-by-cols with
%  random integers between min and max.
%  truerand(n,min,max) returns an n by 1 vector
%  truerand(n, m) uses the default values min = 1, max = 100
%  truerand(n) and truerand
%
% EXAMPLES:
%  y = truerand
%  y = truerand(9)
%  y = truerand(6,6)
%  y = truerand(5,1,20)
%  y = truerand(3, 4, 15, 30)
%
% For more information visit random.org

m = 1;
n = 1;
min = 1;
max = 100;
seed = 'new'; % New implies new random numbers - they will not be the same with each run

if nargin > 0
    m = varargin{1};
end
if nargin == 2 || nargin == 4
    n = varargin{2};
end

if m*n>=10000
    error('The maximum number of values requested must be no greater than 10000')
end

if min < -1e9 || max > 1e9
    error('The range of values must lie in [-1e9, 1e9]')
end

if nargin == 3 || nargin == 4
    [min,max] = varargin{end-1:end};
end

url = sprintf('http://www.random.org/integers/?num=%d&min=%d&max=%d&col=1&base=10&format=plain&rnd=%s',...
              m*n, min, max, seed);
data = urlread(url);
y = reshape(str2num(data),m,n); %#ok<ST2NM>

author: Sable

repo: mcbench-benchmarks


text:
function this = load(this, fileName)
% Loads Geometry info from affine image header (.nii/.hdr/.img) or Philips
% (par/rec) or recon6-ImageData (.mat, Geometry-object)
%
% NOTE: .mat-header files (for 4D niftis) are ignored, since the same voxel
%       position is assumed in each volume for MrImage
%
%   geom = MrAffineTransformation()
%   geom = geom.load(fileName)
%
% This is a method of class MrAffineTransformation.
%
% IN
%
%
% EXAMPLE
%   geom = MrAffineTransformation()
%   geom.load('test.nii')
%   geom.load('test.hdr/img')
%   geom.load('test.par/rec')
%
%   See also MrAffineTransformation MrAffineTransformation.update

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-15
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if exist(fileName, 'file')
    [~, ~, ext] = fileparts(fileName);

    switch ext
        case {'.hdr', '.nii', '.img'}
            V = spm_vol(fileName);
            affineMatrix = V.mat;

            if ~isempty(affineMatrix)
                this.update_from_affine_matrix(affineMatrix);
            end

        case {'.par', '.rec'}
            this = this.load_par(fileName);
        case {'.mat'} % recon 6
            this = this.load_recon6_mat(fileName);
        otherwise
            warning('Only Philips (.par/.rec), nifti (.nii) and analyze (.hdr/.img) files are supported');
    end
else
    fprintf('Geometry data could not be loaded: file %s not found.\n', ...
        fileName);
end

end

author: translationalneuromodeling

repo: tapas


text:
    %This example generates NE elements for a rectangular structure of
    %length = Ly units and width = Lx units with Nx divisions on the x
    %axis using femTriangularMeshGenerator function

    %   Kehinde Orolu
    %   Systems Engineering
    %   University of Lagos, Nigeria
    %   olukeh@yahoo.com

    cla
    Lx=10;
    Ly=10;
    Nx=8;
    NE=144;

    [coords cT nNodes ]=femTriangularMeshGenerator(Lx,Ly,Nx,NE);

    disp(['Number of nodes =  ',num2str(nNodes)])
    disp('Connectivity Table')
    disp(cT)

    z=1;
    for i=1:NE
        figure(1),patch('Vertices',coords(z:z+2,:),'Faces',[1,2,3],'FaceColor','none','EdgeColor','g')
        hold on
        z=z+3;
    end

    figure(1),scatter(coords(:,1),coords(:,2),'MarkerFaceColor','r')

    hold off


text:
function linplus_test525 ( )

%*****************************************************************************80
%
%% TEST525 tests R8PP_FA, R8PP_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST525\n' );
  fprintf ( 1, '  R8PP_FA factors a R8PP system,\n' );
  fprintf ( 1, '  R8PP_SL solves a R8PP system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8pp_random ( n, seed );

  r8pp_print ( n, a, '  The R8PP matrix:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );

  r8vec_print ( n, x, '  The desired solution:' );
%
%  Compute the corresponding right hand side.
%
  b = r8pp_mxv ( n, a, x );

  r8vec_print ( n, b, '  The right hand side:' );
%
%  Factor the matrix.
%
  [ a_lu, info ] = r8pp_fa ( n, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Fatal error!\n' );
    fprintf ( 1, '  R8PP_FA declares the matrix is singular!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The R8PP matrix has been factored.\n' );
%
%  Solve the linear system.
%
  x = r8pp_sl ( n, a_lu, b );

  r8vec_print ( n, x, '  Solution:' );

  return
end


text:
% Change to your downloaded location
clear
addpath('C:\liblinear\matlab')
addpath('../training_code')
addpath('../utilities')

%% load shared definitions and AU data
shared_defs;

% Set up the hyperparameters to be validated
hyperparams.c = 10.^(-7:1:1);
hyperparams.e = 10.^(-3);

hyperparams.validate_params = {'c', 'e'};

% Set the training function
svr_train = @svm_train_linear;

% Set the test function (the first output will be used for validation)
svr_test = @svm_test_linear;

%%
for a=1:numel(aus)

    au = aus(a);

    rest_aus = setdiff(all_aus, au);

    % make sure validation data's labels are balanced
    [users_train, users_valid] = get_balanced_fold(DISFA_dir, users, au, 1/3, 1);

    % need to split the rest
    [train_samples, train_labels, valid_samples, valid_labels, ~, PC, means, scaling, valid_ids, valid_success] = Prepare_HOG_AU_data_generic_dynamic(users_train, users_valid, au, rest_aus, DISFA_dir, hog_data_dir);

    train_labels(train_labels > 1) = 1;
    valid_labels(valid_labels > 1) = 1;

    train_samples = sparse(train_samples);
    valid_samples = sparse(valid_samples);

    %% Validate here
    hyperparams.success = valid_success;

    [ best_params, ~ ] = validate_grid_search_no_par(svr_train, svr_test, false, train_samples, train_labels, valid_samples, valid_labels, hyperparams);

    model = svr_train(train_labels, train_samples, best_params);

    [~, prediction] = svr_test(valid_labels, valid_samples, model);

    name = sprintf('classifiers/AU_%d_dyn.mat', au);

    [ accuracies, F1s, corrs, ccc, rms, classes ] = evaluate_regression_results( prediction, valid_labels );

    save(name, 'model', 'accuracies', 'F1s', 'corrs', 'rms', 'ccc', 'prediction', 'valid_labels');

    name = sprintf('classifiers/AU_%d_dyn.dat', au);

    pos_lbl = model.Label(1);
    neg_lbl = model.Label(2);

    w = model.w(1:end-1)';
    b = model.w(end);

    svs = bsxfun(@times, PC, 1./scaling') * w;
    write_lin_svm(name, means, svs, b, pos_lbl, neg_lbl);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function M = multinomialfactory(n, m)
% Manifold of n-by-m column-stochastic matrices with positive entries.
%
% function M = multinomialfactory(n, m)
%
% The returned structure M is a Manopt manifold structure to optimize over
% the set of n-by-m matrices with (strictly) positive entries and such that
% the entries of each column sum to one.
%
% The metric imposed on the manifold is the Fisher metric such that
% the set of n-by-m column-stochastic matrices (aka the multinomial manifold)
% is a Riemannian submanifold of the space of n-by-m matrices. Also it
% should be noted that the retraction operation that we define
% is first order and as such the checkhessian tool cannot verify
% the slope correctly.
%
% The file is based on developments in the research paper
% Y. Sun, J. Gao, X. Hong, B. Mishra, and B. Yin,
% "Heterogeneous tensor decomposition for clustering via manifold
% optimization", arXiv:1504.01777, 2015.
%
% Link to the paper: http://arxiv.org/abs/1504.01777.
%
% Please cite the Manopt paper as well as the research paper:
%     @Techreport{sun2014multinomial,
%       Title   = {Heterogeneous tensor decomposition for clustering via manifold optimization},
%       Author  = {Sun, Y. and Gao, J. and Hong, X. and Mishra, B. and Yin, B.},
%       Journal = {Arxiv preprint arXiv:1504.01777},
%       Year    = {2014}
%     }

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, April 06, 2015.
% Contributors:
% Change log:

    M.name = @() sprintf('%dx%d column-stochastic matrices with positive entries', n, m);

    M.dim = @() (n-1)*m;

    % We impose the Fisher metric.
    M.inner = @iproduct;
    function ip = iproduct(X, eta, zeta)
        ip = sum((eta(:).*zeta(:))./X(:));
    end

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(X, Y) error('multinomialfactory.dist not implemented yet.');

    M.typicaldist = @() m*pi/2; % This is an approximation.

    % Column vector of ones of length n.
    e = ones(n, 1);

    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)
        lambda = -sum(X.*egrad, 1); % Row vector of length m.
        rgrad = X.*egrad + (e*lambda).*X; % This is in the tangent space.
    end

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, eta)

        % Riemannian gradient computation.
        % lambda is a row vector of length m.
        lambda = - sum(X.*egrad, 1);
        rgrad =  X.*egrad + (e*lambda).*X;

        % Directional derivative of the Riemannian gradient.
        % lambdadot is a row vector of length m.
        lambdadot = -sum(eta.*egrad, 1) - sum(X.*ehess, 1);
        rgraddot = eta.*egrad + X.*ehess + (e*lambdadot).*X + (e*lambda).*eta;

        % Correction term because of the non-constant metric that we
        % impose. The computation of the correction term follows the use of
        % Koszul formula.
        correction_term = - 0.5*(eta.*rgrad)./X;
        rhess = rgraddot + correction_term;

        % Finally, projection onto the tangent space.
        rhess = M.proj(X, rhess);
    end

    % Projection of the vector eta in the ambeint space onto the tangent
    % space.
    M.proj = @projection;
    function etaproj = projection(X, eta)
        alpha = sum(eta, 1); % Row vector of length m.
        etaproj = eta - (e*alpha).*X;
    end

    M.tangent = M.proj;
    M.tangent2ambient = @(X, eta) eta;

    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        % A first-order retraction.
        Y = X.*exp(t*(eta./X)); % Based on mapping for positive scalars.
        Y = Y./(e*(sum(Y, 1))); % Projection onto the constraint set.
        % For numerical reasons, so that we avoid entries going to zero:
        Y = max(Y, eps);
    end

    M.exp = @exponential;
    function Y = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Y = retraction(X, eta, t);
        warning('manopt:multinomialfactory:exp', ...
            ['Exponential for the Multinomial manifold' ...
            'manifold not implemented yet. Used retraction instead.']);
    end

    M.hash = @(X) ['z' hashmd5(X(:))];

    M.rand = @random;
    function X = random()
        % A random point in the ambient space.
        X = rand(n, m); %
        X = X./(e*(sum(X, 1)));
    end

    M.randvec = @randomvec;
    function eta = randomvec(X)
        % A random vector in the tangent space
        eta = randn(n, m);
        eta = M.proj(X, eta); % Projection onto the tangent space.
        nrm = M.norm(X, eta);
        eta = eta / nrm;
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(X) zeros(n, m);

    M.transp = @(X1, X2, d) projection(X2, d);

    % vec and mat are not isometries, because of the scaled metric.
    M.vec = @(X, U) U(:);
    M.mat = @(X, u) reshape(u, n, m);
    M.vecmatareisometries = @() false;
end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
classdef PostProcTestsSuite < handle

    methods

        function obj = PostProcTestsSuite()
            warning('off', 'MATLAB:structOnObject')
            results = runtests("PostProcTests","ProcedureName","testFastDisplacement", 'Debug', true);
            table(results)
            warning('on', 'MATLAB:structOnObject')
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
close all; clear all; clc;
rng('default');
resize_images = true;

target.psnr = 35;
% target mean square error
target.mse = 255^2 * 10^(-target.psnr/10);
image_width = 42;
image_height = 48;
num_pixels = image_width * image_height;
% target sum of squared errors
target.sse = num_pixels * target.mse;
% target residual norm
target.rnorm = sqrt(target.sse);
fprintf('Target PSNR: %.2f dB\n', target.psnr);
fprintf('Target Mean Square Error: %.2f\n', target.mse);
fprintf('Target Root Mean Square Error: %.2f\n', sqrt(target.mse));
fprintf('Target Residual Norm: %.2f\n', target.rnorm);

N =  num_pixels;
redundancy_factor = 8;
D = redundancy_factor * N;

A = spx.dict.simple.gaussian_mtx(N, D);
dc_atom = (1/sqrt(N)) * ones(N, 1);
% throw away one original atom and replace with DC.
A = [dc_atom A(:, 1:D-1)];
dictionary = spx.dict.MatrixOperator(A);

yf = spx.data.image.YaleFaces();
yf.load();
fprintf('\n Resizing images:\n');
tstart = tic;
yf.resize_all(42, 48);
images = yf.ImageData;
[~, S] = size(images);
%images = yf.get_subject_images_resized(1, 42, 48);
elapsed = toc(tstart);
fprintf('Time taken: %.2f seconds \n', elapsed);
representations = zeros(D, S);
elapsed_times = zeros(1, S);
support_sizes = zeros(1, S);
res_norms = zeros(1, S);
psnrs = zeros(1, S);
fprintf('Total images: %d\n', S);
for s=1:S
    Y = images(:, s);
    solver = spx.pursuit.single.OrthogonalMatchingPursuit(dictionary);
    solver.MaxResNorm = target.rnorm;
    solver.Verbose = true;
    fprintf('sparsifying image: %d\n', s);
    tstart = tic;
    result = solver.solve_qr(Y);
    elapsed = toc(tstart);
    fprintf('Time taken: %.2f seconds \n', elapsed);
    fprintf('Support size: %d\n', result.iterations);
    support_sizes(s) = result.iterations;
    elapsed_times(s) = elapsed;
    representations(:, s) = result.z;
    res_norms(s) = result.rnorm;
    sse = result.rnorm^2;
    mse =  sse / num_pixels;
    psnr = 10 * log10(255^2 / mse);
    psnrs(s) = psnr;
    fprintf('PSNR : %.3f dB\n', psnr);
    if mod(s, 100) == 0
        % We wish to save intermediate data too
        save('omp_representations', 'representations');
        save('omp_stats', 'elapsed_times', 'support_sizes', 'res_norms', 'psnrs', 's');
    end
end
% final saving of all data.
save('bin/omp_dictionary', 'dictionary', 'N', 'D');
save('bin/omp_images', 'images');
save('bin/omp_representations', 'representations');
save('bin/omp_stats', 'elapsed_times', 'support_sizes', 'res_norms', 'psnrs', 'S');


text:
function [metabolite_structure,IDsAdded,InchiKeyList,InchiStringList ] = parseChebiIdWebPage(metabolite_structure,startSearch,endSearch,printInchis)

if ~exist('printInchis','var')
    printInchis = 0; % default setting is not to print the InchiKeyList,InchiStringList
end

Mets = fieldnames(metabolite_structure);

if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(Mets);
end

% get inchiKey and smiles from chebiid

warning off;
annotationSource = 'ChEBI website';
annotationType = 'automatic';

IDsAdded = ''; a=1;
InchiStringList = ''; c=1;
InchiKeyList = '';d=1;

fields = fieldnames(metabolite_structure.(Mets{1}));
for i = startSearch : endSearch
    % check that chebiid exists
    if ~isempty(metabolite_structure.(Mets{i}).cheBIId) && isempty(find(isnan(metabolite_structure.(Mets{i}).cheBIId),1))
        % check that smile or inchiKey does not exist
        % go to chebi and parse website for smile
        if isnumeric(metabolite_structure.(Mets{i}).cheBIId)
metabolite_structure.(Mets{i}).cheBIId = num2str(isnumeric(metabolite_structure.(Mets{i}).cheBIId));
        end

        metabolite_structure.(Mets{i}).cheBIId = regexprep(metabolite_structure.(Mets{i}).cheBIId,',',';');
        if contains(metabolite_structure.(Mets{i}).cheBIId,';')
            list = split(metabolite_structure.(Mets{i}).cheBIId,';');
        else
            list{1} = metabolite_structure.(Mets{i}).cheBIId;
        end
        list = regexprep(list,' ','');
        for t = 1 : length(list)
            url=strcat('https://www.ebi.ac.uk/chebi/searchId.do?chebiId=CHEBI:',num2str(list{t}));
            syst = urlread(url);
            if isempty(metabolite_structure.(Mets{i}).smile) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).smile),1))
                startvalue=strfind(syst,'>SMILES	');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+100);
                    [tok,rem] = strtok(string,'>');
                    rem = regexprep(rem,'','');
                    [tok2,rem2] = strtok(rem,'<');
                    tok2 = regexprep(tok2,'>','');
                    idNew = regexprep(tok2,'\s','');
                    idNew = cellstr(idNew);
                    metabolite_structure.(Mets{i}).smile = char(idNew);
                    metabolite_structure.(Mets{i}).smile_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'smile';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).inchiKey) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).inchiKey),1))
                startvalue=strfind(syst,'>InChIKey	');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+100);
                    [tok,rem] = strtok(string,'>');

                    rem = regexprep(rem,'','');
                    [tok2,rem2] = strtok(rem,'<');
                    tok2 = regexprep(tok2,'>','');
                    idNew = regexprep(tok2,'\s','');
                    metabolite_structure.(Mets{i}).inchiKey = idNew;
                    metabolite_structure.(Mets{i}).inchiKey_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'inchiKey';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).inchiString) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).inchiString),1))
                startvalue=strfind(syst,'td>InChI=');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+1000);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).inchiString = idNew;
                    metabolite_structure.(Mets{i}).inchiString_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'inchiString';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end

            if isempty(metabolite_structure.(Mets{i}).wikipedia) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).wikipedia),1))
                startvalue=strfind(syst,'Read full article at Wikipedia');
                % smile is in the next line
                try
                    string = syst(startvalue(1)-100:startvalue(1));
                    [tok,rem] = strtok(string,'wiki');

                    rem = regexprep(rem,'wikipedia.org/wiki/','');
                    rem = regexprep(rem,'target="_blank">R','');
                    idNew = regexprep(rem,'" ','');
                    metabolite_structure.(Mets{i}).wikipedia = idNew;
                    metabolite_structure.(Mets{i}).wikipedia_source =[annotationSource,':',annotationType,':',datestr(now)];

                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'wikipedia';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).knapsack) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).knapsack),1))
                startvalue=strfind(syst,'http://www.knapsackfamily.com/knapsack');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+200);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).knapsack = idNew;
                    metabolite_structure.(Mets{i}).knapsack_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'knapsack';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).keggId) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).keggId),1))
                startvalue=strfind(syst,'http://www.genome.ad.jp');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+200);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).keggId = idNew;
                    metabolite_structure.(Mets{i}).keggId_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'keggId';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).biocyc) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).biocyc),1))
                startvalue=strfind(syst,'http://biocyc.org/');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+200);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).biocyc = idNew;
                    metabolite_structure.(Mets{i}).biocyc_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'biocyc';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).hmdb) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).hmdb),1))
                startvalue=strfind(syst,'http://www.hmdb.ca');
                try
                    % smile is in the next line
                    string = syst(startvalue(1):startvalue(1)+200);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).hmdb = idNew;
                    metabolite_structure.(Mets{i}).hmdb_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'hmdb';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).casRegistry) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).casRegistry),1))
                startvalue=strfind(syst,'https://chem.nlm.nih.gov/chemidplus/');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+200);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).casRegistry = idNew;
                    metabolite_structure.(Mets{i}).casRegistry_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'casRegistry';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if isempty(metabolite_structure.(Mets{i}).lipidmaps) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).lipidmaps),1))
                startvalue=strfind(syst,'http://www.lipidmaps.org');
                % smile is in the next line
                try
                    string = syst(startvalue(1):startvalue(1)+200);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).lipidmaps = idNew;
                    metabolite_structure.(Mets{i}).lipidmaps_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    IDsAdded{a,1} = Mets{i};
                    IDsAdded{a,2} = 'lipidmaps';
                    IDsAdded{a,3} = char(idNew);
                    a=a+1;
                end
            end
            if printInchis
                % inchiKeys and inchiStrings separately printed out
                startvalue=strfind(syst,'>InChIKey	');
                try
                    string = syst(startvalue(1):startvalue(1)+100);
                    [tok,rem] = strtok(string,'>');

                    rem = regexprep(rem,'','');
                    [tok2,rem2] = strtok(rem,'<');
                    tok2 = regexprep(tok2,'>','');
                    idNew = regexprep(tok2,'\s','');
                    metabolite_structure.(Mets{i}).inchiKey = idNew;
                    metabolite_structure.(Mets{i}).inchiKey_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    InchiKeyList{c,1} = Mets{i};
                    InchiKeyList{c,2} = list{t};
                    InchiKeyList{c,3} = idNew;
                    InchiKeyList{c,4} = [annotationSource,':',annotationType,':',datestr(now)];
                    c = c + 1;
                end
                % get Inchistring (again)
                startvalue=strfind(syst,'td>InChI=');

                try
                    string = syst(startvalue(1):startvalue(1)+1000);
                    [tok,rem] = strtok(string,'>');
                    [tok2,rem2] = strtok(rem,'<');
                    idNew = regexprep(tok2,'>','');
                    metabolite_structure.(Mets{i}).inchiString = idNew;
                    metabolite_structure.(Mets{i}).inchiString_source =[annotationSource,':',annotationType,':',datestr(now)];
                    clear startvalue;
                    InchiStringList{d,1} = Mets{i};
                    InchiStringList{d,2} = list{t};
                    InchiStringList{d,3} = idNew;
                    InchiStringList{d,4} = [annotationSource,':',annotationType,':',datestr(now)];
                    d=d+1;
                end

            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function varargout = mrFlatMesh(varargin)
% MRFLATMESH Application M-file for mrFlatMesh.fig (cortical flattening)
%
%    FIG = MRFLATMESH launch mrFlatMesh GUI.
%
%    MRFLATMESH('callback_name', ...) invokes the named callback.
%
% Last Modified by GUIDE v2.5 17-Sep-2009 15:29:36

if (nargin == 0 || iscell(varargin{1}))  % LAUNCH GUI

    fig = openfig(mfilename,'reuse');

    % Generate a structure of handles to pass to callbacks, and store it.
    handles = guihandles(fig);
    if(nargin>0 && iscell(varargin{1}))
        % Set the specified defaults
        param = varargin{1};
        for(ii=1:2:length(param))
            if(strcmp(param{ii},'grayPath'))
                set(handles.inputGrayPath, 'String', param{ii+1});
            end
            if(strcmp(param{ii},'meshPath'))
                set(handles.inputMeshPath, 'String', param{ii+1});
            end
            if(strcmp(param{ii},'savePath'))
                set(handles.inputSavePath, 'String',  param{ii+1});
            end
            if(strcmp(param{ii},'startXYZ'))
                set(handles.editStartX, 'String', num2str(param{ii+1}(1)));
                set(handles.editStartY, 'String', num2str(param{ii+1}(2)));
                set(handles.editStartZ, 'String', num2str(param{ii+1}(3)));
            end
            if(strcmp(param{ii},'unfoldRadiusMM'))
                set(handles.editUnfoldSize, 'String', num2str(param{ii+1}));
            end
            if(strcmp(param{ii},'hemi'))
                hemi = getPopupValue(handles.popHemi,param{ii+1});
                set(handles.popHemi, 'Value',hemi);
            end
        end
    end

    guidata(fig, handles);

    if (nargout > 0)
        varargout{1} = fig;
    end % endif nargout

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

%     try
         [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
%     catch
%         disp(lasterr);
%     end % end try/catch

end % end function

%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and
%| sets objects' callback properties to call them through the FEVAL
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

%-------------------------
function retVar=browseGrayButton_Callback(h,eventData,handles)
%thisFig=gcf;
set (handles.status,'String','Waiting for gray matter segmentation file.');

grayFile = get(handles.inputGrayPath,'String');
if isempty(grayFile), grayFile = meshFindPath(handles); end

[fileName,inputGrayPath]=myUiGetFile(grayFile,{'*.?ray'; '*.mat'},'Get .gray file or coords.mat file');
if fileName == 0, disp('User canceled.'); return; end

set(handles.inputGrayPath,'String',fullfile(inputGrayPath,fileName));
set (handles.status,'String','Gray file set');
return;

%-------------------------
function retVar=browseMeshButton_Callback(h,eventData,handles)
%thisFig=gcf;
set (handles.status,'String','Waiting for mrGray mesh file *.MrM or *.mrm.');

inputMeshPath = get(handles.inputMeshPath,'String');
if isempty(inputMeshPath), inputMeshPath = meshFindPath(handles); end

[fileName,inputMeshPath]=myUiGetFile(inputMeshPath,{'*.mrm'; '*nii.gz'},'Get mesh file *.MrM or *.mrm or *.nii.gz');
if  fileName == 0, disp('User canceled.'); return; end

set(handles.inputMeshPath,'String',fullfile(inputMeshPath,fileName));
set (handles.status,'String','Mesh file set');
return;

%-------------------------
function retVar=saveMeshFile_Callback(h,eventData,handles)
%thisFig=gcf;
set (handles.status,'String','Setting output path.');

savePath = get(handles.inputSavePath,'String');
if isempty(savePath), savePath = meshFindPath(handles); end

[fileName,outputMeshPath]=myUiPutFile(savePath,'*.mat','Select an output file');
if  fileName == 0, disp('User canceled.'); return; end

set(handles.inputSavePath,'String',fullfile(outputMeshPath,fileName));
set (handles.status,'String','Output path set');
return;

%-------------------------
function retVar=go_Callback(h,eventData,handles)
% Press Go (green) button (Go!)

set (handles.status,'String','Running unfold.');
% Get filenames, coords and unfold size
meshFileName=get(handles.inputMeshPath,'String');
grayFileName=get(handles.inputGrayPath,'String');
flatFileName=get(handles.inputSavePath,'String');

xPos=str2num(get(handles.editStartX,'String'));
yPos=str2num(get(handles.editStartY,'String'));
zPos=str2num(get(handles.editStartZ,'String'));
if isempty(xPos) | isempty(yPos) | isempty(zPos)
    errordlg('You must fill in the starting positions.');
    return;
end
startPos=[xPos,yPos,zPos];

classes = get(handles.popHemi,'String');
hemi = classes{get(handles.popHemi,'value')};

% sagMM=get(handles.editScaleSag,'String');
% axMM=get(handles.editScaleAx,'String');
% corMM=get(handles.editScaleCor,'String');
% scaleFactor=[eval(sagMM),eval(axMM),eval(corMM)];
scaleFactor = [];

unfoldSize=eval(get(handles.editUnfoldSize,'String'));

contents = get(handles.popSpacing,'String');
spacingMethod = contents{get(handles.popSpacing,'Value')};
gridSpacing = str2double(get(handles.editGridSpacing,'String'));

showFigures=get(handles.showFiguresCheck,'Value');
saveExtra=get(handles.saveExtraCheck,'Value');
truePerimDist=get(handles.perimDistCheck,'Value');

set(handles.status,'UserData','unfoldMesh v1.7 2001');
% Go, go ,go !!!!
statusHandle=handles.status;
busyHandle=handles.busyBar;
unfoldMeshFromGUI(meshFileName,grayFileName,flatFileName,startPos,scaleFactor,unfoldSize,statusHandle,busyHandle,spacingMethod,gridSpacing,showFigures,saveExtra,truePerimDist,hemi);

return;

%-------------------------
function retVal=cancel_Callback(h,eventData,handles)
% This button actually says Close (in red) on the screen.
close(handles.figure1);
return;

%-------------------------
function retVal=help_Callback(h,eventData,handles)
web1 = 'web([''file:///'' which(''unfoldMeshHelp.html'')],''-browser'')';
web2 = 'web([''file:///'' which(''unfoldMeshHelp.html'')])';

eval(web1,web2);

return;

%-----------------------------------------------
function p = meshFindPath(handles)
% Method for finding the best current path during selection. This
% code is supposed to reduce the amount of clicking around to find files
%

p = get(handles.inputGrayPath,'String');
if ~isempty(p), p=fileparts(p); return;
else p = get(handles.inputMeshPath,'String');
    if ~isempty(p), p=fileparts(p); return;
    else
        p = fileparts(getAnatomyPath('wandell'));
    end
end
return;

% --- Executes during object creation, after setting all properties.
function xaxisX_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
return;

%-------------------------
function xaxisX_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function xAxisY_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
return;

%-------------------------
function xAxisY_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function xAxisZ_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
return;

function xAxisZ_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function xAxisX_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

return;

function xAxisX_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function popSpacing_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
return;

% --- Executes on selection change in popSpacing.
function popSpacing_Callback(hObject, eventdata, handles)
% Don't do anything.  We read the pop-up when we calculate in the go
% routine above.
return;

% --- Executes during object creation, after setting all properties.
function editGridSpacing_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
return;

function editGridSpacing_Callback(hObject, eventdata, handles)
str2double(get(hObject,'String'));
return;

% --------------------------------------------------------------------
function menuFile_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function menuFileSaveParams_Callback(hObject, eventdata, handles)
% File |Save
% This should get the handles from the edit boxes and save them to a
% user-chosen matlab file

% something like:
% params = mrfGet(handles,'editFields');
% fullFileName = mrSelectFile('w');
% if isempty(fullFileName), disp('User canceled.'); end
% save(fullFileName,'params');%
%

disp('Not yet implemented')
return

% --------------------------------------------------------------------
function menuFileLoadParams_Callback(hObject, eventdata, handles)
% File | Load Params
% This should read a matlab file containing the editable fields and set
% them.
%
% Inverse of the save above with a
% handles = mrfSet(handles,'editFields',params);
%
disp('Not yet implemented')
return

% --------------------------------------------------------------------
function menuFileQuit_Callback(hObject, eventdata, handles)
cancel_Callback(hObject,eventdata,handles);
return;

% --------------------------------------------------------------------
function ii = getPopupValue(h,target)
% Given a handle to a popup, h, find the value
str = get(h,'String');
if ~iscell(str)  % Only one term
    if strcmpi(str,target), ii = 1; return; end
else
    for ii=1:length(str)
        if strcmpi(str{ii},target)
            return;
        end
    end
end

% Returns a safe value, but maybe it should be null?
ii = 1;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function result  = omp_spr(Y, K, epsilon, sparse_output)
    if nargin < 4
        sparse_output = 1;
    end
    Y_normalized = spx.norm.normalize_l2(Y);
    result = mex_omp_spr(Y, Y_normalized, K, epsilon, sparse_output);
end


text:
function [S, V, D, Sigma2] = MySVDtau(A, tau)
[m, n] = size(A);
if 2*m < n
    AAT = A*A';
    [S, Sigma2, D] = svd(AAT);
    Sigma2 = diag(Sigma2);
    V = sqrt(Sigma2);
    tol = max(size(A)) * eps(max(V));
    R = sum(V > max(tol, tau));

    %tol = min(size(A)) * eps(max(V));
    % R = sum(V > 1/tau)

%     tol = min(size(A)) * eps(max(V));
%     R = sum(V > tau)

    V = V(1:R);
    S = S(:,1:R);
    D = A'*S*diag(1./V);
    V = diag(V);
    return;
end
if m > 2*n
    [S, V, D, Sigma2] = MySVDtau(A', tau);
    mid = D;
    D = S;
    S = mid;
    return;
end
[S,V,D] = svd(A);
Sigma2 = diag(V).^2;
R = sum(diag(V) > tau);
S = S(:, 1:R);
V = V(1:R, 1:R);
D = D(:, 1:R);

author: andrewssobral


text:
% Evaluates hamming kNN classifier using 1 upto K nearest neighbors for prediction

function [acc acc2] = eval_hammknn(data, W, K, nonlinearity)

Ntraining = data.Ntraining;
Xtraining = data.Xtraining;
Ltraining = double(data.Ltraining);

Ntest = data.Ntest;
Xtest = data.Xtest;
Ltest = double(data.Ltest);

if (~iscell(W))
  B1 = W * Xtraining;
  B2 = W * Xtest;
else
  resp1 = compute_NN_output(W, Xtraining, nonlinearity);
  B1 = resp1{end};
  resp2 = compute_NN_output(W, Xtest, nonlinearity);
  B2 = resp2{end};
end

B1 = logical(single(B1 > 0));
B1 = compactbit(B1);

B2 = logical(single(B2 > 0));
B2 = compactbit(B2);

[nw1 n1] = size(B1);
[nw2 n2] = size(B2);
nb = nw1*8;

if (nw1 ~= nw2)
  error('nw1 ~= nw2\n');
end

if ~exist('K', 'var')
  K = 3;
end

% Assuming there are only 10 labels
% TODO: change if you have more labels
nlabels = 10;

[ret ret2] = hammknn_mex(nlabels, B1, uint32(Ltraining), B2, uint32(Ltest), n1, nb, K);

acc = double(ret) / Ntest;
acc2 = double(ret2) / Ntraining;

author: norouzi

repo: hdml

save_path: github-repos/MATLAB/norouzi-hdml


text:
function alpha_star=getAlpha_star(mask)
% alpha_star=getAlpha_star(mask) get the pri-known alpha values according
% to the mask. See equation (6) in our iccv2009 paper.
%
% Input arguments:
% mask:   MxN matrix specifying scribbles, with 1 foreground, -1 background
%         and 0 otherwise
%
% Output arguments:
% alpha_star:     (MxN) matrix showing the prio-known value of alpha.
%                 The value is 1 for foreground scribble pixels, and 0
%                 otherwise
%
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% @InProceedings{ZhengICCV09,
%   author = {Yuanjie Zheng and Chandra Kambhamettu},
%   title = {Learning Based Digital Matting},
%   booktitle = {The 20th IEEE International Conference on Computer Vision},
%   year = {2009},
%   month = {September--October}
% }
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% Copyright Yuanjie Zheng @ PICSL @ UPenn on May 9th, 2011
% zheng.vision@gmail.com
% http://sites.google.com/site/zhengvision/
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

disp('Computing preknown alpha values ... ...')

alpha_star=zeros(size(mask,1),size(mask,2));
alpha_star(mask>0)=1;
alpha_star(mask<0)=-1;


text:
function h = makescale(varargin)
%MAKESCALE creates a scale for map data.
%
%   MAKESCALE creates a scale on the current axis based on the current axis
%       limits. The scale is made to occupy 1/5th of the map. It is placed
%       in the southeast corner of the map. The units will either be in
%       milimeters, meters or kilometers, depending on the size of the map.
%
%   MAKESCALE(H_AXIS) creates a scale on the axis specificed by the handle
%       H_AXIS based on the its axes limits. H_AXIS must be a scalar.
%
%   MAKESCALE(SCALE) creates a scale made to occupy 1/SCALE of the map.
%       SCALE must be a scalar, and is bounded to be between 1.1 and 10. If
%       a larger value is passed in, 10 will be used. If a smaller value is
%       passed in, 1.1 will be used.
%
%   MAKESCALE(LOCATION) places the scale in the location specified by
%       LOCATION. Acceptable values for location are as follows
%           'northeast'     'ne'
%           'northwest'     'nw'
%           'southeast'     'se'
%           'southwest'     'sw'
%           'north'         'n'
%           'south'         's'
%
%   MAKESCALE('units',UNITS) changes the units systems from SI to imperical
%       units. UNITS should be either 'si' or 'imp.' The units displayed
%       are automatically switched between milimeters, meters, and
%       kilometers for the SI system, or between inches, feet, and statuate
%       miles for the imperical system.
%
%   H = MAKESCALE(...) outputs H, a 3x1 containing the handles of the of
%       box, line, and text.
%
%   Any number of these input sets may be passed in any order.
%
%   The map scale will automatically be updated as the figure is zoomed,
%       panned, resized, or clicked on. It will not, however, be updated
%       upon using the commands "axis", "xlim", or "ylim" as these do not
%       have callback functionality.
%
%   Example:
%       load conus
%       figure
%       plot(uslon,uslat);
%       makeScale
%
%   Example: Placed in the south
%       load conus
%       figure
%       plot(uslon,uslat);
%       makeScale('south')
%
%   Example: Half the size of the Window
%       load conus
%       figure
%       plot(uslon,uslat);
%       makeScale(2,'south')
%
%   Example: Use Imperical Units
%       load conus
%       figure
%       plot(uslon,uslat);
%       makeScale(2,'south','units','imp')
%
%   Example: Zooming In
%       load conus
%       figure
%       plot(uslon,uslat);
%       makeScale(2,'south')
%       zoom(2)
%
%   Note: This assumes axis limits are in degrees. The scale is sized
%       correctly for the center latitude of the map. As the size of
%       degrees longitude change with latitude, the scale becomes invalid
%       with very large maps. Spherical Earth is assumed. Ideally, the map
%       will be proportioned correctly in order to reflect the relationship
%       between a degree latitude and a degree longitude at the center of
%       the map.
%
% By Jonathan Sullian - October 2011

% Check to make sure the correct number of inputs are passed in.
error(nargchk(0,5,nargin,'struct'));

% Parse Inputs
[anum,latlim,lonlim,scale,location,units] = parseInputs(varargin{:});
if ~isreal(scale)
    error('MAKESCALE:ScaleVal','SCALE must be a real number')
end

% Bound the scale
if scale < 1.1 || scale > 10
    warning('MAKESCALE:ScaleVal','SCALE has been capped to be between 1.1 and 10 for readability.')
end
scale = min(max(scale,1.1),10);
earthRadius = 6371000;

% Get the distance of the map
mlat = mean(latlim);
if abs(mlat) > 90;
    d = 0;
else
    d = earthRadius.*cosd(mlat).*deg2rad(diff(lonlim));
end
dmax = d/1.1;
dmin = d/10;
dlat = diff(latlim);
dlon = diff(lonlim);

% Calculate the distance of the scale bar
rnd2 = floor(log10(d/scale))-1;
dscale = round2(d/scale,10^rnd2);

% Cap it
if dscale > dmax;
    rnd2 = rnd2 - 1;
    dscale = round2(dmax,10^rnd2);
end
if dscale < dmin
    rnd2 = rnd2 - 1;
    dscale = round2(dmin,10^rnd2);
end

% Make the text string
if strcmpi(units,'si')
    if d > 1e3*scale
        dst = num2str(dscale/1e3);
        lbl = ' km';
    elseif d > scale
        dst = num2str(dscale);
        lbl = ' m';
    else
        dst = num2str(dscale*1e3);
        lbl = ' mm';
    end
else
    if d > scale/0.000621371192
        rnd2 = floor(log10(d/scale*0.000621371192))-1;
        dscale = round2(d/scale*0.000621371192,10^rnd2);
        dst = num2str(dscale);
        lbl = ' mi';
        dscale = dscale/0.000621371192;
    elseif d > scale*0.3048
        rnd2 = floor(log10(d/scale/0.3048))-1;
        dscale = round2(d/scale/0.30482,10^rnd2);
        dst = num2str(dscale);
        lbl = ' ft';
        dscale = dscale*.3048;
    else
        rnd2 = floor(log10(d/scale/0.3048*12))-1;
        dscale = round2(d/scale/0.30482*12,10^rnd2);
        dst = num2str(dscale);
        lbl = ' in';
        dscale = dscale/12*.3048;
    end
end

% Get the postions
d1 = [-0.02 0.05];
issouth = 0;
iseast = 0;
iswest = 0;
switch lower(location)
    case {'southeast','se'}
        issouth = 1;
        iseast = 1;
    case {'northeast','ne'}
        iseast = 1;
    case {'southwest','sw'}
        issouth = 1;
        iswest = 1;
    case {'northwest','nw'}
        iswest = 1;
    case {'north','n'}
    case {'south','s'}
        issouth = 1;
end

if issouth
    slat = latlim(1)+0.05*diff(latlim);
else
    slat = latlim(end)-0.08*diff(latlim);
end

if iseast
    slon = lonlim(end)-0.05*diff(lonlim);
    slon = [slon slon-rad2deg(dscale./(earthRadius.*cosd(mlat)))];
    slat = [slat slat];
elseif iswest
    slon = lonlim(1)+0.05*diff(lonlim);
    slon = [slon slon+rad2deg(dscale./(earthRadius.*cosd(mlat)))];
    slat = [slat slat];
    slon = fliplr(slon);
else
    slon = mean(lonlim);
    slon = slon + [-rad2deg(dscale./(earthRadius.*cosd(mlat)))/2 rad2deg(dscale./(earthRadius.*cosd(mlat))/2)];
    slat = [slat slat];
    slon = fliplr(slon);
end

% Get the box location
blat = [slat([2 1])+[d1(1)*dlat d1(2)*dlat] slat([1 2])+[d1(2)*dlat d1(1)*dlat]];
blat = blat([2:4 1]);
blon = [slon+[0.02*dlon -0.02*dlon] slon([2 1])+[-0.02*dlon 0.02*dlon]];

% Delete Old Scale
aold = gca;
axes(anum);
ch = get(anum,'Children');
isOldScale = strcmpi(get(ch,'Tag'),'MapScale');
delete(ch(isOldScale));

% Make the scale
washold = ishold;
hold on
hbox = patch(blon,blat,'w');
set(hbox,'Tag','MapScale');
hline = plot(slon,slat,'k','LineWidth',3);
set(hline,'Tag','MapScale');
units_axis = get(gca,'Units');
set(gca,'Units','Inches')
pos = get(gca,'OuterPosition');
sz = mean(pos(4));
htext = text(mean(blon),mean(blat)+.01*dlat,[dst lbl],'HorizontalAlignment','center','FontSize',sz*2.3);
hzoom = zoom;
hpan = pan(gcf);
set(htext,'Tag','MapScale')
set(gca,'Units',units_axis);

% Set Resizer/Zoom/Pan/Click Callbacks
set(gcf,'ResizeFcn',{@ChangeTextSize,gca,htext});
set(hzoom,'ActionPostCallback',{@remakeZoomPanClick,anum,location,scale,units});
set(hpan,'ActionPostCallback',{@remakeZoomPanClick,anum,location,scale,units});
set(anum,'ButtonDownFcn',{@remakeZoomPanClick,anum,location,scale,units});
axes(aold);

% Output Handles
if nargout > 0
    h = [hbox; hline; htext];
end

% Restore Hold Off
if ~washold
    hold off
end

% Change the text font on figure resize.
function ChangeTextSize(~,~,anum,htext)
units = get(anum,'Units');
set(anum,'Units','Inches')
pos = get(anum,'OuterPosition');
sz = mean(pos(4));
set(htext,'FontSize',sz*2.3);
set(anum,'Units',units);

function remakeZoomPanClick(~,~,anum,location,scale,units)
makescale(anum,location,scale,'units',units);

function x = round2(x,base)
x = round(x./base).*base;

function [anum,latlim,lonlim,scale,location,units] = parseInputs(varargin)
% Default Values
anum = gca;
scale = 5;
location = 'se';
units = 'si';

% Loop through number of arguments in
ii = 1;
while ii <= length(varargin)

    % Either a axis number, or a scale value
    if isscalar(varargin{ii}) && isnumeric(varargin{ii})

        % Is it a non-root handle?
        if ishandle(varargin{ii}) && varargin{ii} ~= 0

            % Check if it is an axis number, not a figure number
            pos = get(varargin{ii},'ActivePositionProperty');
            if strcmpi(pos,'outerposition')
                anum = varargin{1};
                ii = ii + 1;
                continue;
            end
        end

        % Scale Value
        scale = varargin{ii};
        ii = ii + 1;

    % Locations
    elseif ischar(varargin{ii})
        if strcmpi(varargin{ii},'units')
            units = varargin{ii+1};
            ii = ii + 2;
        else
            locs = {'northeast','ne','north','n','southeast','se','south',...
                's','southwest','sw','northwest','nw'};
            if ~ismember(lower(varargin{ii}),locs)
                locOut = 'northeast, ne, north, n, southeast, se, south, s, southwest, sw, northwest, nw';
                error('MAKESCALE:LOCS',['LOCATION must be one of the following: ' locOut])
            end
            location = varargin{ii};
            ii = ii + 1;
        end
    end
end

% Get limits
latlim = get(anum,'YLim');
lonlim = get(anum,'XLim');


text:
function test_suite = test_setopt
initTestSuite;

function test_setopt_all_defaults
  x            = [];
  default_opts = [5 6 7 8];
  z = setopt(x, default_opts);
  z_corr       = [5 6 7 8];
assertVectorsAlmostEqual(z, z_corr, 'relative', 0.0001);

function test_setopt_nonzero_becomes_zero
  x            = [1 0 3];
  default_opts = [5 6 7 8];
  z = setopt(x, default_opts);
  z_corr       = [1 6 3 8];
  %z_corr       = [1 0 3 8];   % This would be more intuitive
assertVectorsAlmostEqual(z, z_corr, 'relative', 0.0001);

author: ricedsp


text:
function w1cw = compute_w1cw(TR, Pulse)
%compute_w1cw Compute the constant wave equivalent power over a period TR for a given pulse

Trf = Pulse.Trf;
omega2 = Pulse.omega2;
if moxunit_util_platform_is_octave
    int = quad(omega2, 0, Trf);
else
    int = integral(omega2, 0, Trf);
end
w1cw = sqrt( int / TR );

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [molStr,molTable] = getMolFormula(model,metabolite)
% getMolFormula takes a metabolite in a specified model and returns the
% molecular formula.

% USAGE:
%
%    [molStr,molTable] = getMolFormula(model,'glu[c]')
%
% INPUTS:
%    model:          COBRA model structure with minimal fields:
%                      * .mets
%                      * .metFormulas
%   metabolite:       a specified metabolite with in the model
%
% OUTPUT:
%   molStr:  the molecular formula for the specified metabolite
%
%   molTable: A cellular array containing the breakdown of the metabolite
%             by element
%
% .. Authors: - Bronson R. Weston   2022

%Get the molStr from metFormulas
ind=find(strcmp(model.mets,metabolite));
if isempty(ind)
    error('Invalid input metabolite!')
end

try
    molStr=model.metFormulas{ind};
catch
    error('A "metFormulas" field is necessary in the input model to use getMolFormula.')
end

%Keep track of each element in the molStr
uC=isstrprop(molStr,'upper');
tracker={};
ind=0;
for i=1:length(uC)
    if uC(i)==1
        ind=ind+1;
        tracker{ind}=molStr(i);

    else
        tracker{ind}=[tracker{ind},molStr(i)];
    end
end

%Create the molTable
molTable={};
for i=1:length(tracker)
    unit=tracker{i};
    num=isstrprop(unit,'digit');
    if ~any(num==1)
        molTable(i,:)={tracker{i} , 1};
    else
        molTable{i,1}=unit(num~=1);
        molTable{i,2}=str2num(unit(num==1));
    end
end
end


text:
function updateScattergroup(obj,scatterIndex)

%check: http://undocumentedmatlab.com/blog/undocumented-scatter-plot-behavior

%----SCATTER FIELDS----%

% x - [DONE]
% y - [DONE]
% r - [HANDLED BY SCATTER]
% t - [HANDLED BY SCATTER]
% mode - [DONE]
% name - [DONE]
% text - [NOT SUPPORTED IN MATLAB]
% error_y - [HANDLED BY ERRORBAR]
% error_x - [NOT SUPPORTED IN MATLAB]
% textfont - [NOT SUPPORTED IN MATLAB]
% textposition - [NOT SUPPORTED IN MATLAB]
% xaxis [DONE]
% yaxis [DONE]
% showlegend [DONE]
% stream - [HANDLED BY PLOTLYSTREAM]
% visible [DONE]
% type [DONE]
% opacity ---[TODO]

% MARKER
% marler.color - [DONE]
% marker.size - [DONE]
% marker.opacity - [NOT SUPPORTED IN MATLAB]
% marker.colorscale - [NOT SUPPORTED IN MATLAB]
% marker.sizemode - [DONE]
% marker.sizeref - [DONE]
% marker.maxdisplayed - [NOT SUPPORTED IN MATLAB]

% MARKER LINE
% marker.line.color - [DONE]
% marker.line.width - [DONE]
% marker.line.dash - [NOT SUPPORTED IN MATLAB]
% marker.line.opacity - [DONE]
% marker.line.smoothing - [NOT SUPPORTED IN MATLAB]
% marker.line.shape - [NOT SUPPORTED IN MATLAB]

% LINE
% line.color - [NA]
% line.width - [NA]
% line.dash - [NA]
% line.opacity [NA]
% line.smoothing - [NOT SUPPORTED IN MATLAB]
% line.shape - [NOT SUPPORTED IN MATLAB]
% connectgaps - [NOT SUPPORTED IN MATLAB]
% fill - [HANDLED BY AREA]
% fillcolor - [HANDLED BY AREA]

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(scatterIndex).AssociatedAxis);

%-SCATTER DATA STRUCTURE- %
scatter_data = get(obj.State.Plot(scatterIndex).Handle);

%-SCATTER CHILDREN-%
scatter_child = get(obj.State.Plot(scatterIndex).Handle,'Children');

%-SCATTER CHILDREN DATA-%
scatter_child_data = get(scatter_child);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

%-------------------------------------------------------------------------%

%-scatter xaxis-%
obj.data{scatterIndex}.xaxis = ['x' num2str(xsource)];

%-------------------------------------------------------------------------%

%-scatter yaxis-%
obj.data{scatterIndex}.yaxis = ['y' num2str(ysource)];

%-------------------------------------------------------------------------%

%-scatter type-%
if any(nonzeros(scatter_data.ZData))
    obj.data{scatterIndex}.type = 'scatter3d';
else
    obj.data{scatterIndex}.type = 'scatter';
end

%-------------------------------------------------------------------------%

%-scatter mode-%
obj.data{scatterIndex}.mode = 'markers';

%-------------------------------------------------------------------------%

%-scatter visible-%
obj.data{scatterIndex}.visible = strcmp(scatter_data.Visible,'on');

%-------------------------------------------------------------------------%

%-scatter name-%
obj.data{scatterIndex}.name = scatter_data.DisplayName;

%-------------------------------------------------------------------------%

%-scatter patch data-%
for m = 1:length(scatter_child_data)

    %reverse counter
    n = length(scatter_child_data) - m + 1;

    %---------------------------------------------------------------------%

    %-scatter x-%
    if length(scatter_child_data) > 1
        obj.data{scatterIndex}.x(m) = scatter_child_data(n).XData;
    else
        obj.data{scatterIndex}.x = scatter_child_data.XData;
    end

    %---------------------------------------------------------------------%

    %-scatter y-%
    if length(scatter_child_data) > 1
        obj.data{scatterIndex}.y(m) = scatter_child_data(n).YData;
    else
        obj.data{scatterIndex}.y = scatter_child_data.YData;
    end

    %---------------------------------------------------------------------%

    %-scatter z-%
    if any(nonzeros(scatter_data.ZData))
        if length(scatter_child_data) > 1
            obj.data{scatterIndex}.z(m) = scatter_child_data(n).ZData;
        else
            obj.data{scatterIndex}.z = scatter_child_data.ZData;
        end
    end

    %---------------------------------------------------------------------%

    %-scatter showlegend-%
    leg = get(scatter_data.Annotation);
    legInfo = get(leg.LegendInformation);

    switch legInfo.IconDisplayStyle
        case 'on'
            showleg = true;
        case 'off'
            showleg = false;
    end

    obj.data{scatterIndex}.showlegend = showleg;

    %---------------------------------------------------------------------%

    %-scatter marker-%
    childmarker = extractPatchMarker(scatter_child_data(n));

    %---------------------------------------------------------------------%

    %-line color-%
    if length(scatter_child_data) > 1
        if iscell(childmarker.line.color)
            obj.data{scatterIndex}.marker.line.color{m} = childmarker.line.color{1};
        else
            obj.data{scatterIndex}.marker.line.color{m} = childmarker.line.color;
        end
    else
        obj.data{scatterIndex}.marker.line.color = childmarker.line.color;
    end

    %---------------------------------------------------------------------%

    %-marker color-%
    if length(scatter_child_data) > 1
        if iscell(childmarker.color)
            obj.data{scatterIndex}.marker.color{m} = childmarker.color{1};
        else
            obj.data{scatterIndex}.marker.color{m} = childmarker.color;
        end
    else
        obj.data{scatterIndex}.marker.color = childmarker.color;
    end

    %---------------------------------------------------------------------%

    %-sizeref-%
    obj.data{scatterIndex}.marker.sizeref = childmarker.sizeref;

    %---------------------------------------------------------------------%

    %-sizemode-%
    obj.data{scatterIndex}.marker.sizemode = childmarker.sizemode;

    %---------------------------------------------------------------------%

    %-symbol-%
    obj.data{scatterIndex}.marker.symbol{m} = childmarker.symbol;

    %---------------------------------------------------------------------%

    %-size-%
    if length(scatter_child_data) > 1 || ischar(childmarker.color)
        obj.data{scatterIndex}.marker.size(m) = childmarker.size;
    else
        obj.data{scatterIndex}.marker.size(1:length(childmarker.color)) = childmarker.size;
    end

    %---------------------------------------------------------------------%

    %-line width-%

    if length(scatter_child_data) > 1 || ischar(childmarker.line.color)
        obj.data{scatterIndex}.marker.line.width(m) = childmarker.line.width;
    else
        obj.data{scatterIndex}.marker.line.width(1:length(childmarker.line.color)) = childmarker.line.width;
    end

    %---------------------------------------------------------------------%

end
end

author: plotly


text:
% StackExchange Signal Processing Q80994
% https://dsp.stackexchange.com/questions/80994
% Building a Pipeline for Image Classification / Clustering Tasks with Features Extractor (Example on MNIST Data)
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     29/12/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

dataSetBaseUrl      = 'http://yann.lecun.com/exdb/mnist/';
imageSetUrl         = 'train-images-idx3-ubyte.gz'; %<! http://yann.lecun.com/exdb/mnist/train-images-idx3-ubyte.gz
annotationSetUrl    = 'train-labels-idx1-ubyte.gz'; %<! http://yann.lecun.com/exdb/mnist/train-labels-idx1-ubyte.gz

imageSetFileName        = 'ImageSet';
annotationSetFileName   = 'AnnotationSet';

fileExt = '.gz';

% Parsing Parameters
trimImg     = 1;
sclaeImg    = 1;

% Train / Test Partition
numTestImg = 5000;

%% Generate / Load Data

if(~exist(strcat(imageSetFileName, fileExt), 'file'))
    websave(strcat(imageSetFileName, fileExt), strcat(dataSetBaseUrl, imageSetUrl));
end
if(~exist(strcat(annotationSetFileName, fileExt), 'file'))
    websave(strcat(annotationSetFileName, fileExt), strcat(dataSetBaseUrl, annotationSetUrl));
end

if(~isfile(imageSetFileName))
    gunzip(strcat(imageSetFileName, fileExt));
end
if(~isfile(annotationSetFileName))
    gunzip(strcat(annotationSetFileName, fileExt));
end

[tImg, vLabel] = ParseMnistData(imageSetFileName, annotationSetFileName, trimImg, sclaeImg);

numRows = size(tImg, 1);
numCols = size(tImg, 2);
numImg = size(tImg, 3);

%% Train / Test Partition

vTestSet    = randperm(numImg, numTestImg);
vTrainSet   = 1:numImg;
vTrainSet(vTestSet) = [];

tTrainSetImg    = tImg(:, :, vTrainSet);
vTrainSetLabel  = vLabel(vTrainSet);
tTestSetImg     = tImg(:, :, vTestSet);
vTestSetLabel   = vLabel(vTestSet);

save('MNISTDataSet', 'tTrainSetImg', 'vTrainSetLabel', 'tTestSetImg', 'vTestSetLabel');

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [dirmap,exy,jxy] = anaskel(wsk);

% Performs some analysis on the skeleton provided,
% which should be thinned to a single pixel width.
%
% dirmap returns a cell array with masks of the skeleton
% pixels that are north, northeast, east, and southeast
% or opposite.
%
% exy and jxy return the pixel coordinates of points that
% are endpoints and junctions of the skeleton, respectively.
%
% Example:
% skg = skeleton(bw);
% wsk = bwmorph((skg > 50),'skel',inf);
% [dmap,exy,jxy] = anaskel(wsk);
%
% See also SKELETON

disp('Compiling...  (You may need to execute ''mex -setup'' first.)');
mex 'G:\Projects\Matlab\anaskel.cpp'


text:
function Xsol = nonlinear_eigenspace(L, k, alpha)
% Example of nonlinear eigenvalue problem: total energy minimization.
%
% function Xsol = nonlinear_eigenspace(L, k, alpha)
%
% L is a discrete Laplacian operator,
% alpha is a given constant, and
% k corresponds to the dimension of the least eigenspace sought.
%
% This example demonstrates how to use the Grassmann geometry factory
% to solve the nonlinear eigenvalue problem as the optimization problem:
%
% minimize 0.5*trace(X'*L*X) + (alpha/4)*(rho(X)*L\(rho(X)))
% over X such that X'*X = Identity,
%
% where L is of size n-by-n,
% X is an n-by-k matrix, and
% rho(X) is the diagonal part of X*X'.
%
% This example is motivated in the paper
% "A Riemannian Newton Algorithm for Nonlinear Eigenvalue Problems",
% Zhi Zhao, Zheng-Jian Bai, and Xiao-Qing Jin,
% SIAM Journal on Matrix Analysis and Applications, 36(2), 752-774, 2015.
%

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Bamdev Mishra, June 19, 2015.
% Contributors:
%
% Change log:

    % If no inputs are provided, generate a  discrete Laplacian operator.
    % This is for illustration purposes only.
    % The default example corresponds to Case (c) of Example 6.2 of the
    % above referenced paper.

    if ~exist('L', 'var') || isempty(L)
        n = 100;
        L = gallery('tridiag', n, -1, 2, -1);
    end

    n = size(L, 1);
    assert(size(L, 2) == n, 'L must be square.');

    if ~exist('k', 'var') || isempty(k) || k > n
        k = 10;
    end

    if ~exist('alpha', 'var') || isempty(alpha)
        alpha = 1;
    end

    % Grassmann manifold description
    Gr = grassmannfactory(n, k);
    problem.M = Gr;

    % Cost function evaluation
    problem.cost =  @cost;
    function val = cost(X)
        rhoX = sum(X.^2, 2); % diag(X*X');
        val = 0.5*trace(X'*(L*X)) + (alpha/4)*(rhoX'*(L\rhoX));
    end

    % Euclidean gradient evaluation
    % Note: Manopt automatically converts it to the Riemannian counterpart.
    problem.egrad = @egrad;
    function g = egrad(X)
        rhoX = sum(X.^2, 2); % diag(X*X');
        g = L*X + alpha*diag(L\rhoX)*X;
    end

    % Euclidean Hessian evaluation
    % Note: Manopt automatically converts it to the Riemannian counterpart.
    problem.ehess = @ehess;
    function h = ehess(X, U)
        rhoX = sum(X.^2, 2); %diag(X*X');
        rhoXdot = 2*sum(X.*U, 2);
        h = L*U + alpha*diag(L\rhoXdot)*X + alpha*diag(L\rhoX)*U;
    end

    % Check whether gradient and Hessian computations are correct.
    % checkgradient(problem);
    % pause;
    % checkhessian(problem);
    % pause;

    % Initialization as suggested in above referenced paper.
    X = randn(n, k);
    [U, S, V] = svd(X, 0); %#ok<ASGLU>
    X = U*V';
    [U0, S0, V0] = eigs(L + alpha*diag(L\(sum(X.^2, 2))), k,'sm'); %#ok<NASGU,ASGLU>
    X0 = U0;

    % Call manoptsolve to automatically call an appropriate solver.
    % Note: it calls the trust regions solver as we have all the required
    % ingredients, namely, gradient and Hessian, information.
    Xsol = manoptsolve(problem, X0);

end