text:
function varargout = mcint(h, n, opt)

%MCINT Monte Carlo integrator for multiple functions
%
% Satisfy your wildest integration fantasies!
%
%
%
% written by Lee Ferchoff, 2006
% Do not modify this file in any way unless for your own personal
% use. Do not re-release any part of the code contained in this file.
% Do not attach your own name to any part of the code contained in this
% file. Please reference me and e-mail to inform me if you use any part
% of this code, or original ideas from the algorithm. Thanks! Hope you find
% it helpful.
% E-mail: umferch1 at cc DOT umanitoba DOT ca
%
%
%
%
% A real-valued integral is calculated in arbitrary n-dimensional
% coordinates by averaging the value of a
% function over a large number of randomly selected points within the
% hypervolume to be integrated. Multiple functions can be integrated
% simultaneously over the same domain. The hypervolume can be of
% arbitrary shape, as long as it can be expressed as a series of
% logical conditions on the coordinates.
%
% Examples of the usage of mcint can be found in learnmcint.m
%
% Throughout, x(i,j) contains the matrix of "darts", or random
% points used to evaluate the function to be integrated.
% The rows i, refer to the ith coordinate of a
% point, where each column j, is a separate point. In standard
% notation, the coordinate vectors of all points 'n' are given by:
%
%       1-coordinate vector is     x(1,:)
%       2-coordinate vector is     x(2,:)
%       3-coordinate vector is     x(3,:)
%       ...
%       mth-coordinate is   x(m,:)
%
% In Cartesian coordinates, you might define the 1-dimension to be x,
% the 2-dimension to be y, or if you are using spherical coordinates,
% you could define the 1-dimension to be r, the 2-dimension to be
% theta, and the 3-dimension to be phi.
%
% This is the format that must be used when passing functions
% of coordinates to mcint. Operations on these coordinates must
% be written as vector operations (i.e. .* and .^, not * and ^).
%
% To integrate in a non-Cartesian coordinate system, in your
% function definition, multiply by the Jacobian from
% the Cartesian to 'system' transformation, for whatever coordinate
% system you are using. Jacobians for common coordinate systems
% can be obtained from jacobian.m. For example, to integrate
% f(r,theta,phi) in spherical coordinates, on paper you would write
%
%   triple integral ( f(r,theta,phi) * r^2 * sin(theta) dr d(theta) d(phi) )
%
% So the integrand should be defined as:
%
%   "f(x)" .* x(1,:).^2 + sin(x(2,:))
%       or making use of jacobian.m;
%   "f(x)" .* jacobian(x,'spherical')
%
% You can write your own Jacobian in, if you choose not to call jacobian.m.
% This flexibility allows you to mix coordinate systems. For example, you
% can do a 6-dimensional integral where three of the dimensions are in
% spherical coordinates and the other three are in elliptical. Furthermore,
% you can choose whatever order you want for the coordinate
% defintions, as long as you write the Jacobian accordingly. For example,
% x(1,:) could be theta, x(2,:) r, x(3,:) the angular variable in the
% elliptic coordinates, and so on.
%
%
%INPUT
%
% h -- data structure defining a hypervolume integration problem
%
%   h.b(i,j) -- bounding parallepiped of the hypervolume to be
%       integrated.
%       size(b) = [dim 2] where dim is the number of dimensions
%       of the hypervolume. (i,1) is the minimum bounding coordinate
%       in the ith dimension and (i,2) is the maximum bounding coordinate
%       in the ith dimension.
%
%       Example: To integrate x from -1 to 1, and y from 0 to Inf,
%                b(:,1) = [-1; 0]
%                b(:,2) = [1; Inf]
%
%   h.cond{1,j} (OPTIONAL) -- cell array which contains the functions
%       expressing additional logical conditional statements that define
%       a hypervolume within its bounding box. These conditional
%       statements, indexed by j, are functions.
%       h.cond may be empty, if one wishes to integrate
%       over a rectangular parallepiped and doesn't
%       require any further logical conditions on the
%       hypervolume.
%
%       Example: cond = {@outCircle @aboveLine}
%                where the functions are defined as:
%
%                function tf = outcircle(x, par)
%                   tf = x(1,:).^2 + x(2,:).^2 >= 1;
%
%                function tf = aboveLine(x, par)
%                   tf = x(2,:) > x(1,:);
%
%                In this example, our hypervolume of integration
%                is restricted to the region outside of circle of
%                radius 1 at the origin, and to the upper left
%                of the line y=x.
%
%       All h.cond function definitions should take a structure par as
%       their second input, even if no parameters from 'par' are used in
%       the funciton.
%
%   h.funcs{1,j} -- cell array which contains the functions to be
%       integrated. When the jth function is evaluated, it returns
%       the evaluation of the jth integrand. The usual usage is
%       to enter a function that returns an row vector for the integrand
%       f(1,k), where the kth column is the integrand evaluated at the
%       kth dart point. mcint tests the conditions in the order you enter
%       them in the cell array, therefore it is most efficient to list
%       them from most likely to fail to least likely to fail.
%
%       Example: function f = simpleExample(x, par)
%                   f = x(1,:).^2 .* 2*exp(-x(2,:));
%
%       In this case, the integral returned will be a scalar, I(1,1)
%       containing the value of the integral.
%
%       If h.funcs contains several functions, I(1,m) will be the
%       value of the integral for the mth function in h.funs.
%
%       Example: function f = simpleExample(x, par)
%                   f = x(1,:).^2 .* 2*exp(-x(2,:));
%                function f = Id(x)
%                   f = ones( 1, size(x,2));
%
%       In this example, I(1,1) will be the integral of @simpleExample
%       and I(1,2) will be the integral of @Id. Note that @Id is the
%       identity function; to obtain the hypervolume itself, simply
%       integrate this identity function.
%
%       To allow greater flexibility, a function handle may contain
%       more than one function. If h.funcs contains one function handle
%       which returns a function with m rows, then I(m,1) will be
%       the integral for the integrand defined by the mth row returned
%       by the function.
%
%       Example: function f = energy(x, par)
%                   p = calculatePsi(x);
%                   H = getHamiltonian(x);
%                   f(1,:) = conj(p).*H.*p;
%                   f(2,:) = conj(p).*p;
%
%       This example shows you how you might calculate the expectation
%       value of energy for a quantum system. Here we've calculated
%       p (the wavefunction psi) and H (the Hamiltonian) on the random
%       dart matrix (x). Row one of the function will be the
%       evaluation of the integrand for the expectation value of the
%       Hamiltonian, and row two will be the probability density.
%       mcint will give I(1,1) which is the expectation value of the
%       Hamiltonian, and I(2,1) which is the total probability. Then
%       the user can simply find the expectation value of energy as
%       I(1,1)/I(2,1).
%
%       This method of putting multiple integrands in one function
%       handle is preferable when they are calculated with the same
%       ingredients. It would be more time consuming to make a separate
%       function handle for both expectation values calculated above
%       since the wavefunction psi would have to be calculated in
%       BOTH of the function handles with a call to @calculatePsi:
%                   p = calculatePsi(x);
%
%       This inefficiency can thus be avoided. In Quantum Energy Solver
%       (Lee Ferchoff), I use this method to calculate the expectation
%       value of the Hamiltonian, and the identity operator for ALL of
%       the individuals in the population at at once under one function
%       handle.
%
%       NOTE: As of now, if multiple function handles are used, each one
%             must contain the same number of rows. Support will be added
%             for varying numbers of rows in a future version.
%
%       In general, if h.funcs contains j function handles, then mcint
%       returns I, where I(i,j) is the integral of the function in the ith
%       row of the jth function handle.
%
%       All h.funcs function definitions should take a structure 'par' as
%       their second input, even if no parameters from 'par' are used in
%       the funciton.
%
%
% n -- the number of randomly chosen sampling points
%
%
% opt (OPTIONAL) -- integration options
%
%   opt.maxArray -- adjustable max array size parameter. mcint has
%       built-in memory management. This is the largest number
%       of elements that will be stored in memory for
%       use by mcint at any given time.
%
%       NOTE: The default max array size of 1e4 has been tested
%       extensively, this gives the shortest run time for a fixed amount
%       of points. This has been tested for integrals of dimensions 1 and
%       3 and n of 1e5 and 5e6, and it was found that run time increased
%       if the max array size was either increased or decreased from 1e4
%       in all cases.
%       If it is decreased from 1e4, the run time sharply increases. This
%       is because we lose most of the vectorization of the code,
%       since the function evaluations are broken up into vectorized
%       groups of the maximum array size. If the max array size is
%       increased from 1e4, the run time gradually increases by about
%       25% until memory management is not being used and the whole
%       integration is done in one shot. I suppose this is because
%       doing stuff with really large matricies is unwieldy and in
%       this case breaking up a bit of the vectorization at the benefit
%       of using smaller matricies seems to be a good thing.
%
%       1e-4 is optimum for the linux boxes in Room 523.
%       Anything above 1e-3 seems to be handled equally well by
%       hactar.
%   opt.noerror -- if true, statistical errors are not calculated and
%       only the integral value is returned. Note that this mode cannot
%       be used when a tolerance is requested with opt.tol.
%   opt.tol -- TOLERANCE STOPPING CONDITION
%       absolute tolerance requested (integration will stop
%       when reached)
%       opt.tol(i,j) is the absolute tolerance
%         for the ith row of the jth function in h.funcs. If a tolerance
%         is set to zero, then it is taken to mean that no tolerance is
%         requested for that particular integral.
%   opt.time -- TIME STOPPING CONDITION
%       maximum time allowed in seconds. Integration
%       will end normally after this time.
%   opt.warnOff -- if true, turn off printed warnings, such as not meeting a
%       requested tolerance
%   opt.mcintPar -- a struture containing parameters that is passed to the
%       functions in h.funcs, which can be used in function evaluation
%
%OUTPUT
%
% varargout(1) is
% I(i,j,k) -- integral (with optional error)
%     I(i,j,1) -- value of the integral for the integrand defined
%     in the ith row of the jth function handle in h.funcs.
%     I(i,j,2) -- value of the statistical error of the integral
%     for the integrand defined in the ith row of the jth function
%     handle in h.funcs. The error is calculated from random fluctations
%     of the integrand about its mean value.
%
%     Errors are calculated by default, so k=1:2 by default.
%     If user passes opt.noerror=1 to mcint, then no errors are
%     calculated, and so k=1 only.
%
% varargout(2) is
% info (OPTIONAL) -- technical data about the integration. The data
%     can be used for error checking.
%
%     info.maxArray -- maximum array size used
%     info.memoryManaged -- if true, the requested number of points
%         n was large enough, that the integration had to be
%         performed with chunks of size maxArray
%     info.noerror -- if true, no statistical error calculated
%     info.extraParameters -- if true, extra parameters were passed
%         to the functions to be integrated in opt.mcintPar
%
%     info.dim -- number of dimensions
%     info.n -- number of random points requested
%     info.points -- number of random points used
%     info.pointdiff -- info.n-info.points
%
%     info.tolReq(i,j) -- tolerance requested by user. The abs tolerance
%         for the ith row of the jth function in h.funcs. If a tolerance
%         is set to zero, then it is taken to mean that no tolerance is
%         requested for that particular integral. This is what the user
%         entered in opt.tol.
%     info.someTolReq -- if true, at least one integral tolerance has
%         been requested
%     info.tolMet(i,j) -- true in the indicies for integrals which
%         have had their requested tolerance met. (i,j) labels the ith
%         row of the jth function in h.funcs
%         info.tolMet=-1 in the indicies representing integrals for which
%         a specific tolerance was not requested. If no tolerance was
%         required for any integral, info.tolMet is a scalar -1.
%     info.allTolMet -- if true, the requested tolerance for every
%         tolerance was met
%     info.timeReq -- max time requested by the user. info.timeReq=-1
%         if a specific time limit was not requested
%     info.timeMet -- if true, the integration was within the requested
%         time. info.timeMet=-1 if a specific time limit was not
%         requested
%     info.totalTime -- total time required for integration
%
%
%FEATURES TO BE ADDED IN THE NEXT VERSION:
%
% - relative tolerance stopping condition
% - differening amounts of rows in the functions in h.funcs
% - protection against overflow from a chance hit near a singularity
%
%Lee Ferchoff

