text:
function msm_to_mm_coordinate_real_general ( output_filename, a )

%*****************************************************************************80
%
%% MSM_TO_MM_COORDINATE_REAL_GENERAL writes a "matrix coordinate real general" Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
  [ nrow, ncol ] = size ( a );
  nnzeros = nnz ( a );

  fid = fopen ( output_filename, 'wt+' );

  if ( fid < 0 );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_COORDINATE_REAL_GENERAL - Fatal error!\n' );
    fprintf ( 1, '  Cannot open the output file.\n' );
    error ( 'MSM_TO_MM_COORDINATE_REAL_GENERAL - Fatal error!');
  end;

  fprintf ( fid, '%%%%MatrixMarket matrix coordinate real general\n');
  fprintf ( fid, '%%%%  Created by MSM_TO_MM_COORDINATE_REAL_GENERAL.M\n' );
  fprintf ( fid, '  %d  %d  %d\n', nrow, ncol, nnzeros );

  for j = 1 : ncol
    [ rows, temp, vals ] = find ( a(:,j) );
    sz = size ( rows );
    for k2 = 1 : sz
      fprintf ( fid, '  %d  %d  %f\n', rows(k2), j, vals(k2) );
    end

  end

  fclose ( fid );

  return
end


text:
function [stats, handles] = canlab_force_directed_graph(activationdata, varargin)
% Creates a force-directed graph from a set of variables, and plots
% clusters on 3-D brain as well if entered. Requires matlab BGL toolbox.
%
% :Usage:
% ::
%
%    canlab_force_directed_graph(activationdata OR connection matrix, ['cl', cl])
%
%
% ..
%     Author and copyright information:
%     -------------------------------------------------------------------------
%     Copyright (C) 2013  Tor Wager
%
%     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.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **activationdata:**
%        observations x variables matrix of data to be
%        inter-correlated (variables will be inter-correlated)
%
%        OR
%
%        signed, thresholded connection matrix to be used (e.g.,
%        average time-series correlation matrix across individuals
%        it's also possible to use thresholded t-values from multi-subject group analysis
%
%        If working with multi-subject time series data, best to enter connection matrix
%        as this function does not run models with participants as random
%        effects.
%
% :Optional Inputs: Enter keyword followed by variable with values
%
%   **'cl':**
%        followed by clusters or region structure with brain clusters
%
%   **'threshtype':**
%        followed by threshold type; 'bonf' [default] or 'fdr'
%
%   **'connectmetric':**
%        followed by node connection metric 'corr' or 'partial_corr'
%
%   **'sizescale':**
%        Followed by values to use in sizing of nodes on graph
%        'linear' 'sigmoidal' [default] or 'custom'
%
%   ** 'sizes' **
%       Followed by vector of point sizes for all points - works ONLY if sizescale
%       is set to 'custom'
%
%   **'rset', 'partitions':**
%        Cell of vectors, with indices (integers) of member
%        elements in each group, [1 x g] cell
%        rset can ALSO be a vector of integers, i.e., output
%        from clusterdata
%
%   **'setcolors', 'partitioncolors':**
%        Cell array of colors for each group, [1 x g]
%
%   **'names':**
%       followed by cell array of names for each region/object
%
%   **'namesfield':**
%       followed by name of field to extract region/object names from in
%       clusters structure
%
%   **'linewidth':**
%       followed by line width
%
%   **'linestyle':**
%       followed by line style. Default is 'curved', anything else is
%       'straight'
%
% :Output:
%
%   **stats:**
%        structure with descriptive statistics, including
%        betweenness-centrality, degree of each node
%
% :Examples:
% ::
%
%    [stats, handles] = canlab_force_directed_graph(activationdata, 'cl', cl, 'namesfield', 'shorttitle');
%    [stats, handles] = canlab_force_directed_graph(activationdata, 'cl', cl, 'namesfield', 'shorttitle', 'degree');
%    [stats, handles] = canlab_force_directed_graph(activationdata, 'cl', cl, 'namesfield', 'shorttitle', 'degree', 'rset', rset, 'setcolors', setcolors);
%
%    [stats, handles] = canlab_force_directed_graph(G, 'sizescale', 'custom', 'sizes', [6*ones(17, 1); 12*ones(7, 1)]);

% ..
%    DEFAULTS AND INPUTS
% ..

shan = [];              % outputs
spherehan = [];
handles = [];

cl = [];
threshtype = 'bonf';
connectmetric = 'corr';  % or partial_corr
sizescale = 'sigmoid';
sizes = ones(min(size(activationdata)), 1);
linestyle = 'curved';

setcolors = [];         % control of color subgroups
rset = [];

ptsizetype = 'bc';
names = [];
namesfield = [];        % enter field name, e.g., 'shorttitle'
linewidth = 1;

dofigure = true;

% Variable arg inputs
% -----------------------------------

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch lower(varargin{i})
            % reserved keywords
            case 'degree', ptsizetype = 'degree';
                %case 'design'

                % functional commands
            case 'cl', cl = varargin{i + 1}; varargin{i + 1} = [];

            case 'linewidth', linewidth = varargin{i + 1}; varargin{i + 1} = [];

            case 'linestyle', linestyle = varargin{i + 1}; varargin{i + 1} = [];

            case 'partitions', rset = varargin{i + 1}; varargin{i + 1} = [];

            case 'partitioncolors', setcolors = varargin{i + 1}; varargin{i + 1} = [];

            case {'threshtype' 'connectmetric' 'sizescale' 'setcolors' 'rset' 'names' 'namesfield' 'sizes' 'fdr'}
                eval([varargin{i} ' = varargin{i + 1}; varargin{i + 1} = [];'])

            case {'nofigure', 'nofig'}
                dofigure = false;

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% Check vars, etc.
% ---------------------------------------------------------------------
% if isdir('/Users/tor/Documents/matlab_code_external/matlab_bgl')
%     g = genpath('/Users/tor/Documents/matlab_code_external/matlab_bgl');
%     addpath(g);
% end
%
% if ~exist('fruchterman_reingold_force_directed_layout.m', 'file')
%     error('Must have Matlab BGL toolbox on path (external toolbox)');
% end

activationdata = double(activationdata);

if ~isempty(cl) && (size(activationdata, 2) ~= length(cl))
    disp('activationdata must have as many columns as clusters has elements.')
    error('check activationdata and clusters to make sure they match.')
end

if isempty(rset)
    rset = {1:size(activationdata, 2)};
end

if isempty(setcolors)
    setcolors = scn_standard_colors(length(rset));
end

if ~isempty(namesfield)
    names = {cl.(namesfield)}';
end

if ~iscell(rset) % then integer vector
    for i = 1:max(rset)
        rset2{i} = find(rset == i);
    end
    rset = rset2;
end

% -------------------------------------------------------------------------
% CALCULATIONS
% -------------------------------------------------------------------------
if issymmetric(activationdata)
    fprintf('Thresholded association matrix detected.\n')
    [b, C, r, sig] = deal(activationdata);
    thr = NaN;

else
    fprintf('activationdata appears to be raw data. Running inter-correlations. See also xcorr_multisubject.\n')

     switch connectmetric
        case 'partial_corr'

            % Partial correlation slopes (r is actually slope)
            [r, p] = calc_partial_r(activationdata);

            b = r;

         case 'corr' % default

            [r, p] = corr(double(activationdata));

            b = [];

         otherwise error('Unknown connectmetric');
     end

     sz = size(r, 1);

     switch threshtype
         case 'bonf'
             thr = .05 ./ (sz * (sz - 1) / 2);  % bonferroni...

         case 'fdr'

             thr = p_matrix2fdrthresh(p);

     end

     C = r;

     % Threshold
     sig = sparse(double(p < thr & r > 0));
     signeg = sparse(-double(p < thr & r < 0));
     sig = sig + signeg;

     C(~sig) = 0;

     % Threshold based on significant partial regression effects
     % C is connectivity matrix for graph

    % end if issymmetric
end

% Enforce format for graph
C = (C + C') ./ 2;
C = sparse(C);

% Graph Stats
% -----------------------------------------------

% 2020: Use Matlab graph object
G = graph(C);

%bc = betweenness_centrality(abs(C));  % abs if ignoring neg connections
bc = centrality(G, 'betweenness');  % 'Importance', IMP

% shortest paths: used as estimate of connectivity
% [D S] = mean_path_length(r, rset);
% shortest path using POSITIVE edges only
GG = G;
GG.Edges.Weight( G.Edges.Weight < 0) = 0;

[D S] = mean_path_length(GG, rset);

deg = full(sum(C ~= 0))';

stats.names = names;
stats.rset = rset;
stats.C = C;
stats.r = r;
stats.b = b;

if exist('p', 'var')
    stats.p = p;
end

if exist('sig', 'var')
    stats.sig = sig;
end

stats.thr = thr;
stats.betweenness = bc;
%stats.path_length_distance = D;
stats.degree = deg;
%stats.mean_path_by_rset = S;

switch ptsizetype
    case 'bc'
        ptsizemetric = bc;
    case 'degree'
        ptsizemetric = deg;
    case 'none'
        ptsizemetric = ones(size(bc));
    otherwise error('Unknown ptsizetype');
end
if strcmp(sizescale,'custom')
   fprintf('Node size reflects custom scaling\n');
else
    fprintf('Node size reflects %s\n', ptsizetype);
end
% -------------------------------------------------------------------------
% Force-directed graph
% -------------------------------------------------------------------------

% Xc is coordinate matrix
% Xc = fruchterman_reingold_force_directed_layout(C);

han = plot(G, 'Layout', 'force');   % Get layout
Xc = [han.XData; han.YData]';
delete(han);                        % erase

% could rotate here
% ***

if dofigure
    if ~isempty(cl)
        create_figure('graph', 1, 2);
    else
        create_figure('graph');
    end
end

switch lower(linestyle)

    case 'curved'

        han = nmdsfig_tools('drawlines',Xc, sig, [.5 .5 .5; 0 1 1],[{'-'} {':'}], .1);
        handles.lh_pos = han.hhp;
        handles.lh_neg = han.hhn;% Dots
        %set(handles.lh_pos, 'Color', [.5 .5 .5])
        set(handles.lh_neg, 'Color', [0 .5 1], 'LineStyle', ':')

    otherwise

        gplot(sig>0, Xc, 'k');
        axis image
        lh = findobj(gca, 'Type', 'line');
        set(lh, 'Color', [.5 .5 .5])

        handles.lh_pos = lh;

        gplot(sig<0, Xc, 'b');
        lh = findobj(gca, 'Type', 'line', 'color', 'b');
        set(lh, 'Color', [0 .5 1], 'LineStyle', ':')

        handles.lh_neg = lh;

end

switch sizescale
    case 'linear'
        %         sizes = 4 + 15 * intensity.(cnames{i}) ./max(intensity.(cnames{i}));  % LC none of these variables exist in the workspace
        sizes = 4 + 15 * ptsizemetric ./ max(ptsizemetric);

    case 'sigmoid'
        % Rationale: avoids some VERY large points in areas highly
        % focused on as a priori ROIs, with extreme z-scores in
        % intensity relative to other areas.

        sizes = sigmoidscale(ptsizemetric);

    case 'custom'
        % we already have sizes
end

ph = [];
texth = {};

for k = 1:length(rset)

    texth{k} = [];

    for j = 1:length(rset{k})

        ph(end+1) = plot(Xc(rset{k}(j), 1), Xc(rset{k}(j), 2), 'o', 'Color', setcolors{k} ./ 2, 'MarkerSize', sizes(rset{k}(j)), 'MarkerFaceColor', setcolors{k}, 'LineWidth', linewidth);

        if ~isempty(names)
            offset = .02 * range(get(gca, 'XLim'));
            texth{k}(end+1) = text(Xc(rset{k}(j), 1)+offset, Xc(rset{k}(j), 2), strrep(names{rset{k}(j)}, '_', ' '), 'Color', 'k', 'FontSize', 14);
        end
    end

end

handles.ph = ph;
handles.texth = texth;

xlim = get(gca, 'XLim');  rg = range(xlim) * [-.05 .05]; xlim = xlim + rg;
ylim = get(gca, 'YLim'); rg = range(ylim) * [-.05 .05]; ylim = ylim + rg;
set(gca, 'XLim', xlim, 'YLim', ylim);
axis off
drawnow

stats.Xc = Xc;
stats.sizes = sizes;

if isempty(cl), return, end

% -------------------------------------------------------------------------
% Subcortical surface
% Only if cl is entered
% -------------------------------------------------------------------------

% NOTE: you need 3dHeadUtilityLite on your matlab path.
subplot(1, 2, 2)

% regioncenters
xyz = cat(1, cl.mm_center);

DB = struct('xyz', xyz, 'x', xyz(:, 1), 'y', xyz(:, 2), 'z', xyz(:, 3));

% sizes = sigmoidscale(ptsizemetric, 2, 6);
sizes = sizes ./ 2.5; % rescale to avoid huge ones

% Make Brain with Spheres
% ------------------------------------------------------------

[shan, spherehan] = connectivity3dbrain(xyz, rset, sizes, setcolors, names);

% Add lines
% ------------------------------------------------------------

[~, linehandles] = cluster_nmdsfig_glassbrain(cl,ones(length(cl), 1), {[.5 .5 .5]}, sig, [], 'samefig', 'nobrain', 'noblobs', linestyle);
set(linehandles, 'LineWidth', 2)

drawnow
%saveas(gcf, fullfile(savefigdir, ['graph3d_v2_' cnames{i} '.png']));

handles.surfhan = shan;
handles.spherehan = spherehan;
handles.linehandles = linehandles;

end % function

function [b p] = calc_partial_r(activationdata)

X = activationdata;
X = zscore(X);
%X(isnan(X)) = mean(X(~is;

clear b p

for i = 1:size(X, 2)
    y = X(:, i);
    xx = X;
    xx(:, i) = 1;  % intercept; need it, and also placeholder

    [bb, dev, statsglm] = glmfit(xx, y, 'normal', 'constant', 'off');

    b(:, i) = bb;
    p(:, i) = statsglm.p;
    p(i, i) = 1;
    b(i, i) = NaN;

end

end % function

function sizes = sigmoidscale(sizes, varargin)
% sizes = sigmoidscale(sizes, [lower bound], [upper bound])
%
% rescales a vector based on sigmoid function of zcore(input values)

% Rationale: avoids some VERY large points in areas highly
% focused on as a priori ROIs, with extreme z-scores in
% intensity relative to other areas.

% scale size of nodes - intensity, or betweenness
sizes = zscore(sizes);

A = 4; % lower asymptote
K = 15; % upper asymptote

if nargin > 2
    A = varargin{1};
    K = varargin{2};
end

B = 2;   % growth rate
v = .5;  % high-growth asymptote
Q = .5;
M = .5;  % time of max growth, if v = Q

richards = @(x, A, K, B, v, Q, M) A + (K - A) ./ ((1+Q*exp(-B*(x-M))).^(1/v));

%figure; plot(sort(sizes), richards(sort(sizes), A,K,B,v,Q,M));
sizes = richards(sizes, A,K,B,v,Q,M);

end

function [C S] = mean_path_length(G, rset)
% compute matrix of 1/path lengths for n x n matrix of regions, C (connectivity)
% and 1 / mean path length for sets of regions specified by rset
%
% r is a correlation matrix
% rset is a cell array of length k, for k sets, with vectors describing the
% indices of members of each set.
%
% e.g., r = region_r{i};

% 2020: Use Matlab Graph object
%
% r(isnan(r)) = 0;
% r = (r' + r) ./ 2;  % enforce symmetry, just in case
%
%
% % shortest paths: used to calculate connectivity
% rtmp = sparse(r);
% rtmp(rtmp < 0) = 0;

C = distances(G);

% OLDNote: Uses Floyd-Warshall method if 10% non-zero elements or more
% C = all_shortest_paths(rtmp);

% for each set, average elements of D corresponding to pairs of SETS
% return averages in S
%
% e.g., if D is min path length, S is the average min path length for Set 1
% vs. 2, 1 vs. 3, 2 vs. 3, etc.

% for inf values, we must impute some finite value or every average is Inf
% (i.e., if any regions are unconnected)
% impute max:
% but maybe better to work on unthresholded r matrix.
% Or, even better, return 1/S, 1 / mean path length, so unconnected
% regions/sets get 0.
%  mx = max(D(~isinf(D)));
%  D(isinf(D)) = mx;

S = zeros(length(rset));

for m = 1:length(rset)-1
    for n = m:length(rset)

        % set m to set n relationships (e.g., shortest path lengths)
        vals = 1 ./ C(rset{m}, rset{n});

        % average, excluding "self-connections" with value Inf
        S(m, n) = mean(vals(~isinf(vals)));

    end
end

S = S + S';

end

function [shan, spherehan] = connectivity3dbrain(regioncenters, rset, sizes, setcolors, names)

% Colors
% ----------------------------------------------------------

ctxcolors = cell(1, size(regioncenters, 1));

for i = 1:length(rset)

    for j = 1:length(rset{i})

        wh = rset{i}(j);

        ctxcolors(rset{i}) = setcolors(i);

    end
end

% Brain
% ----------------------------------------------------------

%create_figure('subcortex');

shan = addbrain('limbic');
delete(shan(end));
shan = shan(1:end-1);

shan = [shan addbrain('hires left')];
shan = [shan addbrain('brainstem')];

set(shan(end-1), 'FaceColor', [.5 .5 .5], 'FaceAlpha', .15);
set(shan(1:end-2), 'FaceColor', [.5 .5 .5], 'FaceAlpha', .3);
set(shan(end), 'FaceColor', [.5 .5 .5], 'FaceAlpha', .3);

view(89, 1);
lightRestoreSingle;

lighting gouraud;

% Spheres
% ----------------------------------------------------------

% add spheres
% cortex is special, because it involves different colors

% xyz = regioncenters(rset{1}, :);
% sz = sizes(rset{1});
% if ~isempty(names)
%     mynames = names(rset{1});
% end

xyz = regioncenters;
sz = sizes;
if ~isempty(names)
    mynames = names;
end

for i = 1:size(xyz, 1)
    spherehan(i) = cluster_image_sphere(xyz(i, :), 'color', ctxcolors{i}, 'radius', sz(i));

    if ~isempty(names)
        offset = .04 * range(get(gca, 'XLim'));
        text(xyz(i, 1)+offset, xyz(i, 2)+offset, xyz(i, 3)+offset, mynames{i}, 'Color', 'k', 'FontSize', 14);
    end
end
spherehan = {spherehan};

% for i = 2:length(rset)
%
%     xyz = regioncenters(rset{i}, :);
%     sz = sizes(rset{i});
%     if ~isempty(names)
%     mynames = names(rset{i});
%     end
%
%     spherehan{i} = cluster_image_sphere(xyz, 'color', setcolors{i}, 'radius', sz);
%
%     if ~isempty(names)
%         offset = .02 * range(get(gca, 'XLim'));
%         for j = 1:size(xyz, 1)
%             text(xyz(j, 1)+offset, xyz(j, 2)+offset, xyz(j, 3)+offset, mynames{j}, 'Color', 'k', 'FontSize', 14);
%         end
%     end
%
% end

end % function

function thr = p_matrix2fdrthresh(p)

trilp = tril(p, -1);
wh = logical(tril(ones(size(p)), -1)); % Select off-diagonal values (lower triangle)
trilp = double(trilp(wh));             % Vectorize and enforce double
trilp(trilp < 10*eps) = 10*eps;        % Avoid exactly zero values
fdrthr = FDR(trilp, 0.05);

if isempty(fdrthr), fdrthr = -Inf; end
thr = fdrthr;

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [lb, ub] = getMultiFlybyLBUB(lOpenUT, lDuration, bodiesInfo, celBodyData)
%getMultiFlybyLBUB Summary of this function goes here
%   Detailed explanation goes here

    nvars = length(bodiesInfo)+(length(bodiesInfo)-1);
    IntCon = 1:nvars;
    IntCon = IntCon(end-(length(bodiesInfo)-1)+1:end);
    lb = zeros(1,nvars);
    ub = zeros(1,nvars);

    bodyInfo = bodiesInfo{1};
    parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
    xferGmu = parentBodyInfo.gm;

    lb(1) = lOpenUT;
    ub(1) = lOpenUT + lDuration;

    for(i=1:length(bodiesInfo)-1) %#ok<*NO4LP>
        b1 = bodiesInfo{i};
        b2 = bodiesInfo{i+1};

        mSMA = mean([b1.sma, b2.sma]);
        xP = computePeriod(mSMA, xferGmu);

        lb(i+1) = xP/25;
        ub(i+1) = xP*5;
    end

    lb(IntCon) = 1;
    ub(IntCon) = 2;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% TOEPLITZ_BLOCK - Constructs a block-Toeplitz matrix.
%
%   Y = TOEPLITZ_BLOCK(X)
%
% X = [X_1, X_2, X_3, ..., X_N]
%
% Y = [X_1 X_2 X_3 ... X_N
%      X_2 X_1 X_2 ... ...
%      X_3 X_2 X_1 ... ...
%      ... ... ... ... ...
%      X_N ... X_3 X_2 X_1]

function Y = toeplitz_block(X)

[M,N] = size(X);

if mod(N,M) ~= 0
  error('The number of columns must be a multiple of the number of rows.')
end

D = N/M;

if issparse(X)
  t = cputime();

  nonzeros = zeros(D,1);
  i = cell(D,1);
  j = cell(D,1);
  v = cell(D,1);
  [i{1},j{1},v{1}] = find(X(:,1:M));
  nonzeros(1) = D*length(v{1});
  for d=2:D
    k = (d-1)*M+1;
    l = d*M;
    [i{d},j{d},v{d}] = find(X(:,k:l));
    nonzeros(d) = 2 * (D-d+1) * length(v{d});
  end
  nzs = sum(nonzeros);

  %time = cputime() - t, t = cputime();

  I = zeros(nzs,1);
  J = zeros(nzs,1);
  values = zeros(nzs,1);

  %time = cputime() - t, t = cputime();

  z = 1;
  for d2=1:D
    for d1=1:D
      ind = abs(d1-d2) + 1;
      jnd = z:(z+length(v{ind})-1);
      I(jnd) = i{ind} + (d1-1)*M;
      J(jnd) = j{ind} + (d2-1)*M;
      values(jnd) = v{ind};
      z = z + length(v{ind});
    end
  end

  %time = cputime() - t, t = cputime();

  Y = sparse(I,J,values,N,N,nzs);

  %time = cputime() - t, t = cputime();

else

  Y = zeros(N,N);
  for i=1:D
    for j=1:D
      ind = abs(i-j) + 1;
      k = (i-1)*M+1;
      l = (j-1)*M+1;
      m = (ind-1)*M+1;
      Y(k:(k+M-1),l:(l+M-1)) = X(:,m:(m+M-1));
    end
  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function magnifier
% MAGNIFIER - Pixel Magnification for any graphic.
%
% The figure window created by the MAGNIFIER command will display a
% magnification of pixels of the area under the mouse.  To use it,
% click in the MAGNIFIER window, and drag the mouse to the area of
% interest.

% Written by Eric Ludlam <eludlam@mathworks.com>
% Copyright 2002 The MathWorks Inc

  ad.f=figure('renderer','painters','menubar','none',...
              'busyaction','cancel','doublebuffer','off','tag','magnifier',...
              'name','Magnifier','numbertitle','off');
  ad.a=axes('parent',ad.f);

  f=getframe(ad.f, [ get(ad.f,'currentpoint') 100 100 ] );
  ad.i = image(f.cdata,'parent',ad.a,'erasemode','none');
  set(ad.a, 'color','none','xticklabel',[],'yticklabel',[],'box','on',...
             'position',[.05 .05 .9 .9],'layer','bottom');

  title('Click and drag to magnify');
  set(ad.f,'windowbuttondownfcn',@magnifybuttondown);
  set(ad.f,'windowbuttonupfcn',@magnifybuttonup);

  set(ad.f ,'handlevisibility','callback');
  setappdata(ad.f,'magnifier',ad);

function magnifymotion(a,b)

  global inmotion;

  if isempty(inmotion)

    inmotion = 1;

    ad=getappdata(a,'magnifier');
    cp = get(ad.f,'currentpoint');

    f = getframe(ad.f, [ cp-50 100 100 ]);

    set(ad.i,'cdata',f.cdata);

    inmotion = [];

  end

function magnifybuttondown(a,b)
% Button down function in the magnifier

  ad=getappdata(a,'magnifier');
  ad.down = 1;
  set(ad.f,'windowbuttonmotionfcn',@magnifymotion)
  setappdata(ad.f,'magnifier',ad);

function magnifybuttonup(a,b)

  ad=getappdata(a,'magnifier');
  ad.down=0;
  set(ad.f,'windowbuttonmotionfcn',[])
  setappdata(ad.f,'magnifier',ad);

author: Sable

repo: mcbench-benchmarks


text:
%% set up path

addpath ../misc;
addpath ../filehandling

%rootdatadir = '../experiments/wildtype/results';
% if ispc,
%   rootdatadir = 'C:\Data\JAABA\groundtruth_pBDPGAL4U_data';
% else
%   rootdatadir = '/groups/branson/bransonlab/projects/JAABA/data/groundtruth_pBDPGAL4U_data';
% end
rootdatadir='/groups/branson/home/riveraalbam/Desktop/experiments';
%% parameters

% experiment_names = {
%   'EXT_CSMH_None_Rig1Plate15BowlA_20120519T170213'
%   'FCF_attP2_1500062_None_Rig1Plate15BowlA_20120519T172815'
%   'FCF_cantons_1500002_None_Rig1Plate15BowlA_20120519T160453'
%   };
% experiment_names = {
%   'pBDPGAL4U_TrpA_Rig1Plate10BowlB_20110609T091905'
%   'pBDPGAL4U_TrpA_Rig1Plate15BowlD_20120210T152130'
%   'pBDPGAL4U_TrpA_Rig2Plate17BowlA_20120202T145723'
%   'pBDPGAL4U_TrpA_Rig2Plate17BowlC_20110811T133805'
%   'pBDPGAL4U_TrpA_Rig2Plate17BowlD_20111130T094443'
% };
% experiment_names = {
%   'GMR_14C07_AE_01_TrpA_Rig1Plate15BowlB_20120404T141159'
%   'GMR_14C07_AE_01_TrpA_Rig2Plate14BowlD_20110511T152004'
%   'GMR_34H05_AE_01_TrpA_Rig1Plate10BowlC_20110623T090005'
%   'GMR_34H05_AE_01_TrpA_Rig1Plate15BowlA_20120222T145008'
%   'GMR_71G01_AE_01_TrpA_Rig2Plate14BowlA_20110707T154658'
%   'GMR_71G01_AE_01_TrpA_Rig2Plate14BowlD_20110707T154929'
%   'GMR_71G01_AE_01_TrpA_Rig2Plate17BowlD_20110916T155353'
% };
experiment_names={'/CantonS/CantonSA' '/Dvirilis/virilisA' 'Dmercatorum/mercatorumA'}
expdirs = cellfun(@(x) fullfile(rootdatadir,x),experiment_names,'UniformOutput',false);

%I don't have these
% classifierparamsfiles = {'../experiments/JAABA_classifier_params1.txt'
%   '../experiments/JAABA_classifier_params2.txt'
%   '../experiments/JAABA_classifier_params_wing.txt'};
classifierparamsfiles={'/groups/branson/home/riveraalbam/Desktop/experiments/classifierslistNov25.txt'};
nintervals = 20;
nframesperseg = 500;
keyword = 'larvaspecies';
timestamp = datestr(now,'yyyymmddTHHMMSS');
figpos = [100,100,1280,720];
forcecompute = false;
%what is going to be written in the movie
% behaviorname_dict = {'pivot_tail','TailPivot'
%   'pivot_head','CenterPivot'
%   'AttemptedCopulation','AttemptedCop'};
behaviorname_dict = {'Crawling','Crawling' 'Head_casting','HeadCasting' 'Turning','Turning' 'Allburrows','Burrowing'};

% behaviororder = {
%   'WingFlick'
%   'Copulation'
%   'AttemptedCop'
%   'WingExtension'
%   'Chase'
%   'Touch'
%   'Righting'
%   'Jump'
%   'Walk'
%   'Crabwalk'
%   'CenterPivot'
%   'TailPivot'
%   'Backup'
%   'Stop'
%   'WingGrooming'
%   };
MakeJAABAResultsMovie(expdirs,classifierparamsfiles,...
   'nframesperseg',nframesperseg,...
   'targets',targets,...
   'framestarts',framestarts,...
   'aviname',aviname,...
    'behaviorname_dict',behaviorname_dict,...
   'fixedaxis','false');
   %'figpos',figpos,...

%% make sure all the scores are computed, up-to-date

% classifierparams = ReadClassifierParamsFile(classifierparamsfiles);
%
% nbehaviors = numel(classifierparams);
% dobreak = false;
% for i = 1:nbehaviors,
%   if dobreak,
%     break;
%   end
%   classifierfilecurr = classifierparams(i).classifierfile;
%   configfilecurr = classifierparams(i).configfile;
%   behavior = classifierparams(i).behaviors.names;
%   if iscell(behavior),
%     behavior = sprintf('%s_',behavior{:});
%     behavior = behavior(1:end-1);
%   end
%   tmp = load(classifierfilecurr);
%   classifier_timestamp = tmp.classifierTS;
%   for j = 1:numel(expdirs),
%     if dobreak,
%       break;
%     end
%     scoresfile_curr = fullfile(expdirs{j},classifierparams(i).file.scorefilename);
%     docompute = false;
%     if ~exist(scoresfile_curr,'file'),
%       fprintf('Scores file %s does not exist.\n',scoresfile_curr);
%       docompute = true;
%     else
%       tmp = load(scoresfile_curr);
%       scores_timestamp = tmp.timestamp;
%       if scores_timestamp ~= classifier_timestamp,
%         keyboard;
%         if forcecompute,
%           docompute = true;
%         else
%           if scores_timestamp > classifier_timestamp,
%             question = sprintf('Scores are newer than classifier for %s, %s. Recompute scores?',behavior,experiment_names{j});
%             def = 'No';
%           else
%             question = sprintf('Scores are older than classifier for %s, %s. Recompute scores?',behavior,experiment_names{j});
%             def = 'Yes';
%           end
%           res = questdlg(question,'Recompute?','Yes','No','Cancel',def);
%           if strcmpi(res,'Yes'),
%             docompute = true;
%           elseif strcmpi(res,'Cancel'),
%             dobreak = true;
%             break;
%           end
%         end
%       end
%     end
%
%     if docompute,
%       fprintf('Computing %s...\n',scoresfile_curr);
%       JAABADetect(expdirs{j},'classifierfiles',{classifierfilecurr},'configfiles',{configfilecurr});
%     end
%
%   end
%
% end

%% choose what to show

% %behaviorsuse = {'Chase','Jump','WingGrooming','Righting','Backup','pivot_tail','pivot_head','WingExtension','WingFlick','AttenptedCopulation'};
% %behaviorsuse = {'WingExtension','WingFlick','AttemptedCopulation','Copulation','pivot_head'};
% behaviorsuse={'Crawling','Head_casting','Turning', 'Allburrows'};
% weightingscheme = 'sumlogsumscore';
%
% [targets,expdirs_chosen,framestarts,nframesperseg] = ...
%   ChooseExampleIntervals(expdirs,classifierparamsfiles,...
%   'nintervals',nintervals,'nframesperseg',nframesperseg,...
%   'weightingscheme',weightingscheme);
%
% % sort based on expdir, then framestart
% maxframestart = max(framestarts);
% [~,~,expi] = unique(expdirs_chosen);
% tmp = expi'*10*maxframestart + framestarts;
% [~,order] = sort(tmp);
% targets = targets(order);
% expdirs_chosen = expdirs_chosen(order);
% framestarts = framestarts(order);
% [unique_expdirs,~,expi] = unique(expdirs_chosen);
%
% for i = 1:nintervals,
%   fprintf('%d: %s, fly %d, frame %d\n',i,expdirs_chosen{i},targets(i),framestarts(i));
% end

% %% make the movies
%
% avibasename = sprintf('JAABAResultsMovie_%s_%s',keyword,timestamp);
%
% for i = 1:nintervals,
%   [~,experiment_name] = fileparts(expdirs_chosen{i});
%   aviname = sprintf('%s_%02d_%s_%02d.avi',avibasename,i,experiment_name,targets(i));
%   MakeJAABAResultsMovie(expdirs_chosen{i},classifierparamsfiles,...
%     'nframesperseg',nframesperseg,...
%     'targets',targets(i),...
%     'framestarts',framestarts(i),...
%     'aviname',aviname,...
%     'figpos',figpos,...
%     'behaviorname_dict',behaviorname_dict);
% end
%
% %% BDP
%
% keyword = 'larvaspecies';
% expi = 2;
% expdir = expdirs{expi};
% avibasename = sprintf('JAABAResultsMovie_%s_%s',keyword,timestamp);
% targets = [9,18,3,17];
% framestarts = [9646,10500,14201,14650];
%
% [~,experiment_name] = fileparts(expdirs_chosen{expi});
% aviname = sprintf('%s_%02d_%s.avi',avibasename,expi,experiment_name);
% MakeJAABAResultsMovie(expdir,classifierparamsfiles,...
%   'nframesperseg',nframesperseg,...
%   'targets',targets,...
%   'framestarts',framestarts,...
%   'aviname',aviname,...
%   'figpos',figpos,...
%   'printnone',false,...
%   'behaviorname_dict',behaviorname_dict,...
%   'behaviororder',behaviororder,...
%   'textprefix','pBDPGAL4U',...
%   'textsuffix','1/3rd speed');
%
% %%
%
% keyword = 'larvaspecies';
% expi = 2;
% expi = 5;
% expdir = expdirs{expi};
% avibasename = sprintf('JAABAResultsMovie_%s_%s',keyword,timestamp);
% targets = [9,18,3,17];
% targets = [10,20];
% framestarts = [9646,10500,14201,14650];
% framestarts = [4396,7345];
%
% [~,experiment_name] = fileparts(expdir);
% aviname = sprintf('%s_%02d_%s.avi',avibasename,expi,experiment_name);
% MakeJAABAResultsMovie(expdir,classifierparamsfiles,...
%   'nframesperseg',nframesperseg,...
%   'targets',targets,...
%   'framestarts',framestarts,...
%   'aviname',aviname,...
%   'figpos',figpos,...
%   'printnone',false,...
%   'behaviorname_dict',behaviorname_dict,...
%   'behaviororder',behaviororder,...
%   'textprefix','R71G01',...
%   'textsuffix','1/6th speed');

author: kristinbranson

repo: JAABA


text:
%% <../index.html MATLAB xUnit Test Framework>: How to Run a Specific Test
% To run all the test cases in just one M-file, ignoring other test
% cases that might be in other files in the same directory, give
% the name of the file (without the ".m" extension) as an argument
% to |runtests|.
%
% For example

cd example_subfunction_tests

runtests testFliplr

%%
% To run a single test case, add the name of the test case using a
% colon (":"), like this:

runtests testFliplr:testFliplrVector

%%
% <../index.html Back to MATLAB xUnit Test Framework>

%%
% Copyright 2008-2010 The MathWorks, Inc.

author: ricedsp


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

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

% Last Modified by GUIDE v2.5 28-Jul-2017 10:52:50

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

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

% --- Executes just before Initialization is made visible.
function Initialization_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to Initialization (see VARARGIN)

% Choose default command line output for Initialization
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes Initialization wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = Initialization_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1

% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on slider movement.
function slider3_Callback(hObject, eventdata, handles)
% hObject    handle to slider3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function slider3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function slider4_Callback(hObject, eventdata, handles)
% hObject    handle to slider4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function slider4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function slider5_Callback(hObject, eventdata, handles)
% hObject    handle to slider5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function slider5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in radiobutton1.
function radiobutton1_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton1

% --- Executes on button press in radiobutton3.
function radiobutton3_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton3

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [ bin, bin_limit ] = r8vec_bin ( n, x, bin_num, bin_min, bin_max )

%*****************************************************************************80
%
%% R8VEC_BIN computes bins based on a given R8VEC.
%
%  Discussion:
%
%    The user specifies minimum and maximum bin values, BIN_MIN and
%    BIN_MAX, and the number of bins, BIN_NUM.  This determines a
%    "bin width":
%
%      H = ( BIN_MAX - BIN_MIN ) / BIN_NUM
%
%    so that bin I will count all entries X(J) such that
%
%      BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I).
%
%    The array X does NOT have to be sorted.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of X.
%
%    Input, real X(N), an (unsorted) array to be binned.
%
%    Input, integer BIN_NUM, the number of bins.  Two extra bins,
%    #1 and #BIN_NUM+2, count extreme values.
%
%    Input, real BIN_MIN, BIN_MAX, define the range and size
%    of the bins.  BIN_MIN and BIN_MAX must be distinct.
%    Normally, BIN_MIN < BIN_MAX, and the documentation will assume
%    this, but proper results will be computed if BIN_MIN > BIN_MAX.
%
%    Output, integer BIN(1:BIN_NUM+2).
%    BIN(1) counts entries of X less than BIN_MIN.
%    BIN(BIN_NUM+2) counts entries greater than or equal to BIN_MAX.
%    For 2 <= I <= BIN_NUM+1, BIN(I) counts the entries X(J) such that
%      BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I).
%    where H is the bin spacing.
%
%    Output, real BIN_LIMIT(1:BIN_NUM+1), the "limits" of the bins.
%    BIN(I) counts the number of entries X(J) such that
%      BIN_LIMIT(I-1) <= X(J) < BIN_LIMIT(I).
%
  if ( bin_max == bin_min )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_BIN - Fatal error!\n' );
    fprintf ( 1, '  BIN_MIN = BIN_MAX.\n' );
    error ( 'R8VEC_BIN - Fatal error!' )
  end

  bin(1:bin_num+2) = 0;

  for i = 1 : n

    t = ( x(i) - bin_min ) / ( bin_max - bin_min );

    if ( t < 0.0 )
      j = 0;
    elseif ( 1.0 <= t )
      j = bin_num + 1;
    else
      j = 1 + floor ( bin_num * t );
    end

    bin(j+1) = bin(j+1) + 1;

  end
%
%  Compute the bin limits.
%
  for i = 0 : bin_num
    bin_limit(i+1) = (   ( bin_num - i ) * bin_min  ...
                       + (           i ) * bin_max ) ...
                       / ( bin_num    );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function newelem=meshreorient(node,elem)
%
% newelem=meshreorient(node,elem)
%
% reorder nodes in a surface or tetrahedral mesh to ensure all
% elements are oriented consistently
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2010/05/05
%
% input:
%    node: list of nodes
%    elem: list of elements (each row are indices of nodes of each element)
%
% output:
%    newelem: the element list with consistent ordering
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

% calculate the canonical volume of the element (can be a 2D or 3D)
vol=elemvolume(node,elem,'signed');

% make sure all elements are positive in volume
idx=find(vol<0);
elem(idx,[end-1,end])=elem(idx,[end,end-1]);
newelem=elem;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function record=PASreadrectxt(filename)
  [fd,syserrmsg]=fopen(filename,'rt');
  if (fd==-1),
    PASmsg=sprintf('Could not open %s for reading',filename);
    PASerrmsg(PASmsg,syserrmsg);
  end;

  matchstrs=initstrings;
  record=PASemptyrecord;
  notEOF=1;
  while (notEOF),
    line=fgetl(fd);
    notEOF=ischar(line);
    if (notEOF),
      matchnum=match(line,matchstrs);
      switch matchnum,
    case 1, [imgname]=strread(line,matchstrs(matchnum).str);
	        record.imgname=char(imgname);
	case 2, [x,y,c]=strread(line,matchstrs(matchnum).str);
	        record.imgsize=[x y c];
	case 3, [database]=strread(line,matchstrs(matchnum).str);
	        record.database=char(database);
	case 4, [obj,lbl,xmin,ymin,xmax,ymax]=strread(line,matchstrs(matchnum).str);
	        record.objects(obj).label=char(lbl);
		record.objects(obj).bbox=[min(xmin,xmax),min(ymin,ymax),max(xmin,xmax),max(ymin,ymax)];
	case 5, tmp=findstr(line,' : ');
	        [obj,lbl]=strread(line(1:tmp),matchstrs(matchnum).str);
	        record.objects(obj).label=char(lbl);
    		record.objects(obj).polygon=sscanf(line(tmp+3:end),'(%d, %d) ')';
	case 6, [obj,lbl,mask]=strread(line,matchstrs(matchnum).str);
	        record.objects(obj).label=char(lbl);
    		record.objects(obj).mask=char(mask);
	case 7, [obj,lbl,orglbl]=strread(line,matchstrs(matchnum).str);
            lbl=char(lbl);
	        record.objects(obj).label=lbl;
    		record.objects(obj).orglabel=char(orglbl);
            if strcmp(lbl(max(end-8,1):end),'Difficult')
                record.objects(obj).difficult=true;
                lbl(end-8:end)=[];
            else
                record.objects(obj).difficult=false;
            end
            if strcmp(lbl(max(end-4,1):end),'Trunc')
                record.objects(obj).truncated=true;
                lbl(end-4:end)=[];
            else
                record.objects(obj).truncated=false;
            end
            t=find(lbl>='A'&lbl<='Z');
            t=t(t>=4);
            if ~isempty(t)
                record.objects(obj).view=lbl(t(1):end);
                lbl(t(1):end)=[];
            else
                record.objects(obj).view='';
            end
            record.objects(obj).class=lbl(4:end);

	otherwise, %fprintf('Skipping: %s\n',line);
      end;
    end;
  end;
  fclose(fd);
return

function matchnum=match(line,matchstrs)
  for i=1:length(matchstrs),
    matched(i)=strncmp(line,matchstrs(i).str,matchstrs(i).matchlen);
  end;
  matchnum=find(matched);
  if isempty(matchnum), matchnum=0; end;
  if (length(matchnum)~=1),
    PASerrmsg('Multiple matches while parsing','');
  end;
return

function s=initstrings
  s(1).matchlen=14;
  s(1).str='Image filename : %q';

  s(2).matchlen=10;
  s(2).str='Image size (X x Y x C) : %d x %d x %d';

  s(3).matchlen=8;
  s(3).str='Database : %q';

  s(4).matchlen=8;
  s(4).str='Bounding box for object %d %q (Xmin, Ymin) - (Xmax, Ymax) : (%d, %d) - (%d, %d)';

  s(5).matchlen=7;
  s(5).str='Polygon for object %d %q (X, Y)';

  s(6).matchlen=5;
  s(6).str='Pixel mask for object %d %q : %q';

  s(7).matchlen=8;
  s(7).str='Original label for object %d %q : %q';

return

author: quantombone

repo: exemplarsvm


text:
function asa076_test02 ( )

%*****************************************************************************80
%
%% TEST02 demonstrates the use of THA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02:\n' );
  fprintf ( 1, '  THA evaluates Owen''s T function.\n' );
  fprintf ( 1, '  Compare to tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          H            A        ' );
  fprintf ( 1, 'T                         T\n' );
  fprintf ( 1, '                                ' );
  fprintf ( 1, '(Tabulated)               (THA)                   DIFF\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, h, a, t1 ] = owen_values ( n_data );

    if ( n_data == 0 )
      break
    end

    t2 = tha ( h, 1.0, a, 1.0 );

    fprintf ( 1, '  %12.8f  %12.8f  %24.16e  %24.16e  %10.4e\n', ...
    h, a, t1, t2, abs ( t1 - t2 ) );

  end

  return
end


text:
%> @brief struct_permute Generate a set of full factorial permutations
%> based on arguments given as fields in a structure.
%> Meant to be used stand-alone or in conjunction with struct_overlay()
%> @note Derived from optPermute() on Mathworks
%> @author Eric Cousineau <eacousineau@gmail.com>
%> @example ???
%>
%> 	opt_values = struct('mass', {{5, 10, 15}}, 'pos', {{1, 2}}, 'extra',
%> 	struct('x', {{2, 3}}));
%> 	opt_set = struct_permute(opt_values);
function [opt_set] = struct_permute(opt_values)
	opt_set = {};

	%Test overlay?

	% Full factorial experiment
	fields = fieldnames(opt_values);
	nfields = length(fields);

	% Start recursion
	opt_blank = struct();
	do_recursion(opt_blank, 1);

	function do_recursion(opt, index)
		if index <= nfields
			field = fields{index};
			value = opt_values.(field);
			if iscell(value)
				% Loop through, concatenate
				for i = 1:length(value)
					opt.(field) = value{i};
					% Recurse through the rest of the options
					do_recursion(opt, index + 1);
				end
			elseif isstruct(value)
				% Recurse with a new set of options for that given option...
				% Get the set of options for that set
				% Incorporate these into permutations
				opt_subset = struct_permute(value);
				% And like normal
				for i = 1:length(opt_subset)
					opt.(field) = opt_subset{i};
					% And so on
					do_recursion(opt, index + 1);
				end
			end
		else
			% End of permute
			% Add final options to the list of permutations
			opt_set{end + 1} = opt;
		end
	end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [ output ] = FW_T( par)

% this function implements Frank-Wolfe-thresholding method
% 0.5*norm(Omega.*(L+S-M), 'fro')^2 + lambda_1 ||L||_* + lambda_2 ||S||_1
%
%
%%  Cun Mu and John Wright, Mar '14

M = par.M; % data matrix
[m,n] = size(M);

epsilon = 10^-3; compare = 0; % by default (if no specification)
if isfield(par, 'epsilon') epsilon = par.epsilon; end
lambda_1 = par.lambda_1;
lambda_2 = par.lambda_2;
iter = par.iter;
display = par.display;
method = par.method; % power or exact
Omega = par.Omega;
rho = sum(sum(Omega))/m/n;
%% FW-T method

% initialization
L = zeros(m,n); S = zeros(m,n);
t_1 = 0; t_2 = 0;

U_1 = 0.5*norm(M,'fro')^2/lambda_1; % initial rough guess
U_2 = 0.5*norm(M,'fro')^2/lambda_2; % initial rough guess

history = 0; % to store the values of each iteration
tic
%% full observation scenario
%if ~isfield(par, 'Omega')
if rho == 1
    fprintf('full observation. \n');
    history(1) = norm(M,'fro')^2/2;
    temp = L + S - M;  % gradient
    count = 0;
    for k = 1: iter

        if k>=2
            if abs(history(k)-history(k-1))/history(k-1)<epsilon
                count = count+1;
            else
                count=0;
            end
            if count==5
                break;
                fprintf('--------------------------------\n')
                fprintf('total # of iter. used is %d. \n', k);
            end
        end

        fprintf('the current function value is %d \n', history(k));

        if mod(k,1) == 0
            fprintf('--------------------------------\n')
            fprintf('this is the %d th iter. \n', k);
        end

        %------------------linearization subproblem-----------------------%

        if  strcmp(method,'power') % approximate
            fprintf('power method is leveraged \n')
            [U ev V] =  power_method(temp, 5);
        else % exact
            [U,ev,V] = lansvd(temp,1,'l'); % top eigen
        end

        D_L = -U*V';

        if lambda_1 >= ev
            V_L = 0; V_t_1 = 0;
        else
            V_L = U_1*D_L; V_t_1 = U_1;
        end

        [mag ind] = max(vec(abs(temp)));
        j = floor((ind-1)/m)+1; i = mod(ind-1,m)+1;
        sign_ = sign(temp(i,j));
        D_S = zeros(m,n);
        D_S(i,j) = -sign_;

        if lambda_2 >= mag
            V_S = 0; V_t_2 = 0;
        else
            V_S = U_2*D_S; V_t_2 = U_2;
        end

        H = zeros(2,2);
        temp_1 = V_L-L; temp_2 = V_S-S; % temp = L + S - M;
        H(1,1) = norm(temp_1, 'fro')^2;
        H(2,2) = norm(temp_2, 'fro')^2;
        H(1,2) = sum(sum(temp_1.*temp_2));
        H(2,1) = H(1,2);
        f = zeros(1,2);
        f(1) = sum(sum(temp_1.*temp));
        f(2) = sum(sum(temp_2.*temp));
        f = f + [lambda_1*(V_t_1-t_1), lambda_2*(V_t_2-t_2)];

        % using QP solvers
        lb = zeros(2,1);
        ub = [1;1];
        options.Display = 'off';
        options.TolFun = 10^-5;
        x = quadprog(H,f,[],[],[],[],lb,ub,[],options);
        x = quadprog(H,f,[],[],[],[],lb,ub,[],options);

        %----------------------- update L and S --------------------------%
        alpha = x(1);
        beta = x(2);
        L = (1-alpha)*L + alpha*V_L; t_1 = t_1 + alpha*(V_t_1-t_1);
        S = (1-beta)*S + beta*V_S; t_2 = t_2 + beta*(V_t_2-t_2);

        %------------------------ thresholding----------------------------%

        temp_3 = M-L;

        S = max(temp_3 - lambda_2, 0);
        S = S + min(temp_3 + lambda_2, 0);

        t_2 = sum(sum(abs(S)));

        %----------- update U_1 and U_2 to a better esitmate -------------%
        temp = L + S - M;
        history(k+1) = norm(temp,'fro')^2/2+lambda_1*t_1+lambda_2*t_2;
        U_1 = min(history(k+1)/lambda_1,U_1);
        U_2 = min(history(k+1)/lambda_2,U_2);

    end

end

%% partial observation scenario
if rho<1
    fprintf('partial observation. \n');
    temp = Omega.*(L + S - M);  % gradient
    history(1) = norm(temp,'fro')^2/2+lambda_1*t_1+lambda_2*t_2;
    count = 0;

    for k = 1: iter

        if k>=2
            if abs(history(k)-history(k-1))*2/history(1)<epsilon
                count = count+1;
            else
                count=0;
            end
            if count==5
                break;
                fprintf('--------------------------------\n')
                fprintf('total # of iter. used is %d. \n', k);
            end
        end

        % current function value
        fprintf('the current function value is %d \n', history(k));
        fprintf('--------------------------------\n')
        fprintf('this is the %d th iter. \n', k);

        %------------------linearization subproblem-----------------------%

        if  strcmp(method,'power') % approximate
            [U ev V] =  power_method(temp, 5);
        else % exact
            [U,ev,V] = lansvd(temp,1,'l'); % top eigen
        end

        D_L = -U*V';

        if lambda_1 >= ev
            V_L = 0; V_t_1 = 0;
        else
            V_L = U_1*D_L; V_t_1 = U_1;
        end

        [mag ind] = max(vec(abs(temp)));
        j = floor((ind-1)/m)+1; i = mod(ind-1,m)+1;
        sign_ = sign(temp(i,j));
        D_S = zeros(m,n);
        D_S(i,j) = -sign_;

        if lambda_2 >= mag
            V_S = 0; V_t_2 = 0;
        else
            V_S = U_2*D_S; V_t_2 = U_2;
        end

        %--------------------- use QP (exact search) ---------------------%
        H = zeros(2,2);
        temp_1 = Omega.*(V_L-L); temp_2 = Omega.*(V_S-S); % temp = L + S - M;
        H(1,1) = norm(temp_1, 'fro')^2;
        H(2,2) = norm(temp_2, 'fro')^2;
        H(1,2) = sum(sum(temp_1.*temp_2));
        H(2,1) = H(1,2);
        f = zeros(1,2);
        f(1) = sum(sum(temp_1.*temp));
        f(2) = sum(sum(temp_2.*temp));
        f = f + [lambda_1*(V_t_1-t_1), lambda_2*(V_t_2-t_2)];

        % using QP solvers
        lb = zeros(2,1);
        ub = [1;1];
        options.Display = 'off';
        options.TolFun = 10^-5;
        x = quadprog(H,f,[],[],[],[],lb,ub,[],options);
        x = quadprog(H,f,[],[],[],[],lb,ub,[],options);

        %----------------------- update L and S --------------------------%
        alpha = x(1);
        beta = x(2);
        L = (1-alpha)*L + alpha*V_L; t_1 = t_1 + alpha*(V_t_1-t_1);
        S = (1-beta)*S + beta*V_S; t_2 = t_2 + beta*(V_t_2-t_2);

        %------------------------ thresholding----------------------------%
        temp_3 = S-Omega.*(L+S-M);;
        S = max(temp_3 - lambda_2, 0);
        S = S + min(temp_3 + lambda_2, 0);
        t_2 = norm(vec(S),1);

        %----------- update U_1 and U_2 to a better esitmate -------------%
        temp = Omega.*(L + S - M);  % gradient
        history(k+1) = norm(temp,'fro')^2/2+lambda_1*t_1+lambda_2*t_2;
        U_1 = min(history(k)/lambda_1,U_1);
        U_2 = min(history(k)/lambda_2,U_2);

    end

end
toc
%% output
output.L = L;
output.S = S;
output.hist = history;
output.iter = k;
end

author: andrewssobral

repo: lrslibrary


text:
% Set up axis and view

% Axis
ax = gca;
axis equal;
ax.NextPlot = 'add';
ax.Clipping = 'off';
ax.Box = 'off';
xlabel('X'); ylabel('Y'); zlabel('Z');
grid off;
if ~showAxis
	axis off;
end

% View
cameratoolbar('Show') % show camera toolbar
cameratoolbar('SetCoordSys', 'y')
camup(ground_plane_normal);
campos(C(:, 68));
if showAxis
    camva(50);
else
	camva(30);
end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
function [lineData, logFooter, linesFooter] = tapas_physio_read_physlogfiles_siemens_raw(...
    fileNameLog)
% Read in raw data/footer lines of logfiles, without data selection/sorting
%
% [lineData, linesFooter] = tapas_physio_read_physlogfiles_siemens_raw(...
%        fileNameLog)
%
% IN
%   fileNameLog     file name of physiological log, e.g. *.ecg
%
% OUT
%   lineData        all recording amplitude samples are saved in first line
%                   of file
%   logFooter       struct() of read-out meta information from log file, i.e.
%                   LogStart/StopTimeSeconds
%                   ScanStart/StopTimeSeconds
%                   => uses MDH time stamp in log file to sync to DICOMs
%   linesFooter     all meta-information (e.g. sampling start/stop) is
%                   saved in remaining lines of log file
%
% EXAMPLE
%   tapas_physio_read_physlogfiles_siemens_raw
%
%   See also

% Author: Lars Kasper
% Created: 2016-02-29
% Copyright (C) 2016 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO 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/>.

fid             = fopen(fileNameLog);

if verLessThan('matlab', '8.5') % use buffer size for speed-up, as long as it exists
    C               = textscan(fid, '%s', 'Delimiter', '\n', 'bufsize', 1e9);
else
    C               = textscan(fid, '%s', 'Delimiter', '\n');
end

fclose(fid);

linesFooter = C{1}(2:end);
lineData = C{1}{1};

%Get time stamps from footer:

logFooter.LogStartTimeSeconds =   str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,...
    'LogStartMDHTime'))),'\D',''))) / 1000;
