text:
function [post nlZ dnlZ] = infLaplace(hyp, mean, cov, lik, x, y, opt)

% Laplace approximation to the posterior Gaussian process.
% The function takes a specified covariance function (see covFunctions.m) and
% likelihood function (see likFunctions.m), and is designed to be used with
% gp.m. See also infMethods.m.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch 2013-05-02.
%
% See also INFMETHODS.M.

persistent last_alpha                                   % copy of the last alpha
if any(isnan(last_alpha)), last_alpha = zeros(size(last_alpha)); end   % prevent

if nargin<=6, opt = []; end                        % make opt variable available
if isfield(opt,'postL'), postL = opt.postL;        % recompute matrix L for post
else postL = true; end                                           % default value

inf = 'infLaplace';
n = size(x,1);
if isnumeric(cov),  K = cov;                    % use provided covariance matrix
else K = feval(cov{:},  hyp.cov,  x); end       % evaluate the covariance matrix
if isnumeric(mean), m = mean;                         % use provided mean vector
else m = feval(mean{:}, hyp.mean, x); end             % evaluate the mean vector
likfun = @(f) feval(lik{:},hyp.lik,y,f,[],inf);        % log likelihood function

if any(size(last_alpha)~=[n,1])     % find a good starting point for alpha and f
  alpha = zeros(n,1);                      % start at mean if sizes do not match
else
  alpha = last_alpha;                                             % try last one
  if Psi(alpha,m,K,likfun) > -sum(likfun(m))     % default f==m better => use it
    alpha = zeros(n,1);
  end
end

% switch between optimisation methods
alpha = irls(alpha, m,K,likfun, opt);                         % run optimisation

f = K*alpha+m;                                  % compute latent function values
last_alpha = alpha;                                     % remember for next call
[lp,dlp,d2lp,d3lp] = likfun(f); W = -d2lp; isWneg = any(W<0);
post.alpha = alpha;                            % return the posterior parameters
post.sW = sqrt(abs(W)).*sign(W);             % preserve sign in case of negative

% diagnose optimality
err = @(x,y) norm(x-y)/max([norm(x),norm(y),1]);   % we need to have alpha = dlp
% dev = err(alpha,dlp);  if dev>1e-4, warning('Not at optimum %1.2e.',dev), end

if postL || nargout>1
  if isWneg                  % switch between Cholesky and LU decomposition mode
    % For post.L = -inv(K+diag(1./W)), we us the non-default parametrisation.
    [ldA, iA, post.L] = logdetA(K,W);   % A=eye(n)+K*W is as safe as symmetric B
    nlZ = alpha'*(f-m)/2 - sum(lp) + ldA/2;
  else
    sW = post.sW; post.L = chol(eye(n)+sW*sW'.*K);                   % recompute
    nlZ = alpha'*(f-m)/2 + sum(log(diag(post.L))-lp);   % ..(f-m)/2 -lp +ln|B|/2
  end
end

if nargout>2                                           % do we want derivatives?
  dnlZ = hyp;                                   % allocate space for derivatives
  if isWneg                  % switch between Cholesky and LU decomposition mode
    Z = -post.L;                                                 % inv(K+inv(W))
    g = sum(iA.*K,2)/2; % deriv. of ln|B| wrt W; g = diag(inv(inv(K)+diag(W)))/2
  else
    Z = repmat(sW,1,n).*solve_chol(post.L,diag(sW)); %sW*inv(B)*sW=inv(K+inv(W))
    C = post.L'\(repmat(sW,1,n).*K);                     % deriv. of ln|B| wrt W
    g = (diag(K)-sum(C.^2,1)')/2;                    % g = diag(inv(inv(K)+W))/2
  end
  dfhat = g.*d3lp;  % deriv. of nlZ wrt. fhat: dfhat=diag(inv(inv(K)+W)).*d3lp/2
  for i=1:length(hyp.cov)                                    % covariance hypers
    dK = feval(cov{:}, hyp.cov, x, [], i);
    dnlZ.cov(i) = sum(sum(Z.*dK))/2 - alpha'*dK*alpha/2;         % explicit part
    b = dK*dlp;                            % b-K*(Z*b) = inv(eye(n)+K*diag(W))*b
    dnlZ.cov(i) = dnlZ.cov(i) - dfhat'*( b-K*(Z*b) );            % implicit part
  end
  for i=1:length(hyp.lik)                                    % likelihood hypers
    [lp_dhyp,dlp_dhyp,d2lp_dhyp] = feval(lik{:},hyp.lik,y,f,[],inf,i);
    dnlZ.lik(i) = -g'*d2lp_dhyp - sum(lp_dhyp);                  % explicit part
    b = K*dlp_dhyp;                        % b-K*(Z*b) = inv(eye(n)+K*diag(W))*b
    dnlZ.lik(i) = dnlZ.lik(i) - dfhat'*( b-K*(Z*b) );            % implicit part
  end
  for i=1:length(hyp.mean)                                         % mean hypers
    dm = feval(mean{:}, hyp.mean, x, i);
    dnlZ.mean(i) = -alpha'*dm;                                   % explicit part
    dnlZ.mean(i) = dnlZ.mean(i) - dfhat'*(dm-K*(Z*dm));          % implicit part
  end
end

% Evaluate criterion Psi(alpha) = alpha'*K*alpha + likfun(f), where
% f = K*alpha+m, and likfun(f) = feval(lik{:},hyp.lik,y,  f,  [],inf).
function [psi,dpsi,f,alpha,dlp,W] = Psi(alpha,m,K,likfun)
  f = K*alpha+m;
  [lp,dlp,d2lp] = likfun(f); W = -d2lp;
  psi = alpha'*(f-m)/2 - sum(lp);
  if nargout>1, dpsi = K*(alpha-dlp); end

% Run IRLS Newton algorithm to optimise Psi(alpha).
function alpha = irls(alpha, m,K,likfun, opt)
  if isfield(opt,'irls_maxit'), maxit = opt.irls_maxit; % max no of Newton steps
  else maxit = 20; end                                           % default value
  if isfield(opt,'irls_Wmin'),  Wmin = opt.irls_Wmin; % min likelihood curvature
  else Wmin = 0.0; end                                           % default value
  if isfield(opt,'irls_tol'),   tol = opt.irls_tol;     % stop Newton iterations
  else tol = 1e-6; end                                           % default value

  smin_line = 0; smax_line = 2;           % min/max line search steps size range
  nmax_line = 10;                          % maximum number of line search steps
  thr_line = 1e-4;                                       % line search threshold
  Psi_line = @(s,alpha,dalpha) Psi(alpha+s*dalpha, m,K,likfun);    % line search
  pars_line = {smin_line,smax_line,nmax_line,thr_line};  % line seach parameters
  search_line = @(alpha,dalpha) brentmin(pars_line{:},Psi_line,5,alpha,dalpha);

  f = K*alpha+m; [lp,dlp,d2lp] = likfun(f); W = -d2lp; n = size(K,1);
  Psi_new = Psi(alpha,m,K,likfun);
  Psi_old = Inf;  % make sure while loop starts by the largest old objective val
  it = 0;                          % this happens for the Student's t likelihood
  while Psi_old - Psi_new > tol && it<maxit                       % begin Newton
    Psi_old = Psi_new; it = it+1;
    % limit stepsize
    W = max(W,Wmin); % reduce step size by increasing curvature of problematic W
    sW = sqrt(W); L = chol(eye(n)+sW*sW'.*K);            % L'*L=B=eye(n)+sW*K*sW
    b = W.*(f-m) + dlp;
    dalpha = b - sW.*solve_chol(L,sW.*(K*b)) - alpha; % Newton dir + line search
    [s_line,Psi_new,n_line,dPsi_new,f,alpha,dlp,W] = search_line(alpha,dalpha);
  end                                                  % end Newton's iterations

% Compute the log determinant ldA and the inverse iA of a square nxn matrix
% A = eye(n) + K*diag(w) from its LU decomposition; for negative definite A, we
% return ldA = Inf. We also return mwiA = -diag(w)/A.
function [ldA,iA,mwiA] = logdetA(K,w)
  [m,n] = size(K); if m~=n, error('K has to be nxn'), end
  A = eye(n)+K.*repmat(w',n,1);
  [L,U,P] = lu(A); u = diag(U);           % compute LU decomposition, A = P'*L*U
  signU = prod(sign(u));                                             % sign of U
  detP = 1;                 % compute sign (and det) of the permutation matrix P
  p = P*(1:n)';
  for i=1:n                                                       % swap entries
    if i~=p(i), detP = -detP; j = find(p==i); p([i,j]) = p([j,i]); end
  end
  if signU~=detP  % log becomes complex for negative values, encoded by infinity
    ldA = Inf;
  else            % det(L) = 1 and U triangular => det(A) = det(P)*prod(diag(U))
    ldA = sum(log(abs(u)));
  end
  if nargout>1, iA = U\(L\P); end               % return the inverse if required
  if nargout>2, mwiA = -repmat(w,1,n).*iA; end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function polar_to_c8_test ( )

%*****************************************************************************80
%
%% POLAR_TO_C8_TEST tests POLAR_TO_C8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2015
%
%  Author:
%
%    John Burkardt
%
  test_num = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLAR_TO_C8 converts polar form to C8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      (R1,T1)=R8_UNIFORM_01()    C2=POLAR_TO_C8(R1,T1)' );
  fprintf ( 1, '    (R3,T3)=C8_TO_POLAR(C2)\n' );
  fprintf ( 1, '       ----------------------       ----------------------' );
  fprintf ( 1, '      -----------------------\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : test_num

    [ r1, seed ] = r8_uniform_01 ( seed );
    [ t1, seed ] = r8_uniform_01 ( seed );
    t1 = t1 * 2.0 * pi;
    c2 = polar_to_c8 ( r1, t1 );;
    [ r3, t3 ] = c8_to_polar ( c2 );

    fprintf ( 1, '  %12f  %12f  (%12f  %12f)  %12f  %12f\n', ...
      r1, t1, real ( c2 ), imag ( c2 ), r3, t3 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = logGauss(X, mu, sigma)
% Compute log pdf of a Gaussian distribution.
% Input:
%   X: d x n data matrix
%   mu: d x 1 mean vector of Gaussian
%   sigma: d x d covariance matrix of Gaussian
% Output:
%   y: 1 x n probability density in logrithm scale y=log p(x)
% Written by Mo Chen (sth4nth@gmail.com).
d = size(X,1);
X = X-mu;
[U,p]= chol(sigma);
if p ~= 0
    error('ERROR: sigma is not PD.');
end
Q = U'\X;
q = dot(Q,Q,1);  % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(U)));   % normalization constant
y = -(c+q)/2;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function basis_mn_t6_test ( )

%*****************************************************************************80
%
%% BASIS_MN_T6_TEST verifies BASIS_MN_T6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  node_num = 6;

  t = [ ...
    2.0, 0.0; ...
    4.0, 3.0; ...
    0.0, 4.0; ...
    3.0, 1.5; ...
    2.0, 3.5; ...
    1.0, 2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BASIS_MN_T6_TEST\n' );
  fprintf ( 1, '  Verify basis functions for element T6.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Physical Nodes:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       I    X         Y\n' );
  fprintf ( 1, '\n' );

  for j = 1 : node_num
    fprintf ( 1, '  %6d  %7f  %7f\n', j, t(1:2,j) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The basis function values at basis nodes\n' );
  fprintf ( 1, '  should form the identity matrix.\n' );
  fprintf ( 1, '\n' );

  [ phi, dphidx, dphidy ] = basis_mn_t6 ( t, node_num, t );

  for i = 1 : node_num
    for j = 1 : node_num
      fprintf ( 1, '  %7f', phi(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The X and Y derivatives should sum to 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      dPhidX sum      dPhidY sum\n' );
  fprintf ( 1, '\n' );
  for j = 1 : node_num
    sum_x = sum ( dphidx(1:node_num,j) );
    sum_y = sum ( dphidy(1:node_num,j) );
    fprintf ( 1, '  %14f  %14f\n', sum_x, sum_y );
  end

  return
end


text:
function WtW = pairwise_dists(WU, WUinit)

WU = reshape(WU, [], size(WU,3));
WUinit = reshape(WUinit, [], size(WUinit,3));
WtW = WU' * WUinit;

mu = sum(WU.^2,1);
mu = mu(:);

muinit = sum(WUinit.^2,1);
muinit = muinit(:);

mu     = repmat(mu, 1, size(WUinit,2));
muinit = repmat(muinit', size(WU,2), 1);

WtW = 1 - 2*WtW ./ (muinit + mu);

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
% Installer for the Covariance toolbox

HOME = pwd;

if isunix
    path([HOME,'/lib'],path);
    path([HOME,'/lib/distance'],path);
    path([HOME,'/lib/geodesic'],path);
    path([HOME,'/lib/riemann'],path);
    path([HOME,'/lib/visu'],path);
    path([HOME,'/lib/estimation'],path);
    path([HOME,'/lib/mean'],path);
    path([HOME,'/lib/simulation'],path);
    path([HOME,'/lib/jointdiag'],path);
    path([HOME,'/lib/classification'],path);
    path([HOME,'/lib/potato'],path);
else
    path([HOME,'\lib'],path);
    path([HOME,'\lib\distance'],path);
    path([HOME,'\lib\geodesic'],path);
    path([HOME,'\lib\riemann'],path);
    path([HOME,'\lib\visu'],path);
    path([HOME,'\lib\estimation'],path);
    path([HOME,'\lib\mean'],path);
    path([HOME,'\lib\simulation'],path);
    path([HOME,'\lib\jointdiag'],path);
    path([HOME,'\lib\classification'],path);
    path([HOME,'\lib\potato'],path);
end
disp('Covariance toolbox activated');

author: alexandrebarachant

repo: covariancetoolbox


text:
function [ tempimg ] = visualCimage( img, cimages )
%VISUALCIMAGE Summary of this function goes here
%   Detailed explanation goes here

rl=[255  255   0  0    0    255  255];
gl=[255  0   255  255  0    0    255];
bl=[0    0   255  0    255  255  255];

[aa indd]=max(cimages,[],3);

% figure(102);
clear mask_color;
mask_r = rl(indd);
mask_g = gl(indd);
mask_b = bl(indd);
mask_color(:,:,1) = mask_r;
mask_color(:,:,2) = mask_g;
mask_color(:,:,3) = mask_b;

hsvmask=rgb2hsv(mask_color);
hsvmask(:,:,3)=aa*255;
%     hsvmask(:,:,2)=aa;
mask_color=hsv2rgb(hsvmask);

tempimg = double(img)*0.5 + mask_color*0.5;
%         tempimg =  mask_color;
% imshow(uint8(tempimg));

end

author: zouchuhang

repo: LayoutNet


text:
function imsegs = processSuperpixelImage(fn)
% imsegs = processSuperpixelImage(fn)
% Creates the imsegs structure from a segmentation image
%
% INPUT:
% fn - filenames of segmentation images. Use '/' (not '\') to separate directories.
% Segments are denoted by different RGB colors.
%
% OUTPUT:
% imsegs - image segmentation data
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2006

if isstr(fn)
    fn = {fn};
end

imsegs(length(fn)) = struct('imname', '', 'imsize', [0 0]);
for f = 1:length(fn)
    toks = strtokAll(fn{f}, '/');
    imname = toks{end};
    basename = strtok(imname, '.');
    im = imread(fn{f});
    im = double(im);

    imsegs(f).imname = [strtok(imname,'.') '.jpg'];
    imsegs(f).imsize = size(im);
    imsegs(f).imsize = imsegs(f).imsize(1:2);
    im = im(:, :, 1) + im(:, :, 2)*256 + im(:, :, 3)*256^2;
    [gid, gn] = grp2idx(im(:));
    imsegs(f).segimage = uint16(reshape(gid, imsegs(f).imsize));
    imsegs(f).nseg = length(gn);
end
imsegs = APPgetSpStats(imsegs);

author: Cloud-CV

repo: object-proposals


text:
function [ point_num, face_num, face_data_num ] = xyf_example_size (  )

%*****************************************************************************80
%
%% XYF_EXAMPLE_SIZE sizes the data to be created by XYF_EXAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer POINT_NUM, the number of points.
%
%    Output, integer FACE_NUM, the number of faces.
%
%    Output, integer FACE_DATA_NUM, the number of face items.
%
  n_t = 13;
  n_r = 5;

  face_data_num = 4 * ( n_t - 1 ) * ( n_r - 1 );
  face_num = ( n_t - 1 ) * ( n_r - 1 );
  point_num = n_t * n_r;

  return
end


text:
function varargout = mrmMakeMovieGUI_Pause(varargin)
% MRMMAKEMOVIEGUI_PAUSE M-file for mrmMakeMovieGUI_Pause.fig
%      MRMMAKEMOVIEGUI_PAUSE, by itself, creates a new MRMMAKEMOVIEGUI_PAUSE or raises the existing
%      singleton*.
%
%      H = MRMMAKEMOVIEGUI_PAUSE returns the handle to a new MRMMAKEMOVIEGUI_PAUSE or the handle to
%      the existing singleton*.
%
%      MRMMAKEMOVIEGUI_PAUSE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MRMMAKEMOVIEGUI_PAUSE.M with the given input arguments.
%
%      MRMMAKEMOVIEGUI_PAUSE('Property','Value',...) creates a new MRMMAKEMOVIEGUI_PAUSE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before mrmMakeMovieGUI_Pause_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to mrmMakeMovieGUI_Pause_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help mrmMakeMovieGUI_Pause

% Last Modified by GUIDE v2.5 18-Jun-2010 20:12:09

    % Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @mrmMakeMovieGUI_Pause_OpeningFcn, ...
                       'gui_OutputFcn',  @mrmMakeMovieGUI_Pause_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end

    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % End initialization code - DO NOT EDIT
end

function mrmMakeMovieGUI_Pause_OpeningFcn(hObject, eventdata, handles, parentFigure, editData)
    % Choose default command line output for mrmMakeMovieGUI_Pause
    if (~exist('editData', 'var')), editData = []; end

    handles.output      = hObject;
    handles.parent      = parentFigure;
    handles.editData    = editData;

    % Update handles structure
    guidata(hObject, handles);
    Pause_InitFcn(hObject, handles);

end

function Pause_InitFcn(hObject, handles)
    if (~isempty(handles.editData))
        editData = handles.editData;

        set(handles.PauseAddButton, 'String', 'Update');
        set(handles.FrameCountTextField, 'String', num2str(editData.frames));
        set(handles.PauseLabelTextField, 'String', editData.label);
    end
end

function varargout = mrmMakeMovieGUI_Pause_OutputFcn(hObject, eventdata, handles)

    varargout{1} = handles.output;

end

function PauseAddButton_Callback(hObject, eventdata, handles)
    frames = get(handles.FrameCountTextField, 'String');
    if (isempty(frames))
        errordlg('Please specify a frame count.');
        return;
    end
    pause.eventType = 'pause';
    pause.frames    = str2num(frames);
    pause.label     = get(handles.PauseLabelTextField, 'String');
    set(handles.parent, 'UserData', pause);
    close(handles.figure1);

end

function mrmMakeMovieGUI_Pause_DeleteFcn(hObject, eventdata, handles)
    set(handles.parent, 'UserData', []);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [v,beta,xnorm] = hmake1 (x)
%HMAKE1 construct a Householder reflection
% Example:
%   [v,beta,xnorm] = hmake1 (x)
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

n = length (x) ;
if (n == 1)
    v = 1 ;
    xnorm = norm (x) ;
    if (x (1) < 0)
        beta = 2 ;
    else
        beta = 0 ;
    end
    return
end
sigma = x (2:n)'*x(2:n) ;
xnorm = sqrt (x (1)^2 + sigma) ;
v = x ;
if (sigma == 0)
    v (1) = 1 ;
    if (x (1) < 0)
        beta = 2 ;
    else
        beta = 0 ;
    end
else
    if (x (1) <= 0)
        v (1) = x(1) - xnorm ;
    else
        v (1) = -sigma / (x(1) + xnorm) ;
    end
    beta = (2*v(1)^2) / (sigma + v(1)^2) ;
    v = v / v(1) ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [q,N] = quantile2(X,p,dim,method)
% Quantiles of a sample via various methods.
%
%   Q = QUANTILE2(X,P) returns quantiles of the values in X. P is a scalar
%   or a vector of cumulative probability values.  When X is a vector, Q is
%   the same size as P, and Q(i) contains the P(i)-th quantile.  When X is
%   a matrix, the i-th row of Q contains the P(i)-th quantiles of each
%   column of X.  For N-D arrays, QUANTILE2 operates along the first
%   non-singleton dimension.
%
%   Q = QUANTILE2(X,P,DIM) calculates quantiles along dimension DIM.  The
%   DIM'th dimension of Q has length LENGTH(P).
%
%   Q = QUANTILE2(X,P,DIM,METHOD) calculates quantiles using one of the
%   methods described in http://en.wikipedia.org/wiki/Quantile. The method
%   are designated 'R-1'...'R-9'; the default is R-8 as described in
%   http://bit.ly/1kX4NcT, whereas Matlab uses 'R-5'.
%
%   Q = QUANTILE2(X,P,DIM,METHOD) calculates quantiles using one of the
%   methods described in http://en.wikipedia.org/wiki/Quantile. The method
%   are designated 'R-1'...'R-9'; the default is 'R-8' as described in
%   http://bit.ly/1kX4NcT, whereas Matlab uses 'R-5'.
%
%   Q = QUANTILE2(X,P,[],METHOD) uses the specified METHOD, but calculates
%   quantiles along the first non-singleton dimension.
%
%   [Q,N] = QUANTILE2(...) returns an array that is the same size as Q such
%   that N(i) is the number of points used to calculate Q(i).
%
%   Further reading
%
%   Hyndman, R.J.; Fan, Y. (November 1996). "Sample Quantiles in
%     Statistical Packages". The American Statistician 50 (4): 361-365.
%   Frigge, Michael; Hoaglin, David C.; Iglewicz, Boris (February 1989).
%     "Some Implementations of the Boxplot". The American Statistician 43
%     (1): 50-54.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LICENSE FILE:
% -----------------------------------------------------------------------
% Copyright (c) 2015, Christopher Hummersone
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% =========================================================================
% Last changed:     $Date: 2015-06-16 13:50:46 +0100 (Tue, 16 Jun 2015) $
% Last committed:   $Revision: 385 $
% Last changed by:  $Author: ch0022 $
% =========================================================================

    %% Check input and make default assignments

    assert(isnumeric(X),'X must be a numeric');
    assert(isvector(p) & isnumeric(p),'P must be a numeric vector');
    assert(all(p>=0 & p<=1),'Values in P must be in the interval [0,1].')

    if nargin<2
        error('Not enough input arguments.')
    end

    dims = size(X);
    if nargin<3 || isempty(dim)
        dim = find(dims>1,1,'first'); % default dim
    else % validate input
        assert(isnumeric(dim) | isempty(dim),'DIM must be an integer or empty');
        assert(isint(dim) | isempty(dim),'DIM must be an integer or empty');
        assert(dim>0,'DIM must be greater than 0')
    end

    if nargin<4
        method = 'r-8'; % default method
    else % validate input
        assert(ischar(method),'METHOD must be a character array')
    end

    %% choose method

    % See http://en.wikipedia.org/wiki/Quantile#Estimating_the_quantiles_of_a_population

    switch lower(method)
        case 'r-1'
            min_con = @(N,p)(p==0);
            max_con = @(N,p)(false);
            h = @(N,p)((N*p)+.5);
            Qp = @(x,h)(x(ceil(h-.5)));
        case 'r-2'
            min_con = @(N,p)(p==0);
            max_con = @(N,p)(p==1);
            h = @(N,p)((N*p)+.5);
            Qp = @(x,h)((x(ceil(h-.5))+x(floor(h+.5)))/2);
        case 'r-3'
            min_con = @(N,p)(p<=(.5/N));
            max_con = @(N,p)(false);
            h = @(N,p)(N*p);
            Qp = @(x,h)(x(round(h)));
        case 'r-4'
            min_con = @(N,p)(p<(1/N));
            max_con = @(N,p)(p==1);
            h = @(N,p)(N*p);
            Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
        case 'r-5'
            min_con = @(N,p)(p<(.5/N));
            max_con = @(N,p)(p>=((N-.5)/N));
            h = @(N,p)((N*p)+.5);
            Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
        case 'r-6'
            min_con = @(N,p)(p<(1/(N+1)));
            max_con = @(N,p)(p>=(N/(N+1)));
            h = @(N,p)((N+1)*p);
            Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
        case 'r-7'
            min_con = @(N,p)(false);
            max_con = @(N,p)(p==1);
            h = @(N,p)(((N-1)*p)+1);
            Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
        case 'r-8'
            min_con = @(N,p)(p<((2/3)/(N+(1/3))));
            max_con = @(N,p)(p>=((N-(1/3))/(N+(1/3))));
            h = @(N,p)(((N+(1/3))*p)+(1/3));
            Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
        case 'r-9'
            min_con = @(N,p)(p<((5/8)/(N+.25)));
            max_con = @(N,p)(p>=((N-(3/8))/(N+.25)));
            h = @(N,p)(((N+.25)*p)+(3/8));
            Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
        otherwise
            error(['Method ''' method ''' does not exist'])
    end

    %% calculate quartiles

    % reshape data so function works down columns
    order = mod(dim-1:dim+length(dims)-2,length(dims))+1;
    dims_shift = dims(order);
    x = rearrange(X,order,[dims_shift(1) prod(dims_shift(2:end))]);

    % pre-allocate q
    q = zeros([length(p) prod(dims_shift(2:end))]);
    N = zeros([length(p) prod(dims_shift(2:end))]);
    for m = 1:length(p)
        for n = 1:numel(q)/length(p)
            x2 = sort(x(~isnan(x(:,n)),n)); % sort
            N(m,n) = length(x2); % sample size
            switch N(m,n)
                case 0
                    q(m,n) = NaN;
                case 1
                    q(m,n) = x2;
                otherwise
                    if min_con(N(m,n),p(m)) % at lower limit
                        q(m,n) = x2(1);
                    elseif max_con(N(m,n),p(m)) % at upper limit
                        q(m,n) = x2(N(m,n));
                    else % everything else
                        q(m,n) = Qp(x2,h(N(m,n),p(m)));
                    end
            end
        end
    end

    % restore dims of q to equate to those of input
    q = irearrange(q,order,[length(p) dims_shift(2:end)]);
    N = irearrange(N,order,[length(p) dims_shift(2:end)]);

    % if q is a vector, make same shape as p
    if numel(p)==numel(q)
        q=reshape(q,size(p));
        N=reshape(N,size(p));
    end

end

function y = isint(x)
%ISINT check if input is whole number
    y = x==round(x);
end

function y = rearrange(x,order,shape)
%REARRANGE reshape and permute to make target dim column
    y = permute(x,order);
    y = reshape(y,shape);
end

function y = irearrange(x,order,shape)
%IREARRANGE reshape and permute to original size
    y = reshape(x,shape);
    y = ipermute(y,order);
end

author: brainstorm-tools


text:
function t = isupper(c)
%ISUPPER  True for uppercase letters.
%
%   For a string C, ISUPPER(C) is 1 for uppercase letters and 0
%   otherwise.
%
%   See also: ISALNUM, ISALPHA, ISASCII, ISDIGIT, ISLOWER, ISPRTCHR,
%   ISXDIGIT.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:50:48 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 1, nargin));

   t = isletter(c) & (c == upper(c));

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%kxvregister ' '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros xvregister.pane file
%
% Parameters:
% Toggle: dummy 'Dummy Flag Variable', default: 0: 'A dummy variable since the CLUI cannot be empty'
%
% Example: kxvregister( {'dummy',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% xvregister -
%
%  DESCRIPTION
%
% \fBxvregister" is used to gather information about the customer's of
% Khoros so that Khoral Research can continually provide service, support,
% and updates of Khoros and Khoros based products.  Whenever users of
% Khoros execute any interactive Khoros program, they will be prompted
% to register and license their version of Khoros.  Registration involves
% providing contact information as well as answering a few simple questions.
% After filling out the form, if the "Register" button is selected, the
% user will have the opportunity to either email the registration to Khoral
% Research or save it to a file where it could then be later submitted to
% KRI.
%
% Users are not required to register their version.  If the "Do Not Register"
% button is selected, the user will be marked as unregistered and the
% interactive application will run anyway.  Note however that as long as
% a user is unregistered, they will be prompted to register each time a
% Khoros interactive application is run.
%
% After the registration process is completed, the user will then be presented
% with the Khoros license agreement.  Similar to registration, if the license
% is not agreed to, it will be presented each time an Khoros application is
% run until the user agrees to the license.
%
% Note that xvregister can also be run from the command line.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function subset_sum_test ( )

%*****************************************************************************80
%
%% SUBSET_SUM_TEST tests the SUBSET_SUM library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 May 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUBSET_SUM_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SUBSET_SUM library.\n' );
%
%  Problem #1.
%
  w = [ 15, 22, 14, 26, 32, 9, 16, 8 ];
  t = 53;
  [ c, index ] = subset_sum_test01 ( w, t );
%
%  There's a second solution for problem #1.
%
  subset_sum_test01 ( w, t, [ index+1, 2^8-1] );
%
%  Problem #2.
%
  w = [   267,  493,  869,  961, 1000, 1153, 1246, 1598, 1766, 1922 ];
  t = 5842;
  [ c, index ] = subset_sum_test01 ( w, t );
%
%  Problem #3.
%
  w = [  518533, 1037066, 2074132, 1648264, 796528, ...
        1593056,  686112, 1372224,  244448, 488896, ...
         977792, 1955584, 1411168,  322336, 644672, ...
        1289344,   78688,  157376,  314752, 629504, ...
        1259008 ];
  t = 2463098;
  [ c, index ] = subset_sum_test01 ( w, t );
%
%  Problem #4.
%
  w = [ 41, 34, 21, 20,  8,  7,  7,  4,  3,  3 ];
  t = 50;
  [ c, index ] = subset_sum_test01 ( w, t );
%
%  Problem #5.
%
  w = [ 81, 80, 43, 40, 30, 26, 12, 11, 9 ];
  t = 100;
  [ c, index ] = subset_sum_test01 ( w, t );
%
%  Problem #6.
%
  w = [ 1, 2, 4, 8, 16, 32 ];
  t = 22;
  r = [ 0, 2^6 - 1 ];
  [ c, index ] = subset_sum_test01 ( w, t, r );
%
%  Problem #7.
%
  w = [ 25, 27, 3, 12, 6, 15, 9, 30, 21, 19 ];
  t = 50;
  r = [ 0, 2^10 - 1 ];
  [ c, index ] = subset_sum_test01 ( w, t, r );
%
%  Problem #1.
%
  w = [ 15, 22, 14, 26, 32, 9, 16, 8 ];
  t = 53;
  count = subset_sum_test02 ( w, t );
%
%  Problem #2.
%
  w = [   267,  493,  869,  961, 1000, 1153, 1246, 1598, 1766, 1922 ];
  t = 5842;
  count = subset_sum_test02 ( w, t );
%
%  Problem #3.
%
  w = [  518533, 1037066, 2074132, 1648264, 796528, ...
        1593056,  686112, 1372224,  244448, 488896, ...
         977792, 1955584, 1411168,  322336, 644672, ...
        1289344,   78688,  157376,  314752, 629504, ...
        1259008 ];
  t = 2463098;
  count = subset_sum_test02 ( w, t );
%
%  Problem #4.
%
  w = [ 41, 34, 21, 20,  8,  7,  7,  4,  3,  3 ];
  t = 50;
  count = subset_sum_test02 ( w, t );
%
%  Problem #5.
%
  w = [ 81, 80, 43, 40, 30, 26, 12, 11, 9 ];
  t = 100;
  count = subset_sum_test02 ( w, t );
%
%  Problem #6.
%
  w = [ 1, 2, 4, 8, 16, 32 ];
  t = 22;
  r = [ 0, 2^6 - 1 ];
  count = subset_sum_test02 ( w, t );
%
%  Problem #7.
%
  w = [ 25, 27, 3, 12, 6, 15, 9, 30, 21, 19 ];
  t = 50;
  r = [ 0, 2^10 - 1 ];
  count = subset_sum_test02 ( w, t );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUBSET_SUM_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function MDP = DEMO_MDP_maze
% Demo of mixed continuous and discrete state space modelling
%__________________________________________________________________________
%
% This demonstration of active inference focuses on navigation and planning
% in a fairly complicated maze. The idea is to demonstrate how epistemic
% foraging and goal (target) directed behaviour are integrated in the
% minimisation of expected free energy. In this illustration, and 8 x 8
% maze is learned through novelty driven evidence accumulation - to learn
% the likelihood mapping between hidden states (locations in the maze) and
% outcomes (whether the current location is open or closed). This
% accumulated experience is then used to plan a path from a start to an end
% (target location) under a task set specified by prior preferences over
% locations. These priors are based upon a simple diffusion (CF backwards
% induction) heuristic that specifies subgoals. The subgoals (i.e.,
% locations) contain the most paths from the target within the horizon of
% the current policy.
%
% We will first illustrate the novelty driven epistemic foraging that
% efficiently scans the maze to learn its structure. We then simulate
% planning of (shortest path) trajectory to the target under the assumption
% the maze has been previously learned. Finally, we consider exploration
% under prior preferences to simulate behaviour when both epistemic and
% goal directed imperatives are in play. The focus on this demo is on
% behavioural and electrophysiological responses over moves.
%
% A key aspect of this formulation is the  hierarchical decomposition of
% goal directed behaviour into subgoals that are within the horizon of a
% limited policy - here, to moves that correspond to a trial. The prior
% preferences then contextualise each policy or trial to ensure that the
% ultimate goal is achieved.
%
% Empirically, this sort of construction suggests the existence of Path
% cells; namely, cells who report the initial location of any subsequence
% and continue firing until the end of the local path. This is illustrated
% by plotting simulated activity as a function of trajectories during
% exploration.
%
% see also: spm_MPD_VB_X.m
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEMO_MDP_maze.m 7766 2020-01-05 21:37:39Z karl $

% set up and preliminaries: first level
%--------------------------------------------------------------------------
rng('default')

% generative model at the sensory level (DEM): continuous states
%==========================================================================
% the generative model has two outcome modalities; namely, what (open
% versus closed) and where (the current location in a maze). These outcomes
% are generated from a single hidden factor (location), where the structure
% of the maze is encoded in the likelihood of observation mapping (that can
% be learned through experience). Allowable actions  include for moves (up,
% down, left, right) and staying at the current location. These induce five
% transition matrices that play the role of empirical priors. Finally,
% prior preferences are based upon allowable transitions (that are function
% of learned accumulated likelihood), which are used to define attractive
% locations within the horizon of two-move policies. These priors implement
% a task set and are returned by a subfunction below: spm_maze_cost
%--------------------------------------------------------------------------
label.factor     = {'where'};
label.modality   = {'what','where'};
label.outcome{1} = {'open','closed'};

MAZE  = [...
    1 1 1 1 1 1 1 1;
    1 0 0 0 0 0 0 1;
    1 1 1 0 1 1 0 1;
    1 1 0 0 0 1 0 1;
    1 1 0 1 0 0 0 1;
    1 1 0 1 1 1 0 1;
    1 0 0 0 0 0 0 1;
    1 0 1 1 1 1 1 1];
END   = sub2ind(size(MAZE),5,5);                  % goal or target location
START = sub2ind(size(MAZE),8,2);                  % first or start location

% prior beliefs about initial states: D
%--------------------------------------------------------------------------
D{1}  = zeros(numel(MAZE),1);
Ns    = numel(D{1});

% probabilistic mapping from hidden states to outcomes: A
%--------------------------------------------------------------------------
A{1}  = [1 - MAZE(:), MAZE(:)]';                  % what
A{2}  = eye(Ns,Ns);                               % where
Ng    = numel(A);
for g = 1:Ng
    No(g)  = size(A{g},1);
end

% controlled transitions: B (up, down, left, right, stay)
%--------------------------------------------------------------------------
u    = [1 0; -1 0; 0 1; 0 -1; 0 0];               % allowable actions
nu   = size(u,1);                                 % number of actions
B{1} = zeros(Ns,Ns,nu);
[n,m] = size(MAZE);
for i = 1:n
    for j = 1:m

        % allowable transitions from state s to state ss
        %------------------------------------------------------------------
        s     = sub2ind([n,m],i,j);
        for k = 1:nu
            try
                ss = sub2ind([n,m],i + u(k,1),j + u(k,2));
                B{1}(ss,s,k) = 1;
            catch
                B{1}(s, s,k) = 1;
            end
        end
    end
end

% allowable policies (2 moves): V
%--------------------------------------------------------------------------
V     = [];
for i = 1:nu
    for j = 1:nu
        V(:,end + 1) = [i;j];
    end
end

% priors: (negative cost) C:
%--------------------------------------------------------------------------
for g = 1:Ng
    C{g} = zeros(No(g),1);
end

% basic MDP structure
%--------------------------------------------------------------------------
mdp.V = V;                      % allowable policies
mdp.A = A;                      % observation model or likelihood
mdp.B = B;                      % transition probabilities
mdp.C = C;                      % preferred outcomes
mdp.D = D;                      % prior over initial states

mdp.label = label;
mdp       = spm_MDP_check(mdp);

% exploratory (32 trial) sequence (no experience or task set)
%==========================================================================
% These simulations use a subroutine (spm_maze_search) to perform recursive
% variational inversions of the active inference scheme under different
% levels of experience and prior preferences
%--------------------------------------------------------------------------
SDP = spm_maze_search(mdp,32,START,END,0,0);

% show results - behavioural
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1'); clf
spm_maze_plot(SDP,END)

% show results - electrophysiological
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf
spm_MDP_VB_LFP(SDP);

% exploratory (8 trial) sequence (with experience and task set)
%==========================================================================
MDP = spm_maze_search(mdp,8,START,END,128,1);

% show results in terms of path
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 3'); clf
spm_maze_plot(MDP,END)

% and implicit subgoals
%--------------------------------------------------------------------------
for i = 1:8
    subplot(8,2,(i - 1)*2 + 2);
    c = spm_softmax(MDP(i).C{2}(:,1));
    imagesc(spm_unvec(c,MAZE))
    axis image off, title(sprintf('%s %i','Trial',i))
end

% evidence accumulation under task set
%==========================================================================
MDP   = mdp;
for i = 1:4

    % return to start
    %----------------------------------------------------------------------
    MDP = spm_maze_search(MDP(end),8,START,END,0,1);

    % show results in terms of path
    %----------------------------------------------------------------------
    spm_figure('GetWin',sprintf('%s %i','Search',i)); clf
    spm_maze_plot(MDP,END)

    % and electrophysiological responses
    %----------------------------------------------------------------------
    spm_figure('GetWin',sprintf('%s %i','Responses',i)); clf
    spm_MDP_VB_game(MDP(1:8));
    subplot(6,1,3), set(gca,'YLim',[-1 1]/3);
    subplot(6,1,4), set(gca,'YLim',[0 1]);

end

% in silico psychophysical experiment
%==========================================================================
% Here, we return to the exploratory simulation above and probe each level
% of experience by asking the subject to execute a path to target. The
% behaviour is then assessed in terms of the latency with which the target
% will go is required  - and the number of mistakes or exploratory
% excursions into closed locations.
%--------------------------------------------------------------------------
N     = [];
M     = [];
T     = [];
for i = 1:2:numel(SDP)

    % execute path to target with increasing experience
    %----------------------------------------------------------------------
    MDP = spm_maze_search(SDP(i),8,START,END,0,1);

    % record behaviour
    %----------------------------------------------------------------------
    s   = spm_cat({MDP.s}); s(3:3:end) = [];
    s   = [s,END];
    N   = [N,find(s == END,1)];                      % latency
    M   = [M,sum(MAZE(s))];                          % mistakes
    T   = [T,i/2];                                   % exposure to maze

end

% show results
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 4'); clf
subplot(2,1,1), bar([M;N]'), xlabel('Exposure (seconds)'), axis square
title('Performance','FontSize',16),legend({'Mistakes','Latency'})

% place and direction specific responses
%==========================================================================
spm_figure('GetWin','Figure 5'); clf

% illustrate simulated responses (first eight subsequences)
%--------------------------------------------------------------------------
spm_MDP_VB_LFP(SDP(1:8))

% extract simulated responses
%--------------------------------------------------------------------------
[u,v] = spm_MDP_VB_LFP(SDP);               % responses
E     = kron(eye(32*3),pinv(ones(16,1)));  % expectation matrix
v     = E*spm_cat(v);                      % firing rates (hidden states)
u     = E*spm_cat({SDP.un})';              % firing rates (policies)
s     = spm_cat({SDP.s})';                 % location

% Position
%--------------------------------------------------------------------------
[X,Y] = ndgrid(1:8,1:8);
L     = [X(:),Y(:)];
X     = L(s,:);

% accumulate responses at each location
%--------------------------------------------------------------------------
PC    = zeros(8,8,size(v,2));
for i = 1:size(v,1)
    for j = 1:size(v,2)
       PC(X(i,1),X(i,2),j) = PC(X(i,1),X(i,2),j) + v(i,j);
    end
end
QC    = zeros(8,8,size(u,2));
for i = 1:size(u,1)
    for j = 1:size(u,2)
       QC(X(i,1),X(i,2),j) = QC(X(i,1),X(i,2),j) + u(i,j);
    end
end

% identify path and place cells that fire above threshold (U)
%--------------------------------------------------------------------------
U     = .8;
jpath = [];
jplac = [];
for i = 1:64
    if sum(spm_vec(PC(:,:,i)) > U) > 2
        jpath(end + 1) = i;
    end
end
for i = (1:64) + 128;
    if sum(spm_vec(PC(:,:,i)) > U) > 0
        jplac(end + 1) = i;
    end
end

% and plot as a function of trajectory in space
%--------------------------------------------------------------------------
x     = spm_conv(X + randn(size(X))/8,2,0);
for i = 1:64
    col{i} = spm_softmax(randn(3,1));
end

% path cells
%--------------------------------------------------------------------------
subplot(4,2,3)
plot(x(:,2),x(:,1),'r:'), hold on
for i = 1:size(x,1)
    for j = 1:3:length(jpath)
        if v(i,jpath(j)) > .80
            plot(x(i,2),x(i,1),'.','MarkerSize',32,'Color',col{j})
        end
    end
end
axis([0 9 0 9]), axis ij square, hold off
title('Path cell responses','fontsize',16)

% please cells
%--------------------------------------------------------------------------
subplot(4,2,4)
plot(x(:,2),x(:,1),'r:'), hold on
for i = 1:size(x,1)
    for j = 1:3:length(jpath)
        if v(i,jplac(j)) > .80
            plot(x(i,2),x(i,1),'.','MarkerSize',32,'Color',col{j})
        end
    end
end
axis([0 9 0 9]), axis ij square, hold off
title('Place cell responses','fontsize',16)

return

function MDP = spm_maze_search(mdp,N,START,END,alpha,beta)
% FORMAT MDP = spm_maze_search(mdp,N,START,END,alpha,beta)
% mdp   - MDP structure
% N     - number of trials (i.e., policies: default 8)
% START - index of intial state (default 1)
% END   - index of target state (default 1)
% alpha - prior concentration parameter for likelihood (default 128)
% beta  - precision of prior preference (default 0)
%the argument is
% MDP   - MDP structure array

% preliminaries
%--------------------------------------------------------------------------
try, N;     catch, N     = 8;   end
try, START; catch, START = 1;   end
try, END;   catch, END   = 1;   end
try, alpha; catch, alpha = 128; end
try, beta;  catch, beta  = 0;   end

% initialise concentration parameters: a (if unspecified)
%--------------------------------------------------------------------------
if ~isfield(mdp,'a')
    mdp.a{1} = ones(size(mdp.A{1}))/8 + mdp.A{1}*alpha;
    mdp.a{2} = mdp.A{2}*128;
end
if ~isfield(mdp,'o')
    mdp.o = [];
end
if ~isfield(mdp,'u')
    mdp.u = [];
end
mdp.s = START;

% Evaluate a sequence of moves - recomputing prior preferences at each move
%==========================================================================
for i = 1:N

    % Evaluate preferred states (subgoals) on the basis of current beliefs
    %----------------------------------------------------------------------
    mdp.C{2} = spm_maze_cost(mdp,END)*beta;

    % proceed with subsequent trial
    %----------------------------------------------------------------------
    MDP(i)   = spm_MDP_VB_X(mdp);
    mdp      = MDP(i);
    mdp.s    = mdp.s(:,end);
    mdp.D{1} = MDP(i).X{1}(:,end);
    mdp.o    = [];
    mdp.u    = [];

end

return

function C = spm_maze_cost(MDP,END)
% Evaluate subgoals using graph Laplacian
%==========================================================================
START = MDP.s(1);
if isfield(MDP,'a')
    Q = MDP.a{1};
else
    Q = MDP.A{1};
end
Q   = Q/diag(sum(Q));
Q   = Q(1,:);                                % open states
P   = diag(Q)*any(MDP.B{1},3);               % allowable transitions
ns  = length(Q);                             % number of states
X   = zeros(ns,1);X(START) = 1;              % initial state
Y   = zeros(ns,1);Y(END)   = 1;              % target state

% Preclude transitions to closed states and evaluate graph Laplacian
%--------------------------------------------------------------------------
P   = P - diag(diag(P));
P   = P - diag(sum(P));
P   = expm(P);

% evaluate (negative) cost as a path integral conjunctions
%--------------------------------------------------------------------------
for t = 1:size(MDP.V,1)
    X = P*X;
end
X     = X > exp(-3);
C     = log(X.*(P*Y) + exp(-32));

return

function spm_maze_plot(MDP,END)
% illustrate  search graphically
%--------------------------------------------------------------------------
A  = spm_vec(MDP(1).A{1}(1,:));
ns = numel(A);
ni = sqrt(ns);
A  = reshape(A,ni,ni);
subplot(2,2,1), imagesc(A), axis image
title('Scanpath','fontsize',16);

% Cycle of the trials
%--------------------------------------------------------------------------
h     = [];
MS    = {};
MC    = {};
for p = 1:numel(MDP)

    %  current beliefs and preferences: A likelihood
    %----------------------------------------------------------------------
    if isfield(MDP,'a')
        Q = MDP(p).a{1};
    else
        Q = MDP(p).A{1};
    end
    Q     = Q/diag(sum(Q));
    Q     = Q(1,:);
    a     = reshape(Q(:),ni,ni);
    subplot(2,2,2), imagesc(a), axis image
    title('Likelihood','fontsize',16);

    %  current beliefs and preferences: B transitions
    %----------------------------------------------------------------------
    try
        b = diag(Q)*any(MDP(p).B{1},3);
    catch
        b = diag(Q)*any(MDP(p).B{1},3);
    end
    subplot(2,2,4), imagesc(-b), axis image
    title('Allowable transitions','fontsize',16);

    %  current beliefs and preferences: C preferences
    %----------------------------------------------------------------------
    C     = MDP(p).C{2}(:,1);
    C     = spm_softmax(C);
    C     = reshape(C,ni,ni);
    subplot(2,2,3), imagesc(C), axis image
    title('Preferences','fontsize',16);
    try
        [i,j] = ind2sub([ni,ni],MDP(p).s(1)); hold on
        plot(j,i,'.','MarkerSize',32,'Color','g');
        [i,j] = ind2sub([ni,ni],END);
        plot(j,i,'.','MarkerSize',32,'Color','r'); hold off
    end

    % cycle over  short-term searches
    %----------------------------------------------------------------------
    subplot(2,2,1),hold on
    s     = MDP(p).s;
    for t = 1:numel(s)

        % location
        %------------------------------------------------------------------
        [i,j] = ind2sub([ni,ni],s(t));
        h(end + 1) = plot(j,i,'.','MarkerSize',32,'Color','r');
        try
            set(h(end - 1),'Color','m','MarkerSize',16);
            j = [get(h(end - 1),'Xdata'), get(h(end),'Xdata')];
            i = [get(h(end - 1),'Ydata'), get(h(end),'Ydata')];
            plot(j,i,':r');
        end

        % save
        %------------------------------------------------------------------
        if numel(MS)
            MS(end + 1) = getframe(gca);
        else
            MS = getframe(gca);
        end

    end

    % save
    %----------------------------------------------------------------------
    subplot(2,2,3)
    if numel(MC)
        MC(end + 1) = getframe(gca);
    else
        MC = getframe(gca);
    end

end

% save movie
%--------------------------------------------------------------------------
subplot(2,2,1)
xlabel('click axis for movie')
set(gca,'Userdata',{MS,16})
set(gca,'ButtonDownFcn','spm_DEM_ButtonDownFcn')

subplot(2,2,3)
xlabel('click axis for movie')
set(gca,'Userdata',{MC,16})
set(gca,'ButtonDownFcn','spm_DEM_ButtonDownFcn')

author: spm


text:
function Rx = getPrescribedDose(planNum,planC)
% Get prescribed dose (Gy) for selected plan
%
%AI 07/09/21

indexS = planC{end};

%Get referenced beam SOPInstanceUID
DCMheaderS = planC{indexS.dose}(planNum).DICOMHeaders;
ReferencedSOPInstanceUID = DCMheaderS.ReferencedRTPlanSequence.Item_1.ReferencedSOPInstanceUID;

%Identify corresponding beamNum
SOPInstanceUIDv = {planC{indexS.beams}.SOPInstanceUID};
beamNum = strcmpi(ReferencedSOPInstanceUID,SOPInstanceUIDv);

%Get prescribed dose
Rx = planC{indexS.beams}(beamNum).DoseReferenceSequence.Item_1.DeliveryMaximumDose;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function a = rutis3 ( )

%*****************************************************************************80
%
%% RUTIS3 returns the RUTIS3 matrix.
%
%  Example:
%
%    4 -5  0  3
%    0  4 -3 -5
%    5 -3  4  0
%    3  0  5  4
%
%  Properties:
%
%    A is not symmetric: A' /= A.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A has distinct eigenvalues.
%
%    A has a pair of complex eigenvalues.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%

%
%  Note that the matrix entries are listed by row.
%
  a(1:4,1:4) = [ ...
    4.0, -5.0,  0.0,  3.0; ...
    0.0,  4.0, -3.0, -5.0; ...
    5.0, -3.0,  4.0,  0.0; ...
    3.0,  0.0,  5.0,  4.0 ];

  return
end


text:
function SDataOut = DicomImport(sPath, lRecursive, csAdditionalTags)
% DICOMIMPORT Import DICOM series with preview.
%
%   SDATA = DICOMIMPORT(sPATH, lRECURSIVE) Scans the folder specified by
%   sPATH for DICOM files. If lRECURSIVE is true, DICOMIMPORT
%   hierarchically scans all subolders of sPATH. Valid DICOM files are
%   sorted into series according the DICOM tag SeriesInstanceUID and the
%   folder in which the files were found. For each series, a thumbnail
%   image is created and an overview of the obtained series is displayed in
%   a GUI that lets you select the series to load. Select series as if you
%   would in a file browser by klicking and using the SHIFT and CNTL keys
%   to select multiple series Close the GUI by pressing <ENTER> or <ESCAPE>
%   (the latter returns an empty array). Output SDATA is a struct where
%   length(SDATA) equals the amout of selected series. SDATA contains at
%   least the following fields:
%
%       Img:                The image volume
%       SeriesDescription:  The contents of the corresponding DICOM tag
%       Orientation:        The most prominent image volume orientation
%       Aspect:             The pixel spacing in all 3 directions
%       ImageOrientation:   The direction cosines (compare DICOM standard)
%       ImagePosition:      The coordinates of the upper left corner of
%                           each image slice (compare DICOM standard)
%
%   The tags specified in csADDITIONALTAGS (see below) are appended to this
%   structure.
%
%
%   DICOMIMPORT lets you select the base folder in a dialog and asks if the
%   sobfolders are to be scanned.
%
%   DICOMIMPORT(sPATH, lRECURSIVE, csADDITIONALTAGS) Lets you specify
%   additional DICOM tags you can use for further processing in the cell
%   array of strings csADDITIONALTAGS
%
%
% Copyright 2013 Christian Wuerslin, University of Tuebingen and
% University of Stuttgart, Germany.
% Contact: christian.wuerslin@med.uni-tuebingen.de

% =========================================================================
% *** FUNCTION DicomImport
% ***
% *** Main GUI function. See above for description.
% ***
% =========================================================================

% -------------------------------------------------------------------------
% Control the figure's appearence
SAp.sTITLE              = 'Dicom Import [press enter to confirm]';
SAp.iTHUMBNAILSIZE      = 200;
SAp.iINFOBARHEIGHT      = 125;
SAp.dBGCOLOR            = [0.1 0.2 0.3];
SAp.iNROWS              = 2;
SAp.iNCOLS              = 4;
SAp.iTITLEHEIGHT        = 16;
SAp.dOPACITY            = 0.5;
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Define a list of DICOM tags to read (those that I really need)
csTags = {'SeriesInstanceUID',       'Rows',                 'Columns', ...
          'ImageOrientationPatient', 'ImagePositionPatient', 'PixelSpacing', ...
          'SeriesDescription', 'Modality'};
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Parse input arguments
if nargin
    if ~ischar(sPath), error('Input argument must be a string!'); end
else
    sPath = uigetdir;
    if isnumeric(sPath)
        SDataOut = [];
        return
    end
end

if nargin < 2, lRecursive = strcmp(questdlg('Search subdirectories for DICOM files?', 'Stupid question', 'Yup', 'Noop', 'Yup'), 'Yup'); end
if nargin > 2, csTags = [csTags, csAdditionalTags]; end
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Get all DICOM files and their headers
fprintf(1, 'Scanning for DICOM files...');
SData = fGetDicomFiles(sPath, csTags, lRecursive);
fprintf(1, 'done.\n');
if isempty(SData), return, end
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Sort the series into a struct and get a thumbnail of each set
fprintf(1, 'Creating thumbnails...');
iHashs = [SData.iHash];
iSeriesHash = fUnique(iHashs);
iNDatasets = length(iSeriesHash);
iThumbnails = zeros(SAp.iTHUMBNAILSIZE, SAp.iTHUMBNAILSIZE, iNDatasets);
for iDataset = 1:iNDatasets
    iThisSeries = find([SData.iHash] == iSeriesHash(iDataset));
    SSeries(iDataset).SData = SData(iThisSeries);

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Load the center image and create thumbnail
    iInd = iThisSeries(round(length(iThisSeries)./2));
    dImg = double(dicomread(SData(iInd).sFilename));
    if isempty(dImg), continue, end

    dImg = dImg - min(dImg(:));
    dX = linspace(1, length(dImg), SAp.iTHUMBNAILSIZE);
    [dYY, dXX] = meshgrid(dX, dX);
    dImg = interp2(dImg, dYY, dXX, 'linear*', 0);
    if max(dImg(:)), dImg = dImg./max(dImg(:)); end
    iThumbnails(:,:,iDataset) = dImg;
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end
lSelected = false(iNDatasets, 1);
fprintf(1, 'done.\n');
fprintf(1, 'Found %u files in %u datasets!\n', length(SData), iNDatasets);
clear iDataset iThisSeries iInd dImg dX dXX dYY
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Create the figure.
dFigureWidth  = (SAp.iTHUMBNAILSIZE + 1).*SAp.iNCOLS;
dFigureHeight = (SAp.iTHUMBNAILSIZE + SAp.iTITLEHEIGHT + 1).*SAp.iNROWS + SAp.iINFOBARHEIGHT + 2;
dScreenSize = get(0, 'MonitorPositions');
try
    hF = figure(...
        'Position'             , [(dScreenSize(3) - dFigureWidth )./2, ...
        (dScreenSize(4) - dFigureHeight)./2, ...
        dFigureWidth, dFigureHeight], ...
        'Units'                , 'pixels', ...
        'Color'                , SAp.dBGCOLOR, ...
        'Resize'               , 'off', ...
        'DockControls'         , 'off', ...
        'MenuBar'              , 'none', ...
        'Name'                 , SAp.sTITLE, ...
        'NumberTitle'          , 'off', ...
        'BusyAction'           , 'cancel', ...
        'KeyPressFcn'          , @fKeyPressFcn, ...
        'CloseRequestFcn'      , @fCloseGUI, ...
        'WindowButtonDownFcn'  , @fWindowButtonDownFcn, ...
        'WindowButtonMotionFcn', @fWindowMouseMoveFcn, ...
        'WindowScrollWheelFcn' , @fWindowScrollWheelFcn);
catch %#ok<CTCH>
    hF = figure(...
        'Position'             , [(dScreenSize(3) - dFigureWidth )./2, ...
        (dScreenSize(4) - dFigureHeight)./2, ...
        dFigureWidth, dFigureHeight], ...
        'Units'                , 'pixels', ...
        'Color'                , SAp.dBGCOLOR, ...
        'Resize'               , 'off', ...
        'DockControls'         , 'off', ...
        'MenuBar'              , 'none', ...
        'Name'                 , SAp.sTITLE, ...
        'NumberTitle'          , 'off', ...
        'BusyAction'           , 'cancel', ...
        'KeyPressFcn'          , @fKeyPressFcn, ...
        'CloseRequestFcn'      , @fCloseGUI, ...
        'WindowButtonDownFcn'  , @fWindowButtonDownFcn, ...
        'WindowButtonMotionFcn', @fWindowMouseMoveFcn);
end
clear dScreenSize
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Create the axes
hAxes = zeros(SAp.iNCOLS.*SAp.iNROWS, 1);
hImg  = zeros(SAp.iNCOLS.*SAp.iNROWS, 1);
hText = zeros(SAp.iNCOLS.*SAp.iNROWS, 1);
for iM = 1:SAp.iNROWS
    for iN = 1:SAp.iNCOLS
        iLinInd = (iM - 1).*SAp.iNCOLS + iN;
        dVertPos = dFigureHeight - iM.*(SAp.iTHUMBNAILSIZE + SAp.iTITLEHEIGHT + 1) + 1;
        iPosition = [(iN - 1).*(SAp.iTHUMBNAILSIZE + 1) + 2, dVertPos, SAp.iTHUMBNAILSIZE, SAp.iTHUMBNAILSIZE];
        hAxes(iLinInd) = axes('Parent' , hF, 'Units', 'pixels', 'Position', iPosition);
        hImg(iLinInd) = image(zeros(SAp.iTHUMBNAILSIZE, SAp.iTHUMBNAILSIZE, 3), 'Parent', hAxes(iLinInd));
        hText(iLinInd) = uicontrol(...
            'Parent'            , hF, ...
            'Style'             , 'text', ...
            'Units'             , 'pixels', ...
            'Position'          , [iPosition(1), iPosition(2) + iPosition(4), SAp.iTHUMBNAILSIZE, SAp.iTITLEHEIGHT], ...
            'BackgroundColor'   , SAp.dBGCOLOR./2, ...
            'ForegroundColor'   , 'w', ...
            'HorizontalAlign'   , 'left', ...
            'FontUnits'         , 'pixels', ...
            'FontSize'          , 12);
    end
end
axis(hAxes, 'off');
clear iM iN
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Create the info text objects
uicontrol(...
    'Parent'            , hF, ...
    'Style'             , 'text', ...
    'Units'             , 'pixels', ...
    'Position'          , [2, 2, 200, SAp.iINFOBARHEIGHT], ...
    'BackgroundColor'   , SAp.dBGCOLOR./2, ...
    'ForegroundColor'   , 'w', ...
    'HorizontalAlign'   , 'left', ...
    'FontUnits'         , 'pixels', ...
    'FontSize'          , 14, ...
    'String'            , {'Series Description:', 'Path:', 'Modality:', 'Rows:', 'Columns:', 'Images:', 'Pixel Spacing:'});
hInfoText = uicontrol(...
    'Parent'            , hF, ...
    'Style'             , 'text', ...
    'Units'             , 'pixels', ...
    'Position'          , [202, 2, dFigureWidth - 202, SAp.iINFOBARHEIGHT], ...
    'BackgroundColor'   , SAp.dBGCOLOR./2, ...
    'ForegroundColor'   , 'w', ...
    'HorizontalAlign'   , 'left', ...
    'FontUnits'         , 'pixels', ...
    'FontSize'          , 14, ...
    'String'            , '');
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Create a beatuiful mask for the selection highlight
dMask = repmat(permute(SAp.dBGCOLOR, [1 3 2]) , [SAp.iTHUMBNAILSIZE, SAp.iTHUMBNAILSIZE, 1]);
dMask = dMask.*repmat(linspace(1, 0.5, SAp.iTHUMBNAILSIZE)', [1, SAp.iTHUMBNAILSIZE, 3]);
dMask = dMask + 0.05.*rand(size(dMask));
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Fill the figure
iStartSeries = 1;
iLastSeries = 0;
sResult = 'esc';
fFillPanels();
% -------------------------------------------------------------------------

uiwait(hF);
% -------------------------------------------------------------------------
% ~~~~~~~~~~~~ A lot of user-friendly GUI interaction going on ~~~~~~~~~~~~
% -------------------------------------------------------------------------

% The GUI was closed in some way, continue execution

% -------------------------------------------------------------------------
% Dialog was aborted: Seriously, who does that?
if strcmp(sResult, 'esc')
    SDataOut = [];
    try delete(hF); end
    return
end
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% No seriers selected? - Return
iInd = find(lSelected);
if isempty(iInd)
    SDataOut = [];
    delete(hF);
    return
end
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Load the selected series
fprintf(1, 'Loading slected DICOM files...');
for i = 1:length(iInd)
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Load the specified images
    SThisData = SSeries(iInd(i)).SData;
    iNImages = length(SThisData);
    iImg = zeros(SThisData(1).Rows, SThisData(1).Columns, iNImages, 'uint16');
    for j = 1:iNImages
        iImg(:,:,j) = dicomread(SThisData(j).sFilename);
    end
    % Now images are in iImg, headers in SThisData
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Find out the major orientation
    dImageOrientation = reshape(SThisData(1).ImageOrientationPatient, [3, 2])'; % Should be the same for a volume
    dImagePosition = [SThisData.ImagePositionPatient]';

    dOrientIndicator = sum(abs(dImageOrientation));
    [temp, iMinInd] = min(dOrientIndicator); %#ok<ASGLU>
    d3rdDimInd = dImagePosition(:, iMinInd);
    [temp, iSortInd] = sort(d3rdDimInd, 'ascend'); %#ok<ASGLU>
    iImg = iImg(:,:,iSortInd);
    iImageOrientation = round(dImageOrientation);
    switch(iMinInd)
        case 1 % The x coordinate changes least in both directions -> Sagittal
            SDataOut(i).Orientation = 'Sag'; % -> x: a->p=pos, y: h->f=neg
            if iImageOrientation(1, 2) == -1, flipdim(iImg, 2); end % <- y-coordinate was flipped
            if iImageOrientation(2, 3) ==  1, flipdim(iImg, 1); end % <- z-coordinate was flipped

        case 2 % The y coordinate changes least in both directions -> Coronal
            SDataOut(i).Orientation = 'Cor'; % -> x: r->l=pos, y: h->f=neg
            if iImageOrientation(1, 1) == -1, flipdim(iImg, 2); end % <- x-coordinate was flipped
            if iImageOrientation(2, 3) ==  1, flipdim(iImg, 1); end % <- z-coordinate was flipped

        case 3 % The z coordinate changes least in both directions -> Transversal
            SDataOut(i).Orientation = 'Tra'; % -> x: r->l=pos, y: a->p=pos
            if iImageOrientation(1, 1) == -1, flipdim(iImg, 2); end % <- x-coordinate was flipped
            if iImageOrientation(2, 2) == -1, flipdim(iImg, 1); end % <- Y-coordinate was flipped
    end
    SDataOut(i).Img = iImg;
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Save the most important infos to the ouput variable
    SDataOut(i).SeriesDescriptions = SThisData(1).SeriesDescription;
    SDataOut(i).Aspect = zeros(3, 1);
    SDataOut(i).Aspect(1:2) = SThisData(1).PixelSpacing;
    if ~isscalar(d3rdDimInd), SDataOut(i).Aspect(3)   = abs(d3rdDimInd(2) - d3rdDimInd(1)); end
    SDataOut(i).ImageOrientation = dImageOrientation;
    SDataOut(i).ImagePosition = dImagePosition(iSortInd, :);
    if nargin == 3
        for j = 1:length(csAdditionalTags) % The additionally requested infos
            sTag = csAdditionalTags{j};
            eval(['xData = [SThisData.', sTag, ']'';']);
            xData = xData(iSortInd);
            eval(['SDataOut(i).', sTag, ' = xData;']);
        end
    end
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end
% -------------------------------------------------------------------------
fprintf(1, 'done.\n');

delete(hF);
% =========================================================================
% ***
% *** The 'end' of the DICOMIMPORT main function.
% ***
% =========================================================================

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * NESTED FUNCTION fCloseGUI (nested in DicomImport)
    % * *
    % * * Figure callback
    % * *
    % * * Closes the figure.
    % * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function fCloseGUI(hObject, eventdata) %#ok<*INUSD> eventdata is repeatedly unused
        uiresume(hObject);
        delete(hObject);
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fCloseGUI
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * NESTED FUNCTION fWindowMouseMoveFcn (nested in DicomImport)
    % * *
    % * * Figure callback
    % * *
    % * * Displays informations about the series under the mouse cursor in
    % * * the texts at the bottom of the figure.
    % * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function fWindowMouseMoveFcn(hObject, eventdata)
        if ~exist('hImg', 'var'), return, end; % Return if called during GUI startup

        iAxisInd = fGetAxes();
        iSeriesInd = iAxisInd + iStartSeries -1;

        % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        % Return if not over axes or out of data bounds
        if ~iAxisInd, set(hInfoText, 'String', ''); return, end
        if iSeriesInd > length(lSelected),  set(hInfoText, 'String', ''); return, end
        % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        % Print Information
        sText = SSeries(iSeriesInd).SData(1).SeriesDescription;
        if isempty(sText), sText = 'N/A'; end, csText(1) = {sText};

        sText = fileparts(SSeries(iSeriesInd).SData(1).sFilename);
        csText(2) = {sText};

        sText = SSeries(iSeriesInd).SData(1).Modality;
        if isempty(sText), sText = 'N/A'; end, csText(3) = {sText};

        sText = SSeries(iSeriesInd).SData(1).Rows;
        if isempty(sText), sText = 'N/A'; end, csText(4) = {sText};

        sText = SSeries(iSeriesInd).SData(1).Columns;
        if isempty(sText), sText = 'N/A'; end, csText(5) = {sText};

        sText = num2str(length(SSeries(iSeriesInd).SData));
        if isempty(sText), sText = 'N/A'; end, csText(6) = {sText};

        dPixelSpacing =  SSeries(iSeriesInd).SData(1).PixelSpacing;
        if ~isempty(dPixelSpacing)
            sText = sprintf('%2.2f x %2.2f mm', dPixelSpacing(1), dPixelSpacing(2));
        else
            sText = 'N/A';
        end
        csText(7) = {sText};

        set(hInfoText, 'String', csText);
        % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fWindowMouseMoveFcn
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * NESTED FUNCTION fWindowButtonDownFcn (nested in DicomImport)
    % * *
    % * * Figure callback
    % * *
    % * * Starting callback for mouse button actions. Manage the selection
    % * * of the series.
    % * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function fWindowButtonDownFcn(hObject, eventdata)
        iPanelInd = fGetAxes();
        iSeriesInd = iPanelInd + iStartSeries - 1;

        if ~iPanelInd, return, end % Exit if Event didn't occurr in axes
        if iSeriesInd > length(lSelected), return, end

        switch get(hF, 'SelectionType')

            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            % Left mouse button without any key
            case 'normal'
                iNSelected = sum(lSelected);
                lState = lSelected(iSeriesInd);
                lSelected = false(size(lSelected));
                lSelected(iSeriesInd) = ~lState || iNSelected > 1;
                if ~lState, iLastSeries = iSeriesInd; end
            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            % Right mouse button/shift key
            case 'extend'
                if iLastSeries
                    iInd = sort([iLastSeries, iSeriesInd], 'ascend');
                    lSelected(iInd(1):iInd(2)) = true;
                else
                    lSelected(iSeriesInd) = true;
                    iLastSeries = iSeriesInd;
                end
            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            % Middle mouse button/alt key
            case 'alt'
                lState = lSelected(iSeriesInd);
                lSelected(iSeriesInd) = ~lState;
                if ~lState, iLastSeries = iSeriesInd; end
            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        end
        fFillPanels;
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fWindowButtonDownFcn
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * NESTED FUNCTION fKeyPressFcn (nested in DicomImport)
    % * *
    % * * Figure callback
    % * *
    % * * Callback for keyboard actions.
    % * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function fKeyPressFcn(hObject, eventdata) %#ok<INUSL>

        switch eventdata.Key

            case 'uparrow'
                iStartSeries = max([1, iStartSeries - SAp.iNCOLS]);
                fFillPanels;

            case 'downarrow'
                if iStartSeries + SAp.iNCOLS < length(lSelected),
                    iStartSeries = iStartSeries + SAp.iNCOLS;
                end

            case 'return' % Return and load the selected series
                sResult = 'OK';
                uiresume(hF);

            case 'escape' % Return and do nothing
                uiresume(hF);

        end % switch

        fFillPanels;
        fWindowMouseMoveFcn(hF, []); % Make sure the info at the bottom is valid

    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fKeyPressFcn
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * NESTED FUNCTION fWindowScrollWheelFcn (nested in DicomImport)
    % * *
    % * * Figure callback
    % * *
    % * * Callback for keyboard actions.
    % * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function fWindowScrollWheelFcn(hObject, eventdata)
        if eventdata.VerticalScrollCount > 0
            SEData.Key = 'downarrow';
        else
            SEData.Key = 'uparrow';
        end
        fKeyPressFcn(hObject, SEData); % This is the laszy way: just call the other callback
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fWindowScrollWheelFcn
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * *
    % * * NESTED FUNCTION fFillPanels (nested in DicomImport)
    % * *
    % * * Display the current data in all axes.
	% * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function fFillPanels()
        for iI = 1:length(hImg);
            iSeriesInd = iI + iStartSeries - 1;
            if iSeriesInd > length(lSelected)
                set(hImg(iI), 'CData', dMask./3);
                set(hText(iI), 'String', '<no data>');
            else
                dImg = repmat(iThumbnails(:,:,iSeriesInd), [1 1 3]);
                if lSelected(iSeriesInd),  dImg = 1 - ((1 - dImg).*(1 - SAp.dOPACITY.*dMask)); end
                set(hImg(iI), 'CData', dImg);
                sTitle = SSeries(iSeriesInd).SData(1).SeriesDescription;
                if length(sTitle) > 22, sTitle = [sTitle(1:21), '...']; end
                if isempty(sTitle), sTitle = '<no description>'; end
                set(hText(iI), 'String', sprintf('[%u]: %s', iSeriesInd, sTitle));
            end
        end
        drawnow expose
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fFillPanels
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * *
    % * * NESTED FUNCTION fGetAxes (nested in DicomImport)
    % * *
    % * * Determine the axes number under the mouse cursor. Returns 0 if
    % * * not over an axes at all.
	% * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function iPanelInd = fGetAxes()
        iCursorPos = get(hF, 'CurrentPoint');
        iPanelInd = uint8(0);
        for iI = 1:length(hAxes)
            dPos = get(hAxes(iI), 'Position');
            if ((iCursorPos(1) >= dPos(1)) && (iCursorPos(1) < dPos(1) + dPos(3)) && ...
                (iCursorPos(2) >= dPos(2)) && (iCursorPos(2) < dPos(2) + dPos(4)))
                iPanelInd = uint8(iI);
            end
        end
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fGetAxes
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * *
    % * * NESTED FUNCTION fGetDicomFiles (nested in DicomImport)
    % * *
    % * * Finds and creates a database of DICOM files
	% * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function SD = fGetDicomFiles(sFolder, csTags, lRecursive)
        SD = [];
        SDir = dir(sFolder);
        for iI = 1:length(SDir)
            sName = SDir(iI).name;

            if (SDir(iI).isdir) && (lRecursive)
                if sName(1) == '.', continue, end
                SD = [SD, fGetDicomFiles([sFolder, filesep, SDir(iI).name], csTags, lRecursive)];
            else
                sFilename = [sFolder, filesep, SDir(iI).name];
                try
                    SHeader = dicominfo(sFilename);
                catch %#ok<CTCH>
                    continue
                end
                SRecord = [];
                for iJ = 1:length(csTags)
                    sTag = csTags{iJ};
                    try
                        eval(['SRecord.', sTag, ' = SHeader.', sTag, ';']);
                    catch %#ok<CTCH>
                        eval(['SRecord.', sTag, ' = [];']);
                    end
                end
                if isfield(SHeader, 'SeriesInstanceUID')
                    SRecord.iHash = fHash([sFolder, SHeader.SeriesInstanceUID]);
                    SRecord.sFilename = sFilename;
                    SD = [SD, SRecord];
                end
            end
        end
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fGetDicomFiles
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * *
    % * * NESTED FUNCTION fHash (nested in DicomImport)
    % * *
    % * * Return a hash number
	% * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function iHashVal = fHash(sString)
        dString = double(sString);
        dHashVal = 5381;
        for iI = 1:length(dString)
            dHashVal = mod(dHashVal * 33 + dString(iI), 2.^32 - 1);
        end
        iHashVal = uint32(dHashVal);
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fHash
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * *
    % * * NESTED FUNCTION fUnique (nested in DicomImport)
    % * *
    % * * Returns unique values in order of appearance
	% * *
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    function xUniqueVals = fUnique(xArray)
        xUniqueVals = zeros(length(xArray), 1, class(xArray));
        iInd = 0;
        while ~isempty(xArray)
            iInd = iInd + 1;
            xUniqueVals(iInd) = xArray(1);
            xArray(xArray == xArray(1)) = [];
        end
        xUniqueVals = xUniqueVals(1:iInd);
    end
    % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
    % * * END NESTED FUNCTION fUnique
	% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

end

author: thomaskuestner


text:
function [AUC_FV,AUC_ges,AUC_GT,AUC_trans]=finalexp(method)
%This script run different methods on different sets of features to fill the final table
% method may be 'NN'(network),'tree','randforest','bayes','log','boost', 'linearsvm','kernelsvm'
%csv_path = '../Annotation/saved_anno.mat';
%load(csv_path);

disp(['Classifier is ', method])
% read CV
fea_path = '../dataset_trial/FVs';
fea_types = {'mbhx','mbhy'};

% Load all FV features
if ~exist('feaMap.mat')
    feaMap = containers.Map;
    for fea_i = 1:length(fea_types)
        fea_ext = [fea_types{fea_i},'.fv.txt'];
        fea_list = dir([fea_path,'/*',fea_ext]);
        for vid_i = 1:length(fea_list)
            tmpfea = importdata(fullfile(fea_path,fea_list(vid_i).name));
            feaMap(fea_list(vid_i).name) = tmpfea;
        end
    end
    save('feaMap.mat','feaMap');
else
    load('feaMap.mat');
end

prob_lf = cell(1, 10);
prob_FV = cell(1, 10);  %Fisher Vector
prob_ges = cell(1,10);  %pred gesture
prob_GT = cell(1, 10);  %GT gesture
prob_MFCC = cell(1,10); %MFCC feature
prob_trans = cell(1,10);    %Transcript feature
lab_vec = cell(1, 10);  %label vector
AUC_gt = zeros(1,10);

%fine grid search on gesture features on different weights (0 0.5 1)
for cv = 0:9
    disp(['CROSS VALIDATION:', num2str(cv)])
    if ~exist(sprintf(['grid_search/final%d_',method,'.mat'],cv+1),'file')
        disp('IDTFV');
        tmpprob_FV = test_dtfv_CV(cv, feaMap, method);  %IDT+FV
        disp('Pred');
        [tmplab,tmpprob, w] = test_gestureSearch(cv, method); %Pred Gesture
        disp('GT')
        tmpprob_GT = test_GTgestureSearch(cv, method);  %GT Gesture, currently use 5 GT
        disp('trans');
        tmpprob_trans = test_Trans(cv, method);
        save(sprintf(['grid_search/final%d_',method,'.mat'],cv+1), 'tmpprob_FV','tmplab','tmpprob','tmpprob_GT','tmpprob_trans');
        fprintf('Save search results.\n');
    else
        disp(['Load from ',sprintf(['grid_search/final%d_',method,'.mat'],cv+1)]);
        load(sprintf(['grid_search/final%d_',method,'.mat'],cv+1));
    end
%     [tmplab,tmpprob, w] = test_gestureSearch(cv, method); %Pred Gesture
    lab_vec{cv+1} = tmplab;
    prob_ges{cv+1} = tmpprob;
    prob_FV{cv+1} = tmpprob_FV;
    prob_GT{cv+1} = tmpprob_GT;
    prob_trans{cv+1} = tmpprob_trans;
    if ~exist(sprintf(['grid_search/MFCC%d_',method,'.mat'],cv+1),'file')
        disp('MFCC');
        tmpprob_MFCC = test_MFCC(cv, method);
        save(sprintf(['grid_search/MFCC%d_',method,'.mat'],cv+1), 'tmpprob_MFCC');
        fprintf('Save MFCC.\n');
    else
        disp(['Load from ',sprintf(['grid_search/MFCC%d_',method,'.mat'],cv+1)]);
        load(sprintf(['grid_search/MFCC%d_',method,'.mat'],cv+1));
    end
    prob_MFCC{cv+1} = tmpprob_MFCC;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
% Compute individual feature AUC
for cv = 0:9
    validID = ~isnan(prob_ges{cv+1}{length(prob_ges{1})});
    % IDTFV only
    [~,~,~,AUC_FV(cv+1)] = perfcurve(lab_vec{cv+1}(validID),prob_FV{cv+1}(validID),0);
    % Gesture only
    for comb_i = 1: length(prob_ges{1})
        [~,~,~,AUC_ges(cv+1, comb_i)] = perfcurve(lab_vec{cv+1}(validID),prob_ges{cv+1}{comb_i}(validID),0);
        [~,~,~,AUC_GT(cv+1, comb_i)] = perfcurve(lab_vec{cv+1}(validID),prob_GT{cv+1}{comb_i}(validID),0);
    end
    % Trans only
    [~,~,~,AUC_trans(cv+1)] = perfcurve(lab_vec{cv+1}(validID),prob_trans{cv+1}(validID),0);
    % MFCC only
    [~,~,~,AUC_MFCC(cv+1)] = perfcurve(lab_vec{cv+1}(validID),prob_MFCC{cv+1}(validID),0);

end

fprintf('IDTFV only: %f.\n', mean(AUC_FV));
fprintf('Pred Gesture only: %f.\n', max(mean(AUC_ges,1)));
fprintf('GT Gesture only: %f.\n', max(mean(AUC_GT,1)));
fprintf('Trans only: %f.\n', mean(AUC_trans));
fprintf('MFCC only: %f.\n', mean(AUC_MFCC));

tmpAUC = mean(AUC_ges,1); comb_maxi = find(tmpAUC==max(tmpAUC)); comb_maxi = comb_maxi(1); %for gesture and GT
tmpAUC = mean(AUC_GT,1); comb_maxi_GT = find(tmpAUC==max(tmpAUC)); comb_maxi_GT = comb_maxi_GT(1);

% Do a grid search to find best weights for ges_id
combs = unique(combnk([0 0 0 0 0 .5 .5 .5 .5 .5 1 1 1 1 1], 5), 'rows');
comb_i = 1;
for i = 1:size(combs,1)
    perW = unique(perms(combs(i,:)), 'rows');
    for j = 1:size(perW, 1)
        %tmpw = repmat(perW(j,:), num_C, 1);
        w{comb_i} = perW(j,:);
        comb_i = comb_i +1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find best ratio for idt and Gesture
disp('---------------------------------------------------------------------------------')
disp('IDT+ Gesture...............')
for comb_i = 1: length(prob_ges{1})
    for r  = 1:11
        ratio = -0.1+ 0.1*r;
        AUC = zeros(1, 10);
        AUC_GT = zeros(1, 10);
        for cv = 0:9
            % pred late fusion
            prob_lf= ratio*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_i})) +(1-ratio)* prob_ges{cv+1}{comb_i}(~isnan(prob_ges{cv+1}{comb_i}));
            % GT late fusion
            if ~isinf(prob_GT{cv+1}{comb_i}(~isnan(prob_ges{cv+1}{comb_i})))
            %prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_i}))
            %prob_GT{cv+1}{comb_i}(~isnan(prob_ges{cv+1}{comb_i}))
                prob_lf_GT = ratio*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_i})) +(1-ratio)* prob_GT{cv+1}{comb_i}(~isnan(prob_ges{cv+1}{comb_i}));
            else
                prob_lf_GT = ratio*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_i}));
            end
            [~,~,~,AUC(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_i})),prob_lf,0);
            [~,~,~,AUC_GT(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_i})),prob_lf_GT,0);
        end
        AUC_r(comb_i, r) = mean(AUC);
        AUC_r_GT(comb_i,r) = mean(AUC_GT);
    end
end
% find best ratio and gesture weights value, only select one if more than one are applicable
[comb_s, r_s] = find(AUC_r == max(AUC_r(:)));
[comb_s_GT, r_s_GT] = find(AUC_r_GT == max(AUC_r_GT(:)));
combs = unique(combnk([0 0 0 0 0 .5 .5 .5 .5 .5 1 1 1 1 1], 5), 'rows');
comb_i = 1;
for i = 1:size(combs,1)
    perW = unique(perms(combs(i,:)), 'rows');
    for j = 1:size(perW, 1)
        %tmpw = repmat(perW(j,:), num_C, 1);
        w{comb_i} = perW(j,:);
        comb_i = comb_i +1;
    end
end
comb_s = comb_s(1);
r_s = r_s(1);
comb_s_GT = comb_s_GT(1);
r_s_GT = r_s_GT(1);
fprintf('IDTFV+Pred: %f.\n', max(AUC_r(:)));

ratio1 = -0.1+0.1*r_s; % ratio for idt and gesture
ratio1_GT = -0.1+0.1*r_s_GT; % ratio for idt and GT gesture

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Also consider MFCC feature
disp('---------------------------------------------------------------------------------')
disp('With MFCC Feature......................')
for r = 1:11
    ratio2 = -0.1+0.1*r;
    AUC = zeros(1, 10);
    AUC_GT = zeros(1,10);
    for cv  = 0:9
        % pred+idt probs and GT+idt probs
        prob_lf0 = ratio1*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_s})) +(1-ratio1)* prob_ges{cv+1}{comb_s}(~isnan(prob_ges{cv+1}{comb_s}));
        prob_lf0_GT = ratio1_GT*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT})) +(1-ratio1_GT)* prob_GT{cv+1}{comb_s_GT}(~isnan(prob_ges{cv+1}{comb_s_GT}));
        % compute fusion
        prob_lf= ratio2*prob_lf0 + (1-ratio2)*prob_MFCC{cv+1}(~isnan(prob_ges{cv+1}{comb_s}));
        [~,~,~,AUC(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_s})),prob_lf,0);
        prob_lf_GT= ratio2*prob_lf0_GT + (1-ratio2)*prob_MFCC{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT}));
        [~,~,~,AUC_GT(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT})),prob_lf_GT,0);
    end
    AUC_lf(r) = mean(AUC(~isnan(AUC)));
    AUC_lf_GT(r) = mean(AUC_GT(~isnan(AUC_GT)));
end

fprintf('IDTFV+Pred+MFCC: %f.\n', max(AUC_lf(:)));
fprintf('IDTFV+GT+MFCC: %f.\n', max(AUC_lf_GT(:)));

r_s = find(AUC_lf == max(AUC_lf));
r_s_GT = find(AUC_lf_GT == max(AUC_lf_GT));
r_s_GT = r_s_GT(1);
r_s = r_s(1);

ratio_M = -0.1+0.1*r_s; % ratio for idt+ges and MFCC
ratio_M_GT = -0.1+0.1*r_s_GT; % ratio for idt+GT and MFCC

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Also consider Trans feature
disp('---------------------------------------------------------------------------------')
disp('With Trans Feature......................')
for r = 1:11
    ratio2 = -0.1+0.1*r;
    AUC = zeros(1, 10);
    AUC_GT = zeros(1,10);
    for cv  = 0:9
        % pred+idt probs and GT+idt probs
        prob_lf0 = ratio1*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_s})) +(1-ratio1)* prob_ges{cv+1}{comb_s}(~isnan(prob_ges{cv+1}{comb_s}));
        prob_lf0_GT = ratio1_GT*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT})) +(1-ratio1_GT)* prob_GT{cv+1}{comb_s_GT}(~isnan(prob_ges{cv+1}{comb_s_GT}));
        % compute fusion
        prob_lf= ratio2*prob_lf0 + (1-ratio2)*prob_trans{cv+1}(~isnan(prob_ges{cv+1}{comb_s}));
        [~,~,~,AUC(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_s})),prob_lf,0);
        prob_lf_GT= ratio2*prob_lf0_GT + (1-ratio2)*prob_trans{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT}));
        [~,~,~,AUC_GT(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT})),prob_lf_GT,0);
    end
    AUC_lf(r) = mean(AUC(~isnan(AUC)));
    AUC_lf_GT(r) = mean(AUC_GT(~isnan(AUC_GT)));
end

fprintf('IDTFV+Pred+trans: %f.\n', max(AUC_lf(:)));
fprintf('IDTFV+GT+trans: %f.\n', max(AUC_lf_GT(:)));

r_s = find(AUC_lf == max(AUC_lf));
r_s_GT = find(AUC_lf_GT == max(AUC_lf_GT));
r_s_GT = r_s_GT(1);
r_s = r_s(1);

ratio_T = -0.1+0.1*r_s; % ratio for idt+ges and MFCC
ratio_T_GT = -0.1+0.1*r_s_GT; % ratio for idt+GT and MFCC

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All features
disp('---------------------------------------------------------------------------------')
disp('With All Features...................')
for r = 1:11
    ratio2 = -0.1+0.1*r;
    AUC = zeros(1, 10);
    AUC_GT = zeros(1,10);
    for cv  = 0:9
        % pred gesture with
        prob_lf0 = ratio1*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_s})) +(1-ratio1)* prob_ges{cv+1}{comb_s}(~isnan(prob_ges{cv+1}{comb_s}));
        % GT gesture with trans
        prob_lf0_GT = ratio1_GT*prob_FV{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT})) +(1-ratio1_GT)* prob_GT{cv+1}{comb_s_GT}(~isnan(prob_ges{cv+1}{comb_s_GT}));
        % compute fusion with MFCC
        prob_M0 = ratio_M*prob_lf0 + (1-ratio_M)*prob_MFCC{cv+1}(~isnan(prob_ges{cv+1}{comb_s}));
        prob_M0_GT = ratio_M_GT*prob_lf0_GT + (1-ratio_M_GT)*prob_MFCC{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT}));
        % compute fusion with Trans
        prob_T= ratio2*prob_M0 + (1-ratio2)*prob_trans{cv+1}(~isnan(prob_ges{cv+1}{comb_s}));
        [~,~,~,AUC(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_s})),prob_T,0);
        prob_T_GT= ratio2*prob_M0_GT + (1-ratio2)*prob_trans{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT}));
        [~,~,~,AUC_GT(cv+1)] = perfcurve(lab_vec{cv+1}(~isnan(prob_ges{cv+1}{comb_s_GT})),prob_T_GT,0);
    end
    AUC_final(r) = mean(AUC(~isnan(AUC)));
    AUC_final_GT(r) = mean(AUC_GT(~isnan(AUC_GT)));
end

fprintf('IDTFV+Pred+MFCC+Trans: %f.\n', max(AUC_final(:)));
fprintf('IDTFV+GT+MFCC+Trans: %f.\n', max(AUC_final_GT(:)));

r_f = find(AUC_final == max(AUC_final));
r_f_GT = find(AUC_final_GT == max(AUC_final_GT));
r_f_GT = r_f_GT(1);
r_f = r_f(1);

ratio_All = -0.1+0.1*r_f; % ratio for idt+ges and MFCC
ratio_All_GT = -0.1+0.1*r_f_GT; % ratio for idt+GT and MFCC

author: Doubaibai

repo: DARE

save_path: github-repos/MATLAB/Doubaibai-DARE


text:
%%
% Test for colorspace histograms.

addpath('../toolbox/');
addpath('../colors_functions/');
addpath('../image_blur/');
addpath('../blur_functions//');
addpath('../convolutional_wasserstein/');
% addpath('../../data/images/colors/'); % low-res images
addpath('../../data/images/colors-big/'); % high-res images

% bins for histograms
N = 120;

colsp = 'XYZ';
colsp = 'LCH';  % CIE L*C*H* (CIELCH)
colsp = 'RGB';
colsp = 'LAB';

name = 'castle';
name = 'river';
name = 'room';
name = 'flowers-2';
name = 'blue-7';
f = rescale( load_image(name) );

fC = colorspace(['RGB->' colsp], f);

arange = range(fC(:,:,2));
brange = range(fC(:,:,3));

% L in [0,100], A in [-85,100], B in [-107,95]
myhist = @(x)hist(x(:), 100);
mmin = @(x)min(x(:));
mmax = @(x)max(x(:));

% compute a 2D histogram
H = compute_histogram_2d(fC(:,:,2),fC(:,:,3), N, arange,brange);

H1 = render_lab_histogram(H,arange,brange);
imageplot(H1);

author: gpeyre


text:
function [EB] = MB2EB(MB)
% Convert computery things from megabytes to exabytes.
% Chad A. Greene 2012
EB = MB*2^-40 ;

author: Sable

repo: mcbench-benchmarks