t0 = clock; % record initial time

if ~isfield(h,'cond')
    h.cond = {};
end

% ===========================
% SET OPTIONS
% ===========================

% --------------------------------------
% mcint is memory managed, and will not
% use an array that has more elements
% than the number maxArray
p.info.maxArray = 1e4; % default max array size
% will be changed if opt.maxArray is a field
% --------------------------------------

% -1 means no request
p.info.tolReq = -1;
p.info.someTolReq = -1;
p.info.tolMet = -1;
p.info.allTolMet = -1;
p.info.timeReq = -1;
p.info.timeMet = -1;

p.info.noerror = 0; % if true, statistical errors are not calculated
p.info.extraParameters = 0; % if true, the user has provided opt.mcintPar

warnOff = 0; % if true, no warnings will be printed to screen
p.par = [];

if nargin==3 % optional parameter structure opt has been provided
    if isfield(opt,'maxArray')
        p.info.maxArray = opt.maxArray;
    end
    if isfield(opt,'noerror')
        p.info.noerror = opt.noerror;
    end
    if isfield(opt,'mcintPar') % user has provided extra parameters
        p.info.extraParameters = 1;
        p.par = opt.mcintPar;
    end
    if isfield(opt,'tol') % user requested tolerances
        p.info.tolReq = abs(opt.tol);
        p.info.someTolReq = 1;
        p.info.tolMet = 0;
        p.info.allTolMet = 0;
    end
    if isfield(opt,'time') % user requested max integration time
        p.info.timeReq = opt.time;
        p.info.timeMet = 0;
    end
    if isfield(opt,'warnOff')
        warnOff = opt.warnOff;
    end
end
% error: can't not calculate error and still request an error tolerance
if p.info.noerror==1 && p.info.someTolReq==1
    if warnOff==0
        disp('WARNING: IF a tolerance is requested, error must be calculated.');
        disp('Integration will continue WITH error calculation.');
    end
    p.info.noerror = 0;
end

% ===========================
% SETUP INTEGRATION PROBLEM
% ===========================

p.dim = size(h.b, 1);
p.nf = length(h.funcs); % # functions
points = n; % running counter that keeps track of how many of the n darts have been thrown

% Transform any improper integral bounds:

% true in rows of b having at least one improper integral limit
p.infBound = sum(isinf(h.b),2) ~= 0;
% perform inverse tangent coordinate transformation on affected rows.
% This will allow us to calculate the improper integral in the old
% coordinate system as a proper integral in the new coordinate system.
h.b(p.infBound,:) = atan(h.b(p.infBound,:));
% p.infDims isthe numbers of the dimensions with improper limits
p.infDims = find(p.infBound==1);
p.improperProblem = ~isempty(p.infDims);  % if true, at least one infinite dim
V = prod(h.b(:,2) - h.b(:,1)); % bounding volume (after transformation)

% ===========================
% INTEGRATE
% ===========================

if p.dim*n <= p.info.maxArray

    p.info.memoryManaged = 0;

    p.cols = n;
    p.xMin=repmat(h.b(:,1),1,p.cols);
    p.xMax=repmat(h.b(:,2),1,p.cols);
    s = updateSums(h,p);

else % memory management needed

    p.info.memoryManaged = 1;

    p.cols = floor(p.info.maxArray/p.dim);
    p.xMin=repmat(h.b(:,1),1,p.cols);
    p.xMax=repmat(h.b(:,2),1,p.cols);

    % # of runs which use the full maxArray size
    fullruns = ceil(p.dim*n/p.info.maxArray)-1;
    for j=1:fullruns
        newS = updateSums(h,p);
        if ~exist('s') % first run through loop; initialize s, s2
            s = zeros(size(newS));
        end
        s = s + newS;
        points = j*p.cols;

        % Check for stopping conditions
        if (p.info.timeReq~=-1 && etime(clock,t0)>p.info.timeReq)
            % executes if max time has elapsed
            p.info.timeMet = 1;
            break;
        elseif p.info.someTolReq==1
            % checks if all tolerances have been met
            E = getError(s,points,V);
            tolMatrix = zeros(size(E));
            tolMatrix(:) = Inf;
            tolMatrix = p.info.tolReq(find(p.info.tolReq~=0));
            p.info.tolMet = E < tolMatrix;
            p.info.allTolMet = prod(+p.info.tolMet(:));
            if p.info.allTolMet==1
                break;
            end
        end
    end

    % if integration wasn't stopped by some condition, finish the leftover
    % points
    if ~(p.info.timeMet==1 || p.info.allTolMet==1)
        % The last update covers the leftover points on the exact
        % number of points the user requested

        p.cols = n - points;
        p.xMin=repmat(h.b(:,1),1,p.cols);
        p.xMax=repmat(h.b(:,2),1,p.cols);

        newS = updateSums(h,p);
        s = s + newS;

        points = points + p.cols;
    end

end

% ===========================
% REPORT RESULT
% ===========================

% return integral and possibly error
I(:,:,1) = V/points * s(:,:,1); % integral
if p.info.noerror==0
    I(:,:,2) = getError(s,points,V); % error
end
varargout(1) = {I};

% some tolerance was requested and wasn't met
if p.info.allTolMet==0 && warnOff==0
    disp('***');
    disp('Max number of requested darts thrown.');
    disp('Requested tolerance was not achieved.');
    disp('***');
end

 % build info structure to pass to user
if nargout==2
    p.info.dim = p.dim;
    p.info.n = n;
    p.info.points = points;
    p.info.pointdiff = n - points;
    p.info.totalTime = etime(clock,t0);

    varargout(2) = {p.info};
end

function E = getError(s,points,V)
%GETERROR Find the current statistical error of the integral
%
%INPUT
%   s -- s(:,:,1) function sums
%     -- s(:,:,2) function squared sums
%   points -- number of points in sum
%   V -- hypervolume
%
%OUTPUT
%
%   E(i,j) -- the current error of the ith row and jth column using the
%             function sums in s

E = V/points * sqrt(s(:,:,2) - s(:,:,1).*s(:,:,1)/points);

function s = updateSums(h, p)

%UPDATESUMS Calculate function sums for mcint
%
%INPUT
%   h -- data structure defining the hypervolume integration problem
%   p -- data structure containing the current parameters of the
%        integration
%
%OUTPUT
%   s(i,j,k)  -- k=1, sum of the functions over all darts
%                k=2, sum of the functions squared over all darts
%                   for the ith row of the jth function handle in h.funcs
%
%   if p.info.noerror=1, then the k=2 matrix is not returned

% p.xMin, p.xMax passed into function rather than being calculated
% in it so that p.xMin, p.xMax only need to be calculated once
% for many memory blocks of size info.maxArray

% ===========================
% GENERATE RANDOM DARTS
% ===========================

% x(i,j) is the ith Cartesian coordinate of the jth random point
x = p.xMin + (p.xMax-p.xMin) .* rand(p.dim,p.cols); % random dart matrix

% ===========================
% COORDINATE TRANSFORMATIONS
% ===========================

% perform coordinate transformations if there is at least one improper
% integral. Transform only the infinite dimensions, since the
% transformation costs time and is not needed for the finite dimensions
if p.improperProblem
    x(p.infDims,:) = tan( x(p.infDims,:) );
end

% remove all points from x which aren't in hypervolume
% by using the extra restrictions given on the hypervolume
for i=1:length(h.cond)
    x = x( :, feval(h.cond{i},x));
end

% we calculate the Jacobian of the arctan transformation AFTER removing
% the points that aren't in the hypervolume so we don't waste time
% calculating the Jacobian for points that aren't in the integration
if p.improperProblem
    jacobian = prod( x(p.infDims,:).^2 + 1, 1 );
end

% ===========================
% INTEGRAND EVALUATION
% ===========================

for i = p.nf: -1: 1
    f = feval(h.funcs{i}, x, p.par);
    if p.improperProblem % multiply by the Jacobian of the arctan transform
        jacobian = repmat(jacobian(1,:), [size(f,1) 1]);
        f = f .* jacobian;
    end
    s(:,i,1) = sum(f,2);
    if p.info.noerror==0
        s(:,i,2) = sum(f.*f,2);
    end
end


text:
function zero_rc_test ( )

%*****************************************************************************80
%
%% ZERO_RC_TEST tests the ZERO_RC library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ZERO_RC_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ZERO_RC library.\n' );

  zero_rc_test01 ( );
  zero_rc_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ZERO_RC_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function node_rhs = rhs ( node_num, node_xy )

%*****************************************************************************80
%
%% RHS gives the right-hand side of the differential equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM),
%    the coordinates of the points.
%
%    Output, real NODE_RHS(NODE_NUM,1), the value of the
%    right hand side function at the points.
%
  node_rhs = zeros ( node_num, 1 );

  for j = 1 : node_num
    z = 4.0 - sqrt ( ( node_xy(1,j) - 2.0 ).^2 + ( node_xy(2,j) - 2.0 ).^2 );
    node_rhs(j,1) = max ( z, 0.0 );
  end

  return
end


text:
function cvx_optval = det_rootn( X )