logFooter.LogStopTimeSeconds =    str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,...
    'LogStopMDHTime'))),'\D',''))) / 1000;

% MPCU  = Computer who controls the scanner => physio logging happens here
% MDH   = Compute who is the host; console => DICOM time stamp here!
%
% - according to Chris Rorden, PART
% (http://www.mccauslandcenter.sc.edu/crnl/tools/part) - MDH is the time we
% should use for phys logging synchronization, since DICOM conversion uses
% this clock

% This is just a different time-scale (of the phys log computer), it does
% definitely NOT match with the Acquisition time in the DICOM-headers
logFooter.ScanStartTimeSeconds = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,...
    'LogStartMPCUTime'))),'\D','')));
logFooter.ScanStopTimeSeconds = str2num(char(regexprep(linesFooter(~cellfun(@isempty,strfind(linesFooter,...
    'LogStopMPCUTime'))),'\D','')));

author: translationalneuromodeling

repo: tapas


text:
function [im1, im2, tu, tv] = read_flow_file(seqName, iSeq)

% read_image_flow reads in two pairs of frames and possibly the ground truth optical flow
% specific to the file location in Deqing Sun's computer
%
% isColor~=0 : readin color image if possible; isColor==0: turn color into
% gray if necessary
%
% called by: read_image_flow_tune_para, compute_flow_test,
%            generate_middelbury, compute_flow_tune_para
%

% Authors: Deqing Sun, Department of Computer Science, Brown University
% Contact: dqsun@cs.brown.edu
% $Date: $
% $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

filePath = ['data' filesep]; % or the folder that you save the data

switch seqName

    case 'middle-eval'

        filePrefex = [filePath 'eval-data/'];

        flowFolder = {'Army',  'Mequon', 'Schefflera', 'Wooden',  'Grove', 'Urban', ...
                      'Yosemite',  'Teddy', 'Basketball',  'Evergreen',  'Backyard',  'Dumptruck'};

        imfilename1 = 'frame10.png';
        imfilename2 = 'frame11.png';

        im1=double(imread([filePrefex flowFolder{iSeq} filesep imfilename1]));
        im2=double(imread([filePrefex flowFolder{iSeq} filesep imfilename2]));

        % GT not available
        tu = nan(size(im1));
        tv = tu;

    case 'middle-other'
        % other sequence with GT provided

        imgFilePath     = [filePath 'other-data/'];
        flowFilePath    = [filePath 'other-gt-flow/'];

        subPath = {'Venus', 'Dimetrodon',   'Hydrangea',    'RubberWhale',...
                    'Grove2', 'Grove3', 'Urban2', 'Urban3', ...
                    'Walking', 'Beanbags',     'DogDance',     'MiniCooper'};

        img1Filename = [imgFilePath subPath{iSeq} filesep 'frame10.png'];
        img2Filename = [imgFilePath subPath{iSeq} filesep 'frame11.png'];

        im1=double(imread(img1Filename));
        im2=double(imread(img2Filename));

        if iSeq <=8

            % Read in ground truth flow fields
            flowFilename = [flowFilePath subPath{iSeq} filesep 'flow10.flo'];
            flow = readFlowFile(flowFilename);
            tu = flow(:,:,1);
            tv = flow(:,:,2);

            % Set unknown values to nan
            UNKNOWN_FLOW_THRESH = 1e9;
            tu (tu>UNKNOWN_FLOW_THRESH) = NaN;
            tv (tv>UNKNOWN_FLOW_THRESH) = NaN;

        else

            % GT not available
            tu = nan(size(im1));
            tv = tu;

        end;

    otherwise
        error('unknown image sequence!');
end;


text:
function svgComposite(s, source1, source2, operator, result, k)
% Adds a feComposite SVG filter
% PRELIMINARY IMPLEMENTATION (Parameters may change)
%
% svgComposite(s, source1, source2, operator, result)
% Parameters:
%   s : Array of plot object handles
%   source1 : Any previous defined filter result string, 'SourceGraphic',
%             or 'SourceAlpha'.
%   source2 : Any previous defined filter result string, 'SourceGraphic',
%             or 'SourceAlpha'.
%   operator : Operator 'over','in','out','atop','xor','arithmetic'
%              -> see SVG 1.1 specification. 'arithmetic' is not yet
%              supported.
%   result : String that identifies the filter result for following filter
%            stages.
for i = 1:length(s)
    userdata = get(s(i),'UserData');
    if isfield(userdata, 'svg') && isfield(userdata.svg, 'Filter')
        next = length(userdata.svg.Filter) + 1;
    else
        next = 1;
    end
    userdata.svg.Filter(next).Subfilter.Type = 'feComposite';
    userdata.svg.Filter(next).Subfilter.Source1 = source1;
    userdata.svg.Filter(next).Subfilter.Source2 = source2;
    userdata.svg.Filter(next).Subfilter.Operator = operator;
    if nargin > 5
        userdata.svg.Filter(next).Subfilter.k = k;
    end
    userdata.svg.Filter(next).Subfilter.Result = result;
    set(s(i),'UserData', userdata);
end


text:
function T = sldlda(X, nums, varargin)
%SLDLDA Performs Direct Linear Discriminant Analysis
%
% $ Syntax $
%   - T = sldlda(X, nums)
%   - T = sldlda(X, nums, ...)
%
% $ Arguments $
%   - X:        the training sample matrix
%   - nums:     the numbers of samples in all classes
%   - T:        the solved transform matrix
%
% $ Description $
%   - T = sldlda(X, nums) performs direct LDA on the samples X using
%     default settings.
%
%   - T = sldlda(X, nums, ...) performs direct LDA on the samples X
%     with the specified properties.
%     \*
%     \t   Table 1.  The properties of Fisher Discriminant Analysis   \\
%     \h     name     &     description                                \\
%           'pdimset' &  The cell containing the arguments for determining
%                        the range space of Sb. They will be input to
%                        slrangespace for dimension determination.
%           'whiten'  &  The cell containing the arguments for computing
%                        the whitening transform in 2nd stage. They will
%                        input to slwhiten_from_cov.
%                        default = {}.       \\
%           'Sb'      &  The pre-computed between-class scattering matrix
%                        or the cell containing the arguments for
%                        computing the scatter matrix in the form
%                        {type, ...}, which is input to slscatter.     \\
%           'Sw'      &  The pre-computed within-class scattering matrix
%                        or the cell containing the arguments for
%                        computing the scatter matrix in the form
%                        {type, ...}, which is input to slscatter.     \\
%         'weights'   &  The sample weights. default = [].             \\
%     \*
%
% $ Remarks $
%   -# The function solves the transform in mainly following stages:
%      First solves the range space of between-class scattering,
%      projecting all samples onto it. Then, solve the whitening transform
%      of the projected within-class scattering.
%
%   -# If Sb or its computing rule is given, the range space is  directly
%      solved from Sb, otherwise the null space is solved from class
%      centers. If both Sb and Sw are given, then the samples are not
%      used in the function. In this cases, you can simply input an empty X.
%
%   -# If both Sb and Sw are given, the pre-pca step will not be conducted.
%      no matter whether prepca is true or false.
%
% $ History $
%   - Created by Dahua Lin on May 1st, 2006
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('slfld', 2);
end

% check size

if ~isempty(X)
    if ndims(X) ~= 2
        error('sltoolbox:invaliddims', ...
            'The sample matrix X should be a 2D matrix');
    end
    [d, n] = size(X);

    k = length(nums);
    if ~isequal(size(nums), [1, k]);
        error('sltoolbox:invaliddims', ...
            'The nums vector should be a row vector');
    end
    if sum(nums) ~= n
        error('sltoolbox:sizmismatch', ...
            'The total number in nums is not consistent with that in X');
    end
end

% check options
opts.pdimset = {};
opts.whiten = {};
opts.Sb = {'Sb'};
opts.Sw = {'Sw'};
opts.weights = [];
opts = slparseprops(opts, varargin{:});

has_Sb = ~isempty(opts.Sb) && isnumeric(opts.Sb);
has_Sw = ~isempty(opts.Sw) && isnumeric(opts.Sw);
if has_Sb && has_Sw
    d = size(opts.Sw, 1);

    if ~isequal(size(opts.Sb), [d, d]) || ~isequal(size(opts.Sw), [d, d])
        error('sltoolbox:sizmismatch', ...
            'Size consistency in Sb and Sw');
    end

else
    if isempty(X)
        error('sltoolbox:invalidargs', ...
            'The samples cannot be empty when Sb or Sw is not pre-computed');
    end
    if (has_Sb && ~isequal(size(opts.Sb), [d, d])) || (has_Sw && ~isequal(size(opts.Sw), [d, d]))
        error('sltoolbox:sizmismatch', ...
            'Size consistency in Sb and Sw');
    end

end
w = opts.weights;

%% Step 1: Compute range space of Sb

if has_Sb
    T1 = slrangespace({'cov', opts.Sb}, opts.pdimset{:});
elseif ~isempty(opts.Sb) && ~isequal(opts.Sb, {'Sb'})
    Sb = slscatter({'cov', X}, opts.Sb{:}, 'sweights', w, 'nums', nums);
    T1 = slrangespace(Sb, opts.pdimset{:});
    clear Sb;
else
    Xc = get_weighted_centers(X, w, nums);
    T1 = slrangespace(Xc, opts.pdimset{:});
    clear Xc wc;
end

%% Step 2: Compute the whiten transform for Sw on range space

if has_Sw
    PSw = T1' * opts.Sw * T1;
else
    X = T1' * X;
    PSw = slscatter(X, opts.Sw{:}, 'sweights', w, 'nums', nums);
end
T2 = slwhiten_from_cov(PSw, opts.whiten{:});
T2 = flipdim(T2, 2);

%% Integrate the transforms

T = T1 * T2;

%% The function for computing weighted centers
function [Xc, wc] = get_weighted_centers(X, w, nums)

Xc = slmeans(X, w, nums);
if isempty(w)
    wc = nums;
else
    k = length(nums);
    [sp, ep] = slnums2bounds(nums);
    wc = zeros(1, k);
    for i = 1 : k
        wc(i) = sum(w(sp(i):ep(i)));
    end
end
wc = sqrt(max(wc, 0));
Xc = slmulvec(Xc, wc, 2);

author: lmthang


text:
function [score, prob, pred] = svmMulticlassTest(model, X)
%function [score, prob, pred] = svmMulticlassTest(model, X)
% Input:
	x = prepareDataSVM(X, model.fNorm, model.iksvmN, model.isSparse, model.mapNanToZero);
	for i = 1:model.numClass,
		[score(i,:), gr] = testSVM(x, model.svmModel(i));
	end
	%Test the Multiclass logistic here
	[prob, pred] = testMCL(sparse(score), model.mlrModel);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function f_x = ParFor26(in1)
%PARFOR26
%    F_X = PARFOR26(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.1.
%    16-Jul-2019 11:58:24

f_x = 0.0;

author: dynamicslab

repo: SINDy-PI


text:
classdef CEC2013_F7 < PROBLEM
% <single> <real> <large>
% 7-nonseparable, 1-separable shifted and rotated Schwefel's function

%------------------------------- Reference --------------------------------
% X. Li, K. Tang, M. N. Omidvar, Z. Yang, and K. Qin, Benchmark functions
% for the CEC'2013 special session and competition on large-scale global
% optimization, RMIT University, Australia, 2013.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        Xopt;	% Optimal decision vector
        R25;    % Rotation matrices
        R50;
        R100;
        p;      % Rank of decision variables
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2013.mat'),'Data');
            obj.Xopt = Data{7}.xopt;
            obj.R25  = Data{7}.R25;
            obj.R50  = Data{7}.R50;
            obj.R100 = Data{7}.R100;
            obj.p    = Data{7}.p;
            obj.M    = 1;
            obj.D    = 1000;
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            S = [50 25 25 100 50 25 25 700];
            W = [6.80e2 9.32e-1 2.12e3 5.06e-1 4.35e2 3.34e4 2.57e0 1];
            PopDec = PopDec - repmat(obj.Xopt,size(PopDec,1),1);
            PopObj = zeros(size(PopDec,1),1);
            for i = 1 : length(S)
                loc = obj.p(sum(S(1:i-1))+1:sum(S(1:i)));
                switch S(i)
                    case 25
                        PopDec(:,loc) = PopDec(:,loc)*obj.R25;
                    case 50
                        PopDec(:,loc) = PopDec(:,loc)*obj.R50;
                    case 100
                        PopDec(:,loc) = PopDec(:,loc)*obj.R100;
                end
                if i < length(S)
                    PopObj = PopObj + W(i)*Schwefel(Tasy(Tosz(PopDec(:,loc)),0.2));
                else
                    PopObj = PopObj + W(i)*Sphere(Tasy(Tosz(PopDec(:,loc)),0.2));
                end
            end
        end
    end
end

function F = Schwefel(X)
    F = sum(cumsum(X,2).^2,2);
end

function F = Sphere(X)
    F = sum(X.^2,2);
end

function Z = Tosz(X)
    X1 = zeros(size(X));
    X1(X~=0) = log(abs(X(X~=0)));
    C1 = zeros(size(X)) + 5.5;
    C1(X>0) = 10;
    C2 = zeros(size(X)) + 3.1;
    C2(X>0) = 7.9;
    Z = sign(X).*exp(X1+0.049*(sin(C1.*X1)+sin(C2.*X1)));
end

function Z = Tasy(X,beta)
    Z = X.^(1+repmat(beta*linspace(0,1,size(X,2)),size(X,1),1).*sqrt(X));
    Z(X<=0) = X(X<=0);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


%% REW_campussquare
% %{
imfolder = 'images\REW_campussquare';
im_n = 2;
imfile = cell(im_n,1);
imfile{1} = [imfolder '\' 'campussquare_01.jpg'];
imfile{2} = [imfolder '\' 'campussquare_02.jpg'];

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = [1,2];

% imsize = zeros(im_n,3);
%
% for ii = 1:im_n
%     imsize(ii,:) = size(im{ii});
%     if imsize(ii,1) > 720
%         scale = 720/size(im{ii}, 1);
%         im{ii} = imresize(im{ii}, scale);
%
%         imsize(ii,:) = size(im{ii});
%     end
% end

mosaic = REW_mosaic( im, edge_list, 1, 'persp', 0, imfolder );
%}


text:
function [outputCONTOUR,varargout] = TRACE_MooreNeighbourhood(data2D,varargin)
% TRACE_MooreNeighbourhood  2D boundary tracing using the Moore neighbourhood
%==========================================================================
% FILENAME:     TRACE_MooreNeighbourhood.m
% AUTHOR:       Adam H. Aitkenhead
% DATE:         12th April 2010
% PURPOSE:      2D boundary tracing using the Moore neighbourhood
%
% USAGE:        [listCONTOUR,listNORMALS] = TRACE_MooreNeighbourhood(data2D,pixelFIRST)
%%
% INPUT PARAMETERS:  data2D      - Mandatory.  An LxM array containing the
%                                  2D data.  All pixels within the region
%                                  to be traced must equal 1, and all
%                                  pixels outside the region must equal 0.
%                    pixelFIRST  - Optional.  A (2x1) array containing the
%                                  [x;y] coordinates of any one pixel
%                                  located on the edge of the region to be
%                                  traced.
%
% OUTPUT PARAMETERS: listCONTOUR - Mandatory.  An Nx2 array containing the
%                                  coordinates of the pixels located around
%                                  the edge of the region to be traced.
%                                  The pixels are listed in a clockwise
%                                  direction around the contour.
%                    listNORMALS - Optional.  An Nx2 array containing the
%                                  (approximate) normal vectors pointing
%                                  out of the region of interest.  The
%                                  vectors are listed in a clockwise
%                                  direction around the contour, and
%                                  correspond to each pixel listed in
%                                  listCONTOUR.
%
% EXAMPLE:      data2D = [0,0,0,0,0,0; 0,1,0,1,0,0; 0,1,1,1,1,0; 0,1,1,1,0,0; 0,1,1,1,1,0; 0,0,0,0,0,0];
%               [listCONTOUR,listNORMALS] = TRACE_MooreNeighbourhood(data2D);
%               imagesc(data2D)
%               colormap('gray')
%               hold on
%                 plot(listCONTOUR(:,2),listCONTOUR(:,1),'ro-','LineWidth',2)
%                 plot([listCONTOUR(:,2),listCONTOUR(:,2)+listNORMALS(:,2)]',[listCONTOUR(:,1),listCONTOUR(:,1)+listNORMALS(:,1)]','b','LineWidth',2)
%               hold off
%
% NOTES:      - The algorithm ends when the *second* pixel in the loop is
%               revisited, entered from the same direction as it was
%               entered on its first visit.  The second pixel is used
%               rather than the first pixel for two reasons:  1. The
%               first pixel may be encountered more than once while going
%               around the contour;  2. It is possible that the first pixel
%               will never be entered from the same direction as it was on
%               its first visit.
%             - The code does not identify any internal holes which exist
%               in the region to be traced.
%             - If more than one separate region exists in the image, the
%               code only traces one of these regions.  The user can define
%               which region by choosing an appropriate starting pixel
%               <pixelFIRST>.  Alternatively, the first region encountered
%               by the algorithm will be traced.
%
% REFERENCES: - For a guide to Contour Tracing, including a description of
%               Moore-Neighbour Tracing, see the tutorial by Abeer George
%               Ghuneim at:
%                http://www.imageprocessingplace.com/downloads_V3/
%                  root_downloads/tutorials/
%                  contour_tracing_Abeer_George_Ghuneim/index.html
%==========================================================================

%==========================================================================
% VERSION:  USER:  CHANGES:
% -------   -----  --------
% 090512    AHA    Original version
% 091021    AHA    Improved the calculation of the normals
% 100412    AHA    General housekeeping
%==========================================================================

%================================================
%  READ INPUT PARAMETERS
%================================================

% Configure the array which contains the region to be contour-traced:
if nargin<1
  error(' ERROR:  A 2D array is required as input.');
end

% Define pixelFIRST, the start coordinate (x,y) for the Moore-neighbour tracing algorithm.  This must be an edge pixel on the region to be contour-traced.
if nargin==2
  pixelFIRST = varargin{1};
else %If pixelFIRST hasn't been defined by the user, define it now:
  [tempR,tempC] = find(data2D==1);
  pixelFIRST = [tempR(1);tempC(1)];
end

%Make sure pixelFIRST is a column vector, rather than a row vector
if size(pixelFIRST,2)==2
  pixelFIRST = pixelFIRST';
end

%Makes sure data2D is a binary image
if numel(unique(data2D))~=2  ||  min(data2D(:))~=0  ||  max(data2D(:))~=1
  error(' Input image must be an array composed only of 0s and 1s.')
end

% Add a one-pixel border around the array before doing the analysis.
% This prevents the code from failing if the region to be traced extends to the edge of the array.
data2Doriginal = data2D;
data2D         = [zeros(1,size(data2Doriginal,2)+2);zeros(size(data2Doriginal,1),1),data2Doriginal,zeros(size(data2Doriginal,1),1);zeros(1,size(data2Doriginal,2)+2)];
pixelFIRST     = pixelFIRST + 1;

%================================================
%  INITIAL DEFINITIONS
%================================================

%Define the initial direction of entry (x,y) into the first pixel.  This
%direction must point from the outside to the inside.
tempI = find(data2D(pixelFIRST(1)-1:pixelFIRST(1)+1,pixelFIRST(2))==0);
if isempty(tempI)==0
  directionFIRST = [2-tempI(1);0];
else
  tempI = find(data2D(pixelFIRST(1),pixelFIRST(2)-1:pixelFIRST(2)+1)==0);
  if isempty(tempI)==0
    directionFIRST = [0;2-tempI(1)];
  else
    imagesc(data2D)
    save('temp.txt','data2D','-ascii','-tabs')
    error(' ERROR:  Unable to identify a vector pointing into the first pixel on the contour.');
  end
end

%Create an array for the output.  It will hold the (x,y) coordinates of
%every pixel around the border of the region to be contour-traced.  The
%directions are also recorded, and will be used as part of the test
%criterion to end the algorithm:
listCONTOUR   = pixelFIRST;
listDIRECTION = directionFIRST;

%Create a variable which will be used to tell the algorithm when to stop:
checkLOOPEND = 0;

%================================================
%  PERFORM THE MOORE NEIGHBOURHOOD TRACING
%================================================

ROTangle     = -pi/2;   %radians
ROTmatrix90  = [ cos(ROTangle) , -sin(ROTangle) ; sin(ROTangle) , cos(ROTangle) ];
ROTangle     = -pi;   %radians
ROTmatrix180 = [ cos(ROTangle) , -sin(ROTangle) ; sin(ROTangle) , cos(ROTangle) ];

tempPIX = pixelFIRST;
tempDIR = directionFIRST;

while checkLOOPEND==0

  tempDIR = round(ROTmatrix180*tempDIR);
  tempPIX = tempPIX + tempDIR;
  if data2D(tempPIX(1),tempPIX(2)) == 1
    listCONTOUR   = [listCONTOUR,tempPIX];
    listDIRECTION = [listDIRECTION,tempDIR];
  else
    tempDIR = round(ROTmatrix90*tempDIR);
    tempPIX = tempPIX + tempDIR;
    if data2D(tempPIX(1),tempPIX(2)) == 1
      listCONTOUR   = [listCONTOUR,tempPIX];
      listDIRECTION = [listDIRECTION,tempDIR];
    else
      tempDIR = round(ROTmatrix90*tempDIR);
      tempPIX = tempPIX + tempDIR;
      if data2D(tempPIX(1),tempPIX(2)) == 1
        listCONTOUR   = [listCONTOUR,tempPIX];
        listDIRECTION = [listDIRECTION,tempDIR];
      else
        tempPIX = tempPIX + tempDIR;
        if data2D(tempPIX(1),tempPIX(2)) == 1
          listCONTOUR   = [listCONTOUR,tempPIX];
          listDIRECTION = [listDIRECTION,tempDIR];
        else
          tempDIR = round(ROTmatrix90*tempDIR);
          tempPIX = tempPIX + tempDIR;
          if data2D(tempPIX(1),tempPIX(2)) == 1
            listCONTOUR   = [listCONTOUR,tempPIX];
            listDIRECTION = [listDIRECTION,tempDIR];
          else
            tempPIX = tempPIX + tempDIR;
            if data2D(tempPIX(1),tempPIX(2)) == 1
              listCONTOUR   = [listCONTOUR,tempPIX];
              listDIRECTION = [listDIRECTION,tempDIR];
            else
              tempDIR = round(ROTmatrix90*tempDIR);
              tempPIX = tempPIX + tempDIR;
              if data2D(tempPIX(1),tempPIX(2)) == 1
                listCONTOUR   = [listCONTOUR,tempPIX];
                listDIRECTION = [listDIRECTION,tempDIR];
              else
                tempPIX = tempPIX + tempDIR;
                if data2D(tempPIX(1),tempPIX(2)) == 1
                  listCONTOUR   = [listCONTOUR,tempPIX];
                  listDIRECTION = [listDIRECTION,tempDIR];
                else
                  tempDIR = round(ROTmatrix90*tempDIR);
                  tempPIX = tempPIX + tempDIR;
                  if data2D(tempPIX(1),tempPIX(2)) == 1
                    listCONTOUR   = [listCONTOUR,tempPIX];
                    listDIRECTION = [listDIRECTION,tempDIR];
                  else
                    tempDIR = round(ROTmatrix90*tempDIR);
                    tempPIX = tempPIX + tempDIR;
                    if data2D(tempPIX(1),tempPIX(2)) == 1
                      listCONTOUR   = [listCONTOUR,tempPIX];
                      listDIRECTION = [listDIRECTION,tempDIR];
                    end %if

                  end %if
                end %if
              end %if
            end %if
          end %if
        end %if
      end %if
    end %if
  end %if

  if size(listCONTOUR,2)>2  &&  min(isequal(tempPIX,listCONTOUR(:,2)))  &&  min(isequal(tempDIR,listDIRECTION(:,2)))
    %The stop criterion has been met, so tell the while-loop to exit:
    checkLOOPEND=1;
    %Remove the extra col from the results corresponding to the second
    %visit to the second pixel.
    listCONTOUR         = listCONTOUR(:,1:end-1);
    listDIRECTION       = listDIRECTION(:,1:end-1);
    listENTRYDIRECTIONS = -listDIRECTION;
    %Make sure the normal vector for the first pixel is the same on its
    %first visit as on its last vist:
    listENTRYDIRECTIONS(:,1) = listENTRYDIRECTIONS(:,end);
  end %if

end %while

% Adjust the results to remove the effect of the one-pixel border that was
% added around the array before the contour was traced:
listCONTOUR = listCONTOUR - 1;

%================================================
%  COMPUTE THE NORMALS
%  A combination of two methods is used to compute the normals:
%  -1-  The gradient of the two surrounding contour points is found, and
%       then rotated 90 degrees to get the gradient at the point of
%       interest.
%  -2-  The midpoint of the two surrounding points is found, and the vector
%       from this point to the point of interest is found.
%  Method 1 handles concave edges and diagonal edges better, but gets
%  confused in places where the region of interest ends in a one-pixel-wide
%  line, since at that point the two surrounding points are coincident.
%================================================

% Remove last point, which is a duplicate of the first point.
listCONTOURtemp = listCONTOUR(:,1:end-1);

% Method 1:  Gradient:
listGRADIENTS = (circshift(listCONTOURtemp,[0,-1])-circshift(listCONTOURtemp,[0,1]))/2;  % Calculate the gradient of the two surrounding contour-points.
ROTangle      = pi/2;   %radians                                                         % Rotate the gradient by +90 degree to produce the normals.
ROTmatrix     = [ cos(ROTangle) , -sin(ROTangle) ; sin(ROTangle) , cos(ROTangle) ];
listNORMALS_G = ROTmatrix * listGRADIENTS;

% Method 2:  Midpoint:
listMIDPOINT  = (circshift(listCONTOURtemp,[0,-1])+circshift(listCONTOURtemp,[0,1]))/2;  %Find the midpoint of the two surrounding contour-points.
listNORMALS_M = (listCONTOURtemp - listMIDPOINT)./2;                                     % Calculate the vector to the point of interest

% Combine results from the two methods:
method1sign = sign(round(10*listNORMALS_G));
method2sign = sign(round(10*listNORMALS_M));
method1sign(method1sign==0) = method2sign(method1sign==0);
listNORMALS_M = method1sign.*abs(listNORMALS_M);
listNORMALS        = listNORMALS_M;
listNORMALS(:,:,2) = listNORMALS_G;
listNORMALS = mean(listNORMALS,3);

% For cases where the ROI is larger than a single pixel, ensure the normals are units vectors
if size(listNORMALS,2)>1
  normalLENGTHS = sqrt(sum(listNORMALS.^2,1));
  listNORMALS   = listNORMALS./([1;1]*normalLENGTHS);
else
  listNORMALS = [-1;0];
end

% Replace last point, which is a duplicate of the first point.
listNORMALS   = [listNORMALS,listNORMALS(:,1)];

%================================================
%  DEFINE THE OUTPUT ARGUMENTS
%================================================

outputCONTOUR = listCONTOUR';
listNORMALS   = listNORMALS';

if nargout == 2
  varargout(1) = {listNORMALS};
end


text:
function c = getColor(cidxs)

% http://colorbrewer2.org/

color = cell(8,4);

% qualitative
color{1}{1} = [55,126,184];
color{2}{1} = [255,127,0];
color{3}{1} = [255,255,51];%[178,223,138];
color{4}{1} = [0,0,0];
color{5}{1} = [77,175,74];
color{6}{1} = [228,26,28];
% color{6}{1} = [255,255,51];
color{7}{1} = [152,78,163];
color{8}{1} = [247,129,191];

% sequential
color{1}{2} = [253,141,60];
color{1}{3} = [254,204,92];
color{1}{4} = [255,255,178];
color{1}{5} = [254,240,217];

color{2}{2} = [107,174,214];
color{2}{3} = [189,215,231];
color{2}{4} = [117,107,177];

color{3}{2} = [194,230,153];

color{4}{2} = [140,150,198];
color{4}{3} = [179,205,227];
color{4}{4} = [237,248,251];

c = color{cidxs(1)}{cidxs(2)};

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function fields = utl_timeseries_fields(signal)
% Get the time-series fields of te given signal.
% function Fields = utl_timeseries_fields(Signal)
%
% This function returns the field names in the given signal that are
% carrying time-series information (and which therefore should be filtered,
% buffered, etc.). The result is a combination of a set of hard-coded fields
% and the names listed in the field .tracking.timeseries_fields. Only fields
% that are present in the signal are returned.
%
% In:
%   Signal : a signal for which time-series field names shall be looked up
%
% Out:
%   Fields : Cell array of field names (row vector). It is assumed that the second dimension of
%            these fields is the time axis along which buffering and/or filtering should happen; any
%            number of other dimensions may be present for any field returned by this function.
%
% See also:
%   utl_register_field
%
%                                       Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                       2013-08-16

% define a database of cached queries and associated results
% (this approach is more than 10x as fast as using unique/intersect every time)
persistent keys;
persistent values;

try
    % generate the query (assuming that the signal usually has the right fields)
    try
        query = [fieldnames(signal)' {'|'} signal.tracking.timeseries_fields];
    catch
        query = fieldnames(signal)';
    end

    % turn into a string
    query = [query{:}];

    try
        % look up from the database (assuming that it is already contained)
        fields = values{strcmp(keys,query)};
    catch
        % not yet in DB: actually determine the timeseries fields
        if isfield(signal,'tracking') && isfield(signal.tracking,'timeseries_fields')
            fields = get_timeseries_fields(fieldnames(signal),signal.tracking.timeseries_fields);
        else
            fields = get_timeseries_fields(fieldnames(signal),[]);
        end
        % initialize the DB if necessary, otherwise append
        if ~iscell(keys)
            keys = {query};
            values = {fields};
        else
            keys{end+1} = query;
            values{end+1} = fields;
        end
    end
catch e
    % diagnose the error
    if ~isstruct(signal)
        error('The given argument must be a signal struct.'); end
    if ~isscalar(signal)
        error('The given argument must be a 1x1 struct.'); end
    if ~isfield(signal,'tracking')
        if ~any(isfield(signal,{'data','srate'}))
            if all(isfield(signal,{'head','parts'}))
                error('The given input is an unevaluated expression but this function expects a signal data structure. You may have to first evaluate the data using exp_eval.');
            else
                error('The given data structure does not appear to be a valid signal struct.');
            end
        else
            error('The given signal is lacking the required .tracking field.');
        end
    else
        rethrow(e);
    end
end

function fields = get_timeseries_fields(field_names,ts_fields)
% This function performs the actual computation

% generate initial list of candidates
candidates = {'data','icaact','srcpot','stamps'};

% append whatever is registered in the signal's .tracking.timeseries_fields
candidates = unique([candidates(:); ts_fields(:)]);

% restrict to those that are actually present in the signal
fields = intersect(candidates(:),field_names);
fields = fields(:)';

repo: REST