text:
function [p,t] = meshpoly(node,edge,qtree,p,options)

%*****************************************************************************80
%
% MESHPOLY: Core meshing routine called by mesh2d and meshfaces.
%
% Do not call this routine directly, use mesh2d or meshfaces instead!
%
% Inputs:
%
%  NODE     : Nx2 array of geometry XY co-ordinates
%  EDGE     : Mx2 array of connections between NODE, defining geometry
%             edges
%  QTREE    : Quadtree data structure, defining background mesh and element
%             size function
%  P        : Qx2 array of potential boundary nodes
%  OPTIONS  : Meshing options data structure
%  WBAR     : Handle to progress bar
%
% Outputs:
%
%  P        : Nx2 array of triangle nodes
%  T        : Mx3 array of triangles as indices into P
%
% Mesh2d is a delaunay based algorithm with a "Laplacian-like" smoothing
% operation built into the mesh generation process.
%
% An unbalanced quadtree decomposition is used to evaluate the element size
% distribution required to resolve the geometry. The quadtree is
% triangulated and used as a backgorund mesh to store the element size
% data.
%
% The main method attempts to optimise the node location and mesh topology
% through an iterative process. In each step a constrained delaunay
% triangulation is generated with a series of "Laplacian-like" smoothing
% operations used to improve triangle quality. Nodes are added or removed
% from the mesh to ensure the required element size distribution is
% approximated.
%
% The optimisation process generally returns well shaped meshes with no
% small angles and smooth element size variations. Mesh2d shares some
% similarities with the Distmesh code:
%
%   [1] P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB.
%       SIAM Review, Volume 46 (2), pp. 329-345, June 2004
%
%  Author:
%
%    Darren Engwirda
%

shortedge   = 0.75;
longedge    = 1.5;
smalltri    = 0.25;
largetri    = 4.0;
qlimit      = 0.5;
dt          = 0.2;

stats = struct('t_init',0.0,'t_tri',0.0,'t_inpoly',0.0,'t_edge',0.0, ...
                  't_sparse',0.0,'t_search',0.0,'t_smooth',0.0,'t_density',0.0, ...
                     'n_tri',0);

% Initialise mesh
%  P     : Initial nodes
%  T     : Initial triangulation
%  TNDX  : Enclosing triangle for each node as indices into TH
%  FIX   : Indices of FIXED nodes in P
tic
if options.output
   fprintf('Initialising Mesh \n');
end
[p,fix,tndx] = initmesh(p,qtree.p,qtree.t,qtree.h,node,edge);
stats.t_init = toc;

% Main loop
if options.output
   fprintf('Iteration   Convergence (%%)\n');
end
for iter = 1:options.maxit

   [p,i,j] = unique(p,'rows');                                             % Ensure unique node list
   fix = j(fix);
   tndx = tndx(i);

   tic
   [p,t] = cdt(p,node,edge);                                               % Constrained Delaunay triangulation
   stats.n_tri = stats.n_tri+1;
   stats.t_tri = stats.t_tri+toc;

   tic
   e = getedges(t,size(p,1));                                    % Unique edges
   stats.t_edge = stats.t_edge+toc;

   % Sparse node-to-edge connectivity matrix
   tic
   nume = size(e,1);
   S = sparse(e(:),[1:nume,1:nume],[ones(nume,1); -ones(nume,1)],size(p,1),nume);
   stats.t_sparse = stats.t_sparse+toc;

   tic
%
% Find enclosing triangle in background mesh for nodes
%
  opt = 4;
  if ( opt == 0 )
    tndx = mytsearch ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2), tndx );
  elseif ( opt == 1 )
    tndx = tsearch ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2) );
  elseif ( opt == 2 )
    tndx = tsearch_mex ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2) );
  elseif ( opt == 3 )
    dtri = DelaunayTri ( qtree.p );
    tndx = pointLocation ( dtri, p );
  elseif ( opt == 4 )
    tndx = tsearchn ( qtree.p, qtree.t, p );
  end

   hn = tinterp(qtree.p,qtree.t,qtree.h,p,tndx);  % Size function at nodes via linear interpolation
   h = 0.5*(hn(e(:,1))+hn(e(:,2)));                                        % from the background mesh. Average to edge midpoints.
   stats.t_search = stats.t_search+toc;

   edgev = p(e(:,1),:)-p(e(:,2),:);
   L = max(sqrt(sum((edgev).^2,2)),eps);                                   % Edge lengths

   % Inner smoothing sub-iterations
   time = cputime;
   move = 1.0;
   done = false;
   for subiter = 1:(iter-1)

      moveold = move;

      % Spring based smoothing
      L0 = h*sqrt(sum(L.^2)/sum(h.^2));
      F = max(L0./L-1.0,-0.1);
      F = S*(edgev.*[F,F]);
      F(fix,:) = 0.0;
      p = p+dt*F;

      % Measure convergence
      edgev = p(e(:,1),:)-p(e(:,2),:);
      L0 = max(sqrt(sum((edgev).^2,2)),eps);                               % Edge lengths
      move = norm((L0-L)./L,'inf');                                        % Percentage change
      L = L0;

      if move<options.mlim                                                 % Test convergence
         done = true;
         break
      end

   end
   stats.t_smooth = stats.t_smooth+(cputime-time);

   if options.output
      fprintf('%2i           %2.1f\n',[iter, 100.0*min(1.0, options.mlim/max(move,eps))]);
   end

   tic
   [p,t] = cdt(p,node,edge);                                               % Constrained Delaunay triangulation
   stats.n_tri = stats.n_tri+1;
   stats.t_tri = stats.t_tri+toc;

   tic
   e = getedges(t,size(p,1));                                              % Unique edges
   stats.t_edge = stats.t_edge+toc;

   edgev = p(e(:,1),:)-p(e(:,2),:);
   L = max(sqrt(sum((edgev).^2,2)),eps);                                   % Edge lengths

   tic
%
% Find enclosing triangle in background mesh for nodes
%
  opt = 4;
  if ( opt == 0 )
    tndx = mytsearch ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2), tndx );
  elseif ( opt == 1 )
    tndx = tsearch ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2) );
  elseif ( opt == 2 )
    tndx = tsearch_mex ( qtree.p(:,1), qtree.p(:,2), qtree.t, p(:,1), p(:,2) );
  elseif ( opt == 3 )
    dtri = DelaunayTri ( qtree.p );
    tndx = pointLocation ( dtri, p );
  elseif ( opt == 4 )
    tndx = tsearchn ( qtree.p, qtree.t, p );
  end

   hn = tinterp(qtree.p,qtree.t,qtree.h,p,tndx);                           % Size function at nodes via linear interpolation
   h = 0.5*(hn(e(:,1))+hn(e(:,2)));                                        % from the background mesh. Average to edge midpoints.
   stats.t_search = stats.t_search+toc;

   r = L./h;
   if done && (max(r)<3.0)                                                 % Main loop convergence
      break
   else
      if (iter==options.maxit)
         fprintf('WARNING: Maximum number of iterations reached. Solution did not converge! \n');
         fprintf('Please email the geometry and settings to d_engwirda@hotmail.com \n');
      end
   end

   % Nodal density control
   tic
   if iter<options.maxit
      % Estimate required triangle area from size function
      Ah = 0.5*tricentre(t,hn).^2;
      % Remove nodes
      i = find(abs(triarea(p,t))<smalltri*Ah);                             % Remove all nodes in triangles with small area
      k = find(sum(abs(S),2)<2);                                           % Nodes with less than 2 edge connections
      j = find(r<shortedge);                                               % Remove both nodes for short edges
      if ~isempty(j) || ~isempty(k) || ~isempty(i)
         prob = false(size(p,1),1);                                        % True for nodes to be removed
         prob(e(j,:)) = true;                                              % Edges with insufficient length
         prob(t(i,:)) = true;                                              % Triangles with insufficient area
         prob(k) = true;                                                   % Remove nodes with less than 2 edge connections
         prob(fix) = false;                                                % Don't remove fixed nodes
         pnew = p(~prob,:);                                                % New node list
         tndx = tndx(~prob);
         j = zeros(size(p,1),1);                                           % Re-index FIX to keep consistent
         j(~prob) = 1;
         j = cumsum(j);
         fix = j(fix);
      else
         pnew = p;
      end
      % Add new nodes
      i = abs(triarea(p,t))>largetri*Ah;                                   % Large triangles
      r = longest(p,t)./tricentre(t,hn);
      k = (r>longedge) & (quality(p,t)<qlimit);                            % Low quality triangles
      if any(i|k)

         k = find(k & ~i);
         i = find(i);

         % Add new nodes at circumcentres
         cc = circumcircle(p,[t(i,:); t(k,:)]);

         % Don't add multiple points in one circumcircle
         ok = [true(size(i)); false(size(k))];
         for ii = (length(i)+1):size(cc,1)
            % Current point
            x = cc(ii,1);
            y = cc(ii,2);
            % Check if inside any accepted circumcircles
            in = false;
            j = find(ok);
            for jj = 1:length(j)
               kk = j(jj);
               dx = (x-cc(kk,1))^2;
               if dx<cc(kk,3) && (dx+(y-cc(kk,2))^2)<cc(kk,3)
                  in = true;
                  break;
               end
            end
            if ~in
               ok(ii) = true;
            end
         end
         cc = cc(ok,:);
         cc = cc(inpoly(cc(:,1:2),node,edge),:);                           % Only take internal points

         % New arrays
         pnew = [pnew; cc(:,1:2)];
         tndx = [tndx; zeros(size(cc,1),1)];
      end
      p = pnew;
   end
   stats.t_density = stats.t_density+toc;

end

if options.debug
   disp(stats);
end

end      % meshpoly()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p,t] = cdt(p,node,edge)

% Approximate geometry-constrained Delaunay triangulation.

t = mydelaunayn(p);                                                        % Delaunay triangulation via QHULL

% Impose geometry constraints
i = inpoly(tricentre(t,p),node,edge);                                      % Take triangles with internal centroids
t = t(i,:);

end      % [p,t] = cdt(p,node,edge)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p,fix,tndx] = initmesh(p,ph,th,hh,node,edge)

% Initialise the mesh nodes

% Boundary nodes for all geometry edges have been passed in. Only take
% those in the current face
i = findedge(p,node,edge,1.0e-08);
p = p(i>0,:);
fix = (1:size(p,1))';

% Initial nodes taken as fixed boundary nodes + internal nodes from
% quadtree.
[i,j] = inpoly(ph,node,edge);
p = [p; ph(i&~j,:)];
tndx = zeros(size(p,1),1);

end      % initmesh()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function e = getedges(t,n)

% Get the unique edges and boundary nodes in a triangulation.

e = sortrows( sort([t(:,[1,2]); t(:,[1,3]); t(:,[2,3])],2) );
idx = all(diff(e,1)==0,2);                                                 % Find shared edges
idx = [idx;false]|[false;idx];                                             % True for all shared edges
bnd = e(~idx,:);                                                           % Boundary edges
e = e(idx,:);                                                              % Internal edges
e = [bnd; e(1:2:end-1,:)];                                                 % Unique edges and bnd edges for tri's

end      % getedges()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fc = tricentre(t,f)

% Interpolate nodal F to the centroid of the triangles T.

fc = (f(t(:,1),:)+f(t(:,2),:)+f(t(:,3),:))/3.0;

end      % tricentre()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = longest(p,t)

% Return the length of the longest edge in each triangle.

d1 = sum((p(t(:,2),:)-p(t(:,1),:)).^2,2);
d2 = sum((p(t(:,3),:)-p(t(:,2),:)).^2,2);
d3 = sum((p(t(:,1),:)-p(t(:,3),:)).^2,2);

d = sqrt(max([d1,d2,d3],[],2));

end      % longest()


text:
% function fval = f6(x)
% % Unimodal function f_1
%
% Bound = [-500 500];
%
% if nargin == 0
%     fval = Bound;
% else
%     fval = sum(-x.*sin(abs(x).^.5));
% end
function fval=f6(x)
Bound=[-100 100];