%DET_ROOTN nth-root of the determinant of an SPD matrix.
%   For a square matrix X, DET_ROOTN(X) returns
%       POW(DET(X),1/(size(X,1))
%   if X is symmetric (real) or Hermitian (complex) and positive semidefinite,
%   and -Inf otherwise.
%
%   This function can be used in many convex optimization problems that call for
%   LOG(DET(X)) instead. For example, if the objective function contains nothing
%   but LOG(DET(X)), it can be replaced with DET_ROOTN(X), and the same optimal
%   point will be produced.
%
%   Disciplined convex programming information:
%       DET_ROOTN is concave and nonmonotonic; therefore, when used in
%       CVX specifications, its argument must be affine.

error( nargchk( 1, 1, nargin ) );
if ndims( X ) > 2 || size( X, 1 ) ~= size( X, 2 ),
    error( 'Second argument must be a square matrix.' );
end
err = X - X';
X   = 0.5 * ( X + X' );
if norm( err, 'fro' ) > 8 * eps * norm( X, 'fro' ),
    cvx_optval = -Inf;
else
    [R,p] = chol(X);
    if p > 0,
        cvx_optval = geo_mean(eig(full(X)));
    else
        cvx_optval = geo_mean(diag(R)).^2;
    end
end

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

repo: REST


text:
function os_seg_test(varargin)
% note: this function runs multiple setups (feature combinations) at the
% same time

% load stuff
[opts, imdb] = os_setup(varargin{:}) ;
classes = find(imdb.meta.inUse) ;

figure(200) ; clf ;
cmap = plot_legend(imdb) ;
for s = 1:numel(opts.suffix)
  vl_xmkdir(opts.segPublishDir{s});
  print(fullfile(opts.segPublishDir{s}, [opts.dataset '-legend.pdf']), '-dpdf') ;
end

info = struct() ;
encoders = struct() ;
models = struct() ;
for s = 1:numel(opts.suffix)
  vl_xmkdir(fileparts(opts.segResultPath{s})) ;
  setupName = opts.suffix{s} ;
  fprintf('preparing setup %s\n', setupName) ;

  % make a directory to contain the segment scores for each image
  [setupLoc,setupBaseName,~] = fileparts(opts.segResultPath{s}) ;
  vl_xmkdir(fullfile(setupLoc, setupBaseName)) ;

  % source all encoders for this setup (could be the same encoder of
  % another setup)
  for j = 1:numel(opts.encoders{s})
    name = opts.encoders{s}{j}.name ;
    if ~isfield(encoders, 'name')
      fprintf('loading encoder %s\n', name) ;
      encoders.(name) = load(opts.encoders{s}{j}.path) ;
      if isfield(encoders.(name), 'net')
        if opts.useGpu
          encoders.(name).net = vl_simplenn_move(encoders.(name).net, 'gpu') ;
          encoders.(name).net.useGpu = true ;
        else
          encoders.(name).net = vl_simplenn_move(encoders.(name).net, 'cpu') ;
          encoders.(name).net.useGpu = false ;
        end
      end
    end
  end

  % load the SVM classifiers for this setup
  tmp = load(opts.resultPath{s}) ;
  models(s).w = tmp.w ;
  models(s).b = tmp.b ;

  info.(setupName) = init_info(classes) ;
end

% -------------------------------------------------------------------------
%                                                            Process images
% -------------------------------------------------------------------------

rng(0) ;
printThisImage = false(1, numel(imdb.images.id)) ;
test = find(imdb.images.set <= 3) ;
printThisImage(vl_colsubset(test, 60)) = true ;
confusion = cell(1, numel(imdb.images.id)) ;

%fprintf('num labs: %d\n', matlabpool('size')) ;
parfor i = 1:numel(imdb.images.id)
  confusion{i} = do_one_image(opts, encoders, models, imdb, classes, cmap, printThisImage, i) ;
end

for s = 1:numel(opts.suffix)
  setupName = opts.suffix{s} ;
  for i = 1:numel(imdb.images.id)
    if ismember(imdb.images.set(i), [1 2]), setName = 'train' ; else setName = 'test' ; end
    info.(setupName).(setName).confusion_ = info.(setupName).(setName).confusion_ + confusion{i}{s} ;
  end
  for setName = {'train', 'test'}
    setName = char(setName) ;
    [info.(setupName).(setName).confusion, ...
      info.(setupName).(setName).acc, ...
      info.(setupName).(setName).msrcAcc] = ...
      nconfusion(info.(setupName).(setName).confusion_) ;
  end

  % save
  [setupLoc,setupBaseName,~] = fileparts(opts.segResultPath{s}) ;
  info_ = info.(setupName) ;
  save(opts.segResultPath{s}, '-struct', 'info_') ;

  % print the confusion
  figure(1) ; clf ; t = 0 ;
  for setName = {'train', 'test'}
    t = t + 1 ;
    setName = char(setName) ;
    subplot(1,2,t) ; imagesc(info.(setupName).(setName).confusion) ; axis equal ;
    title(sprintf('%s %s acc:%.2f%% (msrc acc: %.2f%%)', ...
      setupName, setName, ...
      info.(setupName).(setName).acc*100, ...
      info.(setupName).(setName).msrcAcc*100)) ;
  end
  vl_printsize(1) ;
  print(fullfile(setupLoc, [setupBaseName '.pdf']), '-dpdf') ;
end

function saveIntermediate(scoresPath, scores, confusion_)
save(scoresPath, 'scores', 'confusion_') ;

function confusion_ = loadIntermediate(scoresPath)
load(scoresPath, 'confusion_') ;

% -------------------------------------------------------------------------
function confusion = do_one_image(opts, encoders, models, imdb, classes, cmap, printThisImage, i)
% -------------------------------------------------------------------------
if ismember(imdb.images.set(i), [1 2]), setName = 'train' ; else setName = 'test' ; end
task = getCurrentTask() ;
if ~isempty(task), tid = task.ID ; else tid = 1 ; end
fprintf('Task %5d: segmenting image %d of %d (%s)\n', tid, i, numel(imdb.images.id), setName) ;
[~,baseName,~] = fileparts(imdb.images.name{i}) ;
features = struct() ;
for s = 1:numel(opts.suffix)
  setupName = opts.suffix{s} ;
  [setupLoc,setupBaseName,~] = fileparts(opts.segResultPath{s}) ;
  scoresPath = fullfile(setupLoc, setupBaseName, [baseName '-scores.mat']) ;
  nClasses = numel(classes);

  seg_labels = imdb.segments.label(:, imdb.segments.imageId == imdb.images.id(i));
  if (max(max(seg_labels)) == 0)
    confusion_ = zeros(nClasses, nClasses + 1);
  elseif exist(scoresPath)
    % load previous results
    confusion_ = loadIntermediate(scoresPath) ;
  else
    % get the region proposals
    im = imread(fullfile(imdb.imageDir, imdb.images.name{i})) ;
    switch opts.segProposalType
      case 'crisp'
        regions = read_crisp_regions(...
          fullfile(imdb.segmDir, 'crisp', sprintf('%s.png', baseName)), ...
          [size(im,1) size(im,2)]) ;
      case 'scg'
        regions = read_scg_regions(...
          fullfile(imdb.segmDir, 'scg_proposals', sprintf('%s.mat', baseName))) ;
      case 'scg200'
        regions = read_scg_regions(...
          fullfile(imdb.segmDir, 'scg_proposals', sprintf('%s.mat', baseName)), ...
          'maxNumRegions', 200) ;
    end

    % get the region features
    psi = {} ;
    for j = 1:numel(opts.encoders{s})
      n = opts.encoders{s}{j}.name ;
      if ~isfield(features, n)
        features.(n) = encoder_extract_for_regions(encoders.(n), im, regions) ;
      end
      psi{j} = features.(n) ;
    end
    psi = cat(1, psi{:}) ;

    % evaluate classifiers
    scores = bsxfun(@plus, models(s).w' * psi, models(s).b') ;

    % compute confusion and save results
    [~,~,gt] = os_get_gt_regions(imdb, imdb.images.id(i)) ;
    [confusion_, gt_, pred_] = evaluate_seg(classes, gt, regions, scores) ;
    saveIntermediate(scoresPath, scores, confusion_) ;

    % optionally print the image
    if printThisImage(i)
      err = (gt_ ~= pred_) & (gt_ > 0) ;
      good = (gt_ == pred_) & (gt_ > 0) ;
      name = sprintf('%s-%s', baseName, setupName) ;
      pdir = opts.segPublishDir{s} ;
      imwrite(im, fullfile(pdir, [baseName '.jpg'])) ;
      imwrite(gt_+1, colormap(cmap), fullfile(pdir, [name '-gt.png'])) ;
      imwrite(pred_+1, colormap(cmap), fullfile(pdir, [name '-pred.png'])) ;
      imwrite(1 + err + 2*good, [[1 1 1];[1 0 0];[0 1 0]], fullfile(pdir, [name '-err.png'])) ;
    end
  end
  confusion{s} = confusion_ ;
end

% -------------------------------------------------------------------------
function [c, gt_, pred_] = evaluate_seg(classes, gt, regions, scores)
% -------------------------------------------------------------------------
numClasses = numel(classes) ;
if (numel(size(gt)) == 2)
    pred = zeros(size(gt)) ;
else
    pred = zeros(size(gt, 2), size(gt, 3));
end
%scores = scores(1:23,:) ;
[~, labels] = max(scores,[],1) ;
for r=1:numel(labels)
  mask = ismember(regions.basis, regions.labels{r}) ;
  pred(mask) = classes(labels(r)) ;
end

if (numel(size(gt)) == 2)
sel = find(gt(:) > 0) ;
[~, gt_] = ismember(gt, classes) ;
[~, pred_] = ismember(pred, classes) ;
pred_(pred_ == 0) = numClasses + 1 ;
c = accumarray([gt_(sel), pred_(sel)], 1, [numClasses numClasses+1]) ;
else

    [~, pred_] = ismember(pred, classes) ;
    pred_(pred_ == 0) = numClasses + 1 ;
    gt_ = zeros(size(gt, 2), size(gt, 3));
    % fill all labels; we do not care of overlaps for now;
    for ll = 1 : size(gt, 1)
        gt0 = reshape(gt(ll, :, :), size(gt, 2), size(gt, 3));
        gt_(gt0 == 1) = ll;
    end

    % now make sure ground truth matches the predictions
    for ll = 1 : size(gt, 1)
        gt0 = reshape(gt(ll, :, :), size(gt, 2), size(gt, 3));
        sel_gt = gt0 == 1;
        sel_pred = pred_ == ll;
        mask_ = min(sel_pred, sel_gt);
        gt_(mask_) = ll;
    end
    sel = find(gt_(:) > 0);
    c = accumarray([gt_(sel), pred_(sel)], 1, [numClasses numClasses+1]) ;
end
if 0
  cmap = [[1 1 1] ; distinguishable_colors(numel(classes))] ;
  figure(100) ; clf ;
  subplot(1,2,1); image(gt_+1); axis equal ; title('(partial) gt') ;
  subplot(1,2,2); image(pred_+1); axis equal ; title('predicted') ;
  colormap(cmap) ;
  drawnow ;
  keyboard
end

% -------------------------------------------------------------------------
function [c,acc,msrcAcc]=nconfusion(c)
% -------------------------------------------------------------------------
msrcAcc = sum(diag(c)) / sum(c(:)) ;
c = bsxfun(@times, 1./max(sum(c,2),1e-10), c) ;
acc = mean(diag(c)) ;

% -------------------------------------------------------------------------
function info = init_info(classes)
% -------------------------------------------------------------------------
nc = numel(classes) ;
info.train.confusion_ = zeros(nc,nc+1) ;
info.train.confusion = zeros(nc,nc+1) ;
info.test.confusion_ = zeros(nc,nc+1) ;
info.test.confusion = zeros(nc,nc+1) ;
info.train.acc = 0 ;
info.train.msrcAcc = 0 ;
info.test.acc = 0 ;
info.test.msrcAcc = 0 ;
info.lastProcessedImage = 0 ;

% -------------------------------------------------------------------------
function cmap = plot_legend(imdb)
% -------------------------------------------------------------------------
numClasses = sum(imdb.meta.inUse) ;
classes = find(imdb.meta.inUse) ;
cmap = [[1 1 1] ; distinguishable_colors(numClasses) ;] ;
h = 16 ;
set(gca,'units', 'points', 'position',[0 0 100*8 h*3]) ;
set(gcf,'paperunits', 'points', 'paperposition',[0 0 150*8 h*3]) ;
set(gcf,'paperunits', 'points', 'papersize',[100*8 h*3]) ;

for i=1:3
  for j=1:8
    t = (i-1)*8 + j;
    if t > numClasses, break ; end
    x = (j-1)*100 ;
    y = (4-i-1)*h ;
    text(x+1.5*h, y+0.5*h, imdb.meta.classes{classes(t)}, ...
      'interpreter', 'none', 'verticalalign', 'middle', 'fontsize', 10) ;
    patch([0 h h 0]+x,[0 0 h h]+y,cmap(t+1,:)) ;
  end
end
axis equal off ;
xlim([0 100*8])
ylim([0 h*3]) ;

% -------------------------------------------------------------------------
function code = encoder_extract_for_regions(encoder, im, regions)
% -------------------------------------------------------------------------
switch encoder.type
  case 'rcnn'
    code = get_rcnn_features(encoder.net, ...
      im, regions, ...
      'regionBorder', encoder.regionBorder) ;
  case 'dcnn'
    code = get_dcnn_features(encoder.net, ...
      im, regions, ...
      'encoder', encoder) ;
  case 'dsift'
    code = get_dcnn_features([], im, regions, ...
      'useSIFT', true, ...
      'encoder', encoder) ;
end
code = cat(2, code{:}) ;

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function [inH2O] = mbar2inH2O(mbar)
% Convert pressure from millibars to inches of water.
% Chad Greene 2012
inH2O = mbar*0.401463;

author: Sable

repo: mcbench-benchmarks


text:
function obj=stat_fit(obj,varargin)
% stat_fit() Display a custom fit of the data
%
% Example syntax gramm_object.stat_fit('fun',@(alpha,beta,x)alpha*cos(x-beta),'disp_fit',true)
%
% This fuction uses the curve fitting toolbox function fit() to
% fit a provided anonymous function with arguments
% (param1,param2,...paramN,x) to the data.
% Parameters:
% - 'fun': anonymous function used for the fit
% - 'StartPoint': Array containing starting values for the
% parameter to be fitted [start_param1,start_param2,...start_paramN]
% - 'intopt': Option passed to predint() for the type of bounds
% to compute, 'observation' for bounds of a new observation
% (default), or 'functional' for bounds of the fitted curve.
% - 'geom', 'fullrange', 'disp_fit' options: see stat_glm()
% - 'fullrange': set to true if you want the fits to be
%   displayed over the whole width of each subplot instead of
%   being displayed over the range of x values used for the fit
% - 'disp_fit': set to true to display the fitted parameters

p=inputParser;
my_addParameter(p,'fun',@(a,b,x)a*x+b);
my_addParameter(p,'StartPoint',[]);
my_addParameter(p,'intopt','observation');
my_addParameter(p,'geom','area');
my_addParameter(p,'fullrange',false);
my_addParameter(p,'disp_fit',false);
parse(p,varargin{:});
obj.geom=vertcat(obj.geom,{@(dobj,dd)my_fit(dobj,dd,p.Results)});
obj.results.stat_fit={};
end

function hndl=my_fit(obj,draw_data,params)

combx=comb(draw_data.x);
comby=comb(draw_data.y);

%Remove NaNs
sel=~isnan(combx) & ~isnan(comby);
combx=combx(sel);
comby=comby(sel);

%Do the fit depending on options
if isempty(params.StartPoint)
    mdl=fit(shiftdim(combx),shiftdim(comby),params.fun);
else
    mdl=fit(shiftdim(combx),shiftdim(comby),params.fun,'StartPoint',params.StartPoint);
end

%Create x values for the fit plot
if params.fullrange
    newx=linspace(obj.var_lim.minx,obj.var_lim.maxx,100)';
else
    newx=linspace(min(combx),max(combx),100)';
end
%Get fit value and CI
newy=feval(mdl,newx);
yci=predint(mdl,newx,1-obj.stat_options.alpha,params.intopt);

obj.results.stat_fit{obj.result_ind,1}.x=newx;
obj.results.stat_fit{obj.result_ind,1}.y=newy;
obj.results.stat_fit{obj.result_ind,1}.yci=yci;
obj.results.stat_fit{obj.result_ind,1}.model=mdl;

%Plot fit
hndl=plotci(obj,newx,newy,yci,draw_data,params.geom);

%Store plotted handles
hnames=fieldnames(hndl);
for k=1:length(hnames)
    obj.results.stat_fit{obj.result_ind,1}.(hnames{k})=hndl.(hnames{k});
end

%Do we display the results ?
if params.disp_fit
    %Set Y position of display
    if obj.firstrun(obj.current_row,obj.current_column)
        obj.extra.mdltext(obj.current_row,obj.current_column)=0.05;
    else
        obj.extra.mdltext(obj.current_row,obj.current_column)=obj.extra.mdltext(obj.current_row,obj.current_column)+0.03;
    end
    %Get formula and parameters
    form=formula(mdl);
    cvals=coeffvalues(mdl);
    cnames=coeffnames(mdl);
    %Replace parameter names by their value in the formula
    for c=1:length(cnames)
        form=strrep(form,cnames{c},num2str(cvals(c),2));
    end
    obj.results.stat_fit{obj.result_ind,1}.text_handle=text('Units','normalized','Position',[0.1 obj.extra.mdltext(obj.current_row,obj.current_column)],'color',draw_data.color,...
        'String',form);
end

end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function [m, e, t] = PlotMultipleLaminarProfiles(ss)

% PlotMultipleLaminarProfiles(scanList);
%
% Ress, 6/04

mrGlobals

% Get view:
selectedVOLUME = viewSelected('Volume');
if ~isfield(VOLUME{selectedVOLUME}, 'mmPerVox'), VOLUME{selectedVOLUME} = loadAnat(VOLUME{selectedVOLUME}); end
view = VOLUME{selectedVOLUME};

% Get the list of scans in the present view:
if ~exist('ss', 'var'), ss = chooseScans(view); end

% Calculate the laminar profiles
[m, e, t] = ROIlaminarProfile(ss);
if isempty(m)
  % Something went wrong, so quit:
  return
end
m = m(ss);
e = e(ss);
t = t(ss);
if isempty(m), return, end

% Plot the profiles with appropriate labels
labels = {dataTYPES(view.curDataType).scanParams.annotation};
labels = labels(ss);
rLabel = view.ROIs(view.selectedROI).name;
PlotLaminarProfiles(t, m, e, labels, rLabel);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function f = flops_spmul(a,b)
% FLOPS_SPMUL    Flops for sparse matrix multiplication.
% FLOPS_SPMUL(a,b) returns the number of flops for a*b, where multiplication
% and addition of zero doesn't count.
% For example:
%   flops_spmul(0,4) is 0.
%   flops_spmul([1 0 1], [2;3;4]) is 3.
%   flops_spmul(eye(3), [2;3;4]) is 3.

nza = (a ~= 0);
nzb = (b ~= 0);
f = nza*nzb;
f = 2*f - (f ~= 0);
f = sum(sum(f));

author: andrewssobral

repo: lrslibrary


text:
function pBar = bst_progress(varargin)
% bst_progress: Manage the Brainstorm progress bar
%
% USAGE : pBar = bst_progress('start', title, msg, valStart, valStop) : Create a progress bar with start and stop bounds
%         pBar = bst_progress('start', title, msg)                    : Create a progress bar (unlimited)
%         pBar = bst_progress('stop')        : stop and hide progress bar
%         pBar = bst_progress('inc', valInc) : increment of 'valInc' the position of the progress bar
%         pBar = bst_progress('set', pos)    : set the position
%          pos = bst_progress('get')         : get the position
%         pBar = bst_progress('text', txt)   : set the text
%    isVisible = bst_progress('isvisible')   : return 1 if progress bar is visible, 0 else
%         pBar = bst_progress('show')        : display previously defined progress bar
%         pBar = bst_progress('hide')        : hide progress bar
%         pBar = bst_progress('setimage', imagefile) : display an image in the wait bar
%         pBar = bst_progress('setlink', url)        : clicking on the image opens a browser to display the url
%         pBar = bst_progress('removeimage')         : Remove the image from the wait bar

% NOTES : The window is created once, and then never deleted, just hidden.
%         Progress bar is represented by a structure:
%            |- jWindow
%            |- jLabel
%            |- jProgressBar
%
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2013

% JAVA imports
import org.brainstorm.icon.*;
import java.awt.Dimension;

global GlobalData;

%% ===== PARSE INPUTS =====
if ((nargin >= 1) && ischar(varargin{1}))
    commandName = varargin{1};
else
    error('Usage : bst_progress(commandName, parameters)');
end

%% ===== GET OR CREATE PROGRESS BAR =====
% Do nothing in case of server mode
if ~isempty(GlobalData) && ~isempty(GlobalData.Program) && isfield(GlobalData.Program, 'GuiLevel') && (GlobalData.Program.GuiLevel == -1)
    if ismember(lower(commandName), {'pos','isvisible'})
        pBar = 0;
    else
        pBar = [];
    end
    return;
end
% If running in NOGUI mode: just display the message in the command window
if ~bst_get('isGUI')
    switch lower(commandName)
        case 'start',     disp(['PROGRESS> ' varargin{2} ': ' varargin{3}]); pBar = [];
        case 'text',      disp(['PROGRESS> ' varargin{2}]); pBar = [];
        case 'get',       pBar = 1;
        case 'isvisible', pBar = 0;
    end
    return;
end
% Get progress bar
if ~isempty(GlobalData) && ~isempty(GlobalData.Program) && isfield(GlobalData.Program, 'ProgressBar') && ~isempty(GlobalData.Program.ProgressBar)
    pBar = GlobalData.Program.ProgressBar;
else
    pBar = [];
end
% Get Brainstorm GUI context
jBstFrame = bst_get('BstFrame');
if isempty(jBstFrame)
    return
end

% Default window size
if isempty(pBar) || strcmpi(commandName, 'removeimage') || (strcmpi(commandName, 'stop') && pBar.isImage)
    DefaultSize = java_scaled('dimension', 350, 130);
end

% If progress bar was not created yet : create it
if isempty(pBar)
    % If action=isvisible: no need to create the progress bar
    if strcmpi(commandName, 'isvisible')
        pBar = 0;
        return
    end
    % Create a JDialog, if possible dependent of the main Brainstorm JFrame
    pBar.jWindow = java_create('javax.swing.JDialog');
    % Set icon
    try
        pBar.jWindow.setIconImage(IconLoader.ICON_APP.getImage());
    catch
        % Old matlab... just ignore...
    end
    % Set as always-on-top / non-focusable
    pBar.jWindow.setAlwaysOnTop(1);
    pBar.jWindow.setFocusable(0);
    pBar.jWindow.setFocusableWindowState(0);
    % Non-modal
    pBar.jWindow.setModal(0);

    % Closing callback
%     if bst_iscompiled()
        pBar.jWindow.setDefaultCloseOperation(pBar.jWindow.HIDE_ON_CLOSE);
%     else
%         pBar.jWindow.setDefaultCloseOperation(pBar.jWindow.DO_NOTHING_ON_CLOSE);
%         java_setcb(pBar.jWindow, 'WindowClosingCallback', @(h,ev)CloseCallback);
%     end

    % Configure window
    pBar.jWindow.setPreferredSize(DefaultSize);
    % Main panel
    pBar.jPanel = java_create('javax.swing.JPanel');
    pBar.jPanel.setLayout(java_create('java.awt.GridBagLayout'));

    % Create objects
    pBar.isImage = 0;
    pBar.jImage = java_create('javax.swing.JLabel');
    pBar.jLabel = java_create('javax.swing.JLabel', 'Ljava.lang.String;', '...');
    pBar.jLabel.setFont(bst_get('Font'));
    pBar.jProgressBar = java_create('javax.swing.JProgressBar', 'II', 0, 99);
    % Update constraints
    UpdateConstraints(0);

    % Add the main Panel
    pBar.jWindow.getContentPane.add(pBar.jPanel);
    pBar.jWindow.pack();
    % Set window size and location
    %pBar.jWindow.setLocationRelativeTo(pBar.jWindow.getParent());
    jLoc = jBstFrame.getLocation();
    jSize = jBstFrame.getSize();
    pos = [jLoc.getX() + ((jSize.getWidth() - DefaultSize.getWidth()) / 2), ...
           jLoc.getY() + ((jSize.getHeight() - DefaultSize.getHeight()) / 2)];
    pBar.jWindow.setLocation(pos(1), pos(2));
    % Save progress bar
    GlobalData.Program.ProgressBar = pBar;
end

% Linux: need to print something on the command window (don't know why...)
if strcmpi(commandName, 'stop') && ismember(computer('arch'), {'glnx86', 'glnxa64'})
    drawnow();
    fprintf(' ');
    fprintf('\b');
    drawnow();
    % The dialog needs to be displayed for a short period before being hidden
    pause(0.05);
end

%% ===== SWITCH BETWEEN COMMANDS =====
switch (lower(commandName))
    % ==== START ====
    case 'start'
        % Set as "always on top"
        java_call(pBar.jWindow, 'setAlwaysOnTop', 'Z', 1);
        java_call(pBar.jWindow, 'setFocusable',   'Z', 0);
        java_call(pBar.jWindow, 'setFocusableWindowState', 'Z', 0);
        % Call: bst_progress(''start'', title, msg)
        if (nargin == 3) && ischar(varargin{3}) && ischar(varargin{2})
            % Get window title and comment
            wndTitle    = varargin{2};
            msg = varargin{3};
            % Set Progress bar in inderminate mode
            pBar.jProgressBar.setIndeterminate(1);
            pBar.jProgressBar.setStringPainted(0);
            % Set progress bar bounds
            pBar.jProgressBar.setMinimum(0);
            pBar.jProgressBar.setMaximum(100);
            % Set initial value to start
            pBar.jProgressBar.setValue(0);

        % Call: bst_progress(''start'', title, msg, start, stop)
        elseif ((nargin == 5) && ischar(varargin{2}) && ischar(varargin{3}) && isnumeric(varargin{4}) && isnumeric(varargin{5}))
            % Get window title and comment
            wndTitle = varargin{2};
            msg = varargin{3};
            % Set Progress bar in derminate mode
            pBar.jProgressBar.setIndeterminate(0);
            pBar.jProgressBar.setStringPainted(1);
            % Get progress bar bounds
            valStart = varargin{4};
            valStop  = varargin{5};
            % Test bounds
            if ( (valStart >= valStop) || (valStop <= 0) )
                % Set indeterminate bounds
                pBar.jProgressBar.setIndeterminate(1);
                pBar.jProgressBar.setStringPainted(0);
                valStart = 0;
                valStop  = 100;
            end
            % Set progress bar bounds
            pBar.jProgressBar.setMinimum(valStart);
            pBar.jProgressBar.setMaximum(valStop);
            % Set initial value to start
            pBar.jProgressBar.setValue(valStart);

        else
            error(['Usage : bst_progress(''start'', title, comment) ' 10 '        bst_progress(''start'', title, comment, valStart, valStop)']);
        end
        % Set window title
        pBar.jWindow.setTitle(wndTitle);
        % Set window comment (central label)
        pBar.jLabel.setText(msg);
        % Show window
        java_call(pBar.jWindow, 'setVisible', 'Z', 1);
        % Set watch cursor
        jBstFrame.setCursor(java_create('java.awt.Cursor', 'I', java.awt.Cursor.WAIT_CURSOR));

    % ==== STOP ====
    case 'stop'
        % Remove the "always on top" status
        java_call(pBar.jWindow, 'setAlwaysOnTop', 'Z', 0);
        java_call(pBar.jWindow, 'setFocusable',   'Z', 1);
        java_call(pBar.jWindow, 'setFocusableWindowState', 'Z', 1);
        % Hide window
        java_call(pBar.jWindow, 'setVisible', 'Z', 0);
        % Restore cursor
        jBstFrame.setCursor([]);
        % Remove image
        if pBar.isImage
            GlobalData.Program.ProgressBar.isImage = 0;
            pBar.jImage.setIcon([]);
            java_setcb(pBar.jImage, 'MouseClickedCallback', []);
            UpdateConstraints(0);
            pBar.jWindow.setPreferredSize(DefaultSize);
            pBar.jWindow.pack();
        end

    % ==== INCREMENT ====
    case 'inc'
        % Parse arguments
        if ((nargin == 2) && isnumeric(varargin{2}))
            valInc = varargin{2};
        else
            error('Usage : bst_progress(''inc'', valInc)');
        end
        % Get current value
        curValue = pBar.jProgressBar.getValue();
        % Get the incremented progress bar position
        newVal = min(curValue + valInc, pBar.jProgressBar.getMaximum());
        pBar.jProgressBar.setValue(newVal);

    % ==== SET POSITION ====
    case 'set'
        % Parse arguments
        if ((nargin == 2) && isnumeric(varargin{2}))
            newVal = varargin{2};
        else
            error('Usage : bst_progress(''set'', pos)');
        end
        % Get current value
        curValue = pBar.jProgressBar.getValue();
        % Get the incremented progress bar position
        if (curValue ~= newVal)
            newVal = min(newVal, pBar.jProgressBar.getMaximum());
            pBar.jProgressBar.setValue(newVal);
        end

    % ==== GET POSITION ====
    case 'get'
        % Get the incremented progress bar position
        pBar = pBar.jProgressBar.getValue();

    % ==== SET TEXT ====
    case 'text'
        % Parse arguments
        if ((nargin == 2) && ischar(varargin{2}))
            % Set new label
            pBar.jLabel.setText(varargin{2});
        else
            % Get label
            pBar.jLabel.getText();
        end

    % ==== IS VISIBLE ====
    case 'isvisible'
        pBar = pBar.jWindow.isVisible();

    % ==== SHOW ====
    case 'show'
        % Set as "always on top"
        java_call(pBar.jWindow, 'setAlwaysOnTop', 'Z', 1);
        java_call(pBar.jWindow, 'setFocusable',   'Z', 0);
        java_call(pBar.jWindow, 'setFocusableWindowState', 'Z', 0);
        % Show window
        java_call(pBar.jWindow, 'setVisible', 'Z', 1);
        % Set watch cursor
        jBstFrame.setCursor(java_create('java.awt.Cursor', 'I', java.awt.Cursor.WAIT_CURSOR));

    % ==== HIDE ====
    case 'hide'
        % Remove the "always on top" status
        java_call(pBar.jWindow, 'setAlwaysOnTop', 'Z', 0);
        java_call(pBar.jWindow, 'setFocusable',   'Z', 1);
        java_call(pBar.jWindow, 'setFocusableWindowState', 'Z', 1);
        % Hide window
        java_call(pBar.jWindow, 'setVisible', 'Z', 0);
        % Restore cursor
        jBstFrame.setCursor([]);

    % ==== SET IMAGE ====
    case 'setimage'
        % Get image path
        imagefile = varargin{2};
        if ~file_exist(imagefile)
            imagefile = bst_fullfile(bst_get('BrainstormDocDir'), imagefile);
        end
        if ~file_exist(imagefile)
            warning(['Image not found: ' imagefile]);
            return
        end
        % Image in label
        pBar.jImage.setIcon(javax.swing.ImageIcon(imagefile));
        GlobalData.Program.ProgressBar.isImage = 1;
        % Extend size of the frame
        UpdateConstraints(1);
        pBar.jWindow.setPreferredSize([]);
        pBar.jWindow.pack();

    % ==== SET LINK ====
    case 'setlink'
        url = varargin{2};
        java_setcb(pBar.jImage, 'MouseClickedCallback', @(h,ev)web(url, '-browser'));

    % ==== REMOVE IMAGE ====
    case 'removeimage'
        % Remove image
        GlobalData.Program.ProgressBar.isImage = 0;
        pBar.jImage.setIcon([]);
        java_setcb(pBar.jImage, 'MouseClickedCallback', []);
        UpdateConstraints(0);
        pBar.jWindow.setPreferredSize(DefaultSize);
        pBar.jWindow.pack();
end

%% ===== ADD COMPONENTS =====
    function UpdateConstraints(isImage)
        import java.awt.GridBagConstraints;
        import java.awt.Insets;
        % Remove all components
        pBar.jPanel.removeAll();
        % Generic constraints
        c = GridBagConstraints();
        c.fill = GridBagConstraints.BOTH;
        c.gridx = 1;
        c.weightx = 1;
        % IMAGE
        c.gridy = 1;
        c.weighty = isImage;
        c.insets = Insets(0,0,0,0);
        pBar.jPanel.add(pBar.jImage, c);
        % TEXT
        c.gridy = 2;
        c.weighty = ~isImage;
        c.insets = Insets(5,12,5,12);
        pBar.jPanel.add(pBar.jLabel, c);
        % PROGRESS BAR
        c.gridy = 3;
        c.weighty = 0;
        c.insets = Insets(0,12,9,12);
        pBar.jPanel.add(pBar.jProgressBar, c);
%         % CANCEL BUTTON
%         c.gridy = 4;
%         c.weighty = 0;
%         c.insets = Insets(0,12,9,12);
%         c.weightx = 0;
%         pBar.jPanel.add(pBar.jButtonCancel, c);
    end

%     %% ===== CLOSE CALLBACK =====
%     function CloseCallback()
%         % Hide progress bar
%         %java_call(pBar.jWindow, 'setVisible', 'Z', 0);
%         bst_progress('stop');
%
%         if bst_iscompiled()
%             try
%                 % Get command window
%                 cmdWindow = com.mathworks.mde.cmdwin.CmdWin.getInstance();
%                 cmdWindow.grabFocus();
%                 %2) Wait for focus transfer to complete (up to 2 seconds)
%                 focustransferTimer = tic;
%                 while ~cmdWindow.isFocusOwner
%                     pause(0.1);  %Pause some small interval
%                     if (toc(focustransferTimer) > 2)
%                         error('Error transferring focus for CTRL+C press.')
%                     end
%                 end
%
%                 %3) Use Java robot to execute a CTRL+C in the (now focused) command window.
%
%                 %3.1)  Setup a timer to relase CTRL + C in 0.3 second
%                 %  Try to reuse an existing timer if possible (this would be a holdover
%                 %  from a previous execution)
%                 t_all = timerfindall;
%                 releaseTimer = [];
%                 ix_timer = 1;
%                 while isempty(releaseTimer) && (ix_timer<= length(t_all))
%                     if isequal(t_all(ix_timer).TimerFcn, @releaseCtrl_C)
%                         releaseTimer = t_all(ix_timer);
%                     end
%                     ix_timer = ix_timer+1;
%                 end
%                 if isempty(releaseTimer)
%                     releaseTimer = timer;
%                     releaseTimer.TimerFcn = @(h,ev)releaseCtrl_C;
%                 end
%                 releaseTimer.StartDelay = 0.3;
%                 start(releaseTimer);
%
%                 %3.2)  Press press CTRL+C
%                 pressCtrl_C();
%             catch
%                 disp('BST> Could not post a CTRL+C signal in the command window.');
%             end
%         end
%     end
%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     function pressCtrl_C()
%         SimKey = java.awt.Robot;
%         SimKey.keyPress(java.awt.event.KeyEvent.VK_CONTROL);
%         SimKey.keyPress(java.awt.event.KeyEvent.VK_C);
%     end
%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     function releaseCtrl_C()
%         SimKey = java.awt.Robot;
%         SimKey.keyRelease(java.awt.event.KeyEvent.VK_CONTROL);
%         SimKey.keyRelease(java.awt.event.KeyEvent.VK_C);
%         jBstFrame.setVisible(1);
%     end

end

author: brainstorm-tools


text:
function tests = test_intersectEdgePolygon3d
% Test suite for the file intersectEdgePolygon3d.
%
%   Test suite for the file intersectEdgePolygon3d
%
%   Example
%   test_intersectEdgePolygon3d
%
%   See also
%     intersectEdgePolygon3d

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2021-02-24,    using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2021 INRAE - BIA-BIBS.

tests = functiontests(localfunctions);

function test_Simple(testCase) %#ok<*DEFNU>
% Test computation for one edge and one polygon

% Compute intersection between a 3D edge and a 3D triangle
pts3d = [3 0 0; 0 6 0;0 0 9];
edge1 = [0 0 0 3 6 9];

inter = intersectEdgePolygon3d(edge1, pts3d);

assertEqual(testCase, size(inter), [1 3]);
assertEqual(testCase, inter, [1 2 3]);

function test_TwoEdges(testCase) %#ok<*DEFNU>
% Test computation for two edges

pts3d = [3 0 0; 0 6 0;0 0 9];
edges = [0 0 0 3 6 9;10 0 0 10 2 3];

[inter, inside] = intersectEdgePolygon3d(edges, pts3d);

assertEqual(testCase, size(inter), [2 3]);
assertEqual(testCase, size(inside), [2 1]);
assertEqual(testCase, inside, [true ; false]);
assertEqual(testCase, inter(inside, :), [1 2 3]);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%DEMO_REGRESSION_ADDITIVE1  Regression demonstration with additive model
%
%  Description
%    A regression demonstration with one input variable and one output
%    variable with Gaussian noise. The output is assumed to be
%    realization of two additive functions and Gaussian noise.
%
%    The model constructed is following:
%
%    The observations y are assumed to satisfy
%
%         y = f + g + e,    where e ~ N(0, s^2).
%
%    f and g are underlying latent functions, which we are
%    interested in. We place a zero mean Gaussian process prior for
%    them, which implies that at the observed input locations
%    latent values have prior
%
%         f ~ N(0, Kf) and g ~ N(0,Kg)
%
%    where K is the covariance matrix, whose elements are given as
%    K_ij = k(x_i, x_j | th). The function k(x_i, x_j | th) is
%    covariance function and th its parameters.
%
%    Since both likelihoods and prior are Gaussian, we obtain a
%    Gaussian marginal likelihood
%
%        p(y|th) = N(0, Kf + Kg + I*s^2).
%
%    By placing a prior for parameters, p(th), we can find the
%    maximum a posterior (MAP) estimate for them by maximizing
%
%       argmax   log p(y|th) + log p(th).
%         th
%
%    After finding MAP estimate or posterior samples of parameters,
%    we can use them to make predictions for the latent functions.
%    For example, the posterior predictive distribution of f is:
%
%       p(f | y, th) = N(m, S),
%       m = Kf * (Kf + Kg + s^2I)^(-1) * y
%       S = Kf - Kf * (Kf + Kg + s^2I)^(-1) * Kf
%
%    (We could integrate also over the parameters with, for
%    example, grid integration or MCMC. This is not demonstrated
%    here but it is done exactly similarly as in the
%    demo_regression1.)
%
%    The demo is organised in four parts:
%     1) data analysis with full GP model
%     2) data analysis with FIC approximation
%     3) data analysis with PIC approximation
%     4) data analysis with CS+FIC model
%
%    For more detailed discussion of Gaussian process regression
%    see Rasmussen and Williams (2006) and for a detailed
%    discussion on sparse additive models see Vanhatalo and Vehtari
%    (2008).
%
% References:
%
%    Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian
%    Processes for Machine Learning. The MIT Press.
%
%    Vanhatalo, J. and Vehtari, A. (2008). Modelling local and global
%    phenomena with sparse Gaussian processes. Proceedings of the 24th
%    Conference on Uncertainty in Artificial Intelligence.
%
% See also
%  DEMO_REGRESSION1, DEMO_SPARSEREGRESION
%

% Copyright (c) 2008-2010 Jarno Vanhatalo
% Copyright (c) 2010 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.

%========================================================
% PART 1 data analysis with full GP model
%========================================================

% Load the data
S = which('demo_regression1');
L = strrep(S,'demo_regression1.m','demodata/maunaloa_data.txt');
data=load(L);
y = data(:, 2:13);
y=y';
y=y(:);
x = [1:1:length(y)]';
x = x(y>0);             % Remove contaminated observations
y = y(y>0);
avgy = mean(y);
y = y-avgy;
xt = [0:0.5:565]';

[n,nin] = size(x);
% Now 'x' consist of the inputs and 'y' of the output.
% 'n' and 'nin' are the number of data points and the
% dimensionality of 'x' (the number of inputs).

% ---------------------------
% --- Construct the model ---
%
% First create squared exponential and piecewise polynomial 2
% covariance functions and Gaussian noise structures and set
% priors for their parameters (if SuiteSparse is not
% installed, use gpcf_sexp instead of gpcf_ppcs2)
pl1 = prior_t('s2', 100, 'nu', 10);
pl2 = prior_t('s2', 10, 'nu', 10);
pm1 = prior_sqrtt('s2',300);
pm2 = prior_sqrtt('s2',10);
pn = prior_logunif();
gpcf1 = gpcf_sexp('lengthScale', 100, 'magnSigma2', 200, 'lengthScale_prior', pl1, 'magnSigma2_prior', pm1);
if exist('ldlchol')
  gpcf2 = gpcf_ppcs2('nin', nin, 'lengthScale', 5, 'magnSigma2', 5, 'lengthScale_prior', pl2, 'magnSigma2_prior', pm2);
else
  warning('GPstuff:SuiteSparseMissing',...
  ['SuiteSparse is not properly installed. (in BECS try ''use suitesparse'')\n' ...
   'Using gpcf_sexp (non-compact support) instead of gpcf_ppcs2 (compact support)']);
  gpcf2 = gpcf_sexp('lengthScale', 5, 'magnSigma2', 1, 'lengthScale_prior', pl2, 'magnSigma2_prior', pm);
end
lik = lik_gaussian('sigma2', 0.1, 'sigma2_prior', pn);

% Create the GP structure
gp = gp_set('lik', lik, 'cf', {gpcf1, gpcf2}, 'jitterSigma2', 1e-9)

% -----------------------------
% --- Conduct the inference ---
%
% --- MAP estimate -----------
% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% Make predictions. Below Ef_full is the predictive mean and Varf_full
% the predictive variance.
[Eft_full, Varft_full, lpyt_full, Eyt_full, Varyt_full] = gp_pred(gp, x, y, xt, 'yt', ones(size(xt)));
[Eft_full1, Varft_full1] = gp_pred(gp, x, y, xt, 'predcf', 1);
[Eft_full2, Varft_full2] = gp_pred(gp, x, y, xt, 'predcf', 2);

% Plot the prediction and data
figure
subplot(2,1,1)
hold on
plot(x,y,'.', 'MarkerSize',7)
plot(xt,Eyt_full,'k', 'LineWidth', 2)
plot(xt,Eyt_full-2.*sqrt(Varyt_full),'g--')
plot(xt,Eyt_full+2.*sqrt(Varyt_full),'g--')
axis tight
caption1 = sprintf('Full GP:  l_1= %.2f, s^2_1 = %.2f, \n l_2= %.2f, s^2_2 = %.2f \n s^2_{noise} = %.2f', gp.cf{1}.lengthScale, gp.cf{1}.magnSigma2, gp.cf{2}.lengthScale, gp.cf{2}.magnSigma2, gp.lik.sigma2);
title(caption1)
legend('Data point', 'predicted mean', '2\sigma error',4)

subplot(2,1,2)
[AX, H1, H2] = plotyy(xt, Eft_full2, xt, Eft_full1);
set(H2,'LineStyle','--')
set(H2, 'LineWidth', 2)
%set(H1, 'Color', 'k')
set(H1,'LineStyle','-')
set(H1, 'LineWidth', 0.8)
title('The long and short term trend')

%========================================================
% PART 2 data analysis with FIC approximation
%========================================================

% Here we conduct the same analysis as in part 1, but this time
% using FIC approximation. Notice that both covariance components
% utilize the inducing inputs. This leads to problems since the
% number of inducing inputs is too small to capture the short term
% variation. In CS+FIC (later model) the compact support function
% does not utilize the inducing inputs and for this reason it is
% able to capture also the fast variations.

% Place inducing inputs evenly
Xu = [min(x):24:max(x)+10]';

% Create the FIC GP structure
gp_fic = gp_set('type', 'FIC', 'lik', lik, 'cf', {gpcf1,gpcf2}, 'jitterSigma2', 1e-9, 'X_u', Xu)

% -----------------------------
% --- Conduct the inference ---

% --- MAP estimate using modified Newton algorithm ---

% Now you can choose, if you want to optimize only parameters or
% optimize simultaneously parameters and inducing inputs. Note that
% the inducing inputs are not transformed through logarithm when packed

%gp_fic = gp_set(gp_fic, 'infer_params', 'covariance+likelihood+inducing');  % optimize parameters and inducing inputs
gp_fic = gp_set(gp_fic, 'infer_params', 'covariance+likelihood'); % optimize only parameters

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the scaled conjugate gradient method
gp_fic=gp_optim(gp_fic,x,y,'opt',opt);

% Make the prediction
[Eft_fic, Varft_fic, lpyt_fic, Eyt_fic, Varyt_fic] = gp_pred(gp_fic, x, y, xt, 'yt', ones(size(xt)));

% Plot the solution of FIC
figure
%subplot(4,1,1)
hold on
plot(x,y,'.', 'MarkerSize',7)
plot(xt,Eyt_fic,'k', 'LineWidth', 2)
plot(xt,Eyt_fic-2.*sqrt(Varyt_fic),'g--', 'LineWidth', 2)
plot(gp_fic.X_u, -30, 'rx', 'MarkerSize', 5, 'LineWidth', 2)
plot(xt,Eyt_fic+2.*sqrt(Varyt_fic),'g--', 'LineWidth', 2)
axis tight
caption2 = sprintf('FIC:  l_1= %.2f, s^2_1 = %.2f, \n l_2= %.2f, s^2_2 = %.2f \n s^2_{noise} = %.2f', gp_fic.cf{1}.lengthScale, gp_fic.cf{1}.magnSigma2, gp_fic.cf{2}.lengthScale, gp_fic.cf{2}.magnSigma2, gp_fic.lik.sigma2);
title(caption2)
legend('Data point', 'predicted mean', '2\sigma error', 'inducing input','Location','Northwest')

%========================================================
% PART 3 data analysis with PIC approximation
%========================================================

% set the data points into clusters
edges = linspace(-1,max(x)+1,20);
tot=0;
for i=1:length(edges)-1
    trindex{i} = find(x>edges(i) & x<edges(i+1));
end
% Create the FIC GP structure
gp_pic = gp_set('type', 'PIC', 'lik', lik, 'cf', {gpcf1, gpcf2}, 'jitterSigma2', 1e-6, 'X_u', Xu)
gp_pic = gp_set(gp_pic, 'tr_index', trindex);

% -----------------------------
% --- Conduct the inference ---

% --- MAP estimate using modified Newton algorithm ---

% Now you can choose, if you want to optimize only parameters or
% optimize simultaneously parameters and inducing inputs. Note that
% the inducing inputs are not transformed through logarithm when packed

%gp_pic = gp_set(gp_pic, 'infer_params', 'covariance+likelihood+inducing');  % optimize parameters and inducing inputs
gp_pic = gp_set(gp_pic, 'infer_params', 'covariance+likelihood');           % optimize only parameters

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the scaled conjugate gradient method
gp_pic=gp_optim(gp_pic,x,y,'opt',opt);
gp_pic=gp_optim(gp_pic,x,y,'opt',opt,'optimf',@fminscg);

% Make the prediction
[Eft_pic, Varft_pic, lpyt_pic, Eyt_pic, Varyt_pic] = gp_pred(gp_pic, x, y, x, 'tstind', trindex, 'yt', y);

% Plot the solution of PIC
figure
%subplot(4,1,1)
hold on
plot(x,y,'.', 'MarkerSize',7)
plot(x,Eft_pic,'k', 'LineWidth', 2)
plot(x,Eft_pic-2.*sqrt(Varyt_pic),'g--', 'LineWidth', 2)
plot(gp_pic.X_u, -30, 'rx', 'MarkerSize', 5, 'LineWidth', 2)
plot(x,Eft_pic+2.*sqrt(Varyt_pic),'g--', 'LineWidth', 2)
for i = 1:length(edges)
    plot([edges(i) edges(i)],[-30 35], 'k:')
end
axis tight
caption2 = sprintf('PIC:  l_1= %.2f, s^2_1 = %.2f, \n l_2= %.2f, s^2_2 = %.2f \n s^2_{noise} = %.2f', gp_pic.cf{1}.lengthScale, gp_pic.cf{1}.magnSigma2, gp_pic.cf{2}.lengthScale, gp_pic.cf{2}.magnSigma2, gp_pic.lik.sigma2);
title(caption2)
legend('Data point', 'predicted mean', '2\sigma error', 'inducing input','Location','Northwest')

%========================================================
% PART 4 data analysis with CS+FIC model
%========================================================

% Here we conduct the same analysis as in part 1, but this time we
% use CS+FIC approximation

% Create the CS+FIC GP structure
if ~exist('ldlchol')
  error('GPstuff:SuiteSparseMissing',...
        ['SuiteSparse is not properly installed. (in BECS try ''use suitesparse'')\n' ...
         'Can not use CS+FIC without SuiteSparse']);
end
gp_csfic = gp_set('type','CS+FIC', 'lik', lik, 'cf', {gpcf1, gpcf2}, 'jitterSigma2', 1e-9, 'X_u', Xu)

% -----------------------------
% --- Conduct the inference ---

% --- MAP estimate using modified Newton algorithm ---

% Now you can choose, if you want to optimize only parameters or
% optimize simultaneously parameters and inducing inputs. Note that
% the inducing inputs are not transformed through logarithm when
% packed

% optimize parameters and inducing inputs
%gp_csfic = gp_set(gp_csfic, 'infer_params', 'covariance+likelihood+inducing');
% optimize only parameters (default)
%gp_csfic = gp_set(gp_csfic, 'infer_params', 'covariance+likelihood');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the scaled conjugate gradient method
gp_csfic=gp_optim(gp_csfic,x,y,'opt',opt);

% Make the prediction
[Eft_csfic, Varft_csfic, lpyt_csfic, Eyt_csfic, Varyt_csfic] = gp_pred(gp_csfic, x, y, x, 'yt', y);

% Plot the solution of FIC
figure
%subplot(4,1,1)
hold on
plot(x,y,'.', 'MarkerSize',7)
plot(x,Eft_csfic,'k', 'LineWidth', 2)
plot(x,Eft_csfic-2.*sqrt(Varyt_csfic),'g--', 'LineWidth', 1)
plot(gp_csfic.X_u, -30, 'rx', 'MarkerSize', 5, 'LineWidth', 2)
plot(x,Eft_csfic+2.*sqrt(Varyt_csfic),'g--', 'LineWidth', 1)
axis tight
caption2 = sprintf('CS+FIC:  l_1= %.2f, s^2_1 = %.2f, \n l_2= %.2f, s^2_2 = %.2f \n s^2_{noise} = %.2f', gp_csfic.cf{1}.lengthScale, gp_csfic.cf{1}.magnSigma2, gp_csfic.cf{2}.lengthScale, gp_csfic.cf{2}.magnSigma2, gp_csfic.lik.sigma2);
title(caption2)
legend('Data point', 'predicted mean', '2\sigma error', 'inducing input','Location','Northwest')

[Eft, Varft, lpyt, Eyt, Varyt] = gp_pred(gp_csfic, x, y, x, 'yt', y);
[Eft1, Varft1] = gp_pred(gp_csfic, x, y, x, 'predcf', 1);
[Eft2, Varft2] = gp_pred(gp_csfic, x, y, x, 'predcf', 2);

figure
set(gcf,'units','centimeters');
set(gcf,'DefaultAxesPosition',[0.08  0.13   0.84   0.85]);
set(gcf,'DefaultAxesFontSize',16)   %6 8
set(gcf,'DefaultTextFontSize',16)   %6 8
hold on
[AX, H1, H2] = plotyy(x, Eft2, x, Eft1+avgy);
set(H2,'LineStyle','--')
set(H2, 'LineWidth', 3)
set(H1,'LineStyle','-')
set(H1, 'LineWidth', 1)

set(AX(2), 'XLim', [-1 559])
set(AX(1), 'XLim', [-1 559])
set(AX(2), 'YLim', [310 380])
set(AX(1), 'YLim', [-5 5])
set(AX(2), 'XTick' ,[0 276 557])
set(AX(2), 'XTicklabel' ,[1958 1981 2004])
set(AX(1), 'XTick' ,[0 276 557])
set(AX(1), 'XTicklabel' ,[1958 1981 2004])
set(AX(2),'YTick',[310 350 380])
set(AX(2),'YTicklabel',[310 350 380])
set(AX(1),'YTick',[-5 0 5])
set(AX(1),'YTicklabel',[-5 0 5])
%set(get(AX(2),'Ylabel'),'String','ppmv')
%set(get(AX(1),'Ylabel'),'String','ppmv')
set(get(AX(2),'Xlabel'),'String','year')
set(get(AX(1),'Xlabel'),'String','year')

set(gcf,'pos',[5    3   18  10.7])
set(gcf,'paperunits',get(gcf,'units'))
set(gcf,'paperpos',get(gcf,'pos'))

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cnn_epoch_evaluate(work_info, opts, imdb, net_config, work_info_epoch)

if ~work_info.ref.run_eva
    return;
end

if opts.use_dummy_gt
    return;
end

eva_param=opts.eva_param;
class_info=eva_param.class_info;

eva_result_epoch=work_info_epoch.ref.eva_result;
valid_task_flags=work_info_epoch.ref.valid_task_flags;

seg_eva_result_imgs_coarse=work_info_epoch.ref.seg_eva_result_imgs_coarse;
if ~isempty(valid_task_flags)
    seg_eva_result_imgs_coarse=seg_eva_result_imgs_coarse(valid_task_flags);
end
assert(~isempty(seg_eva_result_imgs_coarse));
seg_eva_result_ds_coarse=gen_predict_result(seg_eva_result_imgs_coarse, class_info);

eva_result_epoch.seg_accuracy_global_coarse=seg_eva_result_ds_coarse.accuracy_global;
eva_result_epoch.seg_accuracy_per_class_coarse=seg_eva_result_ds_coarse.accuracy_per_class;
eva_result_epoch.seg_inter_union_score_per_class_coarse=seg_eva_result_ds_coarse.inter_union_score_per_class;
eva_result_epoch.exclude_class_idxes_coarse=seg_eva_result_ds_coarse.exclude_class_idxes;

seg_eva_result_imgs=work_info_epoch.ref.seg_eva_result_imgs;
if ~isempty(valid_task_flags)
    seg_eva_result_imgs=seg_eva_result_imgs(valid_task_flags);
end
assert(~isempty(seg_eva_result_imgs));
seg_eva_result_ds=gen_predict_result(seg_eva_result_imgs, class_info);

eva_result_epoch.seg_accuracy_global=seg_eva_result_ds.accuracy_global;
eva_result_epoch.seg_accuracy_per_class=seg_eva_result_ds.accuracy_per_class;
eva_result_epoch.seg_inter_union_score_per_class=seg_eva_result_ds.inter_union_score_per_class;
eva_result_epoch.exclude_class_idxes=seg_eva_result_ds.exclude_class_idxes;

seg_eva_result_ds_densecrf=[];
if eva_param.eva_densecrf_postprocess

    seg_eva_result_imgs_densecrf=work_info_epoch.ref.seg_eva_result_imgs_densecrf;
    if ~isempty(valid_task_flags)
        seg_eva_result_imgs_densecrf=seg_eva_result_imgs_densecrf(valid_task_flags);
    end
    assert(~isempty(seg_eva_result_imgs_densecrf));
    seg_eva_result_ds_densecrf=gen_predict_result(seg_eva_result_imgs_densecrf, class_info);

    eva_result_epoch.seg_accuracy_global_densecrf=seg_eva_result_ds_densecrf.accuracy_global;
    eva_result_epoch.seg_accuracy_per_class_densecrf=seg_eva_result_ds_densecrf.accuracy_per_class;
    eva_result_epoch.seg_inter_union_score_per_class_densecrf=seg_eva_result_ds_densecrf.inter_union_score_per_class;
    eva_result_epoch.exclude_class_idxes_densecrf=seg_eva_result_ds_densecrf.exclude_class_idxes;
end

eva_names=cell(0, 1);
eva_names_disp=cell(0, 1);

% eva_names{end+1, 1}='seg_accuracy_global_coarse';
% eva_names_disp{end+1, 1}='seg accuracy global (coarse)';

% eva_names{end+1, 1}='seg_accuracy_per_class_coarse';
% eva_names_disp{end+1, 1}='seg accuracy perclass (coarse)';

eva_names{end+1, 1}='seg_inter_union_score_per_class_coarse';
eva_names_disp{end+1, 1}='IoU (coarse)';

eva_names{end+1, 1}='seg_accuracy_global';
eva_names_disp{end+1, 1}='accuracy global';

eva_names{end+1, 1}='seg_accuracy_per_class';
eva_names_disp{end+1, 1}='accuracy perclass';

eva_names{end+1, 1}='seg_inter_union_score_per_class';
eva_names_disp{end+1, 1}='IoU';

eva_result_epoch.eva_names=cat(1, eva_result_epoch.eva_names, eva_names);
eva_result_epoch.eva_names_disp=cat(1, eva_result_epoch.eva_names_disp, eva_names_disp);

work_info_epoch.ref.eva_result=eva_result_epoch;

fprintf('\n\n------------------------------------------------------------\n');

fprintf('\n\nevaluate result (coarse):\n');
disp(seg_eva_result_ds_coarse);

fprintf('\n\n------------------------------------------------------------\n');
disp_eva_result(seg_eva_result_ds, class_info)

fprintf('\n\nevaluate result:\n');
disp(seg_eva_result_ds);

if ~isempty(seg_eva_result_ds_densecrf)
    fprintf('\n\n------------------------------------------------------------\n');
    disp_eva_result(seg_eva_result_ds_densecrf, class_info)

    fprintf('\n\nevaluate result (densecrf):\n');
    disp(seg_eva_result_ds_densecrf);

end

end

function disp_eva_result(pixel_result, class_info)

fprintf('\n\nCategory result:\n');
class_num=class_info.class_num;
class_label_values=class_info.class_label_values;

for c_idx=1:class_num

	one_label_value=class_label_values(c_idx);
    class_name=class_info.class_names{c_idx};

    one_iou=pixel_result.inter_union_score_classes(c_idx);
    one_acc=pixel_result.accuracy_classes(c_idx);

    fprintf('class_idx:%d, label_vale:%d, class_name:%s, pixel_acc:%.4f, IoU:%.4f\n', c_idx, one_label_value, class_name, one_acc, one_iou);
end

end

function pixel_result=gen_predict_result(pixel_results, class_info)

    exclude_class_idxes=class_info.void_class_idxes;
    class_num=class_info.class_num;
    img_num=length(pixel_results);

    pixel_con_mat=zeros(class_num, class_num);

    for img_idx_idx=1:img_num
        one_pixel_result=pixel_results{img_idx_idx};
        pixel_con_mat=pixel_con_mat+one_pixel_result.confusion_mat;
    end

    pixel_result=seg_eva_gen_result_from_con_mat(pixel_con_mat, exclude_class_idxes);

    predict_result=[];
    predict_result.class_num=class_num;
    predict_result.pixel_result=pixel_result;
    predict_result.exclude_class_idxes=exclude_class_idxes;

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function SO3F = ctranspose(SO3F)
% overloads ctranspose
%
% Syntax
%   SO3F'
%

SO3F = conj(SO3F.');

end

author: mtex-toolbox

repo: mtex

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


text:
function [kcal] = Nm2kcal(Nm)
% Convert energy or work from newton-meters to kilocalories.
% Chad A. Greene 2012
kcal = Nm*0.00023884589663;


text:
function tensor_product_test ( )

%*****************************************************************************80
%
%% TENSOR_PRODUCT_TEST tests TENSOR_PRODUCT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TENSOR_PRODUCT_TEST:\n' );
  fprintf ( 1, '  Given a sequence of 1D quadrature rules, construct the\n' );
  fprintf ( 1, '  tensor product rule.\n' );
%
%  1D rule.
%
  n1D = { [ -1.0; +1.0 ] };
  w1D = { [ 1.0; 1.0 ] };

  [ x, w ] = tensor_product ( n1D, w1D );

  x = x';

  [ m, n ] = size ( x );

  quad_rule_print ( m, n, x, w, '  A 1D rule over [-1,+1]:' );
%
%  2D rule.
%
  n1D{2} = [ 2.0; 2.5; 3.0 ];
  w1D{2} = [ 0.25; 0.50; 0.25 ];

  [ x, w ] = tensor_product ( n1D, w1D );
  x = x';

  [ m, n ] = size ( x );

  quad_rule_print ( m, n, x, w, '  A 2D rule over [-1,+1] x [2.0,3.0]:' );
%
%  3D rule.
%
  n1D{3} = [ 10.0; 15.0 ];
  w1D{3} = [ 2.50; 2.50 ];

  [ x, w ] = tensor_product ( n1D, w1D );
  x = x';
  [ m, n ] = size ( x );

  quad_rule_print ( m, n, x, w, '  A 3D rule over [-1,+1] x [2.0,3.0] x [10.0,15.0]:' );

  return
end


text:
function obj = create(obj, varargin)
% Create an object from an empty obj structure, assigning fieldname/value
% pairs as optional arguments.
%
% :Usage:
% ::
%
%     [obj = create(obj, varargin)
%
% Used in fmri_data.m class constructor.

			  % if 'noverbose' is entered, suppress output
  verbose = isempty(strmatch('noverbose', varargin(cellfun(@ischar, varargin))));

  N = fieldnames(obj);

  for i = 1:length(varargin)
    if ischar(varargin{i})

      % Look for a field (attribute) with the input name
      wh = strmatch(varargin{i}, N, 'exact');

      if ~isempty(wh)

	obj.(varargin{i}) = varargin{i + 1};

	% special methods for specific fields
	switch varargin{i}
          case 'dat'
            xx = isnan(obj.(varargin{i}));
            if any(xx(:))
              if verbose
		fprintf('fmri_data.create: Converting %3.0f NaNs to 0s.', sum(xx(:)));
	      end
              obj.dat(xx) = 0;
            end
	end %Switch
      end % ~isempty
    end
  end
end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef TestImageTemplates < CoreTest
    % TestImageTemplates. Tests for the MimImageTemplates class.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods
        function obj = TestImageTemplates

            mock_reporting = CoreMockReporting;
            mock_datset_disk_cache = MockDatasetDiskCache;
            mock_dataset_results = MockDatasetResults;
            context_def = PTKContextDef;
            image_templates = MimImageTemplates(PTKFrameworkAppDef, mock_dataset_results, context_def, mock_datset_disk_cache, mock_reporting);
            mock_dataset_results.ImageTemplates = image_templates;
            null_dataset_stack = [];

            % By default, all contexts should be enabled
            obj.Assert(image_templates.IsContextEnabled(PTKContext.OriginalImage), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LungROI), 'ROI context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightMiddleLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLowerLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftLowerLobe), 'Context is enabled');

            % Check that code fails where template could not be generated
            % When asking for a template, the code should fail if the template
            % could not be generated
            mock_dataset_results.AddMockResult('PTKOriginalImage', PTKContext.OriginalImage, [], [], [], true);
            mock_reporting.AddExpectation('CoreMockReporting.Error', 'MimImageTemplates:NoContext');
            try
                template_image = image_templates.GetTemplateImage(PTKContext.OriginalImage, null_dataset_stack, mock_reporting);
            catch ex
                if ~isa(ex, 'CoreTestException')
                    CoreErrorUtilities.ThrowException('TestImageTemplates:WrongException', 'Test failure: test exception expected');

                end
            end

            % Due to the failure, the original image context should now be
            % disabled
            obj.Assert(~image_templates.IsContextEnabled(PTKContext.OriginalImage), 'Context is disabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LungROI), 'ROI context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightMiddleLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLowerLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftLowerLobe), 'Context is enabled');

            % Now force a failure for two more contexts (LungROI and LeftLung,
            % leaving RightLung enabled)

            mock_dataset_results.AddMockResult('PTKLungROI', PTKContext.LungROI, [], [], [], true);
            mock_reporting.AddExpectation('CoreMockReporting.Error', 'MimImageTemplates:NoContext');
            try
                template_image = image_templates.GetTemplateImage(PTKContext.LungROI, null_dataset_stack, mock_reporting);
            catch ex
                if ~isa(ex, 'CoreTestException')
                    CoreErrorUtilities.ThrowException('TestImageTemplates:WrongException', 'Test failure: test exception expected');

                end
            end
            mock_dataset_results.AddMockResult('PTKGetContextForSingleLung', PTKContext.LeftLung, [], [], [], true);
            mock_reporting.AddExpectation('CoreMockReporting.Error', 'MimImageTemplates:NoContext');
            try
                template_image = image_templates.GetTemplateImage(PTKContext.LeftLung, null_dataset_stack, mock_reporting);
            catch ex
                if ~isa(ex, 'CoreTestException')
                    CoreErrorUtilities.ThrowException('TestImageTemplates:WrongException', 'Test failure: test exception expected');

                end
            end

            % Now three of the contexts should be disbaled
            obj.Assert(~image_templates.IsContextEnabled(PTKContext.OriginalImage), 'Context is disabled');
            obj.Assert(~image_templates.IsContextEnabled(PTKContext.LungROI), 'ROI context is enabled');
            obj.Assert(~image_templates.IsContextEnabled(PTKContext.LeftLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightMiddleLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLowerLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftLowerLobe), 'Context is enabled');

            % Check normal accessing of templates

            % Add template image results
            mock_roi_image = PTKImage;
            mock_roi_image.Title = 'roi';
            mock_full_image = PTKImage;
            mock_full_image.Title = 'full';
            mock_left_image = PTKImage;
            mock_left_image.Title = 'left';
            mock_right_image = PTKImage;
            mock_right_image.Title = 'right';
            mock_ru_image = PTKImage;
            mock_ru_image.Title = 'ru';
            mock_rm_image = PTKImage;
            mock_rm_image.Title = 'rm';
            mock_rl_image = PTKImage;
            mock_rl_image.Title = 'rl';
            mock_lu_image = PTKImage;
            mock_lu_image.Title = 'lu';
            mock_ll_image = PTKImage;
            mock_ll_image.Title = 'll';

            lung_roi_cache_item = MimDatasetStackItem(MimDependency('PTKLungROI', PTKContext.LungROI, '1.2.3.4.5', '1', []), MimDependencyList, false, false, mock_reporting);
            oi_cache_item = MimDatasetStackItem(MimDependency('PTKOriginalImage', PTKContext.OriginalImage, '1.2.3.4.6', '1', []), MimDependencyList, false, false, mock_reporting);
            ll_cache_item = MimDatasetStackItem(MimDependency('PTKGetContextForSingleLung', PTKContext.LeftLung, '1.2.3.4.7', '1', []), MimDependencyList, false, false, mock_reporting);
            rl_cache_item = MimDatasetStackItem(MimDependency('PTKGetContextForSingleLung', PTKContext.RightLung, '1.2.3.4.8', '1', []), MimDependencyList, false, false, mock_reporting);
            rul_cache_item = MimDatasetStackItem(MimDependency('PTKGetContextForLobe', PTKContext.RightUpperLobe, '1.2.3.4.9', '1', []), MimDependencyList, false, false, mock_reporting);
            rml_cache_item = MimDatasetStackItem(MimDependency('PTKGetContextForLobe', PTKContext.RightMiddleLobe, '1.2.3.4.10', '1', []), MimDependencyList, false, false, mock_reporting);
            rll_cache_item = MimDatasetStackItem(MimDependency('PTKGetContextForLobe', PTKContext.RightLowerLobe, '1.2.3.4.11', '1', []), MimDependencyList, false, false, mock_reporting);
            lul_cache_item = MimDatasetStackItem(MimDependency('PTKGetContextForLobe', PTKContext.LeftUpperLobe, '1.2.3.4.12', '1', []), MimDependencyList, false, false, mock_reporting);
            lll_cache_item = MimDatasetStackItem(MimDependency('PTKGetContextForLobe', PTKContext.LeftLowerLobe, '1.2.3.4.13', '1', []), MimDependencyList, false, false, mock_reporting);

            % Check fetching the template images
            mock_dataset_results.AddMockResult('PTKLungROI', PTKContext.LungROI, mock_roi_image, lung_roi_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKOriginalImage', PTKContext.OriginalImage, mock_full_image, oi_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKGetContextForSingleLung', PTKContext.LeftLung, mock_left_image, ll_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKGetContextForSingleLung', PTKContext.RightLung, mock_right_image, rl_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKGetContextForLobe', PTKContext.RightUpperLobe, mock_ru_image, rul_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKGetContextForLobe', PTKContext.RightMiddleLobe, mock_rm_image, rml_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKGetContextForLobe', PTKContext.RightLowerLobe, mock_rl_image, rll_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKGetContextForLobe', PTKContext.LeftUpperLobe, mock_lu_image, lul_cache_item, [], true);
            mock_dataset_results.AddMockResult('PTKGetContextForLobe', PTKContext.LeftLowerLobe, mock_ll_image, lll_cache_item, [], true);

            roi_template_image = image_templates.GetTemplateImage(PTKContext.LungROI, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(roi_template_image.Title, 'roi'), 'Correct template image returned');

            original_template_image = image_templates.GetTemplateImage(PTKContext.OriginalImage, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(original_template_image.Title, 'full'), 'Correct template image returned');

            left_template_image = image_templates.GetTemplateImage(PTKContext.LeftLung, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(left_template_image.Title, 'left'), 'Correct template image returned');

            right_template_image = image_templates.GetTemplateImage(PTKContext.RightLung, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(right_template_image.Title, 'right'), 'Correct template image returned');

            ru_template_image = image_templates.GetTemplateImage(PTKContext.RightUpperLobe, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(ru_template_image.Title, 'ru'), 'Correct template image returned');

            rm_template_image = image_templates.GetTemplateImage(PTKContext.RightMiddleLobe, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(rm_template_image.Title, 'rm'), 'Correct template image returned');

            rl_template_image = image_templates.GetTemplateImage(PTKContext.RightLowerLobe, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(rl_template_image.Title, 'rl'), 'Correct template image returned');

            lu_template_image = image_templates.GetTemplateImage(PTKContext.LeftUpperLobe, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(lu_template_image.Title, 'lu'), 'Correct template image returned');

            ll_template_image = image_templates.GetTemplateImage(PTKContext.LeftLowerLobe, null_dataset_stack, mock_reporting);
            obj.Assert(strcmp(ll_template_image.Title, 'll'), 'Correct template image returned');

            % Contexts should now all be enabled again
            obj.Assert(image_templates.IsContextEnabled(PTKContext.OriginalImage), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LungROI), 'ROI context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLung), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightMiddleLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.RightLowerLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftUpperLobe), 'Context is enabled');
            obj.Assert(image_templates.IsContextEnabled(PTKContext.LeftLowerLobe), 'Context is enabled');
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function exact = p13_exact ( )

%*****************************************************************************80
%
%% P13_EXACT returns the estimated integral for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = pi / 2.0;

  return
end


text:
function x = discrete_sample(p, n)
% Samples from a discrete distribution
%
%   x = discretesample(p, n)
%       independently draws n samples (with replacement) from the
%       distribution specified by p, where p is a probability array
%       whose elements sum to 1.
%
%       Suppose the sample space comprises K distinct objects, then
%       p should be an array with K elements. In the output, x(i) = k
%       means that the k-th object is drawn at the i-th trial.
%
%   Remarks
%   -------
%       - This function is mainly for efficient sampling in non-uniform
%         distribution, which can be either parametric or non-parametric.
%
%       - The function is implemented based on histc, which has been
%         highly optimized by mathworks. The basic idea is to divide
%         the range [0, 1] into K bins, with the length of each bin
%         proportional to the probability mass. And then, n values are
%         drawn from a uniform distribution in [0, 1], and the bins that
%         these values fall into are picked as results.
%
%       - This function can also be employed for continuous distribution
%         in 1D/2D dimensional space, where the distribution can be
%         effectively discretized.
%
%       - This function can also be useful for sampling from distributions
%         which can be considered as weighted sum of "modes".
%         In this type of applications, you can first randomly choose
%         a mode, and then sample from that mode. The process of choosing
%         a mode according to the weights can be accomplished with this
%         function.
%
%   Examples
%   --------
%       % sample from a uniform distribution for K objects.
%       p = ones(1, K) / K;
%       x = discretesample(p, n);
%
%       % sample from a non-uniform distribution given by user
%       x = discretesample([0.6 0.3 0.1], n);
%
%       % sample from a parametric discrete distribution with
%       % probability mass function given by f.
%       p = f(1:K);
%       x = discretesample(p, n);
%

%   Created by Dahua Lin, On Oct 27, 2008
%

%% parse and verify input arguments

assert(isfloat(p), 'discretesample:invalidarg', ...
    'p should be an array with floating-point value type.');

assert(isnumeric(n) && isscalar(n) && n >= 0 && n == fix(n), ...
    'discretesample:invalidarg', ...
    'n should be a nonnegative integer scalar.');

%% main

% process p if necessary

K = numel(p);
if ~isequal(size(p), [1, K])
    p = reshape(p, [1, K]);
end

% construct the bins

edges = [0, cumsum(p)];
s = edges(end);
if abs(s - 1) > eps
    edges = edges * (1 / s);
end

% draw bins

rv = rand(1, n);
c = histc(rv, edges);
ce = c(end);
c = c(1:end-1);
c(end) = c(end) + ce;

% extract samples

xv = find(c);

if numel(xv) == n  % each value is sampled at most once
    x = xv;
else                % some values are sampled more than once
    xc = c(xv);
    d = zeros(1, n);
    dv = [xv(1), diff(xv)];
    dp = [1, 1 + cumsum(xc(1:end-1))];
    d(dp) = dv;
    x = cumsum(d);
end

% randomly permute the sample's order
x = x(randperm(n));

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS