text:
function t_om_solve_nleqs(quiet)
%T_OM_SOLVE_NLEQS  Tests of NLEQ solvers via OM.SOLVE().

%   MP-Opt-Model
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

if nargin < 1
    quiet = 0;
end

core_sp_newton = struct( ...
    'alg',              'NEWTON', ...
    'name',             'Newton''s', ...
    'default_max_it',   10, ...
    'need_jac',         1, ...
    'update_fcn',       @(x, f, J)newton_update_fcn(x, f, J, '')  );
core_sp_gs = struct( ...
    'alg',              'GS', ...
    'name',             'Gauss-Seidel', ...
    'default_max_it',   1000, ...
    'need_jac',         0, ...
    'update_fcn',       @(x, f)x_update_fcn2(x, f)  );

if have_feature('matlab')
    %%  alg         name        check       opts
    cfg = {
        {'DEFAULT', 'default',  []          []  },
        {'NEWTON',  'Newton',   [],         []  },
        {'FD',      'fast-decoupled Newton',[],[]  },
        {'FSOLVE',  'fsolve-1', 'fsolve',   []  },
        {'FSOLVE',  'fsolve-2', 'fsolve',   struct('Algorithm', 'trust-region-dogleg')  },
        {'FSOLVE',  'fsolve-3', 'fsolve',   struct('Algorithm', 'trust-region-reflective')  },
        {'FSOLVE',  'fsolve-4', 'fsolve',   struct('Algorithm', 'levenberg-marquardt', 'TolX', 1e-11) },
        {'GS',      'Gauss-Seidel',[],      []  },
        {'CORE-N',  'Newton-CORE',   [],    core_sp_newton  },
        {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs  },
    };
    if have_feature('matlab', 'vnum') <= 7.010
        cfg([6]) = [];  %% MATLAB 7.10 does not work w/ fsolve alg 3
    end
else    %% octave
    %%  alg         name        check       opts
    cfg = {
        {'DEFAULT', 'default',  []          []  },
        {'NEWTON',  'Newton',   [],         []  },
        {'FD',      'fast-decoupled Newton',[],[]  },
        {'FSOLVE',  'fsolve', 'fsolve',     struct('TolX', 1e-11)  },
        {'GS',      'Gauss-Seidel',[],      []  },
        {'CORE-N',  'Newton-CORE',   [],    core_sp_newton  },
        {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs  },
    };
end

n = 18;

t_begin(14+n*length(cfg), quiet);

for k = 1:length(cfg)
    alg   = cfg{k}{1};
    name  = cfg{k}{2};
    check = cfg{k}{3};
    opts  = cfg{k}{4};
    if ~isempty(check) && ~have_feature(check)
        t_skip(n, sprintf('%s not installed', name));
    else
        opt = struct('verbose', 0, 'alg', alg, 'tol', 1e-11);
        switch alg
            case {'DEFAULT', 'NEWTON'}
            case {'FD'}
                opt.fd_opt.jac_approx_fcn = @jac_approx_fcn2;
            case 'FSOLVE'
                opt.fsolve_opt = opts;
            case 'GS'
                opt.gs_opt.x_update_fcn = @(x, f)x_update_fcn2(x, f);
            case {'CORE-N', 'CORE-GS'}
                opt.core_sp = opts;
                opt.alg = 'CORE';
        end

        switch alg
            case {'DEFAULT', 'NEWTON', 'FSOLVE', 'CORE-N'}
                t = sprintf('%s - 2-d function : ', name);
                x0 = [-1;0];
                om = opt_model;
                om.add_var('x', 2, x0);
                om.add_nln_constraint('f', 2, 1, @f1, []);
                [x, f, e, out, jac] = om.solve(opt);
                t_is(e, 1, 12, [t 'success']);
                t_is(x, [-3; 4], 8, [t 'x']);
                t_is(f, 0, 10, [t 'f']);
                switch alg
                    case {'DEFAULT', 'CORE-N'}
                        out_alg = 'NEWTON';
                    otherwise
                        out_alg = alg;
                end
                t_str_match(out.alg, out_alg, [t 'out.alg']);
                eJ = [1 1; 6 1];
                t_is(jac, eJ, 5.8, [t 'jac']);

                t = sprintf('%s - 2-d function (max_it) : ', name);
                opt.max_it = 3;
                [x, f, e, out] = om.solve(opt);
                t_is(e, 0, 12, [t 'no success']);
                t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
                opt = rmfield(opt, 'max_it');

                t = sprintf('%s - 5-d lin/nonlin function : ', name);
                x0_1 = [-1;0];
                x0_2 = [0;0;0];
                A2 = sparse([2 -1 0; -3 1 -2; 0 5 -4]);
                b2 = [-5; 1; -7];
                x2 = [-2; 1; 3];
                om = opt_model;
                om.add_var('x1', 2, x0_1);
                om.add_var('x2', 3, x0_2);
                om.add_nln_constraint('f', 2, 1, @f1, [], {'x1'});
                om.add_lin_constraint('Ax_b', A2, b2, b2, {'x2'});
                [x, f, e, out, jac] = om.solve(opt);
                t_is(e, 1, 12, [t 'success']);
                t_is(x, [-3; 4; -2; 1; 3], 8, [t 'x']);
                t_is(f, 0, 10, [t 'f']);
                t_str_match(out.alg, out_alg, [t 'out.alg']);
                eJ = [[1 1; 6 1] zeros(2, 3);
                      zeros(3, 2) A2 ];
                t_is(jac, eJ, 5.8, [t 'jac']);
            otherwise
                t_skip(12, sprintf('not implemented for solver ''%s''', alg));
        end

        t = sprintf('%s - 2-d function2 (struct) : ', name);
        x0 = [1;2];
        om = opt_model;
        om.add_var('x', 2, x0);
        om.add_nln_constraint('f', 2, 1, @f2, []);
        [x, f, e] = om.solve(opt);
        t_is(e, 1, 12, [t 'success']);
        t_is(x, [2; 3], 8, [t 'x']);
        t_is(f, 0, 10, [t 'f']);
        t_ok(~isfield(om.soln, 'var'), [t 'no parse_soln() outputs']);

        opt.max_it = 3;
        t = sprintf('%s - 2-d function2 (max_it) : ', name);
        [x, f, e, out] = om.solve(opt);
        t_is(e, 0, 12, [t 'no success']);
        t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
        opt = rmfield(opt, 'max_it');
    end
end

t = 'om.soln.';
opt.alg = 'DEFAULT';
x0_1 = [-1;0];
x0_2 = [0;0;0];
A2 = sparse([2 -1 0; -3 1 -2; 0 5 -4]);
b2 = [-5; 1; -7];
x2 = [-2; 1; 3];
om = opt_model;
om.add_var('x1', 2, x0_1);
om.add_var('x2', 3, x0_2);
om.add_nln_constraint('f', 2, 1, @f1, [], {'x1'});
om.add_lin_constraint('Ax_b', A2, b2, b2, {'x2'});
opt.parse_soln = 1;
[x, f, e, out, jac] = om.solve(opt);
t_is(om.soln.x, x, 14, [t 'x']);
t_is(om.soln.f, f, 14, [t 'f']);
t_is(om.soln.eflag, e, 14, [t 'eflag']);
t_str_match(om.soln.output.alg, out.alg, [t 'output.alg']);
t_is(om.soln.jac, jac, 14, [t 'jac']);

t = 'om.get_soln(''var'', ''x1'') : ';
t_is(om.get_soln('var', 'x1'), x(1:2), 14, [t 'x1']);

t = 'om.get_soln(''var'', ''x'', ''x2'') : ';
t_is(om.get_soln('var', 'x', 'x2'), x(3:5), 14, [t 'x2']);

t = 'om.get_soln(''lin'', ''g'', ''Ax_b'') : ';
g = om.get_soln('lin', 'g', 'Ax_b');
t_is(g{1}, f(3:5), 14, [t 'A * x - u']);
t_is(g{2}, f(3:5), 14, [t 'l - A * x']);

t = 'om.get_soln(''nle'', ''f'') : ';
g = om.get_soln('nle', 'f');
t_is(g, f(1:2), 14, [t 'f']);

t = 'om.get_soln(''nle'', {''g'', ''dg''}, ''f'') : ';
[g, dg] = om.get_soln('nle', {'g', 'dg'}, 'f');
t_is(g, f(1:2), 14, [t 'f']);
t_is(dg, jac(1:2, 1:2), 14, [t 'jac']);

t = 'parse_soln : ';
t_is(om.soln.var.val.x1, om.get_soln('var', 'x1'), 14, [t 'var.val.x1']);
t_is(om.soln.var.val.x2, om.get_soln('var', 'x2'), 14, [t 'var.val.x2']);

t_end;

%% 2-d problem with 2 solutions
%% from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/
function [f, J] = f1(x)
if iscell(x)    %% in case it was part of a varset
    x = x{1};
end
f = [  x(1)   + x(2) - 1;
      -x(1)^2 + x(2) + 5    ];
if nargout > 1
    J = sparse([1 1; -2*x(1) 1]);
end

%% another 2-d problem
%% from Christi Patton Luks, https://www.youtube.com/watch?v=pJG4yhtgerg
function [f, J] = f2(x)
f = [  x(1)^2 +   x(1)*x(2)   - 10;
       x(2)   + 3*x(1)*x(2)^2 - 57  ];
if nargout > 1
    J = sparse([    2*x(1)+x(2)    x(1);
                    3*x(2)^2       6*x(1)*x(2)+1    ]);
end

function JJ = jac_approx_fcn2()
%% for use with fast-decoupled Newton's method
J = [7 2; 27 37];
JJ = {J(1,1), J(2,2)};

function x = x_update_fcn2(x, f)
%% for use with Gauss-Seidel method
x(1) = sqrt(10 - x(1)*x(2));
x(2) = sqrt((57-x(2))/3/x(1));

function x = newton_update_fcn(x, f, J, lin_solver)
dx = mplinsolve(J, -f, lin_solver);     %% compute update step
x = x + dx;                             %% update x

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function varargout = drawTorus(varargin)
%DRAWTORUS Draw a torus (3D ring).
%
%   drawTorus(TORUS)
%   Draws the torus on the current axis. TORUS is given by:
%   [XC YC ZY  R1 R2  THETA PHI]
%   where (XC YZ ZC) is the center of the torus, R1 is the main radius, R2
%   is the radius of the torus section, and (THETA PHI) is the angle of the
%   torus normal vector (both in degrees).
%
%   drawTorus(..., PNAME, PVALUE)
%   Specifies a set of parameter name-value pairs. Parameter names include
%   plitting options ('facecolor', 'linestyle'...), or options specific to
%   torus:
%   'nPhi'      number of meridians used to draw the torus (default is 60).
%   'nTheta'    number of parallels used to draw the torus (default is 60).
%
%
%   Example
%     % draw sample torus
%     figure;
%     drawTorus([50 50 50 30 10 30 45]);
%     axis equal; view([95 10]); light;
%
%   See also
%   drawEllipsoid, revolutionSurface, torusMesh
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2011-06-22, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

%% Default values

% number of meridians
nPhi    = 60;

% number of parallels
nTheta  = 60;

%% Extract input arguments

% Check if axes handle is specified
if isAxisHandle(varargin{1})
    hAx = varargin{1};
    varargin(1) = [];
else
    hAx = gca;
end

torus = varargin{1};
varargin(1) = [];

center = torus(1:3);
r1 = torus(4);
r2 = torus(5);

normal = [0 0];
if size(torus, 2) >= 7
    normal = torus(6:7);
end

% default set of options for drawing meshes
options = {'FaceColor', 'g', 'linestyle', 'none'};

while length(varargin) > 1
    switch lower(varargin{1})
        case 'nphi'
            nPhi = varargin{2};

        case 'ntheta'
            nTheta = varargin{2};

        otherwise
            % assumes this is drawing option
            options = [options varargin(1:2)]; %#ok<AGROW>
    end

    varargin(1:2) = [];
end

%% Draw the torus

% create base torus
circle = circleToPolygon([r1 0 r2], nTheta);
[x, y, z] = revolutionSurface(circle, linspace(0, 2*pi, nPhi));

% transform torus
trans = localToGlobal3d([center normal]);
[x, y, z] = transformPoint3d(x, y, z, trans);

% draw the surface
hs = surf(hAx, x, y, z, options{:});

%% Process output arguments

if nargout > 0
    varargout = {hs};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function rename_dicm(files, fmt)
% Rename dicom files so the names are meaningful to human.
%
% RENAME_DICM(files, outputNameFormat)
%
% The first input is the dicom file(s) or a folder containing dicom files.
% The second input is the format for the result file names. Support format
% include:
%
% 1: run1_00001.dcm (SeriesDescription_instance). It is the shortest name, but
%    if there is MoCo series, or users did not change run names, there will be
%    name conflict.
%
% 2: Smith^John-0004-0001-00001.dcm (PatientName-Series-Acquisition-Instance).
%    This is the BrainVoyager format. It has no name confict, but it is long and
%    less descriptive. Note that BrainVoyager itself has problem to distinguish
%    the two series of images for Siemens fieldmap, while this code can avoid
%    this problem.
%
% 3: run1_004_00001.dcm (SeriesDescription_Series_Instance). This gives short
%    names, while it is descriptive and there is no name conflict most of time.
%
% 4: 2334ZL_run1_00001.dcm (PatientName_SeriesDescription_instance). This may be
%    useful if files for different subjects are in the same folder.
%
% 5: run1_003_001_00001.dcm (SeriesDescription_Series_acquisition_instance).
%    This ensures no name conflict, and is the default.
%
% Whenever there is name confict, you will see red warning and the remaining
% files won't be renamed.
%
% If the first input is not provided or empty, you will be asked to pick up
% a folder.
%
% See also DICM_HDR, SORT_DICM, ANONYMIZE_DICM

% History (yymmdd):
% 0710?? Write it (Xiangrui Li)
% 1304?? Add more options for output format
% 1306?? Exclude PhoenixZIPReport files to avoid error
% 1306?? Fix problem if illegal char in ProtocolName
% 1309?? Use dicm_hdr to replace dicominfo, so it runs much faster
% 1309?? Use 5-digit InstanceNumber, so works better for GE/Philips
% 1402?? Add Manufacturer to flds (bug caused by dicm_hdr update)
% 140506 Use SeriesDescription to replace ProtocolName non-Siemens
% 151001 Avoid cd so it works if m file is at pwd but path not set
% 171211 Make AcquisitionNumer and Manufacturer not mandidate.

if nargin<1 || isempty(files)
    folder = uigetdir(pwd, 'Select a folder containing DICOM files');
    if folder==0, return; end
    files = dir(folder);
    files([files.isdir]) = [];
    files = {files.name};

    str = sprintf(['Choose Output format: \n\n' ...
                   '1: run1_00001.dcm (SeriesDescription_instance)\n' ...
                   '2: BrainVoyager format (subj-series-acquisiton-instance)\n' ...
                   '3: run1_001_00001.dcm (SeriesDescription_series_instance)\n' ...
                   '4: subj_run1_00001.dcm (subj_SeriesDescription_instance)\n' ...
                   '5: run1_001_001_00001.dcm (SeriesDescription_series_acquisition_instance)\n']);
    fmt = inputdlg(str, 'Rename Dicom', 1, {'5'});
    if isempty(fmt), return; end
    fmt = str2double(fmt{1});
else
    if exist(files, 'dir') % input is folder
        folder = files;
        files = dir(folder);
        files([files.isdir]) = [];
        files = {files.name};
    else % files
        if ~iscell(files), files = {files}; end
        folder = fileparts(files{1});
        if ~isempty(folder), folder = pwd; end
    end
    if nargin<2 || isempty(fmt), fmt = 5; end
end

flds = {'InstanceNumber' 'AcquisitionNumber' 'SeriesNumber' 'EchoNumber' 'ProtocolName' ...
        'SeriesDescription' 'PatientName' 'PatientID' 'Manufacturer'};
dict = dicm_dict('', flds);
tryGetField = dicm2nii('', 'tryGetField', 'func_handle');

nFile = numel(files);
if nFile<1, return; end
if ~strcmp(folder(end), filesep), folder(end+1) = filesep; end
err = '';
str = sprintf('%g/%g', 1, nFile);
more off;
fprintf(' Renaming DICOM files: %s', str);

for i = 1:nFile
    fprintf(repmat('\b', [1 numel(str)]));
    str = sprintf('%g/%g', i, nFile);
    fprintf('%s', str);
    s = dicm_hdr([folder files{i}], dict);
    vendor = tryGetField(s, 'Manufacturer', '');
    try % skip if no these fields
        sN = s.SeriesNumber;
        aN = tryGetField(s, 'AcquisitionNumber', 1);
        iN = s.InstanceNumber;
        if fmt ~= 2
            if strncmp(vendor, 'SIEMENS', 7)
                pName = strtrim(s.ProtocolName);
            else
                pName = strtrim(s.SeriesDescription);
            end
            pName(~isstrprop(pName, 'alphanum')) = '_'; % valid file name
            pName = regexprep(pName, '_{2,}', '_');
        end
        if fmt==2 || fmt==4
            sName = tryGetField(s, 'PatientName');
            if isempty(sName), sName = tryGetField(s, 'PatientID'); end
            sName(~isstrprop(sName, 'alphanum')) = '_';
            sName = regexprep(sName, '_{2,}', '_');
        end
    catch me %#ok
        continue;
    end

    if strncmpi(vendor, 'Philips', 7) % SeriesNumber is useless
        sN = aN;
    elseif strncmpi(vendor, 'SIEMENS', 7) && tryGetField(s, 'EchoNumber', 1)>1
        aN = s.EchoNumber; % fieldmap phase image
    end

    if fmt == 1 % pN_001
        name = sprintf('%s_%05g.dcm', pName, iN);
    elseif fmt == 2 % BrainVoyager
        name = sprintf('%s-%04g-%04g-%05g.dcm', sName, sN, aN, iN);
    elseif fmt == 3 % pN_03_00001
        name = sprintf('%s_%02g_%05g.dcm', pName, s.SeriesNumber, iN);
    elseif fmt == 4 % 2322ZL_pN_001
        name = sprintf('%s_%s_%05g.dcm', sName, pName, iN);
    elseif fmt == 5 % pN_003_001_001
        name = sprintf('%s_%03g_%03g_%05g.dcm', pName, sN, aN, iN);
    else
        error('Invalid format.');
    end

    if strcmpi(files{i}, name), continue; end % done already

    if ispc, cmd = ['rename "' folder files{i} '" ' name];
    else, cmd = ['mv "' folder files{i} '" "' folder name '"'];
    end % matlab movefile is too slow

    [er, foo] = system(cmd);
    if er, err = [err files{i} ': ' foo]; end %#ok
end
fprintf('\n');
if ~isempty(err), fprintf(2, '\n%s\n', err); end

author: xiangruili


text:
function drawCameraPly(PLYfilename,Rt,scale,w,h,f)

if ~exist('w','var')
    w = 640;
end
if ~exist('h','var')
    h = 480;
end
if ~exist('f','var')
    f = 570.3422090067766703214147128164768218994140625;
end
if ~exist('scale','var')
    scale = 0.01;
end

camera=[...
0 -w/2 +w/2 +w/2 -w/2
0 -h/2 -h/2 +h/2 +h/2
0   f    f    f    f];

camera = camera * scale;
camera = Rt(:,1:3) * camera + repmat(Rt(:,4),1,5);
edges = [1 2; 1 3; 1 4; 1 5; 2 3;3 4;4 5;5 2]';

graph2ply(PLYfilename, camera, edges);

%{
plot3(camera(1,[2 3 4 5 2]),camera(2,[2 3 4 5 2]),camera(3,[2 3 4 5 2]),'-k'); hold on;
for i=2:5
    plot3(camera(1,[1 i]),camera(2,[1 i]),camera(3,[1 i]),'-k'); hold on;
end
axis equal
xlabel('x');
ylabel('y');
zlabel('z');
%}

author: jianxiongxiao

repo: ProfXkit


text:
function colors = colorcube_colors(varargin)

wh = [10 11 13 17 18 21 22 23 28 30 32 24 36];

nofigs = isempty(findobj(0, 'Type', 'Figure'));

cm = colormap('colorcube');
cm = [cm(wh, :); cm];

% Maybe different versions of mat2cell? Not so stable...
nc = ones(size(cm, 1), 1);
try
    cm = mat2cell(cm, nc, 3)';
catch
    cm = mat2cell(cm, ones(nc, 1))';
end

colors = cm;

if nofigs
    % we created one for colormap, so close
    close
end

if ~isempty(varargin)
    n = varargin{1};

    while length(colors) < n, colors = [colors colors]; end

    colors = colors(1:n);
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function Y = ttm(X,V,varargin)
%TTM Tensor times matrix.
%
%   Y = TTM(X,A,N) computes the n-mode product of tensor X with a
%   matrix A; i.e., X x_N A.  The integer N specifies the dimension
%   (or mode) of X along which A should be multiplied.  If size(A) =
%   [J,I], then X must have size(X,N) = I.  The result will be the
%   same order and size as X except that size(Y,N) = J.
%
%   Y = TTM(X,{A,B,C,...}) computes the n-mode product of tensor X
%   with a sequence of matrices in the cell array.  The n-mode
%   products are computed sequentially along all dimensions (or modes)
%   of X. The cell array contains ndims(X) matrices.
%
%   Y = TTM(X,{A,B,C,...},DIMS) computes the sequence tensor-matrix
%   products along the dimensions specified by DIMS.
%
%   Y = TTM(...,'t') performs the same computations as above except
%   the matrices are transposed.
%
%   Examples
%   X = tensor(rand(5,3,4,2));
%   A = rand(4,5); B = rand(4,3); C = rand(3,4); D = rand(3,2);
%   Y = ttm(X, A, 1)         %<-- computes X times A in mode-1
%   Y = ttm(X, {A,B,C,D}, 1) %<-- same as above
%   Y = ttm(X, A', 1, 't')   %<-- same as above
%   Y = ttm(X, {A,B,C,D}, [1 2 3 4]) %<-- 4-way multiply
%   Y = ttm(X, {D,C,B,A}, [4 3 2 1]) %<-- same as above
%   Y = ttm(X, {A,B,C,D})            %<-- same as above
%   Y = ttm(X, {A',B',C',D'}, 't')   %<-- same as above
%   Y = ttm(X, {C,D}, [3 4])     %<-- X times C in mode-3 & D in mode-4
%   Y = ttm(X, {A,B,C,D}, [3 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, [1 2 4])   %<-- 3-way multiply
%   Y = ttm(X, {A,B,C,D}, [1 2 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, -3)        %<-- same as above
%   Y = ttm(X, {A,B,C,D}, -3)      %<-- same as above
%
%   See also TENSOR, TENSOR/TTT, TENSOR/TTV.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Check the number of arguments
if (nargin < 2)
    error('TTM requires at least two arguments.');
end

%% Create 'n' and 'tflag' arguments from varargin
n = 1:ndims(X);
tflag = '';
if numel(varargin) == 1
    if ischar(varargin{1})
        tflag = varargin{1};
    else
        n = varargin{1};
    end
elseif numel(varargin) == 2
    n = varargin{1};
    tflag = varargin{2};
end

%% Handle cell array
if iscell(V)
    % Copy n into dims
    dims = n;
    % Check that the dimensions are valid
    [dims,vidx] = tt_dimscheck(dims,ndims(X),numel(V));
    % Calculate individual products
    Y = ttm(X, V{vidx(1)}, dims(1), tflag);
    for i = 2 : numel(dims)
        Y = ttm(Y, V{vidx(i)}, dims(i), tflag);
    end
    % All done
    return;
end

%% Check the second argument
if ndims(V) ~= 2
    error('tensor/ttm: 2nd argument must be a matrix.');
end

%% Check n
if (numel(n) ~= 1 || (n < 0) || n > ndims(X))
    error('Dimension N must be between 1 and NDIMS(X).');
end

%% COMPUTE THE PRODUCT

N = ndims(X);
sz = size(X);
order = [n,1:n-1,n+1:N];
newdata = double(permute(X,order));
newdata = reshape(newdata,sz(n),prod(sz([1:n-1,n+1:N])));
if tflag == 't'
    newdata = V'*newdata;
    p = size(V,2);
else
    newdata = V*newdata;
    p = size(V,1);
end
newsz = [p,sz(1:n-1),sz(n+1:N)];
Y = tensor(newdata,newsz);
Y = ipermute(Y,order);

return;

end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function f = compose(f, op)
%COMPOSE   Compose command for SPHEREFUNV objects.
%   F = COMPOSE(F, G) returns the composition G(F) of the SPHEREFUNV object
%   F and a CHEBFUN3 or CHEBFUN3V with three components G.

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

% Deal with empty CHEBFUN objects:
if ( isempty(f) || isempty(op) )
    if ( isa(op, 'chebfun3v') )
        f = spherefunv();
    else
        f = spherefun();
    end
end

% F is real and has three components, so we can compose with a CHEBFUN3 or
% CHEBFUN3V.
% TODO? compose when op is a spherefun object?

if ( isa(op, 'chebfun3') )
    % Get components of F:
    f1 = f.components{1};
    f2 = f.components{2};
    f3 = f.components{3};

    % Check that image(f) is in domain(op):
    rangef = minandmax2est(f);              % Estimate of image(f).
    tol = 100 * chebfun2eps * max(vscale(f), vscale(op)) * ...
            norm(f1.domain, inf);           % Tolerance.
    if ( ~isSubset(rangef, op.domain, tol) )
        error('CHEBFUN:SPHEREFUNV:COMPOSE:DomainMismatch3', ...
            'OP(F) is not defined, since image(F) is not contained in domain(OP).')
    end

    % Call constructor:
    f = spherefun(@(x,y,z) op(feval(f1, x, y, z), feval(f2, x, y, z), ...
        feval(f3, x, y, z)));

elseif ( isa(op, 'chebfun3v') )
    % Check that OP has three components:
    if ( op.nComponents ~= 3 )
        error('CHEBFUN:SPHEREFUNV:compose:nComponents', ...
            'CHEBFUN3V(SPHEREFUNV) is defined for CHEBFUN3V objects with 3 components.')
    end

    % Get components of op:
    op1 = op(1);
    op2 = op(2);
    op3 = op(3);

    % Get components of f:
    f1 = f.components{1};
    f2 = f.components{2};
    f3 = f.components{3};

    % Check that image(f) is in domain(op):
    rangef = minandmax2est(f);              % Estimate of image(f).
    tol = 100 * chebfun2eps * max(vscale(f), vscale(op)) * ...
            norm(f1.domain, inf);           % Tolerance.
    if ( ~isSubset(rangef, op1.domain, tol) )
        error('CHEBFUN:SPHEREFUNV:COMPOSE:DomainMismatch3v', ...
            'OP(F) is not defined, since image(F) is not contained in domain(OP).')
    end

    % Calling directly the constructor for [ op1(f), op2(f), op3(f) ] as below
    % seems to be slightly faster than building a spherefun for each component
    % and then calling the spherefunv constructor on the three spherefuns.

    % Call constructor:
    f = spherefunv(@(x,y,z) op1(feval(f1, x, y, z), feval(f2, x, y, z), ...
        feval(f3, x, y, z)), ...
        @(x,y,z) op2(feval(f1, x, y, z), feval(f2, x, y, z), ...
        feval(f3, x, y, z)), ...
        @(x,y,z) op3(feval(f1, x, y, z), feval(f2, x, y, z), ...
        feval(f3, x, y, z)));

else
    error('CHEBFUN:SPHEREFUNV:COMPOSE:OP', ...
        'Can compose a SPHEREFUNV with a CHEBFUN3 or CHEBFUN3V.')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%example_accel_grad_mtl

clear all;
clc

% add path
addpath(genpath('../../SLEP/'));

% load data and set regularization parameter
load('../../data/dmoz.mat');
lambda = 10^-4;

% center data
for i = 1:length(Xtrain)
        Xtrain{i} = CenterRowData(Xtrain{i});
        Ytrain{i} = CenterRowData(Ytrain{i});
end

% call the main function
[Wp,fval_vec,itr_counter] = accel_grad_mtl(Xtrain,Ytrain,lambda);

author: zzstefan

repo: BrainNetClass


text:
function handle = xyzankurVisualise2(pos,fid)

% XYZANKURVISUALISE2
%
% COPYRIGHT : Carl Henrik Ek and Neil Lawrence, 2008

% MOCAP

if(nargin<2)
  fid = 1;
end

figure(fid);

joint = xyzankur2joint(pos);
handle{1} = xyzankurDraw(joint);axis equal;view([1 0 0]);
set(gca,'XLim',[-20 20],'YLim',[-20 20],'ZLim',[-5 70]);

figure(fid+2);
joint = xyzankur2joint(pos);
handle{2} = xyzankurDraw(joint);axis equal;view([1 1 1]);
set(gca,'XLim',[-20 20],'YLim',[-20 20],'ZLim',[-5 70]);

return;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef StiffnessDescriptor < FieldNameDescriptor

    methods (Access = protected)

        function loadFieldVariable(obj)
            obj.fieldName = 'stiffness';
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function v = sparse(i,j,v,varargin)

v.x = sparse(i,j,v.x,varargin{:});
v.y = sparse(i,j,v.y,varargin{:});
v.z = sparse(i,j,v.z,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function segment_length = p05_boundary_segment_length ( segment_index, h )

%*****************************************************************************80
%
%% P05_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
%
%    Input, real H, the suggested spacing between points.
%
%    Output, integer SEGMENT_LENGTH, the number of points in the segment.
%
  center1 = [  0.0, 0.0 ];
  center2 = [ -0.4, 0.0 ];
  r1 = 1.00;
  r2 = 0.55;

  if ( h <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P05_BOUNDARY_SEGMENT_LENGTH - Fatal error!\n' );
    fprintf ( 1, '  Nonpositive H = %f\n', h );
    error ( 'P05_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
  end

  if ( segment_index == 1 )

    n = round ( ( pi * ( r1  + r2 ) ...
      + ( ( center2(1) - r2 ) - ( center1(1) - r1 ) ) ...
      + ( ( center1(1) + r1 ) - ( center2(1) + r2 ) ) ) / h );
    n = max ( n, 21 );
    segment_length = n;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P05_BOUNDARY_SEGMENT_LENGTH - Fatal error!\n' );
    fprintf ( 1, '  Illegal SEGMENT_INDEX = %d\n', segment_index );
    error ( 'P05_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [B, Y] = compressSELVE(X, SELVEparam)
%function Y = compressSELVE(X, Anchor, lambda, s, sigma, tempResults)

Anchor = SELVEparam.anchor;
s = SELVEparam.s;
lambda = SELVEparam.lambda;
sigma = SELVEparam.sigma;
tempResults = SELVEparam.tempResults;

[n,dim] = size(X);
m = size(Anchor,1);

%% get Z
Z = zeros(n,m);
Dis = sqdist(X',Anchor');
clear X;
clear Anchor;

val = zeros(n,s);
pos = val;
for i = 1:s
    [val(:,i),pos(:,i)] = min(Dis,[],2);
    tep = (pos(:,i)-1)*n+[1:n]';
    Dis(tep) = 1e60;
end
clear Dis;
clear tep;
val = exp(-val/(1/1*sigma^2));
val = repmat(sum(val,2).^-1,1,s).*val; %% normalize
tep = (pos-1)*n+repmat([1:n]',1,s);
Z([tep]) = [val];
Z = sparse(Z);
clear tep;
clear val;
clear pos;
lamda1 = sum(Z);
Z = diag(lamda1.^-0.5)*Z'; %Z : fea * ins
Z = Z';

tempS = tempResults.T*Z' + repmat(tempResults.beta,1,size(Z,1));
Y = double(tempS > repmat(tempResults.Ui,1,n));
Y = Y';
B = compactbit(Y);

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
% GLG
%
% Files
%   cell_to_dwt2 - Divide wavedec2 output to cell
%   dwt2_to_cell - Divide wavedec2 output to cell
%   dwt_to_cell  - Divide wavedec output to cell
%   psnr         - Compute PSNR between images
%   setup        - Add the relevant sub directories

author: Sable

repo: mcbench-benchmarks


text:
function program_05 ( )

%*****************************************************************************80
%
%% PROGRAM_05 displays the points that sample a triangle.
%
%  Discussion:
%
%    This program is similar to program_04, but displays the
%    sample points.
%
%    The program
%    * reads a triangle T (defined by three points),
%    * reads a random number seed;
%    * reads "1" for bad scheme, "2" for good scheme.
%    * reads N, the number of random values to generate;
%    * it then computes N random points in the triangle;
%    * it determines N1, N2, and N3, the number of points
%      that ended up in subtriangles mbc, amc, and abm,
%      where "m" is the centroid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROGRAM_05 - The Eyeball Test on Sampling\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define a triangle T:\n' );

  t_v1 = input ( '  Enter [ T.v1.x, T.v1.y]: ' );
  t_v2 = input ( '  Enter [ T.v2.x, T.v2.y]: ' );
  t_v3 = input ( '  Enter [ T.v3.x, T.v3.y]: ' );
%
%  Get the random seed.
%
  seed = input ( 'Enter a random number seed:  ' );
  rand ( 'state', seed );
%
%  Get the scheme to use.
%
  scheme = input ( 'Enter 1 for "bad scheme", 2 for "good scheme":  ' );
%
%  Get the number of values to generate.
%
  n = input ( 'Enter the number of samples to generate:  ' );

  n1 = 0;
  n2 = 0;
  n3 = 0;

  p = zeros ( 2, n );

  for i = 1 : n

    if ( scheme == 1 )
      r = rand;
      s = rand;
      t = rand;
      xi1 = r / ( r + s + t );
      xi2 = s / ( r + s + t );
      xi3 = t / ( r + s + t );
    elseif ( scheme == 2 )
      r = rand;
      s = rand;

      xi1 =   1.0       - sqrt ( s );
      xi2 = ( 1.0 - r ) * sqrt ( s );
      xi3 =         r   * sqrt ( s );
    end

    if ( xi1 < xi2 & xi1 < xi3 )
      n1 = n1 + 1;
    elseif ( xi2 < xi1 & xi2 < xi3 )
      n2 = n2 + 1;
    elseif ( xi3 < xi1 & xi3 < xi2 )
      n3 = n3 + 1;
    end

    p(1:2,i) = xi1 * t_v1(1:2) ...
             + xi2 * t_v2(1:2) ...
             + xi3 * t_v3(1:2);

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N1 = %8d  %% = %8.4f\n', n1, 100 * n1 / n );
  fprintf ( 1, '  N2 = %8d  %% = %8.4f\n', n2, 100 * n2 / n );
  fprintf ( 1, '  N3 = %8d  %% = %8.4f\n', n3, 100 * n3 / n );
  fprintf ( 1, '  N  = %8d  %% = %8.4f\n', n, 100 );

  scatter ( p(1,:), p(2,:), 'filled' );
  m = ( t_v1 + t_v2 + t_v3 ) / 3;

  line ( [ t_v1(1), t_v2(1) ], [ t_v1(2), t_v2(2) ], 'Color', 'r', 'LineWidth', 2.0 );
  line ( [ t_v2(1), t_v3(1) ], [ t_v2(2), t_v3(2) ], 'Color', 'r', 'LineWidth', 2.0 );
  line ( [ t_v3(1), t_v1(1) ], [ t_v3(2), t_v1(2) ], 'Color', 'r', 'LineWidth', 2.0 );

  line ( [ t_v1(1), m(1) ], [ t_v1(2), m(2) ], 'Color', 'r', 'LineWidth', 2.0 );
  line ( [ t_v2(1), m(1) ], [ t_v2(2), m(2) ], 'Color', 'r', 'LineWidth', 2.0 );
  line ( [ t_v3(1), m(1) ], [ t_v3(2), m(2) ], 'Color', 'r', 'LineWidth', 2.0 );

  p_min(1) = min ( p(1,:) );
  p_min(1) = min ( p_min(1), t_v1(1) );
  p_min(1) = min ( p_min(1), t_v2(1) );
  p_min(1) = min ( p_min(1), t_v3(1) );

  p_max(1) = max ( p(1,:) );
  p_max(1) = max ( p_max(1), t_v1(1) );
  p_max(1) = max ( p_max(1), t_v2(1) );
  p_max(1) = max ( p_max(1), t_v3(1) );

  p_min(2) = min ( p(2,:) );
  p_min(2) = min ( p_min(2), t_v1(2) );
  p_min(2) = min ( p_min(2), t_v2(2) );
  p_min(2) = min ( p_min(2), t_v3(2) );

  p_max(2) = max ( p(2,:) );
  p_max(2) = max ( p_max(2), t_v1(2) );
  p_max(2) = max ( p_max(2), t_v2(2) );
  p_max(2) = max ( p_max(2), t_v3(2) );

  p_range(1) = p_max(1) - p_min(1);
  p_range(2) = p_max(2) - p_min(2);

  margin = 0.025 * max ( p_range(1), p_range(2) );

  x_min = p_min(1) - margin;
  x_max = p_max(1) + margin;
  y_min = p_min(2) - margin;
  y_max = p_max(2) + margin;
%
%  The TITLE function will interpret underscores in the title.
%  We need to unescape such escape sequences!
%
  title_string = 'Sample points in triangle';
  title ( title_string )

  axis ( [ x_min, x_max, y_min, y_max ] );
  axis equal
%
%  Save data to file.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving data to file "triangle_samples.txt".\n' );

  save 'triangle_samples.txt' p -ASCII
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROGRAM_05\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function varargout = spm_rmpath(d)
% Recursively removes SPM paths from the MATLAB path
%   SPM_RMPATH checks if the file spm.m is found and removes the
%   path to that file and any subdirectories below it from the MATLAB
%   path.
%
%   P = SPM_RMPATH performs the same function as above and returns the
%   cleaned path string in P.
%
%   SPM_RMPATH(D) strips the path string D from the MATLAB path.
%
%   P = SPM_RMPATH(D) strips the path string D from the MATLAB path and
%   returns the cleaned path string in P.
%
%   See also PATH, ADDPATH, RMPATH, GENPATH, PATHTOOL, SAVEPATH.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Darren Gitelman & Guillaume Flandin
% $Id: spm_rmpath.m 4025 2010-07-29 11:10:15Z guillaume $

varargout = {};
if ~nargin
    % Get the actual SPM directory
    try, d = spm('dir');
    catch, return; end
end

% Recursively remove directories in the MATLAB path
p = textscan(path,'%s','delimiter',pathsep); p = p{1};
i = strncmp(d,p,length(d)); P = p(i); p(i) = [];
if ~nargin && ~isempty(P)
    fprintf('Removed %s paths starting from base path: "%s"\n',spm('ver','',1),d);
elseif ~isempty(P)
    fprintf('Removed paths starting from base path: "%s" from:\n',d);
else
    fprintf('No matching path strings found to remove\n')
end
if numel(P), fprintf('\t%s\n',P{:}); end

% Set the new MATLAB path
p = strcat(p,pathsep);
path(strcat(p{:}));

% Return the cleaned path if requested
if nargout
    varargout{1} = p;
end

author: spm


text:
function value = year_is_leap_bahai ( y )

%*****************************************************************************80
%
%% YEAR_IS_LEAP_BAHAI returns TRUE if the Bahai year was a leap year.
%
%  Discussion:
%
%    The leap year rules are the same as those used in the Gregorian
%    calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, logical VALUE, TRUE if the year was a leap year,
%    FALSE otherwise.
%
  if ( y <= 0 )
    value = 0;
    return
  end

  if ( mod ( y, 400 ) == 0 )
    value = 1;
  elseif ( mod ( y, 100 ) == 0 )
    value = 0;
  elseif ( mod ( y, 4 ) == 0 )
    value = 1;
  else
    value = 0;
  end

  return
end


text:
function out = chebcoeffs(f, N, kind)
%CHEBCOEFFS   Chebyshev polynomial coefficients of a TRIGTECH.
%   A = CHEBCOEFFS(F, N) or A = CHEBCOEFFS(F, N, 1) returns a column vector of
%   the first N coefficients in the expansion of F in a series of Chebyshev
%   polynomials of the first kind, ordered starting with the coefficient of
%   T_0(x).
%
%   A = CHEBCOEFFS(F, N, 2) does the same but for a series expansion in
%   Chebyshev polynomials of second kind, ordered starting with the coefficient
%   of U_0(x).
%
%   If F is array-valued with P columns, then A is an NxP matrix.
%
% See also LEGCOEFFS, TRIGCOEFFS.

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

if ( (nargin < 2) || isempty(N) )
    error('CHEBFUN:TRIGTECH:chebcoeffs:input', ...
        'F does not have a finite Chebyshev series. Please input N.');
end

% Use first-kind Chebyshev polynomials by default.
if ( (nargin < 3) || isempty(kind) )
    kind = 1;
end

% Trivial empty case:
if ( N <= 0 )
    out = [];
    return
end

% [TODO]: Is there a fast transfrom from TRIGTECH to CHEBTECH?
% Since f is a TRIGTECH it is assumed to be smooth and periodic on [-1,1].
% Computing the chebyshev coefficients via innner products requires working with
% non-periodic, but smooth functions on [-1,1]. The right representation for f
% is then a Chebyshev expansion. We therefore convert f to a (happy) Chebyshev
% interpolant and return the coefficients. As an arbitrary choice we will
% convert f to a chebtech1 and then compute the resulting coefficients.
out = chebcoeffs(chebtech1(@(x) f.feval(x)), N, kind);

end


text:
%  Figure 10.60      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% script for disk read/write head control design
% the design is based on 1/s^2 with robustness
% with respect to a resonance of damping zeta = z at wm.
% the design calls for a single lead selected to maximize wc
% the phase margin is set by alpha = a
% the compensator has a second order roll-off filter
% the resonance is 'gain stabilized' with GM
% the time is scaled to milliseconds
wm= 5*pi
z=.05;
GM = 4;
a=0.1;
% input the nominal plant with a gain of 1.
numGo=1;
denGo= [1 0 0];
sysGo=tf(numGo,denGo);
% input the plant with the resonance
numG= [1/(50*pi) 1];
denG=[1/(25*pi^2) 1/(50*pi) 1  0 0];
sysG=tf(numG,denG);
% Define the design constant
A = ((a*2*z)/(GM))^.333;
% PM >50; we have selected alpha = a = 0.1
Tinv= A*wm*sqrt(a);
T=1/Tinv;
wc =  A*wm
% K =  wc^2;
K=  wc^2;
% Input the lead compensation
numD =K*sqrt(a)*[T 1];
denD = [a*T 1];
sysD=tf(numD,denD);
% Now design the roll-off filter
w1=wm*sqrt(A/sqrt(a));
% design the roll-off damping
z1=.3 ;
numF=1;
denF=[1/w1^2  2*z1/w1 1];
sysF=tf(numF,denF);
% now compute the entire compensation
sysDc=sysD*sysF;
sysOLo=sysD*sysGo;
sysOL=sysDc*sysG;
% bode(sysOLo);
% hold on
clf;
sysH=tf(1,1);
sysCL=feedback(sysOL,sysH);
t=0:0.1:6;
[y,t]=step(sysCL,t);
plot(t,y);
xlabel('Time (msec)');
ylabel('Amplitude');
grid on
title('Fig. 10.60 Step response of disk control;lead plus roll-off ')


text:
function plaza = create_plaza(B, L)
global plazalength;
topgap = 5;
bottomgap = 1;
plaza = zeros(plazalength,B+2);
plaza(1:plazalength,[1,2+B]) = -888;
if mod(B-L,2)==0
    for col = 2:B/2 - L/2 + 1
        for row = 1:(plazalength-1)/2 - topgap * (col-1)
            plaza(row,[col, B+3-col]) = -888;
        end
        for row = (plazalength+3)/2 + bottomgap*(col-1):plazalength
            plaza(row,[col, B+3-col]) = -888;
        end
    end
else
    plaza(1:plazalength, B+3) = -888;
    for col = 2:(B+1)/2 - L/2 + 1
        for row = 1:(plazalength-1)/2 - topgap * (col-1)
            plaza(row, [col, B+4-col]) = -888;
        end
        for row = (plazalength+3)/2 + bottomgap*(col-1):plazalength
            plaza(row, [col, B+4-col]) = -888;
        end
    end
end

author: Eurus-Holmes


text:
% Plot the per-pixel histogram intersection value when two distributions
% are compared.

% Zoya Bylinskii, April 2016
% linked to: "What do different evaluation metrics tell us about saliency models?"

% Note: if comparing multiple saliency maps, compute the resMap for each of
% the saliency maps using this function, but then normalize them jointly
% before plotting (to make range of values comparable across plots).

function [resMap,figtitle] = visualize_SIM(salMap, fixMap, toplot)
% salMap is the saliency model prediction (distribution)
% fixMap is the ground truth fixation map (distribution)
% if only plotting one map, set toplot=1; but if comparing multiple maps,
% see note above

if nargin < 3
    toplot = 1;
end

figtitle = 'Histogram intersection';

% format, resize, and normalize maps
map1 = im2double(fixMap);
map2 = im2double(imresize(salMap, size(fixMap)));
map1= (map1-min(map1(:)))/(max(map1(:))-min(map1(:)));
map2= (map2-min(map2(:)))/(max(map2(:))-min(map2(:)));

% compute per-pixel histogram interseciton
map1 = map1/sum(map1(:));
map2 = map2/sum(map2(:));
resMap = min(map1, map2);

% plot the visualization
if toplot
    figure;
    resMap_norm = resMap/max(resMap(:));
    imshow(resMap_norm);
    colormap('parula'); colorbar;
    title(figtitle,'fontsize',14);
end

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
%SELCLASSFixed mapping for selecting a single class from a dataset
%
%	[B,J] = SELCLASS(A,CLASS,NAME)
%	[B,J]  = A*SELCLASS([],CLASS,NAME)
%	[B,J]  = A*SELCLASS(CLASS,NAME)
%
% INPUT
%   A      Dataset
%   CLASS  Integer: Indices of desired classes in CLASSNAMES(A)
%          String array:  Class names
%          Cell array: Indices of desired classes in CLASSNAMES(A)
%          Default C = {}, i.e. return all classes separated out in a cell
%          array.
%   NAME   Integer: Index of desired labeling, see GETLABLISTNAMES
%          String:  Name of desired labeling, see GETLABLISTNAMES
%          Default: actual LABLIST
%
% OUTPUT
%   B   Desired classes of the dataset A. In case CLASS is a cell array, B
%       is a cell array of the desired classes. In case CLASS is empty, B
%       is a cell array of all classes.
%   J   Indices of returned objects in dataset A: B = A(J,:). In case CLASS
%       is a cell array, J is a cell array as well.
%
% DESCRIPTION
% B is a subset of the dataset A defined by the set of classes (CLASS). In
% case of a multi-labeling system (see MULTI_LABELING) the desired CLASS
% should refer to the label list NAME.
%
% In case A is soft labeled or is a target dataset by B = SELCLASS(A,CLASS)
% the entire dataset is returned, but the labels or targets are reduced to
% the selected class (target) CLASS.
%
% In case CLASS is a cell array the outputs are organised as cell arrays.
%
% EXAMPLES
% a = gendatm; b = selclass(a,[2 3 4]); % selects 3 classes
% a = gendatm; b = a*selclass;         % returns every class in a cell
% a = gendatb;
% a = addlabels(a,genlab(25*ones(4,1)),'4class'); % add second label list
% b = a*selclass('4class'); % returns 4 cells, preserves label list.
%
% SEE ALSO (PRTools Guide)
% DATASETS, SELDAT, CLASSNAMES, GENDAT, GETLABLIST, GETCLASSI, REMCLASS,
% GETLABLISTNAMES, MULTI_LABELING

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function [b,J] = selclass(varargin)

  argin = shiftargin(varargin,{'vector','cell','char'});
  argin = shiftargin(argin,'char',2);
  argin = setdefaults(argin,[],{},[]);

  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');

  else			% Evaluate

    [a,clas,lablist] = deal(argin{:});
    isa(a,'prdataset');
    if ~isempty(lablist)
      curn = curlablist(a);
      a = changelablist(a,lablist);
      [b,J] = feval(mfilename,a,clas);
      if iscell(b)
        for n=1:numel(b)
          b{n} = changelablist(b{n},curn);
        end
      else
        b = changelablist(b,curn);
      end
    else
      if iscell(clas)
        if  isempty(clas)
          for n=1:getsize(a,3);
            clas{n} = n;
          end
        end
        if isdataset(a)
          b = prdataset;
        else
          b = prdatafile;
        end
        J = cell(1,numel(clas));
        b = cell(1,numel(clas));
        for n=1:numel(clas)
          [b{n},J{n}] = seldat(a,clas{n});
        end
      else
        [b,J] = seldat(a,clas);
      end
    end

  end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function coded_bits = convolution_encoder(in_bits)
ConvCodeGenPoly=[1 0 1 1 0 1 1;1 1 1 1 0 0 1];
Nrow = size(ConvCodeGenPoly,1);
Nbits=size(ConvCodeGenPoly,2)+length(in_bits)-1;
uncoded_bits=zeros(Nrow,Nbits);
for row=1:Nrow
    uncoded_bits(row,1:Nbits)=rem(conv(in_bits,ConvCodeGenPoly(row,:)),2);
end
coded_bits=uncoded_bits;
end

author: LyricYang


text:
function [au,c] = assess(crit1,crit2,y,z,tt)
%   [AU,C] = ASSESS(CRIT1,CRIT2,Y,Z,TT)
%
%   Given two criteria vectors for two different models (CRIT1 and CRIT2),
%   observed time matrix Y, event indicator matrix z (= 0 if event is experienced before tt and  =1 if not)
%   and time vector tt returns Harrel's C and AUC for each value in TT time vector
%
%   Note: Z dimensions must be size(y,1) X size(tt,2)

% Copyright (C) 2012 Ernesto Ulloa, Aki Vehtari

% 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.

ip=inputParser;
ip.addRequired('crit1',@(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('crit2',@(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) isreal(x) && all(isfinite(x(:))))
ip.addRequired('z', @(x) isreal(x) && all(isfinite(x(:))))
ip.addRequired('tt', @(x) isreal(x) && all(isfinite(x(:))))
ip.parse(crit1,crit2,y,z,tt)

    for i=1:size(tt,2)
        c1(i)=hct(crit1(:,i),y(:,i),z(:,size(tt,2)),tt(i));
        c2(i)=hct(crit2(:,i),y(:,i),z(:,size(tt,2)),tt(i));
        au1(i)=auct(crit1(:,i),y(:,i),z(:,i),tt(i));
        au2(i)=auct(crit2(:,i),y(:,i),z(:,i),tt(i));
    end

   c(:,1)=c1;
   c(:,2)=c2;
   au(:,1)=au1;
   au(:,2)=au2;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = expcone(x)
%EXPCONE Defines a low-level exponential cone constraint x(2)*exp(x(1)/x(2)) <= x(3)
%
% Input
%    x       : Linear 3x1 SDPVAR object
%
% Example
%
% Standard  exponential cone constraint x(2)*exp(x(1)/x(2)) <= x(3)
%    F = expcone(x)
%
% To quickly define several cones, the argument can be a matrix, and the
% command is then short-hand for
% for i = 1:size(x,2);F = [F,expcone(x(:,i))];end
%
% See also  @SDPVAR/CONE, @SDPVAR/PCONE

[n,m] = size(x);
if min([n m])==1
    x = reshape(x,3,1);
end
[n,m] = size(x);
if n ~=3
    error('x must be a vector or matrix of height 3')
end
y = x;
if min([n m])>1
	y.typeflag = 22;
else
	y.typeflag = 21;
end
y = lmi(y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%This Matlab script can be used to reproduce Figure 1.18 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.01 (Last edited: 2019-04-17)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Define the SNR
SNR = 1;

%Define betabar (strength of inter-cell interference)
betabar = 1e-1;

%Define the range of number of UEs
K = 1:20;

%Define range of antenna-UE ratios
c = [1 2 4 8];

%Extract the maximal number of UEs and BS antennas
Kmax = max(K);
Mmax = Kmax*max(c);

%Select number of Monte Carlo realizations for the line-of-sight (LoS)
%angles and of the non-line-of-sight (NLoS) Rayleigh fading
numberOfRealizations = 10000;

%Generate NLoS channels using uncorrelated Rayleigh fading
H_NLoS_desired = sqrt(1/2)*(randn(Mmax,Kmax,numberOfRealizations)+1i*randn(Mmax,Kmax,numberOfRealizations));
H_NLoS_interfering = sqrt(betabar/2)*(randn(Mmax,Kmax,numberOfRealizations)+1i*randn(Mmax,Kmax,numberOfRealizations));

%Preallocate matrices for storing the simulation results
SE_MMSE_NLoS_montecarlo = zeros(length(K),length(c));
SE_MMSE_NLoS_nonlinear = zeros(length(K),length(c));

%% Go through all Monte Carlo realizations
for n = 1:numberOfRealizations

    %Output simulation progress
    disp([num2str(n) ' realizations out of ' num2str(numberOfRealizations)]);

    %Go through the range of number of UEs
    for kindex = 1:length(K)

        %Go through the range of antenna-UE ratios
        for cindex = 1:length(c)

            %Compute the number of antennas
            M = K(kindex)*c(cindex);

            %Compute the SE with non-linear processing under NLoS propagation
            %for one realization of the Rayleigh fading. We use the classic
            %log-det formula for the uplink sum SE, when treating the
            %inter-cell interference as colored noise
            SE_MMSE_NLoS_nonlinear(kindex,cindex) = SE_MMSE_NLoS_nonlinear(kindex,cindex) + real(log2(det(eye(M) + SNR* H_NLoS_desired(1:M,1:K(kindex),n)*H_NLoS_desired(1:M,1:K(kindex),n)' +  SNR*H_NLoS_interfering(1:M,1:K(kindex),n)*H_NLoS_interfering(1:M,1:K(kindex),n)' )) - log2(det(eye(M)+SNR*H_NLoS_interfering(1:M,1:K(kindex),n)*H_NLoS_interfering(1:M,1:K(kindex),n)')))/numberOfRealizations;

            %Compute the SE with M-MMSE under NLoS propagation for one
            %realization of the Rayleigh fading

            %Compute the M-MMSE combining vectors
            MMMSEfilter = ( SNR* H_NLoS_desired(1:M,1:K(kindex),n)*H_NLoS_desired(1:M,1:K(kindex),n)' + SNR* H_NLoS_interfering(1:M,1:K(kindex),n)*H_NLoS_interfering(1:M,1:K(kindex),n)' + eye(M) ) \ (SNR*H_NLoS_desired(1:M,1:K(kindex),n));

            %Compute the intra-cell channel powers after M-MMSE combining
            channelgainsIntracell = abs(MMMSEfilter'*H_NLoS_desired(1:M,1:K(kindex),n)).^2;

            %Extract the desired signal power for each UE
            signalpowers = diag(channelgainsIntracell);

            %Extract and compute interference powers for each UE
            interferencepowers = sum(channelgainsIntracell,2) - signalpowers + sum(abs(MMMSEfilter'*H_NLoS_interfering(1:M,1:K(kindex),n)).^2,2);

            %Compute the effective 1/SNR after noise amplification
            scalednoisepower = (1/SNR)*sum(abs(MMMSEfilter').^2,2);

            %Compute the uplink SE with M-MMSE combining
            SE_MMSE_NLoS_montecarlo(kindex,cindex) = SE_MMSE_NLoS_montecarlo(kindex,cindex) + sum(log2(1 + signalpowers./(interferencepowers+scalednoisepower)))/numberOfRealizations;

        end

    end

end

%% Plot the simulation results
figure(1);
hold on; box on;

plot(K,SE_MMSE_NLoS_montecarlo(:,4)./SE_MMSE_NLoS_nonlinear(:,4),'r-','LineWidth',1);
plot(K,SE_MMSE_NLoS_montecarlo(:,3)./SE_MMSE_NLoS_nonlinear(:,3),'k-.','LineWidth',1);
plot(K,SE_MMSE_NLoS_montecarlo(:,2)./SE_MMSE_NLoS_nonlinear(:,2),'b--','LineWidth',1);
plot(K,SE_MMSE_NLoS_montecarlo(:,1)./SE_MMSE_NLoS_nonlinear(:,1),'k:','LineWidth',1);

xlabel('Number of UEs (K)');
ylabel('Fraction of non-linear performance');

legend('M/K=8','M/K=4','M/K=2','M/K=1','Location','SouthWest');
ylim([0.5 1]);

author: emilbjornson

repo: massivemimobook


text:
function [dv, rp, orbitIn, orbitOut, deltaVVect, vInfDNorm, vInfDepart, vInfArrive, totDV, r1B, r2B, v1Output, v2Output, xferRp] =  multiFlybyObjFunc(x, numRevsArr,bodiesInfo,includeDepartVInf,includeArrivalVInf,celBodyData)
%multiFlybyObjFunc Summary of this function goes here
%   Detailed explanation goes here
    numX = size(x,1);

    if(numX==0)
        dv = 9E99;
        rp = [];
        orbitIn = [];
        orbitOut = [];
        deltaVVect = [];
        vInfDNorm = [];
        vInfDepart = [];
        vInfArrive = [];
        totDV = [];
        r1B = [];
        r2B = [];
        v1Output = [];
        v2Output = [];
        xferRp = [];
        return;
    end

    flybyBodies = bodiesInfo(2:end-1);
    numFB = length(flybyBodies);
    numREVS = length(bodiesInfo) - 1;

    tm = x(:,end-numREVS-numFB:end-numREVS);
    tm = round(tm);
    tm(tm==2) = -1;
    numTM = size(tm,2);

    numRevInds = x(:,end-numREVS+1:end);

    x = x(:,1:end-numTM-numREVS);
    daTimes = cumsum(x,2);

    gmu = zeros(length(bodiesInfo),numX);
    rVectsB = zeros(3,numX,size(x,2));
    vVectsB = zeros(3,numX,size(x,2));
    for(i=1:length(bodiesInfo)) %#ok<*NO4LP>
        bodyInfo = bodiesInfo{i};
        [rVect, vVect] = getStateAtTime(bodyInfo, daTimes(:,i)', getParentGM(bodyInfo, celBodyData));

        rVectsB(:,:,i) = rVect;
        vVectsB(:,:,i) = vVect;

        gmu(i,:) = bodyInfo.gm * ones(1,numX);

        if(i==1)
            vBDepart = vVect;
        elseif(i==length(bodiesInfo))
            vBArrive = vVect;
        end
    end
    gmu(1,:) = [];
    gmu(end,:) = [];
    gmu = reshape(gmu',1,numel(gmu'));

    rVectsB = reshape(rVectsB,3,length(bodiesInfo)*numX);
    r1B = rVectsB(:,1:end-size(x,1));
    r2B = rVectsB(:,size(x,1)+1:end);

    vVectsB = reshape(vVectsB,3,length(bodiesInfo)*numX);
    vB = vVectsB(:,size(x,1)+1:end-size(x,1));

    tempXferGmu = zeros(1,length(bodiesInfo)-1);
    for(i=1:length(bodiesInfo)-1)
        parentBodyInfo = bodiesInfo{i}.getParBodyInfo(celBodyData);
        tempXferGmuScalar = parentBodyInfo.gm;

%         tempXferGmuScalar = getParentGM(bodiesInfo{i}, celBodyData);
        tempXferGmu((i-1)*numX+1:i*numX) = tempXferGmuScalar;
    end
    xferGmuL = tempXferGmu;

%     xferGmuL = xferGmu * ones(1,(length(bodiesInfo)-1)*numX);

    dt = x(:,2:end) .* tm;
    dt = reshape(dt,1,numel(dt));

    popSize = size(r1B,2) / (length(bodiesInfo)-1);

    for(i=1:size(numRevInds,2))
        numRevInd = round(numRevInds(:,i));
        numRev((i-1)*popSize+1 : i*popSize) = numRevsArr{i}(numRevInd);
    end

    bool = numRev == 0;

    v1 = zeros(size(r1B));
    v2 = v1;

    xferRp = zeros(1,size(r1B,2));

    if(any(bool)) %use the faster vectorized solver
        [v1(:,bool), v2(:,bool)] = lambertBattinVector(r1B(:,bool), r2B(:,bool), dt(bool), numRev(bool), xferGmuL(bool));

        [xferSma, xferEcc, ~, ~, ~, ~] = vect_getKeplerFromState(r1B(:,bool),v1(:,bool),xferGmuL(bool), true);
        [~, xferRp(bool)] = computeApogeePerigee(xferSma, xferEcc);
    end
    if(any(~bool))                      %use the slower but multi-rev capable lambert loop
        bool = ~bool;
        inds = find(bool);

        for(i=1:length(inds))
            ind = inds(i);
            try
                [v1L, v2L, ~, exitflag] = orbit.lambert(r1B(:,ind)', r2B(:,ind)', dt(ind)/86400, numRev(ind), xferGmuL(ind));

                [xferSma, xferEcc, ~, ~, ~, ~] = getKeplerFromState(r1B(:,ind),v1L,xferGmuL(ind), true);
                [~, xferRp(ind)] = computeApogeePerigee(xferSma, xferEcc);

                if(any(isnan(v1L)) || any(isnan(v2L)) || exitflag == -1 || exitflag == -2)
                    v1L(isnan(v1L)) = 9E99;
                    v2L(isnan(v2L)) = 9E99;
                end
            catch ME %#ok<NASGU>
                v1L = 9E99 * ones(1,3);
                v2L = 9E99 * ones(1,3);
            end

            v1(:,ind) = v1L';
            v2(:,ind) = v2L';
        end
    end

    xferRp(isnan(xferRp)) = Inf;

    v1Output = v1;
    v2Output = v2;

    v1D = v1(:,1:numX);
    v1 = v1(:,numX+1:end);
    v2A = v2(:,end-numX+1:end);
    v2 = v2(:,1:end-numX);

    vInfIn = v2 - vB;
    vInfOut = v1 - vB;

    vInfDepart = v1D - vBDepart;
    vInfDNorm = sqrt(sum(abs(vInfDepart).^2,1));
    vInfDNorm = (vInfDNorm');

    vInfArrive = v2A - vBArrive;
    vInfANorm = sqrt(sum(abs(vInfArrive).^2,1));
    vInfANorm = (vInfANorm');

    [smaIn, eIn, hHat, sHat, smaOut, eOut, oHat] = computeMultiFlyByParameters(vInfIn, vInfOut, gmu);
    [hSMAIn, hEccIn, hIncIn, hRAANIn, hArgIn, ~, ~, rp] = computeHyperOrbitFromMultiFlybyParams(smaIn, eIn, hHat, sHat, true);
    [hSMAOut, hEccOut, hIncOut, hRAANOut, hArgOut, ~, ~, ~] = computeHyperOrbitFromMultiFlybyParams(smaOut, eOut, hHat, oHat, false);

%     if(any(abs(rp-rp2)>1E-6))
%         warning(num2str(abs(rp-rp2)));
%     end

    orbitIn = [hSMAIn', hEccIn', hIncIn', hRAANIn', hArgIn', zeros(size(hSMAIn')), ((1-eIn).*smaIn)'];
    orbitOut = [hSMAOut', hEccOut', hIncOut', hRAANOut', hArgOut', zeros(size(hSMAOut')), ((1-eOut).*smaOut)'];

    [~,vVectIn]=vect_getStatefromKepler(hSMAIn, hEccIn, hIncIn, hRAANIn, hArgIn, zeros(size(hArgIn)), gmu);
    [~,vVectOut]=vect_getStatefromKepler(hSMAOut, hEccOut, hIncOut, hRAANOut, hArgOut, zeros(size(hArgOut)), gmu);
    deltaVVect = vVectOut-vVectIn;
    deltaV = sqrt(sum(abs(deltaVVect).^2,1));
    totDV = sum(reshape(deltaV,numX,length(flybyBodies)),2);

    dv = totDV;

    if(includeDepartVInf)
        dv = dv + vInfDNorm;
    end

    if(includeArrivalVInf)
        dv = dv + vInfANorm;
    end

%     dv = vInfDNorm + totDV + vInfANorm;

    dv(dv>=1E99) = NaN;
    dv(isnan(dv)) = max(dv)+10;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function status = CreateDirectories(rootDir)
% status = CreateDirectories(rootDir);

% Set up the directory structure for reconning and subsequent use
% by mrLoadRet. Return TRUE if all directories are successfully
% created or already present.
%
% DBR, 4/99

status = 0;
% Create the root directory if it doesn't already exist:
if ~MakeDir(rootDir)
    disp(['Could not create home directory ',rootDir]);
    return
end

% Create the other directory structures if they don't already exist:
if ~MakeDir('Inplane', rootDir); return; end
if ~MakeDir('Original', fullfile(rootDir,'Inplane')); return; end
if ~MakeDir('TSeries', fullfile(rootDir,'Inplane','Original')); return; end

status = 1;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = r4_erf ( x )

%*****************************************************************************80
%
%% R4_ERF evaluates the error function of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the error function of X.
%
  persistent erfcs
  persistent nterf
  persistent sqeps
  persistent sqrtpi
  persistent xbig

  sqrtpi = 1.7724538509055160;

  if ( isempty ( nterf ) )

    erfcs = [ ...
      -0.049046121234691808, ...
      -0.14226120510371364, ...
       0.010035582187599796, ...
      -0.000576876469976748, ...
       0.000027419931252196, ...
      -0.000001104317550734, ...
       0.000000038488755420, ...
      -0.000000001180858253, ...
       0.000000000032334215, ...
      -0.000000000000799101, ...
       0.000000000000017990, ...
      -0.000000000000000371, ...
       0.000000000000000007 ]';

    nterf = r4_inits ( erfcs, 13, 0.1 * r4_mach ( 3 ) );
    xbig = sqrt ( - log ( sqrtpi * r4_mach ( 3 ) ) );
    sqeps = sqrt ( 2.0 * r4_mach ( 3 ) );

  end

  y = abs ( x );

  if ( y <= sqeps )
    value = 2.0 * x / sqrtpi;
  elseif ( y <= 1.0 )
    value = x * ( 1.0 + r4_csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) );
  elseif ( y <= xbig )
    value = 1.0 - r4_erfc ( y );
    if ( x < 0.0 )
      value = - value;
    end
  else
    value = 1.0;
    if ( x < 0.0 )
      value = - value;
    end
  end

  return
end


text:
%  INTERNAL FUNCTION: wishart distribution
%
%  Note:
%     One would rarely need to call distribution functions directly.
%     Refer to <densities.html> file for list of supported distributions.
%
%

author: jmaih


text:
function hh = logbar(edges,vals,width,colors,hax)

if nargin < 3,
  width = .8;
end

edges = edges(:)';
nbins = numel(edges)-1;
sz = size(vals);
if ~any(sz==nbins),
  error('vals should be ndata x nbins');
end

if ~ismatrix(vals),
  error('vals must be 1- or 2-dimensional');
end
if sz(2) ~= nbins && sz(1) == nbins,
  vals = vals';
end
ndata = size(vals,1);

if nargin < 5,
  hax = gca;
end
hfig = get(hax,'Parent');

if nargin < 4 || isempty(colors),
  colors = get(hfig,'Colormap');
  if size(colors,1) > ndata,
    colors = colors(round(linspace(1,size(colors,1),ndata)),:);
  elseif size(colors,1) < ndata,
    colors = colors(modrange(1:ndata,1,size(colors,1)),:);
  end
elseif size(colors,1) < ndata,
  colors = colors(modrange(1:ndata,1,size(colors,1)),:);
end

dedge = diff(edges)*width;
off0 = (1-width)/2;
hh = nan(1,ndata);
holdstate = ishold;
for i = 1:ndata,
  off1 = off0 + (i-1)/ndata;
  off2 = off0 + i/ndata;
  x1 = edges(1:end-1)+off1*dedge;
  x2 = edges(1:end-1)+off2*dedge;
  x = [x1;x1;x2;x2];
  y = [zeros(1,nbins);vals([i,i],:);zeros(1,nbins)];
  hh(i) = patch(x(:),y(:),colors(i,:),'Parent',hax);
  if i == 1,
    hold(hax,'on');
  end
end
if holdstate == 0,
  hold(hax,'off');
end

set(hax,'XScale','log');

author: kristinbranson

repo: JAABA


text:
function [pool,logpaths] = par_getworkers_ssh(varargin)
% Acquire workers on some remote machines and return hostnames and ports of those that are available.
% Pool = par_getworkers_ssh(Hostnames,ProcessorsPerNode,MatlabThreads,StartupCommand,ShareResources,SharingPort,LoggingPath,IdentityFile,BinaryWorker,BinaryNmae,MatlabCompilerRoot,MinFreeMemory,MaxCPULoad,LoadWindow,AvoidProcess,RecruitMachines,VerboseOutput)
%
% This function attempts to ensure that the desired number of worker processes on remote machines
% is available for participating in distributed computations. Such computations are usually being
% submitted in the form of tasks to the respective worker processes, via the function par_schedule.
%
% Alternatively, the workers can also be started manually on the remote machines (this is the
% recommended way if Windows machines are involved or use of SSH is not permitted) and their names
% would then just be listed in the cluster settings (rather than to configure acquisition behavior).
%
% In:
%   --- what shall be acquired ---
%
%   Hostnames : cell array of hostnames on which to acquire workers (duplicates and port assignments
%               will be ignored).
%
%   ProcessorsPerNode : number of processes to acquire per node; 0 means as many as there are cores
%                       on each node, divided by the value of the 'matlabthreads' parameter
%                       (default: 0)
%
%   MatlabThreads : assumed number of threads per MATLAB instance (default: 4)
%
%   --- how shall it be acquired ---
%
%   StartupCommand : startup command, if any; if binary_worker is true, this is a sequence of shell
%                    commands, otherwise it is a sequence of MATLAB statements which contain the
%                    placeholder %d where the TCP port # would be inserted; If the startup_command
%                    is passed as 'autogenerate', a bcilab-specific startup command will be
%                    generated (default: 'autogenerate')
%
%   ShareResources : whether worker that are already running will be returned as part of the worker
%                    set, and are thus potentially shared with other users. If false, the requested
%                    number of workers will always be started and be private to the requester -
%                    however, potentially over-subscribing the node. (default: true)
%
%   SharingPort : the starting port for shared workers; if 0, this is a hash of the startup_command
%                 (default: 0)
%
%   --- associated files ---
%
%   LoggingPath : default logfile path; make sure that this does not conflict with other users' log
%                 paths. (default: '')
%
%   IdentityFile : SSH identity file, if needed (implying the -i option for ssh) (default: '')
%
%   --- optional binary worker support ---
%
%   BinaryWorker : whether to use assume a binary worker implementation rather than the MATLAB
%                  worker (default: false)
%
%   BinaryName : name of the worker binary (default: 'build')
%
%   MatlabCompilerRoot : installation path of the MATLAB compiler runtime; if empty, a few locations
%                        will be searched for an MCR that corresponds to the MATLAB version that is
%                        running this function. If nothing is found, it is assumed that the MCR has
%                        been installed into the system path  and is found automatically. (default: '')
%
%   --- exclusion criteria ---
%
%   MinFreeMemory : minimum memory (in bytes) that needs to be available for an instance to be started
%                   (default: 4GB)
%
%   MaxCPULoad : maximum acceptable average CPU load (summed over all cores) for a node to be
%                considered (default: Inf)
%
%   LoadWindow : CPU load measurement window, in minutes; can be 1, 5 or 15 (default: 15)
%
%   AvoidProcess : avoid nodes with the given process running (default: 'avoidme')
%
%   ShutdownTimeout : if non-zero, the worker will shut itself down if it has not received a
%                     heartbeat signal from a client in the given time, in seconds. (default: 0)
%
%   RecruitMachines : whether to recruit other machines, rather than just list them (default: true)
%
%
%   --- misc ---
%
%   HarvestTimeout : Harvesting timeout. This function will attempt to harvest process id's from the
%                    workers until this timeout has expired (in which case the pid of straggers will
%                    not be known). Note that when multiple MATLABs try to launch simultaneously, it
%                    can take longer than normal. (default: 300)
%
%   VerboseOutput : whether to show verbose outputs during acquisition (default: true)
%
%
% Out:
%   Pool : cell array of 'hostname:port' strings specifying the list of available machines
%          Note: in the case that no return value is requested, the global variable
%                tracking.parallel.pool will receive this result. This is the recommended way to
%                use par_getworkers_ssh, as par_schedule uses this pool by default.
%
%  Logpaths : cell array of file paths of the logfiles corresponding to the workers in pool
%
%
% Notes:
%   If workers are shared between users, code which maintains global state (in particular
%   user-specific configuration options) may show unexpected behavior, as these variables are not
%   duplicated across users.
%
%   This function works only across POSIX-compliant (Unix/Linux/Mac) systems. On standard Windows
%   workstations, the workers must be started manually (unless POSIX emulation software is
%   installed).
%
% Examples:
%   % acquire workers on the three hosts (assuming bcilab is in the same file-system location remotely as here)
%   par_getworkers_ssh({'machine1','machine2','machine3'})
%
%   % try to start fast pre-compiled workers (assuming bcilab is in the same file-system location remotely as here)
%   par_getworkers_ssh('hostnames', {'machine1','machine2','machine3'}, 'binary_worker',true)
%
%   % MATLAB version again, also setting a custom logging path
%   par_getworkers_ssh('hostnames', {'machine1','machine2','machine3'}, 'logging_path','/home/christian/workerlogs/')
%
%   % like before, but this time only consider machines that have at least 8GB of free RAM and at most 10% CPU load
%   par_getworkers_ssh('hostnames', {'machine1','machine2','machine3'}, 'logging_path','/home/christian/workerlogs/', 'min_memory',2^33, 'max_cpuload',0.1)
%
%   % like before, but measuring the CPU load within the last minute (instead to 15 minutes)
%   par_getworkers_ssh('hostnames', {'machine1','machine2','machine3'}, 'logging_path','/home/christian/workerlogs/', 'min_memory',2^33, 'max_cpuload',0.1, 'load_window',1)
%
%   % acquire workers on the three hosts, this time using a custom startup command (the %d is where the port goes)
%   par_getworkers_ssh('hostnames', {'machine1','machine2','machine3'},'startup_command','cd path/to/bcilab/; bcilab(''myconfig'',''worker'',{%d,1}')
%
%   % try to start fast pre-compiled workers (instead of entire MATLAB sessions); note that the
%   % startup command is this time a UNIX shell command which CD's into the directory where the run_build.sh is located
%   par_getworkers_ssh('hostnames', {'machine1','machine2','machine3'}, 'binary_worker',true, 'startup_command','cd path/to/bcilab/build/distrib')
%
% See also:
%   par_worker
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-02-15

% Copyright (C) Christian Kothe, SCCN, 2011, christian@sccn.ucsd.edu
%
% 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

global tracking;

opts = arg_define(varargin,...
    ... % resources to acquire
    arg({'hostnames','Hostnames','Hosts'},{''},[],'Host names to acquire. Host names or IP addresses of machines on which to consider starting workers. If this is left empty, the current global worker pool will be used.','type','expression','shape','row'), ...
    arg({'processors_per_node','WorkersPerNode','ProcessorsPerNode'},4,uint32([1 1024]),'Workers per node. A value of 0 means that as many instances as there are cores on each node shall be used, but divided by the value of the "matlabthreads" (Number of MATLAB threads) parameter.'), ...
    arg({'matlabthreads','MatlabThreads','matlab_threads'},4,uint32([1 32]),'Number of MATLAB threads. This is the number of threads that each worker uses internally.'), ...
    ... % SSH command
    arg({'identity_file','IdentityFile'},'', [], 'SSH identity file. Optional, if needed for passwordless login. This corresponds to the -i option in ssh.'), ...
    ... % startup commands within MATLAB
    arg({'matlab_command','MATLABCommand'},'autogenerate', [], 'Command to start MATLAB. This is the command to run on the command line to start the desired version of MATLAB. If this is set to ''autogenerate'' the same path that runs the par_getworkers_ssh command will be used (useful on clusters with identical file systems).'), ...
    arg({'startup_prefix','StartupPrefix'},'',[], 'Startup prefix. Any startup lines to run before running the main BCILAB command.'),...
    arg({'startup_command','StartupCommand'},'autogenerate',[], 'Startup command. If the binary_worker parameter (Use compiled workers) is true, this is a sequence of shell commands, otherwise it is a sequence of MATLAB statements which contain the placeholder %d where the TCP port # would be inserted; If set to ''autogenerate'', a startup command specific to the current BCILAB environment will be generated.'),...
    ... % binary commands for compiled workers
    arg({'binary_worker','BinaryWorker'},false,[], 'Use compiled workers. Whether to start a given binary worker implementation rather than the MATLAB worker.'), ...
    arg({'mcr_root','MatlabCompilerRoot'},'', [], 'MATLAB Compiler directory. Installation path of the MATLAB compiler runtime; if empty, a few locations will be searched for an MCR that corresponds to the MATLAB version that is running this function (this assumes that the compiler is installed in the same directory remotely as it is locally). If nothing is found, it is assumed that the MCR has been installed into the system path and is found automatically.'), ...
    arg({'binary_name','BinaryName'},'build', [], 'Name of the binary. The name of the BCILAB binary (if using compiled workers).'), ...
    ... % port allocation
    arg({'share_resources','ShareResources'},true,[],'Share resources. Whether workers that are already running (started by other users) will be returned as part of the worker set, and are thus potentially shared with other users. If false, the requested number of workers will always be started and be private to the requester -- however, potentially over-subscribing the capacity of the node.'),...
    arg({'sharing_port','SharingPort'},0, uint32([0 65535]), 'Port for sharing. The lowest port used by shared workers -- if 0, this is computed based on the startup_command (recommended).'), ...
    ... % worker status and logging
    arg({'logging_path','LoggingPath'},env_translatepath('home:/.bcilab/logs/workers'), [], 'Logging path. Make sure that this does not conflict with other users'' log paths.'), ...
    arg({'verbose_output','VerboseOutput'},true, [], 'Verbose output. Whether to display information about the acquisition process.'), ...
    arg({'processors_command','NumProcessorsCommand'},'cat /proc/cpuinfo | grep ^processor | wc -l;', [], 'Processor-check command. This command should return the number of processors on a given worker machine.','guru',true), ...
    arg({'cpuload_command','CPULoadCommand'},'cat /proc/loadavg;', [], 'CPU-load command. This command should generate CPU load statistics.','guru',true), ...
    arg({'memfree_command','MemCheckCommand'},'cat /proc/meminfo | grep ^MemFree;', [], 'Avoidance-check command. This command should return the amount of free RAM on the given worker machine.','guru',true), ...
    arg({'avoidcheck_command','AvoidanceCheckCommand'},'ps -A | grep ''%s'' | wc -l;', [], 'Avoidance-check command. This command should return the number of processes on a given client to avoid (based on a regex pattern).','guru',true), ...
    ... % sanity checks
    arg({'sanity_checks','SanityChecks'},true, [], 'Perform sanity checks. This checks whether the necessary paths are present on the worker machines, etc.'), ...
    arg({'path_exists_command','PathExistsCommand'},'ls %s;', [], 'Path-check command. This command should return success if the path exists, and otherwise yield an error return value.','guru',true), ...
    arg({'binary_exists_command','BinaryExistsCommand'},'which %s;', [], 'Binary-check command. This command should return success if the binary exists, and otherwise yield an error return value.','guru',true), ...
    ... % exclusion criteria
    arg({'min_memory','MinFreeMemory','MinMemory'},2^32,[0 Inf], 'Minimum free memory. Minimum amount of RAM memory (in bytes) that needs to be available for an instance to be started.'),...
    arg({'max_cpuload','MaxCPULoad'},Inf, [0 Inf], 'Maximum CPU load. The maximum acceptable average CPU load (summed over all cores) for a node to be considered.'), ...
    arg({'load_window','LoadWindow'},15,[0 Inf],'Load estimation window. The time window in seconds (into the past) over which the CPU load will be estimated to determine whether a worker shall be started on a given node.'), ...
    arg({'avoid_proc','AvoidProcess'},'avoidme',[],'Process to avoid. Optionally a process name that indicates that no workers shall be spawned if the process is running on some machine (e.g. a time-critical computation).'),...
    arg({'avoid_hosts','AvoidHosts'},{},[],'Hosts to avoid. Optionally a list of hostnames to avoid.'), ...
    ... % shutdown management
    arg({'shutdown_timeout','ShutdownTimeout'},0,[0 Inf],'Shutdown timeout. If non-zero, the worker will be shut down if it has not received a heartbeat signal from a client in the given time frame, in seconds. For this to work, the function env_acquire_cluster should be used as it sets up the heartbeat timer.'),...
    ... % misc
    arg({'harvest_timeout','HarvestTimeout'},300,[],'Harvesting timeout. This function will attempt to harvest process id''s from the workers until this timeout has expired (in which case the pid of straggers will not be known). Note that when multiple MATLABs try to launch simultaneously, it can take longer than normal.'), ...
    arg({'recruit_machines','RecruitMachines'},true,[],'Recruit workers. Whether to actually recruit other machines, rather than just list them.'));

% no arguments were passed? bring up GUI dialog
if isempty(varargin)
    opts = arg_guidialog;
    if isempty(opts)
        return; end % -> user clicked cancel
end

arg_toworkspace(opts);

if ispc
    disp('Note: Acquiring workers from BCILAB will only work if your operating system provides an ''ssh'' command. From a Windows machine, you likely have to have some program installed to use this feature. Also note that your cluster needs to run Linux for this to work.');
    disp('You can always fall back to starting your worker processes (either as MATLAB instances or compiled binaries) by hand and listing their host:port''s in the worker pool (see Cluster Settings GUI).');
end

if isempty(mcr_root)
    mcr_root = ''; end
if isempty(identity_file)
    identity_file = ''; end
if isempty(hostnames) || isequal(hostnames,{''})
    hostnames = par_globalsetting('pool'); end

% sanity checks
if isempty(hostnames)
    disp('The list of hostnames to connect to is empty; exiting.');
    return;
end
if ~iscellstr(hostnames)
    error('The hostnames must be given as a cell array of strings.'); end
if ~isscalar(processors_per_node) || ~isnumeric(processors_per_node)
    error('The processors_per_node parameter must be given as a number.'); end
if ~islogical(share_resources)
    error('The share_resources parameter must be either true or false.'); end
if ~isscalar(sharing_port) || ~isnumeric(sharing_port)
    error('The sharing_port parameter must be given as a number.'); end
if ~ischar(logging_path)
    error('The logging_path parameter must be given as a string.'); end
if ~isempty(startup_prefix) && ~ischar(startup_prefix)
    error('The startup_prefix parameter must be given as a string.'); end
if ~ischar(startup_command)
    error('The startup_command parameter must be given as a string.'); end
if ~islogical(binary_worker)
    error('The binary_worker flag must be either true or false.'); end
if strcmp(matlab_command,'autogenerate')
    matlab_command = [matlabroot filesep 'bin' filesep 'matlab']; end
if ~ischar(mcr_root)
    error('The mcr_root parameter must be given as a string.'); end
if ~ischar(identity_file)
    error('The identity_file parameter must be given as a string.'); end
if ~isscalar(max_cpuload) || ~isnumeric(max_cpuload)
    error('The max_cpuload parameter must be given as a number.'); end
if ~isscalar(min_memory) || ~isnumeric(min_memory)
    error('The min_memory parameter must be given as a number.'); end
if ~isscalar(load_window) || ~isnumeric(load_window)
    error('The load_window parameter must be given as a number.'); end
if ~isscalar(shutdown_timeout) || ~isnumeric(shutdown_timeout)
    error('The shutdown_timeout parameter must be given as a number.'); end
window_remap = [1 1 1 2 2 2 2 2 2 2 3 3 3 3 3];
load_window = window_remap(min(load_window,15));
if ~isempty(identity_file)
    identity_file = env_translatepath(identity_file);
    if ~exist(identity_file,'file')
        disp_once('The given SSH identitity file does not appear to exist: %s',identity_file); end
    identity_file = ['-i ' identity_file];
end
if ~isempty(avoidcheck_command)
    avoidcheck_command = sprintf(avoidcheck_command,avoid_proc); end

% give some recommendations regarding good input ranges
if processors_per_node > 64
    disp_once('The WorkersPerNode setting is unusually large -- make sure that you do not use more workers than you have cores.'); end
if matlabthreads > 8
    disp_once('The MatlabThreads setting is relatively large; typically MATLAB wastes increasingly many CPU cycles with more than 8-16 threads (except for some very streamlined operations).'); end
if ~exist(logging_path,'dir')
    disp_once('The logging path does not exist on the local file system; make sure that it is accessible to remote workers.'); end
if ~isempty(mcr_root) && binary_worker && ~exist(mcr_root,'dir')
    disp_once('The matlab compiler root does not exist on the local file system; make sure that it is accessible to remote workers.'); end
if min_memory < (2^10 * 500)
    disp_once('The MinFreeMemory setting is very small; MATLAB alone often needs more than 500MB free memory to run.'); end
if load_window > 15
    disp_once('The LoadWindow setting is larger than what is typically tracked on Linux systems (15 minutes).'); end
if shutdown_timeout == 0
    disp_once('It is recommended that a nonzero ShutdownTimeout value is used (otherwise workers can linger on forever).'); end
if shutdown_timeout ~= 0 && harvest_timeout+30 > shutdown_timeout
    % (otherwise workers can shut themselves down before the heartbeat is being initialized
    disp_once('The HarvestTimeout is larger than the ShutdownTimeout; increasing ShutdownTimeout.');
    shutdown_timeout = harvest_timeout+30;
end

% generate in the bcilab-specific startup command
if strcmp(startup_command,'autogenerate')
    if binary_worker
        startup_command = sprintf('cd %s',env_translatepath('bcilab:/build/distrib'));
    else
        startup_command = sprintf('%s; cd %s; bcilab %s worker {%s,1,%d} parallel {}',startup_prefix, env_translatepath('bcilab:/'), tracking.configscript,'%d',shutdown_timeout);
    end
end

% use a port that depends on the startup command for sharing resources
if sharing_port == 0
    sharing_port = 10000 + mod(23457+hlp_fingerprint([matlab_command startup_command]),50000); end

% remove any port assignments
for i=1:length(hostnames)
    colons = hostnames{i}==':';
    if any(colons)
        hostnames{i} = hostnames{i}(1:find(colons)-1); end
end
% replace localhost by hlp_hostname
for i=1:length(hostnames)
    if strcmp(hostnames{i},'localhost')
        hostnames{i} = hlp_hostname; end
end
% remove duplicates
hostnames = unique(hostnames);

% remove any host to avoid
if ~isempty(avoid_hosts)
    for h=1:length(avoid_hosts)
        hostnames = hostnames(~strncmp(hostnames,avoid_hosts{h},length(avoid_hosts{h}))); end
end

% filter hostnames by machine availability, and retrieve machine stats
disp('Listing compute servers ...');
stats = [];
for host = hostnames(:)'
    % collect system info
    stats_commands = [processors_command cpuload_command avoidcheck_command memfree_command];
    ssh_command = sprintf('ssh %s -x %s "%s"',identity_file,host{1},stats_commands);
    [errcode,info] = system(ssh_command);
    if ~errcode
        try
            % parse info
            lines = hlp_split(info,10);
            loadavg = hlp_split(lines{2},' ');
            meminfo = hlp_split(lines{4},': ');
            % store stats
            stats(end+1).hostname = host{1};
            try
                stats(end).processors = str2num(lines{1});
            catch
                fprintf('The processor-check command (%s) returns an unexpected format on host %s.\n',processors_command,host{1});
                stats(end).processors = 4;
            end
            try
                stats(end).cpuload = str2num(loadavg{load_window});
            catch
                fprintf('The cpuload command (%s) returns an unexpected format on host %s.\n',cpuload_command,host{1});
                stats(end).cpuload = 0;
            end
            try
                stats(end).avoidance = str2num(lines{3});
                if ~isscalar(stats(end).avoidance)
                    stats(end).avoidance = 0; end
            catch
                fprintf('The avoid-processes command (%s) returns an unexpected format on host %s.\n',avoidcheck_command,host{1});
                stats(end).avoidance = 0;
            end
            try
                stats(end).freemem = 1024*str2num(meminfo{2});
            catch
                fprintf('The memory-check command (%s) returns an unexpected format on host %s.\n',memfree_command,host{1});
                stats(end).freemem = Inf;
            end
        catch
            fprintf('The stats command (%s) gave an error on host %s.\n',ssh_command,host{1});
        end
    end
end

% start workers if necessary
pool = {};
logpaths = {};
batchid = ['ssh-' hlp_hostname '-' strrep(strrep(datestr(now),':','.'),' ','_')];
if ~isempty(stats)
    fprintf('%.0f%% of requested hosts are available.\n',100*length(stats)/length(hostnames));

    % reorder nodes by increasing CPU load
    [dummy,idx] = sort([stats.cpuload]); %#ok<ASGLU>
    stats = stats(idx);
    hostnames = {stats.hostname};

    % start the workers...
    for k = 1:length(hostnames)
        host = hostnames{k};
        % check CPU load criterion
        if stats(k).cpuload > max_cpuload
            continue; end
        % check avoid criterion
        if stats(k).avoidance
            continue; end
        % determine number of ports (=procs) to consider
        if processors_per_node == 0
            % as many as there are cores
            rangelen = ceil(stats(k).processors/matlabthreads);
        else
            % at most 4x as many as (#cores/#matlabthreads)
            if ~isempty(stats(k).processors) && stats(k).processors > 0
                rangelen = min(ceil(stats(k).processors/matlabthreads)*4,processors_per_node);
            else
                rangelen = processors_per_node;
            end
        end
        % determine starting port range to consider
        if share_resources
            % the worker process may be shared between users
            rangestart = sharing_port;
        else
            while true
                % chose a random start port that does not overlap with the share port in a range of
                % 512 successive ports
                rangestart = 10000+50000*rand();
                if isempty(intersect(rangestart:rangestart+512-1,sharing_port:sharing_port+512-1))
                    break; end
            end
        end

        % try to start them...
        for port=rangestart:rangestart+rangelen-1
            if verbose_output
                fprintf('  acquiring worker at %s on port %d...\n', host, port); end
            % construct command
            if isempty(logging_path)
                logportion = '';
            else
                try
                    % ensure that the logging path exists
                    io_mkdirs([logging_path filesep batchid filesep],{'+w','a'});
                catch
                    disp_once(['Warning: the logging path "' logging_path '" does not exist locally and could not be created.\nPlease make sure that it exists on all worker machines, as otherwise the computation will not start.']);
                end
                logpath = sprintf('%s/%s/%s_%d.out',logging_path,batchid,host,port);
                if sanity_checks
                    % check whether the logging path is writable on the worker
                    check_command = sprintf('ssh %s -x %s "touch %s"',identity_file,host,logpath);
                    [errcode,info] = system(check_command);
                    if errcode
                        fprintf('The logging path (%s) does not appear to be writable on host %s.\n',logpath,host); end
                end
                logportion = [' > ' logpath];
            end
            if binary_worker
                % search for MCR root path
                if isempty(mcr_root)
                    [major,minor] = mcrversion;
                    mcr_tail = sprintf([filesep 'MATLAB_Compiler_Runtime' filesep 'v%d%d'],major,minor);
                    possible_locations = {'/opt/MATLAB','/usr/local/MATLAB','/usr/common/MATLAB',hlp_homedir,[hlp_homedir filesep 'MATLAB'],'C:\\Program Files\\MATLAB','C:\\Program Files (x86)\\MATLAB','/Applications/MATLAB/'};
                    for p=1:length(possible_locations)
                        if exist([possible_locations{p} mcr_tail],'dir')
                            mcr_root = [possible_locations{p} mcr_tail]; end
                    end
                    if isempty(mcr_root)
                        fprintf('MATLAB compiler runtime v%d%d not found; if you get an error subsequently, please make sure that it is installed in a recognized location (or pass the mcr_root as an acquire options).\n',major,minor); end
                end
                % determine the file to run
                if isempty(mcr_root)
                    run_file = binary_name;
                else
                    run_file = ['run_' binary_name '.sh'];
                end
                % perform sanity checks
                if sanity_checks
                    % check whether the startup command works
                    check_command = sprintf(['ssh %s -x %s "' startup_command '; ' path_exists_command '"'],identity_file,host,run_file);
                    [errcode,info] = system(check_command);
                    if errcode
                        fprintf('The worker binary path (%s) or file (%s) appear to be missing on host %s.\n',startup_command,run_file,host); end
                end
                if isempty(mcr_root)
                    % assume that the binary finds all necessary MCR files
                    command = ['ssh ' identity_file ' -x ' host '  "' startup_command '; ./' run_file ' worker ''{' num2str(port) ',1,' num2str(shutdown_timeout) '}'' parallel {} ' logportion '" &'];
                else
                    % use the run script to point to the correct MCR install location
                    command = ['ssh ' identity_file ' -x ' host ' "' startup_command '; ./' run_file mcr_root ' worker ''{' num2str(port) ',1,' num2str(shutdown_timeout) '}'' parallel {} ' logportion '" &'];
                end
            else
                % perform sanity checks
                if sanity_checks
                    % check presence of MATLAB binary
                    if any(matlab_command==filesep)
                        test_command = path_exists_command;
                    else
                        test_command = binary_exists_command;
                    end
                    check_command = sprintf(['ssh %s -x %s "' test_command '"'],identity_file,host,matlab_command);
                    [errcode,info] = system(check_command);
                    if errcode
                        fprintf('The MATLAB binary (%s) appears to be missing on host %s.\n',matlab_command,host); end
                end
                % run the MATLAB startup command (and substitute the port into it)
                command = sprintf('ssh %s -x %s "%s -nodisplay -r ''%s''%s" &', ...
                    identity_file, host, matlab_command, sprintf(startup_command,port),logportion);
            end
            % issue startup command, if necessary
            if recruit_machines && port_free(host,port)
                % port free: check for memory criterion
                if stats(k).freemem > min_memory
                    % start, and add to pool
                    system(command);
                    logpaths{end+1} = logpath;
                    stats(k).freemem = stats(k).freemem - min_memory;
                    pool{end+1} = sprintf('%s:%d',host,port);
                else
                    fprintf('Host %s does not have the required amount of free memory.\n',host);
                end
            else
                % worker is already running: add to pool
                pool{end+1} = sprintf('%s:%d',host,port);
            end
        end
    end
    fprintf('%i processor slots acquired.\n',length(pool));

    fprintf('\nWaiting for workers to start up to get process IDs...\n');
    harvested_addresses = par_parse_logfiles(logpaths, harvest_timeout);
    % now merge the pool into the harvested addresses (PIDs are collected opportunistically)
    num_matches = 0;
    for p=1:length(harvested_addresses)
        addr = harvested_addresses{p};
        pos = addr=='@';
        if any(pos)
            hostport = addr(find(pos,1)+1:end);
            match = strcmp(pool,hostport);
            if any(match)
                pool{find(match,1)} = addr;
                num_matches = num_matches+1;
            end %#ok<AGROW>
        end
    end
    if num_matches == 0 && length(harvested_addresses) == length(pool)
        % we couldn't match the reported hostnames with those used to schedule the workers, but we
        % got the same number of results as we scheduled, so we just replace the pool by
        % harvested_addresses
        pool = harvested_addresses;
        fprintf('Note: the format of the hostnames used to schedule workers does not match the format of hostnames reported by the workers themselves. Please consider changing it.\n');
    end
    if num_matches < length(pool)
        % if you have workers without PID this means that they will not be immediately released
        % when you call env_release_cluster or par_clearworkers; instead, they will terminate
        % themselves once they are finished with their current computation and and their shutdown
        % timeout has expired (note that if your shutdown timeout is 0, they will linger until you
        % kill them)
        fprintf('Launched %i workers (%i without PID): %s\n',length(pool),length(pool)-num_matches,hlp_tostring(pool));
        if shutdown_timeout == 0
            fprintf('WARNING: the workers without PID have no shutdown timeout set, so you need to kill them manually.'); end
    else
        fprintf('Launched %i workers: %s\n',length(pool),hlp_tostring(pool));
    end

else
    disp('None of the listed machines is available; please make sure that:');
    disp(' * you can ssh into these machines without a password prompt"');
    disp(' * the files /proc/cpuinfo, /proc/meminfo, and /proc/loadavg are available on them"');
    disp(' * the commands ps, cat, grep, and wc are working on them"');
end

if nargout == 0
    % assign to global variable, if necessary...
    par_globalsetting('pool',pool(:)');
    par_globalsetting('logfiles',logpaths(:)');
end

% check whether the given port on the specified host is free
function result = port_free(host,port)
try
    java.net.Socket(host,port);
    result = false;
catch
    result = true;
end

repo: REST


text:
function [] = mainExp(class,trainId,paramStruct)
if(nargin<4)
    paramStruct = {};
end

mainTrain(class,trainId,paramStruct);

% Test with keypoints
mainTest(class,trainId,'withKps',{'opt.rotationPredNum','1'...
'opt.relaxInit','{}','opt.relaxOpt','{''scale'',''translation'',''rotation'',''kps''}'});

%mainTest(class,trainId,strcat(trainId,'withKpsAllTest'),{'opt.rotationPredNum','1',...
%'opt.relaxInit','{}','opt.relaxOpt','{''scale'',''translation'',''rotation'',''kps''}'});

%  mainTest(class,trainId,strcat(trainId,'withKpsSDS'),{'opt.rotationPredNum','1',...
%  'opt.relaxInit','{''mask''}','opt.relaxOpt','{''scale'',''translation'',''rotation'',''kps''}',...
%  'opt.lambda','[1 1 200 0.2 0.2 5]','opt.truncatedLoss','1'});
%
% % % Pose Hypotheses w/o keypoints
%  mainTest(class,trainId,strcat(trainId,'poseHyp1'),{'opt.rotationPredNum','1',...
%  'opt.relaxInit','{''scale'',''translation'',''rotation'',''kps''}',...
%  'opt.relaxOpt','{''scale'',''translation'',''rotation''}','opt.useOrigTri','1'});
%
%  mainTest(class,trainId,strcat(trainId,'poseHyp1SDS'),{'opt.rotationPredNum','1',...
%  'opt.relaxInit','{''scale'',''translation'',''rotation'',''kps'',''mask''}',...
%  'opt.relaxOpt','{''scale'',''translation'',''rotation''}','opt.useOrigTri','1',...
%  'opt.lambda','[1 1 200 0.2 0.2 5]','opt.truncatedLoss','1'});

% mainTest(class,trainId,strcat(trainId,'poseHyp2'),{'opt.rotationPredNum','2',...
% 'opt.relaxInit','{''scale'',''translation'',''rotation'',''kps''}',...
% 'opt.relaxOpt','{''scale'',''translation'',''rotation''}','opt.useOrigTri','1'});
%
% mainTest(class,trainId,strcat(trainId,'poseHyp3'),{'opt.rotationPredNum','3',...
% 'opt.relaxInit','{''scale'',''translation'',''rotation'',''kps''}',...
% 'opt.relaxOpt','{''scale'',''translation'',''rotation''}','opt.useOrigTri','1'});
%
% mainTest(class,trainId,strcat(trainId,'poseHyp4'),{'opt.rotationPredNum','4',...
% 'opt.relaxInit','{''scale'',''translation'',''rotation'',''kps''}',...
% 'opt.relaxOpt','{''scale'',''translation'',''rotation''}','opt.useOrigTri','1'});
%
% Test with keypoints
%mainTest(class,trainId,strcat(trainId,'withKpsGtSubtype'),{...
%'opt.relaxInit','{}','opt.relaxOpt','{''scale'',''translation'',''rotation'',''kps''}'});
%mainTest(class,trainId,strcat(trainId,'withKpsPredSubtype'),{'opt.useGtSubtypes','0',...
%'opt.relaxInit','{}','opt.relaxOpt','{''scale'',''translation'',''rotation'',''kps''}'});

% Test without keypoints
%mainTest(class,trainId,strcat(trainId,'woKpsGtSubtype'),{...
%'opt.relaxInit','{''scale'',''translation'',''rotation'',''kps''}',...
%'opt.relaxOpt','{''scale'',''translation'',''rotation''}'});
%mainTest(class,trainId,strcat(trainId,'woKpsPredSubtype'),{...
%'opt.relaxInit','{''scale'',''translation'',''rotation'',''kps''}',...
%'opt.relaxOpt','{''scale'',''translation'',''rotation''}'});
end
% Run locally
%evalMeshesAll(classes,trainId,0);
%evalMeshesAll(classes,strcat(trainId,'withKpsgtSubtype'),1);
%evalMeshesAll(classes,strcat(trainId,'predSubtype'),1);
%[scoresMeshes, scoresDmaps] = collectEvals(classes,trainId,0);
%[scoresMeshes, scoresDmaps] = collectEvals(classes,strcat(trainId,'predSubtype'),1);
%[scoresMeshes, scoresDmaps] = collectEvals(classes,strcat(trainId,'gtSubtype'),1);
%And the same for IoU sets

repo: CategoryShapes


text:
function [TorF, vstr, rdate] = have_feature_evalc()
%HAVE_FEATURE_EVALC  Detect availability/version info for EVALC
%
%   Feature detection function implementing 'evalc' tag for HAVE_FEATURE
%   to detect support for EVALC function.
%
%   See also HAVE_FEATURE, EVALC.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

TorF = exist('evalc') ~= 0;
if TorF
    if have_feature('matlab')
        v = have_feature('matlab', 'all');
    else
        v = have_feature('octave', 'all');
        TorF = 0;   %% In Octave, evalc was introduced in 4.2.x, but as of 6.x
                    %% we still choose to ignore evalc in Octave since it does
                    %% not capture the output of all functions (e.g. .mex/.oct
                    %% functions such as ipopt, glpk), which is one of the
                    %% main reasons we'd like to use it. Not only does it not
                    %% capture the output, it also does not prevent it from
                    %% going to the console.
    end
    vstr  = v.vstr;
    rdate = v.date;
else
    vstr = '';
    rdate = '';
end

author: MATPOWER

repo: matpower