if nargin==0
    fval = Bound;
else
    x=abs(floor(x+0.5));
    fval= sum(x);
end

author: sfvsfv

repo: Mathematical-modeling


text:
%   X: data matrix, each row is one observation, each column is one feature
%   d: reduced dimension
%   Y: dimensionanlity-reduced data

%   Copyright by Quan Wang, 2011/05/10
%   Please cite: Quan Wang. Kernel Principal Component Analysis and its
%   Applications in Face Recognition and Active Shape Models.
%   arXiv:1207.3538 [cs.CV], 2012.

function Y=PCA(X,d)

%% eigenvalue analysis
Sx=cov(X);
[V,D]=eig(Sx);
eigValue=diag(D);
[eigValue,IX]=sort(eigValue,'descend');
eigVector=V(:,IX);

%% normailization
norm_eigVector=sqrt(sum(eigVector.^2));
eigVector=eigVector./repmat(norm_eigVector,size(eigVector,1),1);

%% dimensionality reduction
Y=X*eigVector(:,1:d);

author: Sable

repo: mcbench-benchmarks


text:
%
% Polynomial subtraction  h = f - g
%
% Syntax:  >> h = PolynomialPlus(f,g)
%     or:  >> h = pplus(f,g)
%
%   Input:    f --- (string or numeric) polynomial one
%             g --- (string or numeric) polynomial two
%
%  Output:    h --- (string) difference f-g
%
%  Example:  >> PolynomialMinus('5*x^2*y+8','3*x^2*y+2')
%
%            ans =
%
%            6 + 2*x^2*y

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
classdef givens

methods(Static)

function [c,s] = rotation(a,b)
    % Givens rotation computation
    % Determines cosine-sine pair (c,s) so that [c s;-s c]'*[a;b] = [r;0]
    % G = [c  s; -s c]
    % x = [a; b]
    % G' x = [r; 0]
    % GVL4: algorithm 5.1.3
    if b==0
        % No rotation needed
        c = 1; s = 0;
    else
        if abs(b)>abs(a)
            tau = -a/b;
            s = 1/sqrt(1+tau^2);
            c = s*tau;
        else
            tau = -b/a;
            c = 1/sqrt(1+tau^2);
            s = c*tau;
        end
    end
end % function

function G = planerot(x)
    % Mimics the planerot function of MATLAB
    a = x(1);
    b = x(2);
    if b==0
        % No rotation needed
        c = 1; s = 0;
    else
        sa = sign(a);
        sb = sign(b);
        if abs(b)>abs(a)
            tau = -a/b;
            s = -sb/sqrt(1+tau^2);
            c = s*tau;
        else
            tau = -b/a;
            c = sa/sqrt(1+tau^2);
            s = c*tau;
        end
    end
    G = [c, -s; s c];
end % function

function theta = theta(c, s)
    theta = atan(s/c);
end

function [Q, R] = qr(A)
    % Givens method for computing the QR factorization A = QR.
    % A is modified in the algorithm
    % GVL4: algorithm 5.2.4
    import spx.la.givens.rotation;
    [m,n] = size(A);
    % Space for saving the Q factor
    Q = eye(m,m);
    % Iterate over columns from first to last
    for j=1:n
        % Iterate over rows from last pair to first pair
        for i=m:-1:j+1
            % pick two elements from consecutive rows
            a = A(i-1,j);
            b = A(i,j);
            % Compute the needed rotation for making b 0.
            [c,s] = rotation(a,b);
            % form the Givens rotation matrix
            G = [c s;-s c];
            % Rotate the two consecutive rows based submatrix
            A(i-1:i,j:n) = G'*A(i-1:i,j:n);
            % Rotate corresponding columns in Q matrix
            % by postmultiplication with G
            Q(:,i-1:i) = Q(:,i-1:i)*G;
        end
    end
    % Return the triangular form
    R = A;
end % function

end % methods

end % classdef

author: indigits

repo: sparse-plex


text:
function [adm,admu,tdmu] = som_distortion(sM, D, arg1, arg2)

%SOM_DISTORTION Calculate distortion measure for the map.
%
% [adm,admu,tdmu] = som_distortion(sMap, D, [radius], ['prob'])
%
%  adm = som_distortion(sMap,D);
%  [adm,admu] = som_distortion(sMap,D);
%  som_show(sMap,'color',admu);
%
%  Input and output arguments:
%   sMap     (struct) a map struct
%   D        (struct) a data struct
%            (matrix) size dlen x dim, a data matrix
%   [radius] (scalar) neighborhood function radius to be used.
%                     Defaults to the last radius_fin in the
%                     trainhist field of the map struct, or 1 if
%                     that is missing.
%   ['prob'] (string) If given, this argument forces the
%                     neigborhood function values for each map
%                     unit to be normalized so that they sum to 1.
%
%   adm      (scalar) average distortion measure (sum(dm)/dlen)
%   admu     (vector) size munits x 1, average distortion in each unit
%   tdmu     (vector) size munits x 1, total distortion for each unit
%
% The distortion measure is defined as:
%                                           2
%    E = sum sum h(bmu(i),j) ||m(j) - x(i)||
%         i   j
%
% where m(i) is the ith prototype vector of SOM, x(j) is the jth data
% vector, and h(.,.) is the neighborhood function. In case of fixed
% neighborhood and discreet data, the distortion measure can be
% interpreted as the energy function of the SOM. Note, though, that
% the learning rule that follows from the distortion measure is
% different from the SOM training rule, so SOM only minimizes the
% distortion measure approximately.
%
% If the 'prob' argument is given, the distortion measure can be
% interpreted as an expected quantization error when the neighborhood
% function values give the likelyhoods of accidentally assigning
% vector j to unit i. The normal quantization error is a special case
% of this with zero incorrect assignement likelihood.
%
% NOTE: when calculating BMUs and distances, the mask of the given
%       map is used.
%
% See also SOM_QUALITY, SOM_BMUS, SOM_HITS.

% Reference: Kohonen, T., "Self-Organizing Map", 2nd ed.,
%    Springer-Verlag, Berlin, 1995, pp. 120-121.
%
%    Graepel, T., Burger, M. and Obermayer, K.,
%    "Phase Transitions in Stochastic Self-Organizing Maps",
%    Physical Review E, Vol 56, No 4, pp. 3876-3890 (1997).

% Contributed to SOM Toolbox vs2, Feb 3rd, 2000 by Juha Vesanto
% Copyright (c) by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 030200

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

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

% map
M = sM.codebook;
munits = prod(sM.topol.msize);

% data
if isstruct(D), D = D.data; end
[dlen dim] = size(D);

% arg1, arg2
rad = NaN;
normalize = 0;
if nargin>2,
  if isnumeric(arg1), rad = arg1;
  elseif ischar(arg1) && strcmp(arg1,'prob'), normalize = 0;
  end
end
if nargin>3,
  if isnumeric(arg2), rad = arg2;
  elseif ischar(arg2) && strcmp(arg2,'prob'), normalize = 0;
  end
end

% neighborhood radius
if isempty(rad) || isnan(rad),
  if ~isempty(sM.trainhist), rad = sM.trainhist(end).radius_fin;
  else rad = 1;
  end
end
if rad<eps, rad = eps; end

% neighborhood
Ud = som_unit_dists(sM.topol);
switch sM.neigh,
 case 'bubble',   H = (Ud <= rad);
 case 'gaussian', H = exp(-(Ud.^2)/(2*rad*rad));
 case 'cutgauss', H = exp(-(Ud.^2)/(2*rad*rad)) .* (Ud <= rad);
 case 'ep',       H = (1 - (Ud.^2)/rad) .* (Ud <= rad);
end
if normalize,
  for i=1:munits, H(:,i) = H(:,i)/sum(H(:,i)); end
end

% total distortion measure
mu_x_1 = ones(munits,1);
tdmu = zeros(munits,1);
hits = zeros(munits,1);
for i=1:dlen,
  x = D(i,:);                        % data sample
  known = ~isnan(x);                 % its known components
  Dx = M(:,known) - x(mu_x_1,known); % each map unit minus the vector
  dist2 = (Dx.^2)*sM.mask(known);    % squared distances
  [qerr bmu] = min(dist2);           % find BMU
  tdmu = tdmu + dist2.*H(:,bmu);     % add to distortion measure
  hits(bmu) = hits(bmu)+1;           % add to hits
end

% average distortion per unit
admu = tdmu;
ind = find(hits>0);
admu(ind) = admu(ind) ./ hits(ind);

% average distortion measure
adm = sum(tdmu)/dlen;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: ilarinieminen

repo: SOM-Toolbox


text:
function CPD = tree_CPD(varargin)
%DTREE_CPD Make a conditional prob. distrib. which is a decision/regression tree.
%
% CPD =dtree_CPD() will create an empty tree.

if nargin==0
  % This occurs if we are trying to load an object from a file.
  CPD = init_fields;
  clamp = 0;
  CPD = class(CPD, 'tree_CPD', discrete_CPD(clamp, []));
  return;
elseif isa(varargin{1}, 'tree_CPD')
  % This might occur if we are copying an object.
  CPD = varargin{1};
  return;
end

CPD = init_fields;

clamped = 0;
fam_sz = [];
CPD = class(CPD, 'tree_CPD', discrete_CPD(clamped, fam_sz));

%%%%%%%%%%%

function CPD = init_fields()
% This ensures we define the fields in the same order
% no matter whether we load an object from a file,
% or create it from scratch. (Matlab requires this.)

%init the decision tree set the root to null
CPD.tree.num_node = 0;
CPD.tree.root=1;
CPD.tree.nodes=[];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [] = ReSaveInfoWithGroundTruth( jobID, taskIDs )

data = readSeqDataFromPlainText( '../data/mocap6/' );
data = ARSeqData( 1, data);

for tt = 1:length( taskIDs )
   INFO = loadSamplerInfo( jobID, taskIDs(tt) );

   INFO.data = data;

   outpath = getUserSpecifiedPath( 'SimulationResults' );
   outpath = fullfile( outpath, num2str(jobID), num2str(taskIDs(tt)), 'Info.mat'  );
   save( outpath, '-struct', 'INFO' );

end

author: michaelchughes

repo: NPBayesHMM


text:
f = findSubjects('/silver/scr1/dti/trilinear_PPD','',{'...','mb041004','nad040610','vt040717','zs040630'});
n = length(f);
for(ii=1:n)
    dt = load(f{ii}, 'anat', 'b0', 'dt6');
    if(ii==1)
        sz = size(dt.anat.img);
        t1_a = zeros([sz(1) sz(2) n]);
        t1_s = zeros([sz(1) sz(3) n]);
        t1_c = zeros([sz(2) sz(3) n]);
        t1Avg = zeros(sz);
        sz = size(dt.b0);
        b0_a = zeros([sz(1) sz(2) n]);
        b0_s = zeros([sz(1) sz(3) n]);
        b0_c = zeros([sz(2) sz(3) n]);
        b0Avg = zeros(sz);
        dt6Avg = zeros(size(dt.dt6));
    end
    t1_a(:,:,ii) = squeeze(dt.anat.img(:,:,91));
    t1_s(:,:,ii) = squeeze(dt.anat.img(:,109,:));
    t1_c(:,:,ii) = squeeze(dt.anat.img(91,:,:));
    t1Avg = t1Avg + dt.anat.img./max(dt.anat.img(:));
    b0_a(:,:,ii) = squeeze(dt.b0(:,:,46));
    b0_s(:,:,ii) = squeeze(dt.b0(:,55,:));
    b0_c(:,:,ii) = squeeze(dt.b0(46,:,:));
    b0Avg = b0Avg + dt.b0./max(dt.b0(:));
    dt6Avg = dt6Avg + dt.dt6;
    %disp(dt.anat.xformToAcPc);
end

figure; imagesc(makeMontage(t1_a)); axis image; colormap gray;
figure; imagesc(makeMontage(t1_s)); axis image; colormap gray;
figure; imagesc(makeMontage(t1_c)); axis image; colormap gray;

figure; imagesc(makeMontage(b0_a)); axis image; colormap gray;
figure; imagesc(makeMontage(b0_s)); axis image; colormap gray;
figure; imagesc(makeMontage(b0_c)); axis image; colormap gray;

figure; imagesc(makeMontage(t1Avg)); axis image; colormap gray;
figure; imagesc(makeMontage(b0Avg)); axis image; colormap gray;
figure; imagesc(makeMontage(permute(t1Avg,[3 1 2]),[29:size(t1Avg,2)-20])); axis image xy; colormap gray;

m = makeMontage(permute(t1Avg,[3 2 1]),[32:size(t1Avg,1)-30]);
figure; imagesc(m); axis image xy; colormap gray;
m = makeMontage(permute(b0Avg,[3 2 1]),[16:size(b0Avg,1)-15]);
figure; imagesc(m); axis image xy; colormap gray;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_suite = test_confusion_matrix
% tests for cosmo_confusion_matrix
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function classes=test_confusion_matrix_basics()
    nsamples=30;
    ntargets=5;
    delta=10;

    targets=[ceil(ntargets*rand(nsamples,1));randperm(ntargets)']+delta;
    predicted=[ceil(ntargets*rand(nsamples,1));randperm(ntargets)']+delta;

    [mx,classes]=cosmo_confusion_matrix(targets,predicted);

    assertEqual(classes,delta+(1:ntargets)');

    assertEqual(size(mx),[ntargets,ntargets]);

    for k=1:ntargets
        for j=1:ntargets
            count=sum(targets==(k+delta) & predicted==(j+delta));
            assertEqual(count, mx(k,j));
        end
    end

    ds=struct();
    ds.samples=predicted;
    ds.sa.targets=targets;

    [mx2,classes2]=cosmo_confusion_matrix(ds);
    assertEqual(mx,mx2);
    assertEqual(classes,classes2);

    predicted3=predicted(randperm(numel(predicted)));
    mx3=cosmo_confusion_matrix(targets,predicted3);

    ds.samples=[predicted predicted3(:)];
    [mx_both,classes3]=cosmo_confusion_matrix(ds);
    assertEqual(mx_both,cat(3,mx2,mx3));
    assertEqual(classes3,classes);

function test_confusion_matrix_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                            cosmo_confusion_matrix(varargin{:}),'');
    % size mismatch
    aet([1;1],1);

    % missing target
    aet([1;1],[1;2]);

    % no dataset
    aet(struct());
    aet({});

    ds=struct();
    ds.samples=1;
    aet(ds,1);
    ds.sa.targets=1;
    % second argument with dataset
    aet(ds,1);

    % missing argument with numeric
    aet(1)

    % target row vector
    aet([1 1],[1;1])
    aet([1;1],[1 1])

    % no target vector
    aet(eye(2),[1;1]);

    % no target vector
    aet(ones([2 2 2]),[1;1]);
    aet([1;1], ones([2 2 2]));

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function varargout = det_internal(varargin)

switch class(varargin{1})

    case 'double'
        X = varargin{1};
        X = reshape(X,sqrt(length(X)),[]);
        varargout{1} = det(X);

    case 'char' % YALMIP send 'model' when it wants the epigraph or hypograph

        operator = struct('convexity','none','monotonicity','none','definiteness','none','model','callback');
        operator.range = [-inf inf];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
end


text:
function croppedImages = crop_all(this, rangeIndexArray, otherImages)
%Crops all given images to same range as this image
%
%   Y = MrDataNd()
%   Y.crop_all(rangeIndexArray, otherImages)
%
% This is a method of class MrDataNd.
%
% IN
%   rangeIndexArray     cell(1,2*nDimsToCrop) of pairs of
%                       dimLabel, [startSampleIndex endSampleIndex]
%                       referring to dimensions of reference image for
%                       cropping (this), e.g.,
%                       {'x', [10 100], 'y', [50 120], 'z', [10 30]}
%
%   otherImages         cell(nImages,1) of MrDataNd
%                           OR
%                       single MrDataNd object to be cropped
% OUT
%   croppedImages       cell(nImages+1,1) of equally cropped images,
%                       including reference image for cropping (this)
% EXAMPLE
%   crop_all
%
%   See also MrDataNd

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

if ~iscell(otherImages)
    otherImages = {otherImages};
end
nImages = numel(otherImages) + 1;

%% create
croppedImages = cell(nImages,1);

dimLabels = rangeIndexArray(1:2:end);
cropIndexRanges = rangeIndexArray(2:2:end);

%% convert index crop of this image into range and do the crop

nDims = numel(dimLabels);
selectionIndexArray = cell(1, 2*nDims);

%% convert cropping to current image and apply crop
croppedImages = [{this}; reshape(otherImages, [], 1)];
for iImage = 1:nImages

    %% setup selectionIndexArray for cropping, dim by dim
    for iDim = 1:nDims
        cropIndexRange = cropIndexRanges{iDim};
        dimLabel = dimLabels{iDim};
        idxDim = this.dimInfo.get_dim_index(dimLabel);
        crop = this.dimInfo.index2sample(cropIndexRange, idxDim);
        res = croppedImages{iImage}.dimInfo.resolutions;
        selectionIndexArray{2*iDim-1} = dimLabel;
        selectionIndexArray{2*iDim} = crop(1):res(idxDim):crop(2);
    end

    %% crop
    croppedImages{iImage} = croppedImages{iImage}.select('type', 'samples', ...
        selectionIndexArray{:});
end

author: translationalneuromodeling

repo: tapas


text:
function [ out ] = func_featureExtraction_filterbank( dat, varargin )

if iscell(varargin)
    opt=opt_cellToStruct(varargin{:});
elseif isstruct(varargin{:})
    opt=varargin{:};
end

for ii=1:length(dat)
    FT{ii}=func_featureExtraction(dat{ii}, {'feature',opt.feature});
end

out=FT;
end

author: PatternRecognition

repo: OpenBMI


text:
classdef MimGuiDataset < CoreBaseClass
    % MimGuiDataset. Handles the interaction between the GUI and the MIM interfaces
    %
    %
    %     You do not need to modify this file. To add new functionality, create
    %     new plguins in the Plugins and GuiPlugins folders.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        CurrentContext
        GuiDatasetState
    end

    properties (Access = private)
        Dataset
        Gui
        ModeSwitcher
        MimMain
        Reporting
        Settings
        AppDef
        ContextDef
        ManualSegmentationsChangedListener
        MarkersChangedListener
        PreviewImageChangedListener
    end

    methods
        function obj = MimGuiDataset(app_def, gui, viewer_panel, settings, reporting)
            obj.AppDef = app_def;
            obj.ContextDef = app_def.GetContextDef;
            obj.GuiDatasetState = MimGuiDatasetState();
            obj.ModeSwitcher = MimModeSwitcher(viewer_panel, obj, app_def, settings, reporting);

            obj.Gui = gui;
            obj.Reporting = reporting;
            obj.Settings = settings;
            obj.MimMain = MimMain(app_def.GetFrameworkAppDef, reporting);
            obj.AddEventListener(obj.GetImageDatabase, 'SeriesHasBeenDeleted', @obj.SeriesHasBeenDeleted);
            obj.AddEventListener(obj.ModeSwitcher, 'ModeChangedEvent', @obj.ModeHasChanged);
        end

        function ChangeMode(obj, mode)
            obj.ModeSwitcher.SwitchMode(mode, obj.Dataset, obj.GuiDatasetState.CurrentPluginInfo, obj.GuiDatasetState.CurrentPluginName, obj.GuiDatasetState.CurrentVisiblePluginName, obj.CurrentContext, obj.GuiDatasetState.CurrentSegmentationName);
        end

        function ModeAutoSave(obj)
            obj.ModeSwitcher.ModeAutoSave();
        end

        function AutoSaveMarkers(obj)
            obj.Gui.AutoSaveMarkers();
        end

        function mode = GetMode(obj)
            mode = obj.ModeSwitcher.CurrentMode;
            if isempty(mode)
                obj.Reporting.Error('MimGuiDataset::NoMode', 'The operation is not possible in this mode');
            end
        end

        function mode = GetCurrentModeName(obj)
            mode = obj.ModeSwitcher.CurrentModeString;
        end

        function mode = GetCurrentSubModeName(obj)
            mode = obj.ModeSwitcher.GetSubModeName;
        end

        function is_dataset = DatasetIsLoaded(obj)
            is_dataset = ~isempty(obj.Dataset);
        end

        function is_dataset = IsPluginResultLoaded(obj)
            is_dataset = ~isempty(obj.Dataset);
        end

        function image_database = GetImageDatabase(obj)
            image_database = obj.MimMain.GetImageDatabase();
        end

        function linked_recorder = GetLinkedRecorder(obj)
            linked_recorder = obj.MimMain.FrameworkSingleton.GetLinkedDatasetRecorder();
        end

        function [uids, patient_ids] = ImportDataRecursive(obj, folder_path)
            [uids, patient_ids] = obj.MimMain.ImportDataRecursive(folder_path);
        end

        function [sorted_paths, sorted_uids] = GetListOfPaths(obj)
            [sorted_paths, sorted_uids] = obj.MimMain.ImageDatabase.GetListOfPaths();
        end

        function template_image = GetTemplateImage(obj)
            template_image = obj.Dataset.GetTemplateImage(PTKContext.OriginalImage);
        end

        function SaveMarkers(obj, name, markers)
            if obj.DatasetIsLoaded()
                obj.Dataset.SaveMarkerPoints(name, markers);
            end
        end

        function markers = LoadMarkers(obj, name)
            markers = obj.Dataset.LoadMarkerPoints(name);
        end

        function DeleteMarkerSet(obj, name)
            if obj.DatasetIsLoaded()
                obj.Dataset.DeleteMarkerSet(name);
            end
        end

        function SaveManualSegmentation(obj, name, segmentation)
            if obj.DatasetIsLoaded()
                obj.Dataset.SaveManualSegmentation(name, segmentation);
            end
        end

        function segmentation = LoadManualSegmentation(obj, segmentation_name)
            if obj.DatasetIsLoaded()
                segmentation = obj.Dataset.LoadManualSegmentation(segmentation_name, []);
            else
                segmentation = [];
            end
        end

        function DeleteManualSegmentation(obj, name)
            if obj.DatasetIsLoaded()
                obj.Dataset.DeleteManualSegmentation(name);
                if strcmp(name, obj.GuiDatasetState.CurrentSegmentationName)
                    obj.Gui.DeleteOverlays();
                end
            end
        end

        function dataset_cache_path = GetDatasetCachePath(obj)
            if obj.DatasetIsLoaded()
                dataset_cache_path = obj.Dataset.GetDatasetCachePath();
            else
                dataset_cache_path = obj.MimMain.GetDirectories.GetCacheDirectory();
            end
        end

        function dataset_cache_path = GetEditedResultsPath(obj)
            if obj.DatasetIsLoaded()
                dataset_cache_path = obj.Dataset.GetEditedResultsPath();
            else
                dataset_cache_path = obj.MimMain.GetDirectories.GetEditedResultsDirectoryAndCreateIfNecessary();
            end
        end

        function dataset_cache_path = GetOutputPath(obj)
            if obj.DatasetIsLoaded()
                dataset_cache_path = obj.Dataset.GetOutputPath();
            else
                dataset_cache_path = obj.MimMain.GetDirectories.GetOutputDirectoryAndCreateIfNecessary();
            end
        end

        function image_info = GetImageInfo(obj)
            if obj.DatasetIsLoaded()
                image_info = obj.Dataset.GetImageInfo();
            else
                image_info = [];
            end
        end

        function patint_name = GetPatientName(obj)
            if obj.DatasetIsLoaded()
                patint_name = obj.Dataset.GetPatientName();
            else
                patint_name = [];
            end
        end

        function ClearCacheForThisDataset(obj)
            if obj.DatasetIsLoaded()
                obj.Dataset.ClearCacheForThisDataset(false);
                obj.Gui.UpdateGuiForNewDataset([]);
            end
        end

        function RefreshPlugins(obj)
            obj.Gui.RefreshPluginsForDataset(obj.Dataset)
        end

        function currently_loaded_image_UID = GetUidOfCurrentDataset(obj)
            currently_loaded_image_UID = obj.GuiDatasetState.CurrentSeriesUid;
        end

        function ClearDataset(obj)
            try
                obj.ModeSwitcher.UpdateCurrentMode([], [], [], [], []);
                obj.Gui.ClearImages;
                delete(obj.Dataset);

                obj.Dataset = [];

                obj.Settings.SetLastImageInfo([], obj.Reporting);

                obj.SetNoDataset(false);
                obj.Gui.UpdateGuiForNewDataset([]);

            catch exc
                if MimErrors.IsErrorCancel(exc.identifier)
                    obj.Reporting.ShowMessage('MimGuiDataset:LoadingCancelled', 'User cancelled');
                else
                    obj.Reporting.ShowMessage('MimGuiDataset:ClearDatasetFailed', ['Failed to clear dataset due to error: ' exc.message]);
                end
            end
        end

        function ClearDatasetKeepPatient(obj)
            try
                obj.ModeSwitcher.UpdateCurrentMode([], [], [], [], []);
                obj.Gui.ClearImages;
                delete(obj.Dataset);

                obj.Dataset = [];

                obj.Settings.SetLastImageInfo([], obj.Reporting);

                obj.SetNoDataset(true);
                obj.Gui.UpdateGuiForNewDataset([]);

            catch exc
                if MimErrors.IsErrorCancel(exc.identifier)
                    obj.Reporting.ShowMessage('MimGuiDataset:LoadingCancelled', 'User cancelled');
                else
                    obj.Reporting.ShowMessage('MimGuiDataset:ClearDatasetFailed', ['Failed to clear dataset due to error: ' exc.message]);
                end
            end
        end

        function SaveEditedResult(obj)
            obj.ModeSwitcher.SaveEditedResult();
        end

        function DeleteThisImageInfo(obj)
            obj.DeleteDatasets(obj.GetUidOfCurrentDataset());
        end

        function DeleteImageInfo(obj, uid)
            obj.DeleteDatasets(uid);
        end

        function DeleteDatasets(obj, series_uids)
            % Removes a dataset from the database and deletes its disk cache. If the dataset
            % is currently loaded then the callback from the image database will cause the
            % current dataset to be cleared.

            obj.MimMain.DeleteDatasets(series_uids);
            obj.Settings.RemoveLastPatientUid(series_uids);
            if any(strcmp(series_uids, obj.GetUidOfCurrentDataset()))
                obj.ClearDataset();
            end
        end

        function SwitchPatient(obj, patient_id)
            if ~strcmp(patient_id, obj.GuiDatasetState.CurrentPatientId)
                obj.ModeSwitcher.UpdateCurrentMode([], [], [], [], []);
                obj.Gui.SetTab(obj.AppDef.DefaultModeOnNewDataset);
                obj.Gui.ClearImages;
                obj.DeleteListeners;
                delete(obj.Dataset);
                obj.Dataset = [];
                obj.GuiDatasetState.SetPatientClearSeries(patient_id, []);
            end
        end

        function InternalLoadImages(obj, image_info_or_uid)

            % Set this to empty in case an exception is thrown before it is
            % set
            series_uid = [];

            delete_image_info = false;
            patient_id = obj.GuiDatasetState.CurrentPatientId;
            patient_visible_name = obj.GuiDatasetState.CurrentPatientVisibleName;

            try
                if isa(image_info_or_uid, 'PTKImageInfo')
                    series_uid = image_info_or_uid.ImageUid;
                    new_dataset = obj.MimMain.CreateDatasetFromInfo(image_info_or_uid);
                elseif ischar(image_info_or_uid)
                    series_uid = image_info_or_uid;
                    new_dataset = obj.MimMain.CreateDatasetFromUid(image_info_or_uid);
                else
                   new_dataset = [];
                end

                obj.ModeSwitcher.UpdateCurrentMode([], [], [], [], []);
                obj.Gui.SetTab(obj.AppDef.DefaultModeOnNewDataset);

                obj.Gui.ClearImages();
                delete(obj.Dataset);

                obj.Dataset = new_dataset;

                obj.ReplaceListeners(new_dataset);

                image_info = obj.Dataset.GetImageInfo();
                series_uid = image_info.ImageUid;
                modality = image_info.Modality;

                [preferred_context, plugin_to_use] = obj.AppDef.GetPreferredContext(modality);

                % If the modality is not CT then we load the full dataset
                load_full_data = isempty(preferred_context);

                % Attempt to obtain the region of interest
                if ~load_full_data
                    if obj.Dataset.IsContextEnabled(preferred_context)
                        try
                            new_image = obj.Dataset.GetResult(plugin_to_use);
                        catch exc
                            if MimErrors.IsErrorCancel(exc.identifier)
                                obj.Reporting.Log('LoadImages cancelled by user');
                                load_full_data = false;
                                rethrow(exc)
                            else
                                obj.Reporting.ShowMessage('MimGuiDataset:CannotGetROI', ['Unable to extract region of interest from this dataset. Error: ' exc.message]);
                                load_full_data = true;
                            end
                        end
                    else
                        load_full_data = true;
                    end
                end

                % If we couldn't obtain the ROI, we load the full dataset
                if load_full_data
                    % Force the image to be saved so that it doesn't have to be
                    % reloaded each time
                    new_image = obj.Dataset.GetResult('PTKOriginalImage', PTKContext.OriginalImage, [], true);
                end

                % The modality may not have been set if it is not contained
                % in the file metadata, but may have been set after
                % loading, in which case we update it here
                if isempty(modality) && ~isempty(new_image.Modality)
                    modality = new_image.Modality;
                    image_info.Modality = new_image.Modality;
                end

                if isfield(new_image.MetaHeader, 'PatientID')
                    patient_id = new_image.MetaHeader.PatientID;
                else
                    patient_id = series_uid;
                end

                % Update and save settings if anything has changed
                obj.Settings.SetLastImageInfo(image_info, obj.Reporting);
                obj.Settings.AddLastPatientUid(patient_id, series_uid);

                if isempty(image_info)
                    patient_visible_name = [];
                    series_name = [];
                else
                    series_info = obj.GetImageDatabase.GetSeries(image_info.ImageUid);
                    patient_info = obj.GetImageDatabase.GetPatient(series_info.PatientId);
                    patient_visible_name = patient_info.ShortVisibleName;
                    series_name = series_info.Name;
                end

                obj.GuiDatasetState.SetPatientAndSeries(patient_id, series_uid, patient_visible_name, series_name, modality);
                obj.GuiDatasetState.ClearPlugin;
                obj.Gui.SetTab(obj.AppDef.DefaultModeOnNewDataset);

                obj.Gui.UpdateGuiForNewDataset(obj.Dataset);

                % Set the image after updating the GuiState. This is necessary because setting
                % the image triggers a GUI resize, and the side panels need to be repopulated
                % first
                if load_full_data
                    obj.SetImage(new_image, PTKContext.OriginalImage);
                else
                    obj.SetImage(new_image, PTKContext.LungROI);
                end

                % Update toolbar again because setting the image will
                % change the visibility of some gui plugins
                obj.Gui.UpdateToolbar();

                obj.Gui.LoadDefaultMarkersIfRequiredWithoutProgressBar;

            catch exc
                if MimErrors.IsErrorCancel(exc.identifier)
                    obj.Reporting.ShowProgress('Cancelling load');
                    obj.ClearDataset;
                    obj.Reporting.ShowMessage('MimGuiDataset:LoadingCancelled', 'User cancelled loading');
                elseif MimErrors.IsErrorUidNotFound(exc.identifier)
                    uiwait(errordlg('This dataset is missing. It will be removed from the patient browser.', [obj.AppDef.GetName ': Cannot find dataset'], 'modal'));
                    obj.Reporting.ShowMessage('MimGuiDataset:FileNotFound', 'The original data is missing. I am removing this dataset.');
                    delete_image_info = true;
                elseif MimErrors.IsErrorFileMissing(exc.identifier)
                    uiwait(errordlg('This dataset is missing. It will be removed from the patient browser.', [obj.AppDef.GetName ': Cannot find dataset'], 'modal'));
                    obj.Reporting.ShowMessage('MimGuiDataset:FileNotFound', 'The original data is missing. I am removing this dataset.');
                    delete_image_info = true;
                elseif MimErrors.IsErrorUnknownFormat(exc.identifier)
                    uiwait(errordlg('This is not an image file or the format is not supported. It will be removed from the Patient Browser.', [obj.AppDef.GetName ': Cannot load this image'], 'modal'));
                    obj.Reporting.ShowMessage('MimGuiDataset:FormatNotSupported', 'This file format is not supported. I am removing this dataset.');
                    delete_image_info = true;
                else
                    uiwait(errordlg(exc.message, [obj.AppDef.GetName ': Cannot load dataset'], 'modal'));
                    obj.Reporting.ShowMessage('MimGuiDataset:LoadingFailed', ['Failed to load dataset due to error: ' exc.message]);
                end

                % We do this outside the catch block, in case it throws another exception
                if delete_image_info
                    try
                        % The series_uid may have been set before the
                        % exception was thrown, in which case we use this
                        % to specify which dataset to delete
                        if isempty(series_uid)
                            obj.DeleteThisImageInfo();
                        else
                            obj.DeleteImageInfo(series_uid);
                        end
                    catch exc
                        obj.Reporting.ShowMessage('MimGuiDataset:DeleteImageInfoFailed', ['Failed to delete dataset due to error: ' exc.message]);
                    end
                end
                obj.ClearDatasetKeepPatient();
            end
        end

        function RunPlugin(obj, plugin_name, context, wait_dialog)
            % Causes the GUI to run the named plugin and display the result

            if ~obj.DatasetIsLoaded()
                return;
            end

            try
                obj.RunPluginTryCatchBlock(plugin_name, context, wait_dialog);
            catch exc
                if MimErrors.IsErrorCancel(exc.identifier)
                    obj.Reporting.ShowMessage('MimGuiDataset:LoadingCancelled', ['The cancel button was clicked while the plugin ' plugin_name ' was running.']);
                else
                    obj.Reporting.ShowMessageFromException('MimGuiDataset:PluginFailed', ['The plugin ' plugin_name ' failed with the following error: ' exc.message], exc);
                    show_error_dialog = true;

                    if isa(exc, 'MimSuggestEditException')
                        default_edited_result = obj.Dataset.GetDefaultEditedResult(exc.PluginToEdit, exc.PluginContext);
                        if ~isempty(default_edited_result)
                            show_error_dialog = false;

                            choice = questdlg(['Segmentation for ' exc.PluginVisibleName ' could not be performed automatically because the algorithm was not able to process this dataset. Do you wish to create the segmentation manually using the editing tools? If you are unsure, click Cancel.'], ...
                                ['Segmentation failed for ' exc.PluginVisibleName '.'], ...
                                'Create segmentation manually', 'Cancel', 'Cancel');
                            switch choice
                                case 'Create segmentation manually'
                                    % Save the new edit
                                    obj.Dataset.SaveEditedResult(exc.PluginToEdit, default_edited_result, exc.PluginContext);

                                    % Run the plugin to load the edit
                                    % into the viewer
                                    obj.RunPluginTryCatchBlock(exc.PluginToEdit, context, wait_dialog);

                                    % Switch to edit mode
                                    obj.Gui.ChangeMode(MimModes.EditMode);

                                    uiwait(warndlg(['The segmentation created for ' exc.PluginVisibleName ' is incomplete. Please review and corect the segmentation before performing any analysis.'], ['Review and correct  ' exc.PluginVisibleName], 'modal'));
                            end
                        end

                    end

                    if show_error_dialog
                        uiwait(errordlg(['The plugin ' plugin_name ' failed with the following error: ' exc.message], [obj.AppDef.GetName ': Failure in plugin ' plugin_name], 'modal'));
                    end
                end
            end
            wait_dialog.Hide;
        end

        function InvalidateCurrentPluginResult(obj)
            % Indicates that the currently loaded result has been deleted or modified in
            % such a way that it is no longer representative of the plugin

            obj.GuiDatasetState.ClearPlugin();
            obj.UpdateModes();
        end

        function LoadManualSegmentationCallback(obj, segmentation_name)
            obj.Gui.LoadSegmentationCallback(segmentation_name);
        end

        function LoadAndDisplayManualSegmentation(obj, segmentation_name, wait_dialog)
            % Causes the GUI to run the named segmentation and display the result

            if ~obj.DatasetIsLoaded()
                return;
            end

            visible_name = segmentation_name;
            wait_dialog.ShowAndHold(['Loading segmentation ' visible_name]);

            new_image = obj.LoadManualSegmentation(segmentation_name);

            image_title = visible_name;
            image_title = ['MANUAL SEGMENTATION ', image_title];
            if isempty(new_image)
                obj.Reporting.Error('MimGuiDataset:EmptyImage', ['The segmentation ' segmentation_name ' did not return an image when expected. ']);
            end
            obj.Gui.ReplaceOverlayImageCallback(new_image, image_title);
            % We need to reset the image changed callback because we have
            % deliberately changed the image, and otherwise UpdateModes()
            % will trigger a prompt to the user to save changes
            obj.ModeSwitcher.MarkOverlayAsUnchanged();
            obj.GuiDatasetState.SetSegmentation(segmentation_name);
            obj.UpdateModes();

            wait_dialog.Hide;
        end

        function OverlayImageChanged(obj)
            obj.ModeSwitcher.OverlayImageChanged;
        end

        function UpdateModeTabControl(obj)
            obj.Gui.UpdateModeTabControl(obj.GuiDatasetState);
        end

        function SetNoDataset(obj, keep_patient)
            if keep_patient
                obj.GuiDatasetState.ClearSeries();
            else
                obj.GuiDatasetState.ClearPatientAndSeries();
            end
            obj.GuiDatasetState.ClearPlugin();
            obj.UpdateModes();
            obj.Gui.UpdateGuiForNewDataset([]);
        end

        function UpdateEditedStatus(obj, is_edited)
            obj.GuiDatasetState.UpdateEditStatus(is_edited);
        end

        function segmentation_list = GetListOfManualSegmentations(obj)
            if isempty(obj.Dataset)
                segmentation_list = CorePair.empty;
            else
                segmentation_list = obj.Dataset.GetListOfManualSegmentations();
            end
        end

        function context_list = GetAllContextsForManualSegmentations(obj)
            if isempty(obj.Dataset)
                context_list = {};
            else
                context_list = obj.Dataset.GetAllContextsForManualSegmentations();
            end
        end

        function segmentation_list = GetListOfMarkerSets(obj)
            if isempty(obj.Dataset)
                segmentation_list = CorePair.empty();
            else
                segmentation_list = obj.Dataset.GetListOfMarkerSets();
            end
        end

        function delete(obj)
            obj.DeleteListeners();
        end

        function plugin_cache = GetPluginCache(obj)
            plugin_cache = obj.MimMain.FrameworkSingleton.GetPluginInfoMemoryCache();
        end

        function preview_image = FetchPreview(obj, plugin_name)
            if isempty(obj.Dataset)
                preview_image = [];
            else
                preview_image = obj.Dataset.GetPluginPreview(plugin_name);
            end
        end

       function is_linked_dataset = IsLinkedDataset(obj, linked_name_or_uid)
            % Returns true if another dataset has been linked to this one, using
            % the name or uid specified

            if isempty(obj.Dataset)
                is_linked_dataset = false;
            else
                is_linked_dataset = obj.Dataset.IsLinkedDataset(linked_name_or_uid);
            end
       end

        function is_gas_mri = IsGasMRI(obj)
            % Check if this is a hyperpolarised gas MRI image

            if isempty(obj.Dataset)
                is_gas_mri = false;
            else
                is_gas_mri = obj.Dataset.IsGasMRI();
            end
        end

        function SaveTableAsCSV(obj, plugin_name, subfolder_name, file_name, description, table, file_dim, row_dim, col_dim, filters)
            if isempty(obj.Dataset)
                obj.Reporting.Error('MimGuiDataset:NoDatasetLoaded', ['Could not save analysis results because no dataset is currently loaded.']);
            else
                obj.Dataset.SaveTableAsCSV(plugin_name, subfolder_name, file_name, description, table, file_dim, row_dim, col_dim, filters);
            end
        end

        function ViewerPanelModeChanged(obj, new_mode)
            % This methods is called when the current mode in the viewer panel has changed

            obj.ModeSwitcher.ViewerPanelModeChanged(new_mode, obj.Dataset, obj.GuiDatasetState.CurrentPluginInfo, obj.GuiDatasetState.CurrentPluginName, obj.GuiDatasetState.CurrentVisiblePluginName, obj.CurrentContext);
        end
    end

    methods (Access = private)

        function DeleteListeners(obj)
            CoreSystemUtilities.DeleteIfValidObject(obj.PreviewImageChangedListener);
            CoreSystemUtilities.DeleteIfValidObject(obj.MarkersChangedListener);
            CoreSystemUtilities.DeleteIfValidObject(obj.ManualSegmentationsChangedListener);
            obj.PreviewImageChangedListener = [];
        end

        function ReplaceListeners(obj, new_dataset)
            obj.DeleteListeners();
            if ~isempty(new_dataset)
                obj.PreviewImageChangedListener = addlistener(new_dataset, 'PreviewImageChanged', @obj.PreviewImageChangedCallback);
                obj.MarkersChangedListener = addlistener(new_dataset, 'MarkersChanged', @obj.MarkersChangedCallback);
                obj.ManualSegmentationsChangedListener = addlistener(new_dataset, 'ManualSegmentationsChanged', @obj.ManualSegmentationsChangedCallback);
            end
        end

        function SeriesHasBeenDeleted(obj, series_uid, ~)
            % If the currently loaded dataset has been removed from the database, then clear
            % and delete
            if strcmp(series_uid, obj.GetUidOfCurrentDataset)
                obj.ClearDataset;
            end
        end

        function ModeHasChanged(obj, ~, mode)
            % Called when the mode changes
            obj.Gui.SetTabMode(mode.Data);
        end

        function result = RunPluginTryCatchBlock(obj, plugin_name, context_to_request, wait_dialog)
            new_plugin = obj.LoadPluginInfoStructure(plugin_name, obj.Reporting);
            visible_name = CoreTextUtilities.RemoveHtml(new_plugin.ButtonText);
            wait_dialog.ShowAndHold(['Computing ' visible_name]);

            if strcmp(new_plugin.PluginType, 'DoNothing')
                % Call with 2 output arguments to prevent fetching of cache
                % info
                [result, ~] = obj.Dataset.GetResultWithCacheInfo(plugin_name, context_to_request);
            else

                % Determine the context we require (full image, lung ROI,
                % etc) if not specified by the caller
                % Normally we keep the last context, but if a context plugin is
                % selected, we switch to the new context
                if isempty(context_to_request)
                    context_to_request = obj.CurrentContext;
                end
                if strcmp(new_plugin.PluginType, 'ReplaceImage')
                    context_to_request = obj.ContextDef.ChooseOutputContext(new_plugin.Context);
                end

                [result, cache_info, new_image] = obj.Dataset.GetResultWithCacheInfo(plugin_name, context_to_request);

                if isa(cache_info, 'MimCompositeResult')
                    cache_info = cache_info.GetFirstResult;
                end

                image_title = visible_name;
                if cache_info.IsEdited
                    image_title = ['EDITED ', image_title];
                end
                if strcmp(new_plugin.PluginType, 'ReplaceOverlay')
                    if isempty(new_image)
                        obj.Reporting.Error('MimGuiDataset:EmptyImage', ['The plugin ' plugin_name ' did not return an image when expected. If this plugin should not return an image, then set its PluginType property to "DoNothing"']);
                    end
                    obj.ModeSwitcher.PrePluginCall();
                    obj.Gui.ReplaceOverlayImageCallback(new_image, image_title);
                    obj.GuiDatasetState.SetPlugin(new_plugin, plugin_name, visible_name, cache_info.IsEdited);
                    obj.UpdateModes();

                elseif strcmp(new_plugin.PluginType, 'ReplaceQuiver')

                    obj.Gui.ReplaceQuiverCallback(new_image);

                elseif strcmp(new_plugin.PluginType, 'ReplaceImage')
                    obj.SetImage(new_image, context_to_request);
                end
            end
        end

        function PreviewImageChangedCallback(obj, ~, event_data)
            notify(obj.GuiDatasetState, 'PreviewImageChanged', CoreEventData(event_data.Data));
        end

        function ManualSegmentationsChangedCallback(obj, ~, event_data)
            notify(obj.GuiDatasetState, 'ManualSegmentationsChanged', CoreEventData(event_data.Data));
        end

        function MarkersChangedCallback(obj, ~, event_data)
            notify(obj.GuiDatasetState, 'MarkersChanged', CoreEventData(event_data.Data));
        end

        function SetImage(obj, new_image, context)
            obj.CurrentContext = context;
            obj.Gui.SetImage(new_image);
        end

        function UpdateModes(obj)
            obj.ModeSwitcher.UpdateCurrentMode(obj.Dataset, obj.GuiDatasetState.CurrentPluginInfo, obj.GuiDatasetState.CurrentPluginName, obj.GuiDatasetState.CurrentVisiblePluginName, obj.CurrentContext);
            obj.Gui.UpdateModeTabControl(obj.GuiDatasetState);
            obj.Gui.UpdateToolbar();
        end

        function new_plugin = LoadPluginInfoStructure(obj, plugin_name, reporting)
            % Obtains a handle to the plugin which can be used to parse its properties
            new_plugin = obj.GetPluginCache.GetPluginInfo(plugin_name, [], reporting);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function excelrecord4fcn(identified, IRFperiods, IRFt, endo, strctident, irf_estimates, n, pref)
% script excelrecord4
% records the information contained in the worksheet 'IRF' of the excel spreadsheet 'results.xls'

% create the cell that will be saved on excel
irfcell={};

% build preliminary elements: space between the tables
horzspace=repmat({''},2,5*identified);
vertspace=repmat({''},IRFperiods+3,1);

% loop over variables (vertical dimension)
for ii=1:n
tempcell={};
   % loop over shocks (horizontal dimension)
   for jj=1:identified
   % create cell of IRF record for variable ii in response to shock jj
      % if a sign restriction identification scheme has been used, use the structural shock labels
      if IRFt==4||IRFt==5||IRFt==6
      temp=['response of ' endo{ii,1} ' to ' strctident.signreslabels_shocks{jj,1} ' shocks'];
      else
      temp=['response of ' endo{ii,1} ' to ' endo{jj,1} ' shocks'];
      end
   irf_ij=[temp {''} {''} {''};{''} {''} {''} {''};{''} {'lw. bound'} {'median'} {'up. bound'};num2cell((1:IRFperiods)') num2cell((irf_estimates{ii,jj})')];
   tempcell=[tempcell irf_ij vertspace];
   end
irfcell=[irfcell;horzspace;tempcell];
end

% trim
irfcell=irfcell(3:end,1:end-1);

% write in excel
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),irfcell,'IRF','B2');

author: european-central-bank

repo: BEAR-toolbox


text:
function op = prox_max( q )

%PROX_MAX    Entry-wise maximum element.
%    OP = PROX_MAX( q ) implements the nonsmooth function
%        OP(X) = q * max( X(:) ).
%    Q is optional; if omitted, Q=1 is assumed. But if Q is supplied,
%    then it must be a positive real scalar.
% Dual: proj_simplex.m (at least if X is a vector)
% See also proj_simplex, prox_linf, prox_linf

if nargin == 0,
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) || numel( q ) ~= 1 || q <= 0,
	error( 'Argument must be positive.' );
end
op = @(varargin)prox_lmax_q( q, varargin{:} );

function [ v, x ] = prox_lmax_q( q, x, t )
if nargin < 2,
    error( 'Not enough arguments.' );
end
% We have two options when input x is a matrix:
%   Does the user want to treat it as x(:),
%   Or does the user want to treat each column separately?
% Most other functions (e.g. l1, linf) treat it as x(:)
% so that will be the default. However, we leave it
% as a hard-coded option so that the user can change it
% if they want.
% (Note: the dual function, proj_simplex, vectorizes it)
VECTORIZE = true;
% right now, we do not have the non-vectorized version implemented.

tau = max( x(:) );
if nargin == 3,
    s   = sort( nonzeros(x), 'descend' );
%     s   = sort( x, 'descend' ); % 'nonzeros' does a x(:) operation
    cs  = cumsum(s);
    ndx = find( cs - (1:numel(s))' .* [s(2:end);0] >= t * q, 1 );
    if ~isempty( ndx ),
        tau = ( cs(ndx) - t * q ) / ndx;
        x   = x .* ( tau ./ max( x, tau ) );
    else
        x(:) = 0;  % adding Oct 21
        tau = 0;
    end
elseif nargout == 2,
    error( 'This function is not differentiable.' );
end
v = q * tau;

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
%%%%%script to generate figure 4 from the paper 'Invariant Scattering Convolution Networks'
%%%%%feb 2012%%%%%%%
%%%%%%Joan Bruna and Stephane Mallat$$$$$$$$$

%%%%%This figure shows two signals with the same first-order
%%%%%scattering coefficients, but different second order.
clear all;

N=256;
foptions.J=6;
foptions.L=8;
soptions.M=2;
[Wop,filters]=wavelet_factory_2d([N N], foptions, soptions);

if 1

raster=1;
%%%%generate a synthetic triangle
tmp=zeros(4*N);
ix=1:4*N;
iix=ones(4*N,1)*ix;
iiy=iix';
rho=1;
tmp=double(iix < rho*iiy);
apert=pi/50;
angles=angle(iix+i*iiy);
radius=sqrt(iix.^2+iiy.^2);

tmp=double(angles < pi/4+apert).*(angles > pi/4-apert).*(radius < 4*N);
tmp=circshift(tmp,[N/2 N/2]);
gg=fspecial('gaussian',[9 9],1);
tmp=imfilter(tmp,gg);
[gr1,gr2]=gradient(tmp);
tmp=sqrt(gr1.^2+gr2.^2);

gg=fspecial('gaussian',[9 9],4);
tmp=imfilter(tmp,gg);
test{raster}=tmp(1:4:end,1:4:end);

[psi,phi,lp]=legacy_reshape_filters(filters, size(test{raster}));

test{raster}=ifft2(fft2(test{raster}).*(sqrt(lp{1})));
test{raster}=test{raster}-mean(test{raster}(:));
test{raster}=test{raster}/norm(test{raster}(:));
tempo=randn(size(test{raster}));
tempo=tempo-mean(tempo(:));
tempo=tempo/norm(tempo(:));

%%%% equalize white noise such that its first-order scattering
%%%% coefficients match the prescribed ones from the first image
[geq,l1f,l1g,l1eq]=equalize_first_order_scattering(test{raster},tempo,psi,phi,lp);

%%%compute scattering transform
[sc{raster}]=scat(test{raster},Wop);
fprintf('scatt done \n')

raster=raster+1;

test{raster}=geq;
%%%compute scattering transform
[sc{raster}]=scat(test{raster},Wop);
fprintf('scatt done \n')

end

dirac=zeros(size(test{1}));
dirac(1)=1;
dirac=fftshift(dirac);
[scdirac]=scat(dirac,Wop);

%%%% construct scattering display
copts.renorm_process=0;
copts.l2_renorm=1;
[out{1}]=scat_display(sc{1},scdirac,copts);
%%%% construct scattering display
[out{2}]=scat_display(sc{2},scdirac,copts);

normvalues(1)=max(max(out{1}{1}(:)),max(out{2}{1}(:)));
normvalues(2)=max(max(out{1}{2}(:)),max(out{2}{2}(:)));

%display results
normvalueseff=1.15*normvalues;
test{1}=max(0,test{1});
show_two_spectra_dlux(test{1},out{1}, .5,normvalueseff);
show_two_spectra_dlux(test{2},out{2}, .5,normvalueseff);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function region = p11_region ( )

%*****************************************************************************80
%
%% P11_REGION returns the name of the integration region for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [N, trilab, xlab] = triconncomp(tri, x)
% TRICONNCOMP  Find connected components in mesh
%
% [N, TRILAB, XLAB] = triconncomp(TRI, X)
%
%   Return the connected components found in the mesh described by TRI, X.
%
%   TRI is a 3-column matrix. Each row represents the indices of the three
%   vertices that form a triangle. TRI as a whole represents the closed
%   surface.
%
%   X is a 3-column matrix. Each row represents the Cartesian coordinates
%   of a vertex on the surface, indexed by TRI values.
%
%   N is the number of connected components found by the algorithm.
%
%   TRILAB and XLAB are vectors that label TRI and X, respectively. The
%   label indicates which connected components they belong to.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 2);
nargoutchk(0, 3);

% save triangle matrix for the labels
if (nargout > 1)
    tri0 = tri;
end

% sort triangles, so that the same vertices are nearby, and we don't need
% to create first lots of many separated connected components, and then
% merge them
tri = sort(tri, 2, 'ascend');
tri = sortrows(tri);

% init output
cc = {unique(tri(1, :))};

% assign vertices to separate connected components
for  I = 2:size(tri, 1)

    % find components that are incident to the current triangle
    idx = find(cellfun(@(x) any(intersect(tri(I, :), x)), cc));

    % if the triangle is in no component, then we have to create a new one
    if isempty(idx)
        cc{end+1} = tri(I, :);

        % if the triangle is incident to one component, we add the rest of
        % the vertices to it
    elseif (length(idx) == 1)
        cc{idx} = union(cc{idx}, tri(I, :));

        % if the triangle is incident to two components, we add the
        % triangle to one of them, and then combine them both
    elseif (length(idx) == 2)
        cc{idx(1)} = union(cc{idx(1)}, tri(I, :));
        cc{idx(1)} = union(cc{idx(1)}, cc{idx(2)});
        cc(idx(2)) = [];

        % if the triangle is incident to three components, we add the
        % triangle to one of them, and then combine all three of them
    elseif (length(idx) == 3)
        cc{idx(1)} = union(cc{idx(1)}, tri(I, :));
        cc{idx(1)} = union(cc{idx(1)}, cc{idx(2)});
        cc{idx(1)} = union(cc{idx(1)}, cc{idx(3)});
        cc(idx(2:3)) = [];

    else
        error('Assert fail: A triangle cannot have more than three vertices')
    end

end

% number of connected components
N = length(cc);

% if the user didn't ask for the labels, we don't need to waste time
% generating them
if (nargout < 3)
    return
end

% create labels for vertices
xlab = zeros(size(x, 1), 1);
for I = 1:N
    xlab(cc{I}) = I;
end

% recover the original ordering of the triangles
tri = tri0;

% create labels for triangles (the triangle will have the same label as its
% first vertex, as by construction all vertices in the triangle have the
% same label)
trilab = xlab(tri(:, 1));

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.7.2b.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, 2010
% University California, San Diego / Delft University of Technology

if strcmpi(get(handles.k,'Visible'),'on')
    if get(handles.listbox1,'Value')~=25
        if get(handles.ka,'Value')
            set(handles.k,'Enable','off');
        else
            set(handles.k,'Enable','on');
        end
    end
end

repo: cluster


text:
function check = gamma_check ( a, b, c )

%*****************************************************************************80
%
%% GAMMA_CHECK checks the parameters of the Gamma PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GAMMA_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= 0.\n' );
    fprintf ( 1, '  B = %f\n', b );
    check = 0;
    return
  end

  if ( c <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GAMMA_CHECK - Fatal error!\n' );
    fprintf ( 1, '  C <= 0.\n' );
    fprintf ( 1, '  C = %f\n', c );
    check = 0;
    return
  end

  check = 1;

  return
end


text:
function [fusedBox] = fusePETMRI(ROIbox_PET,ROIbox_MRI,maskBox_PET,MRIinv,MRIweight,wavelet)
% -------------------------------------------------------------------------
% function [fusedBox] = fusePETMRI(ROIbox_PET,ROIbox_MRI,maskBox_PET,Invert,MRIweight,wavelet)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function fuses the region of interest (ROI) of two registered PET and
% MRI volumes using a technique based on the wavelet transform. See Ref. [1]
% for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - ROIbox_PET: 3D array of the smallest box containing the ROI of the PET
%               volume.
% - ROIbox_MRI: 3D array of the smallest box containing the ROI of the MRI
%               volume.
% - maskBox_PET: 3D array of the mask of the smallest box specifying the
%                ROI of the PET volume. Voxels within the ROI are assigned
%                value of 1, voxels outside a value of 0.
% - MRIinv: String specifying if the intensities of the MRI volume are
%           inverted prior to fusion with PET. Either 'Inv' for inversion,
%           or 'NoInv' for no inversion.
% - MRIweight: Numerical value specifying the weight of the MRI scan in the
%              fusion with PET.
% - wavelet: String specifying the name of the MATLAB wavelet basis used
%            wavelet basis used in the fusion process.
% -------------------------------------------------------------------------
% OUTPUTS:
% - fusedBox: 3D array of the smallest box containing the ROI of the fused
%             PET/MRI volume. Note that the output contains random values
%             outside the ROI. These values are added in the fusion process
%             in order to accurately calculate the mean and the standard
%             deviation of the ROI (only) of the MRI volume when performing
%             Collewet normalization in the wavelet domain. Apply the mask
%             to 'fusedBox' by setting NaNs outside the ROI to accurately
%             visualize the fusion.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2013
% - Revision: May 2015
% -------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% PET PRE_PROCESSING
ROIbox_PET = sqrt(ROIbox_PET);

% RESAMPLING MRI VOLUME TO PET IN-PLANE RESOLUTION (slice spacing
% previously verified to be the same)
szPET = size(ROIbox_PET);
temp=zeros(szPET);
for i = 1:szPET(3)
    temp(:,:,i) = imresize(ROIbox_MRI(:,:,i),[szPET(1),szPET(2)],'Method','cubic','Antialiasing',true);
end
ROIbox_MRI = temp;

% NORMALIZATION (and inversion) OF VOLUMES
ROIOnly_MRI = ROIbox_MRI;
ROIOnly_PET = ROIbox_PET;
ROIOnly_MRI(maskBox_PET==0) = NaN;
ROIboxFill_MR = fillBox(ROIOnly_MRI);
ROIOnly_PET(maskBox_PET==0) = NaN;
minMRI = min(ROIOnly_MRI(:));
ROIboxFill_MR = ROIboxFill_MR - minMRI;
ROIOnly_MRI = ROIOnly_MRI - minMRI;
ROIboxFill_MR = ROIboxFill_MR./max(ROIOnly_MRI(:)).*255;
if strcmp(MRIinv,'Inv')
    ROIboxFill_MR = 255 - ROIboxFill_MR;
end
minPET = min(ROIOnly_PET(:));
ROIbox_PET = ROIbox_PET - minPET;
ROIOnly_PET = ROIOnly_PET - minPET;
ROIbox_PET = ROIbox_PET./max(ROIOnly_PET(:)).*255;

% WAVELET DECOMPOSITION AND FUSION OF COEFFICIENTS
wdecMRI = wavedec3(ROIboxFill_MR,1,wavelet);
wdecPET = wavedec3(ROIbox_PET,1,wavelet);
wdecFUSED = wdecPET;
nbcell = length(wdecMRI.dec);
for i = 1:nbcell
    wdecMRI.dec{i} = CollewetNorm(wdecMRI.dec{i});
    wdecMRI.dec{i}(isnan(wdecMRI.dec{i})) = wdecPET.dec{i}(isnan(wdecMRI.dec{i}));
    wdecFUSED.dec{i} = (1-MRIweight).*wdecPET.dec{i} + MRIweight.*wdecMRI.dec{i};
end

% WAVELET RECONSTRUCTION
fusedBox = waverec3(wdecFUSED);
fusedBox = fusedBox - min(fusedBox(:));
fusedBox = fusedBox./max(fusedBox(:)).*255;

end

author: mvallieres

repo: radiomics


text:
function femIF = genNedIFEM3DFace(mesh,fem,femI)
%% Usage: Generate Quadrature Information on Interface Faces Used in PPIFEM
%         Each interface face (triangle) is cut into 3 small triangles
%            A1                A1
%           D  E      or      E  D         equivalent
%          A2  A3            A3   A2
%
%          p = [A1;D;E;A2;A3];
%          t = [1 2 3; 2 3 4; 3 4 5]
%
% femIF.tL --- global index of left element
% femIF.tR --- global index of right element
% femIF.basL --- 3*nfi-4-4 matrix basis functions of left elem on small tri
% femIF.basR --- 3*nfi-4-4 matrix basis functions of right elem on small tri
% femIF.gx --- Gaussian x nodes on small tri
% femIF.gy --- Gaussian y nodes on small tri
% femIF.gz --- Gaussian z nodes on small tri
% femIF.gw --- Gaussian weights on small tri (3 Gaussian pt for linear IFEM)
% femIF.area --- Areas of all small tri
% femIF.normal --- unit normal vector of each small tri

% Last Modified by Xu Zhang 08/07/2020
%%
nfI = -min(mesh.fLoc); % number of interface faces
intFID = find(mesh.fLoc<0);
nfIB = size(find(mesh.f_t(intFID,2)==0),1); % # of interface faces on Boundary
nfII = nfI - nfIB; % # of internal interface faces: two tetra share it L&R
tL = zeros(3*nfII,6); tR = zeros(3*nfII,6);
tIntfLID = zeros(3*nfII,1); tIntfRID = zeros(3*nfII,1);
basL = zeros(3*nfII,6,6); basR = zeros(3*nfII,6,6);
gx = zeros(3*nfII,3); gy = zeros(3*nfII,3); gz = zeros(3*nfII,3);
A = zeros(3*nfII,1); normal = zeros(3*nfII,3);

tB = zeros(3*nfIB,6); basB = zeros(3*nfIB,6,6);
gxB = zeros(3*nfIB,3); gyB = zeros(3*nfIB,3); gzB = zeros(3*nfIB,3);
AB = zeros(3*nfIB,1); normalB = zeros(3*nfIB,3);

id = 0; idB = 0;
for i = 1:nfI
    %% Form small triangular partition on the interface face
    fID = intFID(i); % face index
    f_e = mesh.f_e(fID,:); % three surrounding edge index
    idE = find(mesh.eLoc(f_e)<0); % find index of two interface edges
    tmp = [mesh.e(f_e(idE(1)),:), mesh.e(f_e(idE(2)),:)];
    nd1 = sum(tmp) - sum(unique(tmp)); % the node of two interface edges
    nd2 = sum(mesh.e(f_e(idE(1)),:)) - nd1;
    nd3 = sum(mesh.e(f_e(idE(2)),:)) - nd1;
    p = [mesh.p(nd1,:); mesh.eIntP(-mesh.eLoc(f_e(idE(1))),:); ...
        mesh.eIntP(-mesh.eLoc(f_e(idE(2))),:); mesh.p(nd2,:); mesh.p(nd3,:)];
    t = [1 2 3; 2 3 4; 3 4 5];

    %% gx gy gz on a triangle with 3 internal point, accurate upto pd = 2
    X1 = p(t(:,1),:);  X2 = p(t(:,2),:);   X3 = p(t(:,3),:);
    G = zeros(3,9);  w1 = 2/3;  w2 = 1/6; % see gaussPtri.m
    G(:,[1,4,7]) = w1*X1 + w2*(X2+X3);
    G(:,[2,5,8]) = w1*X2 + w2*(X1+X3);
    G(:,[3,6,9]) = w1*X3 + w2*(X1+X2);

    %% triangle area on three-dimension.
    x1 = X1(:,1);  y1 = X1(:,2);   z1 = X1(:,3);
    x2 = X2(:,1);  y2 = X2(:,2);   z2 = X2(:,3);
    x3 = X3(:,1);  y3 = X3(:,2);   z3 = X3(:,3);
    AT = 1/2*(((x1-x3).*(y2-y1) - (x1-x2).*(y3-y1)).^2 + ...
        ((y1-y3).*(z2-z1) - (y1-y2).*(z3-z1)).^2 + ...
        ((z1-z3).*(x2-x1) - (z1-z2).*(x3-x1)).^2).^(1/2);

    %% Left and Right Element
    tIDL = mesh.f_t(fID,1); % element index of left element
    tIDR = mesh.f_t(fID,2); % element index of right element
    if tIDR > 0 % Internal Face
        tIDLi = -mesh.tLoc(tIDL); % intf elem index of left element
        tIDRi = -mesh.tLoc(tIDR); % intf elem index of right element

        %% Determine piece
        nd1ID = mesh.pLoc(nd1); tLp = femI.plusPC(tIDLi); tRp = femI.plusPC(tIDRi);
        if (nd1ID < 0 && tLp == 1) || (nd1ID > 0 && tLp == 2)
            basL(id+1,:,:) = femI.bas2(tIDLi,:,:);
            basL(id+2,:,:) = femI.bas1(tIDLi,:,:);
            basL(id+3,:,:) = femI.bas1(tIDLi,:,:);
        elseif (nd1ID < 0 && tLp == 2) || (nd1ID > 0 && tLp == 1)
            basL(id+1,:,:) = femI.bas1(tIDLi,:,:);
            basL(id+2,:,:) = femI.bas2(tIDLi,:,:);
            basL(id+3,:,:) = femI.bas2(tIDLi,:,:);
        end
        if (nd1ID < 0 && tRp == 1) || (nd1ID > 0 && tRp == 2)
            basR(id+1,:,:) = femI.bas2(tIDRi,:,:);
            basR(id+2,:,:) = femI.bas1(tIDRi,:,:);
            basR(id+3,:,:) = femI.bas1(tIDRi,:,:);
        elseif (nd1ID < 0 && tRp == 2) || (nd1ID > 0 && tRp == 1)
            basR(id+1,:,:) = femI.bas1(tIDRi,:,:);
            basR(id+2,:,:) = femI.bas2(tIDRi,:,:);
            basR(id+3,:,:) = femI.bas2(tIDRi,:,:);
        end

        gx(id+1:id+3,:) = G(:,1:3);
        gy(id+1:id+3,:) = G(:,4:6);
        gz(id+1:id+3,:) = G(:,7:9);

        A(id+1:id+3,:) = AT;
        normal(id+1:id+3,:) = repmat(mesh.f_norm(fID,:),3,1);

        %% tL and tR, use locID, b/c index on interface cell is different
        temp = fem.g2ldof(tIDL,:);
        %temp1 = temp(femI.elocIDvec(tIDLi,:));
        tL(id+1:id+3,:) = repmat(temp,3,1);
        tIntfLID(id+1:id+3,:) = tIDL;

        temp = fem.g2ldof(tIDR,:);
        %temp2 = temp(femI.elocIDvec(tIDRi,:));
        tR(id+1:id+3,:) = repmat(temp,3,1);
        tIntfRID(id+1:id+3,:) = tIDR;
        id = id+3;

    elseif tIDR == 0 % Boundary Face
        tIDLi = -mesh.tLoc(tIDL); % intf elem index of left element

        %% Determine piece: only one element.
        nd1ID = mesh.pLoc(nd1); tLp = femI.plusPC(tIDLi);
        if (nd1ID < 0 && tLp == 1) || (nd1ID > 0 && tLp == 2)
            basB(idB+1,:,:) = femI.bas2(tIDLi,:,:);
            basB(idB+2,:,:) = femI.bas1(tIDLi,:,:);
            basB(idB+3,:,:) = femI.bas1(tIDLi,:,:);
        elseif (nd1ID < 0 && tLp == 2) || (nd1ID > 0 && tLp == 1)
            basB(idB+1,:,:) = femI.bas1(tIDLi,:,:);
            basB(idB+2,:,:) = femI.bas2(tIDLi,:,:);
            basB(idB+3,:,:) = femI.bas2(tIDLi,:,:);
        end

        gxB(idB+1:idB+3,:) = G(:,1:3);
        gyB(idB+1:idB+3,:) = G(:,4:6);
        gzB(idB+1:idB+3,:) = G(:,7:9);

        AB(idB+1:idB+3,:) = AT;
        normalB(idB+1:idB+3,:) = repmat(mesh.f_norm(fID,:),3,1);

        %% tB use locID, b/c index on interface cell is different
        temp = fem.t(tIDL,:);
        temp1 = temp(femI.locID(tIDLi,:));
        tB(idB+1:idB+3,:) = repmat(temp1,3,1);
        idB = idB+3;
    end
end

femIF = struct('tL',tL,'tR',tR,'tB',tB,'basL',basL,'basR',basR,'basB',basB, ...
    'gx',gx,'gy',gy,'gz',gz,'gxB',gxB,'gyB',gyB,'gzB',gzB,'area',A,'areaB',AB,...
    'gw',[1/3;1/3;1/3],'normal',normal,'normalB',normalB,'tIntfLID',tIntfLID,...
    'tIntfRID',tIntfRID);

repo: ifem


text:
function [DisS, dis] = DistanceSum_MPP(X, Y, configure)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Define the weighted distance of two point processes' instances
%
% Reference:
% Iwayama, Koji, Yoshito Hirata, and Kazuyuki Aihara.
% "Definition of distance for nonlinear time series analysis
% of marked point process data."
% Physics Letters A 381.4 (2017): 257-262.
%
% Provider:
% Hongteng Xu @ Georgia Tech
% June 13, 2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dis = zeros(length(configure.weight),1);
DisS = 0;
for i = 1:length(configure.weight)
    dis(i) = Distance_MPP(X, Y, configure.M, configure.id{i});
    DisS = DisS + configure.weight(i) * dis(i);
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function kernkernVardistTest(kernType1, kernType2, numData, numIn)

% KERNKERNVARDISTTEST Description

% VARGPLVM

x = randn(numData, numIn);
x2 = randn(numData/2, numIn);

kern1 = kernCreate(x,kernType1);
kern2 = kernCreate(x,kernType2);

params1 = 0.2*randn(1,kern1.nParams)./sqrt(randn(1,kern1.nParams).^2);
kern1 = kernExpandParam(kern1, params1);

params2 = 0.2*randn(1,kern2.nParams)./sqrt(randn(1,kern2.nParams).^2);
kern2 = kernExpandParam(kern2, params2);

vardist = vardistCreate(x, numIn, 'gaussian');

params = randn(1,(numData*numIn*2));
params(1:numData*numIn) = x(:)';
vardist = vardistExpandParam(vardist, params);

Psi2 = kernkernVardistPsi2Compute(kern1, kern2, vardist, x2);
covGrad = ones(size(Psi2));
[gKern1, gKern2, gVarmeans, gVarcovars, gInd] = kernkernVardistPsi2Gradient(kern1, kern2, vardist, x2, covGrad);
%

epsilon = 1e-6;
paramskern1 = kernExtractParam(kern1);
paramskern2 = kernExtractParam(kern2);
paramsvar = vardistExtractParam(vardist);
%xx2 = x2';
params = [paramskern1 paramskern2 paramsvar x2(:)'];
origParams = params;

fprintf('Kernel hyperparameters\n');
if strcmp(kern1.type,'rbfard2')
fprintf('var: %2.6g max and min inputscale: %2.6g, %2.6g.\n',kern1.variance,max(kern1.inputScales),min(kern1.inputScales));
elseif strcmp(kern1.type,'linard2')
fprintf('max and min inputscale: %2.6g, %2.6g.\n',max(kern1.inputScales),min(kern1.inputScales));
else
  % do nothing
end

fprintf('----- Psi2 term ------- \n');
E = eig(Psi2);
fprintf('max and min eigenvalue of Psi2: %2.6g, %2.6g.\n',max(E),min(E));
origParams = params;
for i = 1:length(params);
  params = origParams;
  params(i) = origParams(i) + epsilon;
  kern1 = kernExpandParam(kern1, params(1:kern1.nParams));
  kern2 = kernExpandParam(kern2, params(kern1.nParams+1:kern1.nParams+kern2.nParams));

  vardist = vardistExpandParam(vardist, params(kern1.nParams+kern2.nParams+1:kern1.nParams+kern2.nParams+vardist.nParams));
  xx2 = params(kern1.nParams+kern2.nParams+vardist.nParams+1:end);
  x2 = reshape(xx2,[numData/2 numIn]);
  Lplus(i) = full(sum(sum(kernkernVardistPsi2Compute(kern1, kern2, vardist, x2))));
  params(i) = origParams(i) - epsilon;

  kern1 = kernExpandParam(kern1, params(1:kern1.nParams));
  kern2 = kernExpandParam(kern2, params(kern1.nParams+1:kern1.nParams+kern2.nParams));
  vardist = vardistExpandParam(vardist, params(kern1.nParams+kern2.nParams+1:kern1.nParams+kern2.nParams+vardist.nParams));
  xx2 = params(kern1.nParams+kern2.nParams+vardist.nParams+1:end);
  x2 = reshape(xx2,[numData/2 numIn]);
  Lminus(i) = full(sum(sum(kernkernVardistPsi2Compute(kern1, kern2, vardist, x2))));
end
params = origParams;
gLDiff = .5*(Lplus - Lminus)/epsilon;
g = [gKern1 gKern2 gVarmeans gVarcovars gInd];
% check firstly the kernel hyperparameters
kerndiff1 = abs(g(1:kern1.nParams) - gLDiff(1:kern1.nParams));
[g(1:kern1.nParams); gLDiff(1:kern1.nParams)]
pause

index = [kern1.nParams+1:kern1.nParams+kern2.nParams];
kerndiff2 = abs(g(kern1.nParams+1:kern1.nParams+kern2.nParams) ...
            - gLDiff(kern1.nParams+1:kern1.nParams+kern2.nParams));
[g(kern1.nParams+1:kern1.nParams+kern2.nParams); gLDiff(kern1.nParams+1:kern1.nParams+kern2.nParams)]
pause

index = [kern1.nParams+kern2.nParams+1:kern1.nParams+kern2.nParams+(vardist.nParams/2)];
varmeansdiff =  abs(g(index) - gLDiff(index));
[g(index); gLDiff(index)]
pause

index = [kern1.nParams+kern2.nParams+(vardist.nParams/2)+1:kern1.nParams+kern2.nParams+vardist.nParams];
varcovarsdiff =  abs(g(index) - gLDiff(index));
[g(index); gLDiff(index)]
pause
index = [(kern1.nParams+kern2.nParams+vardist.nParams+1):size(g,2)];
varinddiff =  abs(g(index) - gLDiff(index));
[g(index); gLDiff(index)]
pause
fprintf('Kernel1 hyps max diff: %2.6g.\n', max(kerndiff1));
fprintf('Kernel2 hyps max diff: %2.6g.\n', max(kerndiff2));
fprintf('Variational means max diff: %2.6g.\n', max(varmeansdiff));
fprintf('Variational covars max diff: %2.6g.\n', max(varcovarsdiff));
fprintf('Inducing inputs max diff: %2.6g.\n', max(varinddiff));

author: SheffieldML

repo: vargplvm


text:
function varargout = revolutionSurface(varargin)
%REVOLUTIONSURFACE Create a surface of revolution from a planar curve.
%
%   usage
%   [X Y Z] = revolutionSurface(C1, C2, N);
%   create the surface of revolution of parametrized function (xt, yt),
%   with N+1 equally spaced slices, around the Oz axis.
%   It assumed that C1 corresponds to the x coordinate, and that C2
%   corresponds to the Oz coordinate.
%
%   [X Y Z] = revolutionSurface(CURVE, N);
%   is the same, but generating curve is given in a single parameter CURVE,
%   which is a [Nx2] array of 2D points.
%
%   [X Y Z] = revolutionSurface(..., THETA)
%   where THETA is a vector, uses values of THETA for computing revolution
%   angles.
%
%   [X Y Z] = revolutionSurface(..., LINE);
%   where LINE is a 1x4 array, specifes the revolution axis in the
%   coordinate system of the curve. LINE is a row vector of 4 parameters,
%   containing [x0 y0 dx dy], where (x0,y0) is the origin of the line and
%   (dx,dy) is a direction vector of the line.
%   The resulting revolution surface still has Oz axis as symmetry axis. It
%   can be transformed using transformPoint3d function.
%   Surface can be displayed using :
%   H = surf(X, Y, Z);
%   H is a handle to the created patch.
%
%   revolutionSurface(...);
%   by itself, directly shows the created patch.
%
%   Example
%   % draws a piece of torus
%   circle = circleAsPolygon([10 0 3], 50);
%   [x y z] = revolutionSurface(circle, linspace(0, 4*pi/3, 50));
%   surf(x, y, z);
%   axis equal;
%
%
%
%   See also
%       surf, transformPoint3d, drawSphere, drawTorus, drawEllipsoid
%       surfature (on Matlab File Exchange)
%
%
%   ------
%   Author: David Legland
%   e-mail: david.legland@grignon.inra.fr
%   Created: 2004-04-09
%   Copyright 2005 INRA - CEPIA Nantes - MIAJ Jouy-en-Josas.

%   based on function cylinder from matlab
%   31/06/2006 fix bug when passing 3 parameters
%   20/04/2007 rewrite processing of input parameters, add psb to specify
%       revolution axis
%   24/10/2008 fix angle vector
%   29/07/2010 doc update

%% Initialisations

% default values

% use revolution using the full unit circle, decomposed into 24 angular
% segments (thus, some vertices correspond to particular angles 30°,
% 45°...)
theta = linspace(0, 2*pi, 25);

% use planar vertical axis as default revolution axis
revol = [0 0 0 1];

% extract the generating curve
var = varargin{1};
if size(var, 2)==1
    xt = var;
    yt = varargin{2};
    varargin(1:2) = [];
else
    xt = var(:,1);
    yt = var(:,2);
    varargin(1) = [];
end

% extract optional parameters: angles, axis of revolution
% parameters are identified from their length
while ~isempty(varargin)
    var = varargin{1};

    if length(var) == 4
        % axis of rotation in the base plane
        revol = var;

    elseif length(var) == 1
        % number of points -> create row vector of angles
        theta = linspace(0, 2*pi, var);

    else
        % use all specified angle values
        theta = var(:)';

    end
    varargin(1) = [];
end

%% Create revolution surface

% ensure length is enough
m = length(xt);
if m==1
    xt = [xt xt];
end

% ensure x and y are vertical vectors
xt = xt(:);
yt = yt(:);

% transform xt and yt to replace in the reference of the revolution axis
tra = createTranslation(-revol(1:2));
rot = createRotation(pi/2 - lineAngle(revol));
[xt, yt] = transformPoint(xt, yt, tra*rot);

% compute surface vertices
x = xt * cos(theta);
y = xt * sin(theta);
z = yt * ones(size(theta));

%% Process output arguments

% format output depending on how many output parameters are required
if nargout == 0
    % draw the revolution surface
    surf(x, y, z);

elseif nargout == 1
    % draw the surface and return a handle to the shown structure
    h = surf(x, y, z);
    varargout{1} = h;

elseif nargout == 3
    % return computed mesh
    varargout = {x, y, z};
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function computeAllModelChoice_batchHN(pathExperiments,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% function computeAllModelChoice_batchHN(pathWORK,fSetName,outcomes,freedomMat,maxOrder,nBoot,imbalance,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set selection for all feature set types
% and for all experiments with different degrees of freedom. See ref. [1,2]
% for more details.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to where all experiments need to be
%                     performed.
%                     --> Ex: /myProject/WORKSPACE/RESULTS
% 2. expNumber: Number of the single experiment to perform
%                  --> Ex: 4
% 3. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 4. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 5. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (see ref.[2]).
%               --> Ex: 'IALR'
% 6. nBatch: Number of parallel batch.
%            --> Ex: 8
% 7. matlabPATH: Full path to the MATLAB excutable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% 9. seed: Numerical number to use as seed for bootstrapping experiment
%          --> Ex: 54288
% -------------------------------------------------------------------------
% OUTPUTS: Mutivariable models are saved in a folder named 'MODELS' in the
% corresponding folder of each experiment (e.g. 'Experiment1', 'Experiment2', etc.)
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

% INITIALIZATON
time = 60; % Number of seconds to wait before checking if parallel computations are done
cd(pathExperiments), load('training')
pathSet = fullfile(pwd,'FSET'); mkdir('MODELS'), cd('MODELS'), pathModels = pwd;
mkdir('batchLog_Models'), cd('batchLog_Models'), pathBatch = pwd;
nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
for o = 1:nOutcomes
    outcomes.(nameOutcomes{o}) = training.(nameOutcomes{o}).outcome;
end
setNames = fieldnames(training.(nameOutcomes{1}).text);
[param] = batchExperiments(setNames,outcomes,nBatch); nBatch = length(param);

% PRODUCE BATCH COMPUTATIONS
save('workspace','pathSet','pathModels','training','param','maxOrder','nBoot','imbalance','seed'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    for j = 1:numel(param{i})
        fprintf(fid,['computeAllModelChoice_HN(pathSet,pathModels,training,param{',num2str(i),'}{',num2str(j),'},maxOrder,nBoot,imbalance,seed)\n']);
    end
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

cd(startpath)
end


text:
%UMFPACK_SIMPLE a simple demo
%
% Example:
%   umfpack_simple
%
% Copyright 1995-2007 by Timothy A. Davis.
%
% UMFPACK License:
%
%     Your use or distribution of UMFPACK or any modified version of
%     UMFPACK implies that you agree to this License.  UMFPACK is
%     is free software; you can redistribute it and/or
%     modify it under the terms of the GNU General Public
%     License as published by the Free Software Foundation; either
%     version 2 of the License, or (at your option) any later version.

% Availability: http://www.cise.ufl.edu/research/sparse/umfpack
%
% See also: umfpack, umfpack2, umfpack_details

help umfpack_simple

format short

A = [
 2  3  0  0  0
 3  0  4  0  6
 0 -1 -3  2  0
 0  0  1  0  0
 0  4  2  0  1
] ;
fprintf ('A = \n') ; disp (A) ;

A = sparse (A) ;

b = [8 45 -3 3 19]' ;
fprintf ('b = \n') ; disp (b) ;

fprintf ('Solution to Ax=b via UMFPACK:\n') ;
fprintf ('x1 = umfpack2 (A, ''\\'', b)\n') ;

x1 = umfpack2 (A, '\', b) ;
fprintf ('x1 = \n') ; disp (x1) ;

fprintf ('Solution to Ax=b via MATLAB:\n') ;
fprintf ('x2 = A\\b\n') ;

x2 = A\b ;
fprintf ('x2 = \n') ; disp (x2) ;

fprintf ('norm (x1-x2) should be small: %g\n', norm (x1-x2)) ;

fprintf ('Type ''umfpack_demo'' for a full demo of UMFPACK\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip