text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Gaussian Process Regression 2D Example  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%             1) Load 2D Regression Datasets                 %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Example of Gaussian Process Regression (2D)
%% Generate Random Gaussian distribution
clear all; close all; clc;
K = 90;
a =-50;
b = 50;
x = a + (b-a).*rand(K,1);
y = a + (b-a).*rand(K,1);

a    = 40;
b    = 100;
vars = a + (b-a).*rand(K,1);

Mu    = [x,y]';
Sigma = zeros(2,2,K);

for k=1:K
    Sigma(:,:,k) = eye(2,2) .* vars(k);
end

gmm_x.Priors    = ones(1,K)./K;
gmm_x.Mu        = Mu;
gmm_x.Sigma     = Sigma;

% Generate some training Data
X = ml_gmm_sample(500,gmm_x.Priors,gmm_x.Mu,gmm_x.Sigma )';
f = @(X)ml_gmm_pdf(X',gmm_x.Priors,gmm_x.Mu,gmm_x.Sigma );
y = f(X);
y = y(:);

% Plot Training Data

% Plot Real Function
options           = [];
options.title     = 'Training data';
options.surf_type = 'surf';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_value_func(X,f,[1 2],options);hold on

% Plot Training Data
options.bFigure     = false;
options.surf_type   = 'scatter';
options.points_size = 20;
ml_plot_value_func(X,f,[1 2],options);hold on

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 2) "Train" GPR Model                       %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Train GP (no real training, just give parameters)
epsilon         = 0.1;
dims            = 2;
model.X_train   = X;
model.y_train   = y;
rbf_var         = 25;
gp_f            = @(X)ml_gpr(X,[],model,epsilon,rbf_var);

% Plot Estimated Function
options           = [];
options.title     = 'Estimated y=f(x) from Gaussian Process Regression';
options.surf_type = 'surf';
if exist('h2','var') && isvalid(h2), delete(h2);end
h2 = ml_plot_value_func(X,gp_f,[1 2],options);hold on

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%           3) Test GPR Model on Train points                %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Test GPR
options           = [];
options.bFigure   = true;
options.title     = 'Test data';
options.surf_type = 'pcolor';
dims              = [1,2];

if exist('h2','var') && isvalid(h2), delete(h2);end
h2 = gp_plot(model.X_train,gp_f,dims,options);
colorbar

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%           4) Grid Search for GPR with RBF Kernel           %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% K-fold cross validation
Kfold = 20;

disp('Parameter grid search GP');

rbf_vars = [5,10,20,40,50,100,10000];

test  = cell(length(rbf_var),1);
train = cell(length(rbf_var),1);

for i=1:length(rbf_vars)
    disp(['[' num2str(i) '/' num2str(length(rbf_vars)) ']']);

    f                       = @(X,y,model)ml_gpr(X,y,model,epsilon,rbf_vars(i));
    [test_eval,train_eval]  = ml_kcv(X,y,Kfold,f,'regression');

    test{i}                 = test_eval;
    train{i}                = train_eval;
    disp(' ');
end

%% Get Statistics

[ stats ] = ml_get_cv_grid_states_regression(test,train);

%% Plot Statistics

options             = [];
options.title       = 'GPR k-CV';
options.metrics     = {'nrmse','r'};     % <- you can add many other metrics, see list in next cell box
options.para_name   = 'variance rbf';

if exist('handle','var'), delete(handle); end
[handle,handle_test,handle_train] = ml_plot_cv_grid_states_regression(stats,rbf_vars,options);

%% Full list of evaluation metrics for regression methods

options.metric = {'mse','nmse','rmse','nrmse','mae','mare','r','d','e','me','mre'};

%   '1'  - mean squared error                               (mse)
%   '2'  - normalised mean squared error                    (nmse)
%   '3'  - root mean squared error                          (rmse)
%   '4'  - normalised root mean squared error               (nrmse)
%   '5'  - mean absolute error                              (mae)
%   '6'  - mean  absolute relative error                    (mare)
%   '7'  - coefficient of correlation                       (r)
%   '8'  - coefficient of determination                     (d)
%   '9'  - coefficient of efficiency                        (e)
%  '10'  - maximum absolute error                           (me)
%  '11'  - maximum absolute relative error                  (mre)

author: epfl-lasa


text:
%%  Publish vistasoft functional data on remote data server
%  (Archiva)
%
% These functional data are used as part of the mrvTest protocol. They
% accompany the validation data that are also stored on the remote site.
%
% BW did an svn update on the vistadata directory.  Then ran through this
% script.
%
% BW/ Copyright Vistasoft Team 2016

%% Open the remote data object

% Create remote data toolbox object
rd = RdtClient('vistasoft');

% To write to the archive, you must have a password
rd.credentialsDialog;

rd.crp('/vistadata/functional');

folder = fullfile(vistaRootPath,'local','functional');
a = rd.publishArtifacts(folder,'type','zip','verbose',true);
rdtPrintArtifactTable(a);

%% OLD

%% Open the remote data object

% Create remote data toolbox object
rd = RdtClient('vistasoft');

% To write to the archive, you must have a password
rd.credentialsDialog;

% Act with lots of printouts
rd.configuration.verbosity=0;

% Theses are the paths on the remote data server
p = rd.listRemotePaths;

%% First, copy the elements in the mrBOLD_01 directory
baseDir = fullfile(vistaRootPath,'..','vistadata','functional');
cd(baseDir);

% Note the the remote path and the local path match.
rd.crp('/vistadata/functional/mrBOLD_01');
fullDirectory = fullfile(baseDir,'mrBOLD_01');

% Publish the files in the mrBOLD_01 directory.
rd.publishArtifacts(fullDirectory,'print',true);

% You can check for the files in that remote path
a = rd.listArtifacts('print',true);

% Read one of the ones in the remote path
d = rd.readArtifact(a(12));

%% Now make the base mrBOLD_01 directory and publish the various subdirs

% Sometimes we have to go two deep.
baseDir = fullfile(vistaRootPath,'..','vistadata','functional','mrBOLD_01');
cd(baseDir)

%% Gray
rd.crp('/vistadata/functional/mrBOLD_01/Gray');
fullDirectory = fullfile(baseDir,'Gray');
dir(fullDirectory)

rd.publishArtifacts(fullDirectory,'verbose',true);
a = rd.listArtifacts;
a(:).artifactId

%% Inplane
rd.crp('/functional/mrBOLD_01/Inplane');
fullDirectory = fullfile(baseDir,'Inplane');
rd.publishArtifacts(fullDirectory);
a = rd.listArtifacts;
a(:).artifactId

%% Inplane/Original
rd.crp('/functional/mrBOLD_01/Inplane/Original');
fullDirectory = fullfile(baseDir,'Inplane','Original');
dir(fullDirectory)

rd.publishArtifacts(fullDirectory);
a = rd.listArtifacts;
a(:).artifactId

%% Inplane/ROIs
rd.crp('/functional/mrBOLD_01/Inplane/ROIs');
fullDirectory = fullfile(baseDir,'Inplane','ROIs');
dir(fullDirectory)

rd.publishArtifacts(fullDirectory);
a = rd.listArtifacts;
a(:).artifactId

%% Raw
rd.crp('/functional/mrBOLD_01/Raw');
fullDirectory = fullfile(baseDir,'Raw');
dir(fullDirectory)

rd.publishArtifacts(fullDirectory);
a = rd.listArtifacts;
a(:).artifactId

%% 3DAnatomy
rd.crp('/functional/mrBOLD_01/3DAnatomy');
fullDirectory = fullfile(baseDir,'3DAnatomy');
dir(fullDirectory)

rd.publishArtifacts(fullDirectory);
a = rd.listArtifacts;
a(:).artifactId

%% Volume
rd.crp('/functional/mrBOLD_01/Volume');
fullDirectory = fullfile(baseDir,'Volume');
dir(fullDirectory)

rd.publishArtifacts(fullDirectory);
a = rd.listArtifacts;
a(:).artifactId

%%

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [node,face,elem]=meshunitsphere(tsize,maxvol)
%
% [node,face,elem]=meshunitsphere(tsize,maxvol)
%
% create the surface and/or volumetric mesh of a unit sphere
% centered at [0 0 0] and radius 1
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   tsize: maximum size of the surface triangles (from 0 to 1)
%   maxvol: maximum volume of the tetrahedron; if one wants to return
%           elem without specifying maxvol, maxvol=tsize^3
%
% output:
%   node: node coordinates, 3 columns for x, y and z respectively
%   face: integer array with dimensions of NB x 3, each row represents
%         a surface mesh face element
%   elem: integer array with dimensions of NE x 4, each row represents
%         a tetrahedron. If ignored, this function only produces the surface
%
% example:
%   [node,face]=meshunitsphere(0.05);
%   [node,face,elem]=meshunitsphere(0.05,0.01);
%   plotmesh(node,elem,'x>0'); axis equal;
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

dim=60;
esize=tsize*dim;
thresh=dim/2-1;
[xi,yi,zi]=meshgrid(0:0.5:dim,0:0.5:dim,0:0.5:dim);
dist=thresh-sqrt((xi-30).^2+(yi-30).^2+(zi-30).^2);
dist(dist<0)=0;
clear xi yi zi;

% extract a level-set at v=thresh, being a sphere with R=thresh
% the maximum element size of the surface triangles is tsize*dim

[node,face]=vol2restrictedtri(dist,1,[dim dim dim],dim*dim*dim,30,esize,esize,40000);
node=(node-0.5)*0.5;

[node,face]=removeisolatednode(node,face);

node=(node-30)/28;
r0=sqrt(sum((node.*node)'));
node=node.*repmat(1./r0(:),1,3);

if(nargout==3)
   if(nargin==1)
       maxvol=tsize*tsize*tsize;
   end
   [node,elem,face]=surf2mesh(node,face,[-1 -1 -1]*1.1,[1 1 1]*1.1,1,maxvol,[],[]);
   elem=elem(:,1:4);
end

face=face(:,1:3);

author: fangq


text:
% METAPLOTTOPO - plot concatenated multichannel data epochs in a topographic or
%               rectangular array. Uses a channel location file with the same
%               format as TOPOPLOT, or else plots data on a rectangular grid.
%
% Usage:
%    >> axes = metaplottopo(data, 'key1', 'val1', 'key2', 'val2')
%
% Inputs:
%   data       = data consisting of consecutive epochs of (chans,frames)
%                or (chans,frames,n)
%
% Optional inputs:
%  'chanlocs'  = [struct] channel structure or file plot ERPs at channel
%                locations. See help READLOCS for data channel format.
%  'geom'      = [rows cols] plot ERP in grid (overwrite previous option).
%                Grid size for rectangular matrix. Example: [6 4].
%  'title'     = [string] general plot title {def|'' -> none}
%  'chans'     = vector of channel numbers to plot {def|0 -> all}
%  'axsize'    = [x y] axis size {default [.07 .07]}
%  'plotfunc'  = [string] plot function name. If none is entered, axes
%                are created and returned.
%  'plotargs'  = [cell] plotting function arguments.
%  'datapos'   = [integer] position of data array in the function call.
%                Default is 1.
%  'squeeze'   = ['on'|'off'] squeeze data array after selecting channel.
%                Default is 'off'.
%
% Output:
%  Axes        = [real] array of axes handles of the same length as the
%                number of plotted channels.
%  Channames   = [cell] cell array of channel name for each plot
%
% Author: Arnaud Delorme, Scott Makeig, CERCO, CNRS, 2007-
%
% See also: PLOTTOPO

% Copyright (C) 2007, Arnaud Delorme, CERCO, arno@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [Axes, outchannames ]= metaplottopo(data, varargin);

%
%%%%%%%%%%%%%%%%%%%%% Graphics Settings - can be customized %%%%%%%%%%%%%%%%%%
%
LINEWIDTH     = 0.7;     % data line widths (can be non-integer)
FONTSIZE      = 10;      % font size to use for labels
CHANFONTSIZE  = 7;       % font size to use for channel names
TICKFONTSIZE  = 8;       % font size to use for axis labels
TITLEFONTSIZE = 12;      % font size to use for the plot title
PLOT_WIDTH    = 0.95;     % 0.75, width and height of plot array on figure
PLOT_HEIGHT   = 0.88;    % 0.88
gcapos = get(gca,'Position'); axis off;
PLOT_WIDTH    = gcapos(3)*PLOT_WIDTH; % width and height of gca plot array on gca
PLOT_HEIGHT   = gcapos(4)*PLOT_HEIGHT;
MAXCHANS      = 256;     % can be increased
curfig = gcf;            % learn the current graphic figure number
%
%%%%%%%%%%%%%%%%%%%% Default settings - use commandline to override %%%%%%%%%%%
%
DEFAULT_AXWIDTH  = 0.05; %
DEFAULT_AXHEIGHT = 0.08; %
DEFAULT_SIGN = -1;                        % Default - plot positive-up
ISRECT = 0;                               % default
ISSPEC = 0;                               % Default - not spectral data

outchannames = {};
if nargin < 1
    help metaplottopo
    return
end

if iscell(data), nchans = size(data{1},1);
else             nchans = size(data,1);
end
g = finputcheck(varargin, { 'chanlocs'  ''    []          '';
    'chans'     'integer'               [1 size(data,1)]  [1:nchans];
    'geom'      'integer'               [1 Inf]     [];
    'title'     'string'                []          '';
    'plotfunc'  'string'                []          '';
    'plotargs'  'cell'                  []          {};
    'datapos'   'integer'               []          1;
    'calbar'    'real'                  []          [];
    'axcopycom' 'string'                []          '';
    'squeeze'   'string'                {'on' 'off'} 'off';
    'axcopycom' 'string'                {}          'off';
    'axsize'    'float'                 [0 1]       [nan nan]}, 'metaplottopo' );
if ischar(g), error(g); end
if length(g.chans) == 1 && g.chans(1) ~= 0, error('can not plot a single ERP'); end

[chans,framestotal]=size(data);           % data size

%
%%%%%%%%%%%%%%% Substitute defaults for missing parameters %%%%%%%%%%%%%%%%
%

if (isempty(g.chanlocs) || ~isfield(g.chanlocs, 'theta')) && isempty(g.geom)
    n = ceil(sqrt(length(g.chans)));
    g.geom = [n n];
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Test parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
axwidth  = g.axsize(1);
axheight = g.axsize(2);
if ~isempty(g.geom)
    if isnan(axheight) % if not specified
        axheight = gcapos(4)/(g.geom(1)+1);
        axwidth  = gcapos(3)/(g.geom(2)+1);
    end
else
    axheight = DEFAULT_AXHEIGHT*(gcapos(4)*1.25);
    axwidth =  DEFAULT_AXWIDTH*(gcapos(3)*1.3);
end
fprintf('Plotting data using axis size [%g,%g]\n',axwidth,axheight);

%
%%%%%%%%%%%%% Extend the size of the plotting area in the window %%%%%%%%%%%%
%
curfig = gcf;
h=figure(curfig);
set(h,'PaperUnits','normalized'); % use percentages to avoid US/A4 difference
set(h,'PaperPosition',[0.0235308 0.0272775 0.894169 0.909249]); % equivalent
orient portrait
axis('normal');

axcolor= get(0,'DefaultAxesXcolor'); % find what the default x-axis color is
vertcolor = 'k';
horicolor = vertcolor;

%
%%%%%%%%%%%%%%%%%%%%%%%%% Plot and label specified channels %%%%%%%%%%%%%%%%%%
%

if ~isempty(data)
    if iscell(data)
        data{1} = data{1}(g.chans,:,:,:);
        data{2} = data{2}(g.chans,:,:,:);
    else
        data = data(g.chans,:,:,:);
    end
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%% Print plot info %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(curfig); h=gca;title(g.title, 'interpreter', 'none'); % title plot
hold on
axis('off')

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Read chan_locs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if isempty(g.chanlocs) || ~isfield(g.chanlocs, 'theta') % plot in a rectangular grid
    ISRECT = 1;
    ht = g.geom(1);
    wd = g.geom(2);
    if length(g.chans) > ht*wd
        fprintf('metaplottopo(): (%d) channels to be plotted > grid size [%d %d]\n',...
            length(g.chans),ht,wd);
        return
    end
    halfht = (ht-1)/2;
    halfwd = (wd-1)/2;
    xvals = zeros(ht*wd,1);
    yvals = zeros(ht*wd,1);
    dist  = zeros(ht*wd,1);
    for j=1:ht
        for i=1:wd
            xvals(i+(j-1)*wd) = -halfwd+(i-1);
            yvals(i+(j-1)*wd) =  halfht-(j-1);
            % dist(i+(j-1)*wd) =  sqrt(xvals(j+(i-1)*ht).^2+yvals(j+(i-1)*ht).^2);
        end
    end
    % maxdist = max(dist);
    maxxvals = max(xvals);
    maxyvals = max(yvals);
    for j=1:ht
        for i=1:wd
            xvals(i+(j-1)*wd) = 0.499*xvals(i+(j-1)*wd)/maxxvals;
            yvals(i+(j-1)*wd) = 0.499*yvals(i+(j-1)*wd)/maxyvals;
        end
    end
    if ~exist('channames')
        channames = repmat(' ',ht*wd,4);
        for i=1:ht*wd
            channum = num2str(i);
            channames(i,1:length(channum)) = channum;
        end
    end

else % read chan_locs file
    % read the channel location file
    % ------------------------------
    if isstruct(g.chanlocs)
        nonemptychans = cellfun('isempty', { g.chanlocs.theta });
        nonemptychans = find(~nonemptychans);
        [tmp channames Th Rd] = readlocs(g.chanlocs(nonemptychans));
        channames = strvcat({ g.chanlocs.labels });
    else
        [tmp channames Th Rd] = readlocs(g.chanlocs);
        channames = strvcat(channames);
        nonemptychans = [1:length(channames)];
    end
    Th = pi/180*Th;                 % convert degrees to radians
    Rd = Rd;

    if length(g.chans) > length(g.chanlocs),
        error('metaplottopo(): data channels must be <= ''chanlocs'' channels')
    end

    [yvalstmp,xvalstmp] = pol2cart(Th,Rd); % translate from polar to cart. coordinates
    xvals(nonemptychans) = xvalstmp;
    yvals(nonemptychans) = yvalstmp;

    % find position for other channels
    % --------------------------------
    totalchans = length(g.chanlocs);
    emptychans = setdiff_bc(1:totalchans, nonemptychans);
    totalchans = floor(sqrt(totalchans))+1;
    for index = 1:length(emptychans)
        xvals(emptychans(index)) = 0.7+0.2*floor((index-1)/totalchans);
        yvals(emptychans(index)) = -0.4+mod(index-1,totalchans)/totalchans;
    end
    channames = channames(g.chans,:);
    xvals     = xvals(g.chans);
    yvals     = yvals(g.chans);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% xvals = 0.5+PLOT_WIDTH*xvals;   % controls width of  plot array on page!
% yvals = 0.5+PLOT_HEIGHT*yvals;  % controls height of plot array on page!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(xvals) > 1
    xvals = (xvals-mean([max(xvals) min(xvals)]))/(max(xvals)-min(xvals)); % recenter
    xvals = gcapos(1)+gcapos(3)/2+PLOT_WIDTH*xvals;   % controls width of plot
    % array on current axes
end
yvals = gcapos(2)+gcapos(4)/2+PLOT_HEIGHT*yvals;  % controls height of plot
% array on current axes
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot traces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

Axes = [];
fprintf('Plotting all channel...');
for c=1:length(g.chans), %%%%%%%% for each data channel %%%%%%%%%%%%%%%%%%%%%%%%%%

    xcenter = xvals(c); if isnan(xcenter), xcenter = 0.5; end;
    ycenter = yvals(c); if isnan(ycenter), ycenter = 0.5; end
    Axes = [Axes axes('Units','Normal','Position', ...
        [xcenter-axwidth/2 ycenter-axheight/2 axwidth axheight])];
    hold on;
    axis('off');
    %axes('Units','Normal','Position', [xcenter-axwidth/2 ycenter-axheight/2 axwidth axheight])
    %axis('off');

    %
    %%%%%%%%%%%%%%%%%%%%%%% Plot data traces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ~isempty( g.plotfunc )
        %figure(curfig);
        eval( [ 'func = @' g.plotfunc ';' ] );
        if iscell(data), tmp = { g.plotargs{1:g.datapos(1)-1} data{1}(c,:,:,:,:) g.plotargs{g.datapos(1):g.datapos(2)-1} data{2}(c,:) g.plotargs{g.datapos(2):end}};
        else             tmp = { g.plotargs{1:g.datapos-1}    data(c,:,:,:,:)    g.plotargs{g.datapos:end} };
        end
        if strcmpi(g.squeeze, 'on') tmp{g.datapos} = squeeze(tmp{g.datapos}); end
        tmp = { tmp{:} 'title' channames(c,:) 'plotmode' 'topo'};
        feval(func, tmp{:});
    end
    outchannames{c} = deblank(channames(c,:));
end % c, chans / subplot

fprintf('\n');

%
%%%%%%%%%%%%%%%%%%%%% Make time and amp cal bar %%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~isempty(g.calbar)
    ax = axes('Units','Normal','Position', [0.85 0.1 axwidth axheight]); % FIX!!!!
    axes(ax)
    axis('off');
    if   g.calbar(3) < g.calbar(4), g.ydir = 1;
    else g.calbar(5) = g.calbar(3); g.calbar(3) = []; g.ydir = -1;
    end
    [xmin xmax ymin ymax] = deal(g.calbar(1), g.calbar(2),g.calbar(3), g.calbar(4));
    figure(curfig);p=plot([0 0],[ymin ymax],'color','k'); % draw vert axis at zero
    if g.ydir == -1
        set(gca, 'ydir', 'reverse');
    end
    hold on
    figure(curfig);p=plot([xmin xmax],[0 0],'color','k'); % draw horizontal axis
    xlim([xmin xmax]);
    ylim([ymin ymax]);
    axis off;

    % draw text limits
    % ----------------
    xdiff = xmax - xmin;
    ydiff = ymax - ymin;
    signx = xmin-0.15*xdiff;
    figure(curfig);axis('off');h=text(double(signx),double(ymin),num2str(ymin,3)); % text ymin
    set(h,'FontSize',TICKFONTSIZE);               % choose font size
    set(h,'HorizontalAlignment','right','Clipping','off');

    figure(curfig);axis('off');h=text(double(signx), double(ymax),['+' num2str(ymax,3)]);  % text +ymax
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','right','Clipping','off');

    ytick = g.ydir*(-ymax-0.3*ydiff);
    figure(curfig);tick = [int2str(xmin)]; h=text(double(xmin),double(ytick),tick);
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','center',...
                        'Clipping','off');  % center text

    tick = 'Time (ms)'; figure(curfig);h=text(double(xmin+xdiff/2),double(ytick-0.5*g.ydir*ydiff),tick);
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','center',...
                        'Clipping','off');  % center text

    tick = [int2str(xmax)]; figure(curfig);h=text(double(xmax),double(ytick),tick);
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','center',...
                        'Clipping','off');  % center text

end

%        'set(gcbf, ''''unit'''', ''''pixel'''');' ...
%        'tmp = get(gcbf, ''''position'''');' ...
%        'tmp2 = get(0, ''''screensize'''');' ...
%        'if tmp(2)+500 > tmp2(4), tmp(2) = tmp2(4)-500; end;' ...
%        'set(gcbf, ''''position'''', [ tmp(1) tmp(2) 560   420]);' ...
if ~strcmpi(g.axcopycom, 'off')
    if length(g.chans) > 100
        fprintf('Disabling axis copy feature to speed up plotting');
    else
        if strcmpi(g.axcopycom, 'on')
            g.axcopycom = [ 'axis on;' ...
                    'tmp = get(gca, ''''userdata'''');' ...
                    'if ~isempty(tmp), xlabel(tmp{1});' ...
                    'ylabel(tmp{2});' ...
                    'if length(tmp) == 3 && ~isempty(tmp{3}), legend(tmp{3}{:}); end;' ...
                    'end; clear tmp tmp2;' ];
        end
        axcopy(gcf, g.axcopycom); % turn on popup feature
    end
end
icadefs;
set(gca,'Color',BACKCOLOR);               % set the background color

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% POP_ERPIMAGE - draw an ERP-image plot of a given EEG channel or independent
%                  component. Uses a pop-up window if less than three (or four
%                  in one condition) input arguments are supplied. Calls ERPIMAGE.
%                  For further details see >> help erpimage
% Usage:
%   >> pop_erpimage(EEG, typeplot);          % pop-up a data entry window
%   >> pop_erpimage(EEG, typeplot, lastcom); % pop-up a data entry window
%   >> pop_erpimage(EEG, typeplot, channel); % no pop-up window
%   >> pop_erpimage(EEG, typeplot, channel, projchan, title, ...
%                  smooth, decimate, sortingtype, sortingwin, ...
%                            sortingeventfield, renorm, options...);
% Graphic interface:
%   "Channel or Component" - [edit box] Enter channel number or component
%                 number to plot. ERPIMAGE equivalent: 'channel'
%   "Project to channel #" - [edit box] (for plotting independent components).
%                 Allow reprojecting the component activity
%                 to a given channel or group of channels.
%                 ERPIMAGE equivalent: [none]
%   "Smoothing" - [text box] Smoothing parameter in number of trials.
%                 ERPIMAGE equivalent: 'smooth'
%   "Downsampling" - [edit box] Decimate parameter.
%                 ERPIMAGE equivalent: 'decimate'
%   "Time limits" - [edit box] Enter the time limits in milliseconds.
%                 ERPIMAGE equivalent: the 1st and 2nd parameters of the 'limit' array
%   "Figure title" - [edit box] Enter the figure title here.  If empty, a title
%                 is automatically generated. ERPIMAGE equivalent: 'title'
%   "Plot scalp map" - [checkbox] Setting this option plot a scalp map of the
%                 channel location (or component topography) next to the
%                 erpimage. ERPIMAGE equivalent: 'topo'
%   "plot ERP" - [checkbox] Setting this option plot the channel or component
%                 ERP below the ERP image. ERPIMAGE equivalent: 'erp'
%   "Plot colorbar" - [checkbox] Plot the colorbar on the right of the erpimage.
%                 ERPIMAGE equivalent: 'cbar'
%   "ERP limits" - [edit box] Set the minimum and maximum value for the ERP plot
%                 ERPIMAGE equivalent: 3rd and 4th parameters of the 'limit' array
%   "Color limits" - [edit box] Set the color limits for the ERP image.
%                 ERPIMAGE equivalent: 'caxis'
%   "Epoch sorting field" - [button and edit box] Specify the event field which
%                 values will be used to sort the trials. For instance, if you
%                 select the 'latency' fields, trials will be sorted by the
%                 latency of the selected events.
%                 ERPIMAGE equivalent: 'sortingeventfield'
%   "Event type(s)" - [button and edit box] Specify which event subset to select,
%                 based on the event 'type' field values (to scan for event types, use
%                 menu Edit > Events values and look at the values of the 'type'
%                 field). For instance, entering type 'rt' (if defined) and field
%                 'latency' in option above will sort trials on reaction time latencies.
%                 When several selected events are present in individual trials,
%                 the first event values are used for sorting and a warning is issued.
%                 ERPIMAGE equivalent: 'sortingtype'
%   "Event time range" - [edit box] Specify which event subset to select based on
%                 event latency values. As the option above, this further restrains
%                 the selection of events. For example, entering [200 300] in this
%                 box, 'rt' for the event type (above), and 'latency' for the
%                 epoch sorting field will select trials with reaction-time latencies
%                 in between 200 and 300 ms. Trials with no such event will not be
%                 included in the ERP-image plot. ERPIMAGE equivalent: 'sortingwin'
%   "rescale" - [edit box] 'yes', 'no', or a Matlab formula.
%                 ERPIMAGE equivalent: 'renorm'
%   "align" - [edit box] Set this to 'Inf' to re-align the individual trials
%                 on the median latency of the selected events. Else, enter an epoch time
%                 (in ms) to align the events to (Ex: 0). ERPIMAGE equivalent: 'align'
%   "Don't sort by value" - [checkbox] Check this box if you do not want to
%                 sort the trials but do want to plot the selected event values.
%                 ERPIMAGE equivalent: 'nosort'
%   "Don't plot value" - [checkbox] Check this box if you do not want to
%                 plot the selected event values, but still want to sort
%                 the data trials according to these values.
%                 ERPIMAGE equivalent: 'noplot'
%   "Sort by phase > Frequency" - [edit box] Specify the frequency or frequency
%                 range for sorting trials by phase. ERPIMAGE equivalent:
%                 3rd and 4th inputs to 'phasesort'
%   "Window center (ms)" - [edit box] ERPIMAGE equivalent: 1st 'phasesort' input
%   "Percent low-amp. trials to ignore" - [edit box] ERPIMAGE equivalent:
%                 2nd 'phasesort' input
%   "Wavelet cycles" - [text] Number of wavelet cycles used for spectral decomposition
%                 at the specified latency. To change this, see "More options" {default: 3}
%   "Inter-trial coherence options > Frequency" - [edit box] Frequency at which
%                 to compute coherence. Constrained to be the same as the
%                 "Sort by phase > Frequency" edit box. ERPIMAGE equivalent: 'coher'
%   "Signif. level" - [edit box] Coherence significance cutoff, as a proability
%                  (Ex: .05). ERPIMAGE equivalent: 'signif'
%   "Amplitude limit" - [edit box] Amplitude limits [min max] for the data power
%                 plot at the selected frequency. ERPIMAGE equivalent:
%                 5th and 6th inputs to 'limit'
%   "Coher limits" - [edit box] Upper limit (<=1) for the coherence
%                 plot. ERPIMAGE equivalent: 7th and 8th inputs of 'limit'
%   "Image amps" - [checkbox] Check this box for plotting the spectral amplitude
%                 image at the selected frequency (instead of plotting EEG potential).
%                 ERPIMAGE equivalent: 'plotamp'
%   "Plot spectrum" - [edit box] Plot the channel or component data spectrum in
%                 the top right corner of the ERP image. ERPIMAGE equivalent: 'spec'
%   "Baseline ampl." - [edit box] Baseline amplitude for data power plot at the
%                 selected frequency. ERPIMAGE equivalent: 7th inputs of 'limit'
%   "Mark times" - [edit box] Time(s) in ms to plot vertical lines.
%                 ERPIMAGE equivalent: 'vert'
%   "More options" - [edit box] Enter 'key', 'value' sequences. Other ERPIMAGE
%                 options not handled by this interface, including: 'erpstd' to
%                 plot the ERP standard deviation; 'auxvar' to plot auxiliary
%                 variables; 'ampsort' to sort trials based on amplitude at
%                 the selected frequency, etc.  For further information see
%                 >> help ERPIMAGE
% Inputs:
%   EEG        - dataset structure
%   typeplot   - 1=channel, 0=component {default: 1}
%   lastcom    - string containing previous pop_erpimage command (from LASTCOM)
%                or from the previous function call output.  The values from this
%                function call are used as default in the graphic interface.
%
% Commandline options:
%   channel    - Index of channel or component(s) to plot {default: 1}
%   projchan   - Channel to back-project the selected component(s) to.
%                If plotting channel activity, this argument is ignored.
%                If [], the ICA component activation is plotted {default []}.
%   title      - ['string'] Plot title {default: []}
%   smooth     - Smoothing parameter (number of trials). {Default: 5}
%                ERPIMAGE equivalent: 'avewidth'
%   decimate   - Decimate parameter (i.e. ratio of trials_in/trials_out).
%                ERPAIMGE equivalent: 'decimate' {Default: 0}
%   sortingtype - Sorting event type(s) ([int vector]; []=all). See Notes below.
%                Either a string or an integer.
%   sortingwin - Sorting event window [start, end] in seconds ([]=whole epoch)
%   sortingeventfield - Sorting field name. {default: none}. See Notes below.
%   options    - ERPIMAGE options, separated by commas (Ex: 'erp', 'cbar').
%                {Default: none}. For further details see >> erpimage help
%
% Outputs from pop-up:
%   String containing the command used to evaluate this plotting function
%   (saved by EEGLAB as LASTCOM). Enter it as the 'lastcom' input to restore
%   the previous parameters as defaults in a new POP_ERPIMAGE pop-up window
%
% Outputs from commandline:
%   Same as ERPIMAGE. Note: No outputs are returned when a window pops-up
%   to ask for additional arguments
%
% Notes:
%   1) A new figure is created only when the pop-up window is called,
%      so you may call this command with >3 args to draw in SBPLOT axes.
%   2) To sort epochs, first define the event field to be used with
%      the argument 'sortingeventfield' (for instance 'latency'). Then,
%      because there may be several events with different latencies in a
%      given epoch, it is possible to consider only a subsets of events
%      using the 'sortingtype' and 'sortingwin' arguments. The
%      'sortingtype' argument selects events with definite types. The
%      'sortingwin' argument helps to define a specific time window in the
%      epoch to select events. For instance the epoch time range may be -1
%      to 2 seconds but one may want to select events only in the range 0
%      to 1 second. These three parameters are forwarded to the function
%      EEG_GETEPOCHEVENT, whose help message contains more details.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGLAB, ERPIMAGE, EEG_GETEPOCHEVENT

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad
% 02-12-02 added new event format compatibility -ad
% 02-15-02 text interface editing -sm & ad
% 03-07-02 add the eeglab computeica options -ad
% 02-15-02 modified the function according to the new event/epoch structure -ad
% 03-18-02 added title -ad & sm
% 04-04-02 added outputs -ad & sm

function varargout = pop_erpimage( EEG, typeplot, channel, projchan, titleplot, smooth, decimate, sortingtype, ...
            sortingwin, sortingeventfield, varargin)

varargout{1} = '';
if nargin < 1
   help pop_erpimage;
   return;
end

if typeplot == 0 && isempty(EEG.icasphere)
   error('no ICA data for this set, first run ICA');
end;
if EEG.trials == 1
   error('erpimage of one trial cannot be plotted');
end;

if nargin < 2
	typeplot = 1; %1=signal; 0=component
end
lastcom = [];
if nargin < 3
	popup = 1;
else
	popup = ischar(channel) | isempty(channel);
	if ischar(channel)
		lastcom = channel;
	end
end

opt = [];
if popup
	% get contextual help
	% -------------------
    erpimagefile = which('erpimage.m');
	[txt2 vars2] = gethelpvar(erpimagefile);
	txt  = { txt2{:}};
	vars = { vars2{:}};
	% [txt vars] = gethelpvar('erpimopt.m');
	% txt  = { txt{:} txt2{:}};
	% vars = { vars{:} vars2{:}};

    opt.topo         = getkeyval(lastcom, 'topo', 'present', 1);
    opt.fieldname    = getkeyval(lastcom, 10);
    opt.type         = getkeyval(lastcom, 8);
    opt.renorm       = getkeyval(lastcom, 'renorm','', 'no');
    opt.erp          = fastif(getkeyval(lastcom, '''erp''', 'present', 1), 'on', 'off');
    opt.cbar         = fastif(getkeyval(lastcom, 'cbar', 'present', 1), 'on', 'off');
    opt.nosort       = fastif(getkeyval(lastcom, 'nosort', 'present', 0), 'on', 'off');
    opt.noplot       = fastif(getkeyval(lastcom, 'noplot', 'present', 0), 'on', 'off');
    opt.plotamps     = fastif(getkeyval(lastcom, 'plotamps', 'present', 0), 'on', 'off');
    opt.index        = str2num(getkeyval(lastcom,3,[],'1'));
    opt.smoothing    = str2num(getkeyval(lastcom, 6, [], int2str(min(max(EEG.trials-5,0), 10))));
    opt.downsampling = str2num(getkeyval(lastcom, 7, [], '1'));
    opt.caxis        = str2num(getkeyval(lastcom, 'caxis'));
    opt.eventrange   = str2num(getkeyval(lastcom, 9));
    opt.align        = str2num(getkeyval(lastcom, 'align'));
    opt.phasesort    = str2num(getkeyval(lastcom, 'phasesort'));
    opt.coher        = str2num(getkeyval(lastcom, 'coher'));
    opt.spec         = str2num(getkeyval(lastcom, 'spec'));
    opt.vert         = str2num(getkeyval(lastcom, 'vert'));
    opt.limits       = str2num(getkeyval(lastcom, 'limits'));
    opt.limits       = [ opt.limits NaN NaN NaN NaN NaN NaN NaN NaN NaN ]; opt.limits = opt.limits(1:9);
    opt.coher        = [ opt.coher NaN NaN NaN NaN NaN NaN NaN NaN NaN ];  opt.coher  = opt.coher(1:3);
    opt.phasesort    = [ opt.phasesort NaN NaN NaN NaN NaN NaN NaN NaN NaN ]; opt.phasesort = opt.phasesort(1:4);
    if isnan(opt.limits(1)), opt.limits(1:2) = 1000*[EEG.xmin EEG.xmax]; end

	commandphase = [ 'if ~isempty(get(findobj(''parent'', gcbf, ''tag'', ''phase''),''string'')),' ...
					 '   if ~isempty(get(findobj(''parent'', gcbf, ''tag'', ''coher''),''string'')), ' ...
					 '      set(findobj(''parent'', gcbf, ''tag'', ''coher''), ''string'', get(findobj(''parent'', gcbf, ''tag'', ''phase''),''string''));' ...
					 'end; end;' ];
	commandcoher = [ 'if ~isempty(get(findobj(''parent'', gcbf, ''tag'', ''coher''),''string'')),' ...
					 '   if ~isempty(get(findobj(''parent'', gcbf, ''tag'', ''phase''),''string'')), ' ...
					 '      set(findobj(''parent'', gcbf, ''tag'', ''phase''), ''string'', get(findobj(''parent'', gcbf, ''tag'', ''coher''),''string''));' ...
					 'end; end;' ];

	commandfield = ['tmpEEG = get(gcbf, ''userdata''); if isempty(tmpEEG.event)' ...
				   '   errordlg2(''No events'');' ...
				   'else' ...
				   '   tmpfieldnames = fieldnames(tmpEEG.event);' ...
				   '   [tmps,tmpv] = listdlg2(''PromptString'', ''Select fields'', ''SelectionMode'',''single'',''ListString'', tmpfieldnames);' ...
				   '   if tmpv' ...
				   '       set(findobj(''parent'', gcbf, ''tag'', ''field''), ''string'', tmpfieldnames{tmps});' ...
				   '   end;' ...
				   'end;' ...
                   'if isempty(get(findobj(''parent'', gcbf, ''tag'', ''type''), ''string'')),' ...
                   '   warndlg2(''Do not forget to select an event type in the next edit box'');' ...
                   'end;' ...
				   'clear tmpEEG tmps tmpv tmpfieldnames;' ];
	commandtype = [ 'tmpEEG = get(gcbf, ''userdata''); if ~isfield(tmpEEG.event, ''type'')' ...
				   '   errordlg2(''No type field'');' ...
				   'else' ...
                   '   tmpevent = tmpEEG.event;' ...
                   '   if isnumeric(tmpEEG.event(1).type),' ...
				   '        [tmps,tmpstr] = pop_chansel(unique([ tmpevent.type ]));' ...
				   '   else,' ...
                   '        [tmps,tmpstr] = pop_chansel(unique({ tmpevent.type }));' ...
                   '   end;' ...
				   '   if ~isempty(tmps)' ...
				   '       set(findobj(''parent'', gcbf, ''tag'', ''type''), ''string'', tmpstr);' ...
				   '   end;' ...
				   'end;' ...
                   'if isempty(get(findobj(''parent'', gcbf, ''tag'', ''field''), ''string'')),' ...
                   '   warndlg2(''Do not forget to select an event type in the previous edit box'');' ...
                   'end;' ...
				   'clear tmpEEG tmps tmpv tmpstr tmpevent tmpfieldnames;' ];
    commandchan = 'tmpEEG = get(gcbf, ''userdata''); tmpchanlocs = tmpEEG(1).chanlocs; [tmp tmpval] = pop_chansel({tmpchanlocs.labels}, ''withindex'', ''on'', ''selectionmode'', ''single''); set(findobj(gcbf, ''tag'', ''chan''), ''string'',tmpval); clear tmpEEG tmp tmpchanlocs tmpval';

	geometry = { [1 1 0.4 0.5 2.1] [1 1 1 1 1] [1 1 1 1 1] [1 1 1 1 1] [1] [1] [1 1 1 0.8 0.8 1.2] [1 1 1 0.8 0.8 1.2] [1] [1] ...
				 [1.6 1.7 1.2 1 .5] [1.6 1.7 1.2 1 .5] [1] [1] [1.5 1 1 1 1] [1.5 1 1 1 1] [1] [1] [1.5 1 1 2.2] [1.5 1 1 2.2]};
    uilist = { { 'Style', 'text', 'string', fastif(typeplot, 'Channel', 'Component(s)'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', num2str(opt.index), 'tag', 'chan' } ...
               { 'style' 'pushbutton' 'string'  '...', 'enable' fastif(~isempty(EEG(1).chanlocs) && typeplot, 'on', 'off') 'callback' commandchan } ...
			   { 'Style', 'text', 'string', 'Fig. title', 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', '', 'tag', 'title'  } ...
			   ...
			   { 'Style', 'text', 'string', 'Smoothing', 'fontweight', 'bold', 'tooltipstring', context('avewidth',vars,txt) } ...
			   { 'Style', 'edit', 'string', num2str(opt.smoothing), 'tag', 'smooth' } ...
			   { 'Style', 'checkbox', 'string', 'Plot scalp map', 'tooltipstring', 'plot a 2-d head map (vector) at upper left', ...
				 'value', opt.topo, 'tag', 'plotmap'  } { } { } ...
			   { 'Style', 'text', 'string', 'Downsampling', 'fontweight', 'bold', 'tooltipstring', context('decimate',vars,txt) } ...
			   { 'Style', 'edit', 'string', num2str(opt.downsampling), 'tag', 'decimate' } ...
			   { 'Style', 'checkbox', 'string', 'Plot ERP', 'tooltipstring', context('erp',vars,txt), 'value', fastif(strcmpi(opt.erp, 'on'), 1,0), 'tag', 'erp' } ...
			   { 'Style', 'text', 'string', fastif(typeplot, 'ERP limits (uV)','ERP limits'), 'tooltipstring', [ 'Plotting limits for ERP trace [min_uV max_uV]' 10 '{Default: ERP data limits}']  } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.limits(3)), num2str(opt.limits(3:4)), ''), 'tag', 'limerp'  } ...
			   { 'Style', 'text', 'string', 'Time limits (ms)', 'fontweight', 'bold', 'tooltipstring',  'Select time subset in ms' } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.limits(1)), num2str(opt.limits(1:2)), ''), 'tag', 'limtime' } ...
			   { 'Style', 'checkbox', 'string', 'Plot colorbar','tooltipstring', context('caxis',vars,txt), 'value', fastif(strcmpi(opt.cbar, 'on'), 1,0), 'tag', 'cbar' } ...
			   { 'Style', 'text', 'string', 'Color limits (see Help)','tooltipstring', context('caxis',vars,txt) } ...
			   { 'Style', 'edit', 'string', num2str(opt.caxis), 'tag', 'caxis' } ...
			   {} ...
			   { 'Style', 'text', 'string', 'Sort/align trials by epoch event values', 'fontweight', 'bold'} ...
			   { 'Style', 'pushbutton', 'string', 'Epoch-sorting field', 'callback', commandfield, ...
				 'tooltipstring', 'Epoch-sorting event field name (Ex: latency; default: no sorting):' } ...
			   { 'Style', 'pushbutton', 'string', 'Event type(s)', 'callback', commandtype, 'tooltipstring', ['Event type(s) subset (default: all):' 10 ...
                '(See ''/Edit/Edit event values'' for event types)']  } ...
			   { 'Style', 'text', 'string', 'Event time range', 'tooltipstring', [ 'Sorting event window [start, end] in milliseconds (default: whole epoch):' 10 ...
												  'events are only selected within this time window (can be useful if several' 10 ...
												  'events of the same type are in the same epoch, or for selecting trials with given response time)']} ...
			   { 'Style', 'text', 'string', 'Rescale', 'tooltipstring', 'Rescale sorting variable to plot window (yes|no|a*x+b)(Ex:3*x+2):' } ...
			   { 'Style', 'text', 'string', 'Align', 'tooltipstring',  context('align',vars,txt) } ...
			   { 'Style', 'checkbox', 'string', 'Don''t sort by value', 'tooltipstring', context('nosort',vars,txt), 'value', fastif(strcmpi(opt.nosort, 'on'), 1,0), 'tag', 'nosort' } ...
			   { 'Style', 'edit', 'string', opt.fieldname, 'tag', 'field' } ...
			   { 'Style', 'edit', 'string', opt.type, 'tag', 'type' } ...
			   { 'Style', 'edit', 'string', num2str(opt.eventrange), 'tag', 'eventrange' } ...
			   { 'Style', 'edit', 'string', opt.renorm, 'tag', 'renorm' } ...
			   { 'Style', 'edit', 'string', num2str(opt.align),  'tag', 'align' } ...
			   { 'Style', 'checkbox', 'string', 'Don''t plot values', 'tooltipstring', context('noplot',vars,txt), 'value', fastif(strcmpi(opt.noplot, 'on'), 1,0), 'tag', 'noplot' } ...
			   {} ...
			   { 'Style', 'text', 'string', 'Sort trials by phase', 'fontweight', 'bold'} ...
			   { 'Style', 'text', 'string', 'Frequency (Hz | minHz maxHz)', 'tooltipstring', ['sort by phase at maximum-power frequency' 10 ...
               'in the data within the range [minHz,maxHz]' 10 '(overrides frequency specified in ''coher'' flag)']  } ...
			   { 'Style', 'text', 'string', 'Percent low-amp. trials to ignore', 'tooltipstring', ['percent of trials to reject for low' ...
			    'amplitude. Else,' 10 'if prct is in the range [-100,0] -> percent to reject for high amplitude'] } ...
			   { 'Style', 'text', 'string', 'Window center (ms)', 'tooltipstring', 'Center time of the n-cycle window'  } ...
			   { 'Style', 'text', 'string', 'Wavelet cycles', 'tooltipstring', 'cycles per wavelet window'  } {}...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.phasesort(3)),num2str(opt.phasesort(3:4)),'') 'tag', 'phase', 'callback', commandphase } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.phasesort(2)),num2str(opt.phasesort(2)),''), 'tag', 'phase2' } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.phasesort(1)),num2str(opt.phasesort(1)),''), 'tag', 'phase3' } ...
			   { 'Style', 'text', 'string', '        3'  } {}...
			   {} ...
			   { 'Style', 'text', 'string', 'Inter-trial coherence options', 'fontweight', 'bold'} ...
			   { 'Style', 'text', 'string', 'Frequency (Hz | minHz maxHz)', 'tooltipstring', [ '[freq] -> plot erp plus amp & coher at freq (Hz)' 10 ...
               '[minHz maxHz] -> find max in frequency range' 10 '(or at phase freq above, if specified)']} ...
			   { 'Style', 'text', 'string', 'Signif. level (<0.20)', 'tooltipstring', 'add coher. signif. level line at alpha (alpha range: (0,0.1])' } ...
			   { 'Style', 'text', 'string', 'Amplitude limits (dB)'  } ...
			   { 'Style', 'text', 'string', 'Coher limits (<=1)'  } ...
			   { 'Style', 'checkbox', 'string', 'Image amps', 'tooltipstring',  context('plotamps',vars,txt), 'value', fastif(strcmpi(opt.plotamps, 'on'), 1,0), 'tag', 'plotamps' } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.coher(1)), num2str(opt.coher(1:2)), ''), 'tag', 'coher', 'callback', commandcoher } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.coher(3)), num2str(opt.coher(3)),   ''), 'tag', 'coher2'  } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.limits(5)), num2str(opt.limits(5:6)), ''), 'tag', 'limamp' } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.limits(7)), num2str(opt.limits(7:8)), ''), 'tag', 'limcoher' } ...
               {'style', 'text', 'string', '   (Requires signif.)' } ...
			   {} ...
			   { 'Style', 'text', 'string', 'Other options', 'fontweight', 'bold'} ...
			   { 'Style', 'text', 'string', 'Plot spectrum (minHz maxHz)','tooltipstring',  context('spec',vars,txt)} ...
			   { 'Style', 'text', 'string', 'Baseline ampl. (dB)', 'tooltipstring', 'Use it to fix baseline amplitude' } ...
			   { 'Style', 'text', 'string', 'Mark times (ms)','tooltipstring',  context('vert',vars,txt)} ...
			   { 'Style', 'text', 'string', 'More options (see >> help erpimage)' } ...
			   { 'Style', 'edit', 'string', num2str(opt.spec), 'tag', 'spec' } ...
			   { 'Style', 'edit', 'string', fastif(~isnan(opt.limits(9)), num2str(opt.limits(9)), ''), 'tag', 'limbaseamp' } ...
			   { 'Style', 'edit', 'string', num2str(opt.vert), 'tag', 'vert' } ...
			   { 'Style', 'edit', 'string', '', 'tag', 'others' } ...
			};
	if typeplot == 0 % add extra param for components
		geometry = { [1 1 0.1 0.8 2.1] geometry{:} };
		uilist   = { { } { } { } { } { } uilist{:}};
		uilist{1} = uilist{6};
		uilist{2} = uilist{7};
		uilist{6} = { 'Style', 'text', 'string', 'Project to channel #', 'fontweight', 'bold','tooltipstring', ['Project component(s) to data channel' 10 ...
												  'This allows plotting projected component activity at one channel in microvolts'] };
		uilist{7} = { 'Style', 'edit', 'string', getkeyval(lastcom, 4), 'tag', 'projchan' };
	end
    [oldres a b res] = inputgui( 'geometry', geometry, 'uilist', uilist, 'helpcom', 'pophelp(''pop_erpimage'');', ...
							 'title', fastif( typeplot, 'Channel ERP image -- pop_erpimage()', 'Component ERP image -- pop_erpimage()'), 'userdata', EEG);
	if isempty(oldres), return; end

	% first rows
	% ---------
    channel = eeg_decodechan(EEG.chanlocs, res.chan);
	titleplot    = res.title;
	if isfield(res, 'projchan'),
        if ~isempty(res.projchan)
            if strcmpi(res.projchan(1),'''')
                 projchan = eval( [ '{' res.projchan '}' ]);
            else projchan = parsetxt( res.projchan);
            end
            if ~isempty(projchan) && ~isempty(str2num(projchan{1}))
                projchan = cellfun(@str2num, projchan);
            end
        else
            projchan = [];
        end
    else,
        projchan = [];
    end
    opt = [];
	if ~isempty(res.others)
        try,
            tmpcell = eval( [ '{' res.others '}' ] );
            opt = struct( tmpcell{:} );
        catch, error('Additional options ("More options") requires ''key'', ''val'' arguments');
        end
	end
	if ~typeplot && isempty(projchan)
        opt.yerplabel = '';
    else
        opt.yerplabel = '\muV' ;
	end
	smooth       = eval(res.smooth);
    if res.plotmap
		if isfield(EEG.chanlocs, 'theta')
            if ~isfield(EEG, 'chaninfo'), EEG.chaninfo = []; end
			if typeplot == 0
				     opt.topo = [ ' { mean(EEG.icawinv(:,[' int2str(channel) ']),2) EEG.chanlocs EEG.chaninfo } '];
			else     opt.topo = [ ' { [' int2str(channel) '] EEG.chanlocs EEG.chaninfo } '];
			end;
		end
	end

	decimate     = eval( res.decimate );
	if res.erp
		opt.erp = 'on';
	end

	% finding limits
	% --------------
	limits(1:8)  = NaN;
	if ~isempty(res.limerp)
		limits(3:4) = eval( [ '[' res.limerp ']' ]);
	end
	if ~isempty(res.limtime) % time limits
		if ~strcmp(res.limtime, num2str(1000*[EEG.xmin EEG.xmax]))
			limits(1:2) = eval( [ '[' res.limtime ']' ]);
		end
	end
	if ~isempty(res.limamp)
		limits(5:6) = eval( [ '[' res.limamp ']' ]);
	end
	if ~isempty(res.limcoher)
		limits(7:8) = eval( [ '[' res.limcoher ']' ]);
	end
	if ~isempty(res.limbaseamp)
		limits(9) = eval( res.limbaseamp ); %bamp
	end
	if ~all(isnan(limits))
        opt.limits = limits;
	end

	% color limits
	% --------------
	if res.cbar
		opt.cbar = 'on';
	end
	if res.caxis
		opt.caxis = str2num(res.caxis);
	end

	% event rows
	% ----------
	if res.nosort
		opt.nosort = 'on';
	end
	try, sortingeventfield = eval( res.field ); catch, sortingeventfield = res.field; end
    if ~isempty(res.type)
       if strcmpi(res.type(1),'''')
            sortingtype = eval( [ '{' res.type '}' ] );
       else sortingtype = parsetxt( res.type );
       end
    end
	sortingwin   = eval( [ '[' res.eventrange ']' ] );
	if ~isempty(res.field) && ~strcmp(res.renorm, 'no')
		opt.renorm = res.renorm;
	end
	if ~isempty(res.align)
		opt.align = str2num(res.align);
	end
	if res.noplot
		opt.noplot = 'on';
	end

	% phase rows
	% ----------
	tmpphase = [];
	if ~isempty(res.phase)
		tmpphase = eval( [ '[ 0 0 ' res.phase ']' ]);
	end
	if ~isempty(res.phase2)
		tmpphase(2) = eval( res.phase2 );
	end
	if ~isempty(res.phase3)
		tmpphase(1) = eval( res.phase3 );
	end
	if ~isempty(tmpphase)
		opt.phasesort = tmpphase;
	end

	% coher row
	% ----------
	tmpcoher = [];
	if res.plotamps
		opt.plotamps = 'on';
	end
	if ~isempty(res.coher)
		tmpcoher = eval( [ '[' res.coher ']' ]);
	end
	if ~isempty(res.coher2)
		if length(tmpcoher) == 1
			tmpcoher(2) = tmpcoher(1);
		end
		tmpcoher(3) = eval( res.coher2 );
	end
	if ~isempty(tmpcoher)
		opt.coher = tmpcoher;
	end

	% options row
	% ------------
	if ~isempty(res.spec)
		opt.spec = eval( [ '[' res.spec ']' ]);
	end
	if ~isempty(res.vert)
		opt.vert = eval( [ '[' res.vert ']' ]);
	end
	figure;
    options = '';
else
	options = '';
	if nargin < 4
		projchan = [];
	end
	if nargin < 5
		titleplot = '';
	end
	if nargin < 6
		smooth = 5;
	end
	if nargin < 7
		decimate = 0;
	end
	if nargin < 8
		sortingtype = [];
	end
	if nargin < 9
		sortingwin = [];
	end
	if nargin < 10
		sortingeventfield = [];
	end
    %options = vararg2str(varargin); % NO BECAUSE OF THE CHANNEL LOCATION
    %                                  PROBLEM BELOW
	for i=1:length( varargin )
		if ischar( varargin{ i } )
			options = [ options ', ''' varargin{i} '''' ];
		else
		  if ~iscell( varargin{ i } )
		      options = [ options ',' vararg2str({varargin{i}}) ];
		  else
		      %options = [ options ', { [' num2str(varargin{ i }{1}') ']'' EEG.chanlocs EEG.chaninfo }' ]; % JRI -- why does this ignore value passed as topo option?
        if length(varargin{i})>1
		optchanlocs = varargin{i}{2}; % JRI -- fix
        else
		optchanlocs = EEG.chanlocs;
	end;
  	if length(varargin{i})>2,
            optchaninfo = varargin{i}{3};
        else
            optchaninfo = EEG.chaninfo;
        end
          options = [ options ', { [' num2str(varargin{ i }{1}') ']'' optchanlocs optchaninfo }' ];
		  end;
		end
	end;
end
try, icadefs; set(gcf, 'color', BACKCOLOR,'Name',' erpimage()'); catch, end

% testing inputs
% --------------
if ischar(channel)
    channel = eeg_decodechan(EEG.chanlocs, channel);
end
if typeplot == 0 && length(channel) > 1 && isempty(projchan)
	error('A channel must be selected to plot (the sum of) several component projections');
end

% find sorting latencies
% ---------------------
typetxt = '';
if ~isempty(sortingeventfield)
    %events = eeg_getepochevent( EEG, sortingtype, sortingwin, sortingeventfield);
	events = sprintf('eeg_getepochevent( EEG, %s)', vararg2str({sortingtype, sortingwin, sortingeventfield}));

    % generate text for the command
    % -----------------------------
    for index = 1:length(sortingtype)
        if ischar(sortingtype{index})
            typetxt = [typetxt ' ''' sortingtype{index} '''' ];
        else
            typetxt = [typetxt ' ' num2str(sortingtype{index}) ];
        end
    end;
% $$$ 	% renormalize latencies if necessary
% $$$ 	% ----------------------------------
% $$$ 	switch lower(renorm)
% $$$ 	    case 'yes',
% $$$ 	         disp('Pop_erpimage warning: *** sorting variable renormalized ***');
% $$$ 	         events = (events-min(events)) / (max(events) - min(events)) * ...
% $$$ 	                        0.5 * (EEG.xmax*1000 - EEG.xmin*1000) + EEG.xmin*1000 + 0.4*(EEG.xmax*1000 - EEG.xmin*1000);
% $$$ 	    case 'no',;
% $$$ 	    otherwise,
% $$$ 	        locx = findstr('x', lower(renorm))
% $$$ 	        if length(locx) ~= 1, error('Pop_erpimage error: unrecognize renormalazing formula'); end
% $$$ 	        eval( [ 'events =' renorm(1:locx-1) 'events' renorm(locx+1:end) ';'] );
% $$$ 	end
else
	events = 'ones(1, EEG.trials)*EEG.xmax*1000';
    %events = ones(1, EEG.trials)*EEG.xmax*1000;
    sortingeventfield = '';
end;

if ischar(projchan)
    projchan = { projchan };
end
if iscell(projchan)
    projchannum = std_chaninds(EEG, projchan);
else
    projchannum = projchan;
end

if typeplot == 1
	tmpsig = ['mean(EEG.data([' int2str(channel) '], :),1)'];
else
    % test if ICA was computed or if one has to compute on line
    % ---------------------------------------------------------
    tmpsig = [ 'eeg_getdatact(EEG, ''component'', [' int2str(channel) '], ''projchan'', [' int2str(projchannum) '])' ];
end

% outputs
% -------
outstr = '';
if ~popup
    for io = 1:nargout, outstr = [outstr 'varargout{' int2str(io) '},' ]; end
    if ~isempty(outstr), outstr = [ '[' outstr(1:end-1) '] =' ]; end
end

% plot title
% ----------
if isempty(titleplot)
    if typeplot==1 % if channel plot
        if ~isempty(EEG.chanlocs) % if channel information exist
              titleplot = [ EEG.chanlocs(channel).labels ];
        else, titleplot = [ int2str(channel) ];
        end
    else
        titleplot = [ 'Comp. ' int2str(channel) ];
        if ~isempty(projchan),
            tmpstr = vararg2str({projchan});
            tmpstr(find(tmpstr == '''')) = '"';
            titleplot = [ titleplot ' -> Chan. ' tmpstr ];
        end
    end
end

% plot the data and generate output command
% --------------------------------------------
if isempty( options )
    if isfield(opt, 'topo')
        tmptopo = opt.topo;
        opt = rmfield(opt, 'topo');
    else
        tmptopo = '';
    end
    if ~isempty(opt)
        fields = fieldnames(opt);
        values = struct2cell(opt);
        params = { fields{:}; values{:} };
        options = [ ',' vararg2str( { params{:} } ) ];
        tmpind = find( options == '\' ); options(tmpind(1:2:end)) = [];
        if ~isempty(tmptopo), options = [ options ',''topo'',' tmptopo ]; end
    end
end

% varargout{1} = sprintf('figure; pop_erpimage(%s,%d,%d,''%s'',%d,%d,{%s},[%s],''%s'',''%s''%s);', inputname(1), typeplot, channel, titleplot, smooth, decimate, typetxt, int2str(sortingwin), sortingeventfield, renorm, options);
popcom = sprintf('figure; pop_erpimage(EEG,%d, [%s],[%s],''%s'',%d,%d,{%s},[%s],''%s'' %s);', typeplot, int2str(channel), vararg2str({projchan}), titleplot, smooth, decimate, typetxt, int2str(sortingwin), sortingeventfield, options);

com = sprintf('%s erpimage( %s, %s, linspace(EEG.xmin*1000, EEG.xmax*1000, EEG.pnts), ''%s'', %d, %d %s);', outstr, tmpsig, events, titleplot, smooth, decimate, options);
disp('Command executed by pop_erpimage:');
disp(' '); disp(com); disp(' ');
eval(com)

if popup
	varargout{1} = popcom; % [10 '% Call: ' com];
end

return;

% get contextual help
% -------------------
function txt = context(var, allvars, alltext);
	loc = strmatch( var, allvars);
	if ~isempty(loc)
		txt= alltext{loc(1)};
	else
		disp([ 'warning: variable ''' var ''' not found']);
		txt = '';
	end


text:
function [v,unv] = statmoments(p,n)
%STATMOMENTS Computes statistical central moments of image histogram.
%   [W,UNV] = STATMOMENTS(P,N) computes up to the Nth statistical
%   central moment of a histogram whose components are in vector
%   P. The length of P must equal 256 or 65536.
%
%   The program outputs a vector V with V(1) = mean, V(2) = variance,
%   V(3) = 3rd moment, . . . V(N) = Nth central moment. The random
%   variable values are normalized to the range [0,1], so all
%   moments also are in this range.
%
%   The program also outputs a vector UNV containing the same moments
%   as V, but using un-normalized random variable values (e.g., 0 to
%   255 if length(P) = 2^8). For example, if length(P) = 256 and V(1)
%   = 0.5, then UNV(1) would have the value UNV(1) = 127.5 (half of
%   the [0,255] range).
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

Lp = length(p);
if (Lp ~= 256) && (Lp ~= 65536)
   error('P must be a 256- or 65536-element vector.');
end
G = Lp - 1;

% Make sure the histogram has unit area, and convert it to a
% column vector.
p = p/sum(p); p = p(:);

% Form a vector of all the possible values of the
% random variable.
z = 0:G;

% Now normalize the z's to the range [0,1].
z = z./G;

% The mean.
m = z*p;

% Center random variables about the mean.
z = z - m;

% Compute the central moments.
v = zeros(1,n);
v(1) = m;
for j = 2:n
   v(j) = (z.^j)*p;
end

if nargout > 1
   % Compute the uncentralized moments.
   unv = zeros(1,n);
   unv(1)=m.*G;
   for j = 2:n
      unv(j) = ((z*G).^j)*p;
   end
end

author: dipum

repo: dipum-toolbox

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


text:
function [ abd, ipvt, info ] = zgbfa ( abd, lda, n, ml, mu )

%*****************************************************************************80
%
%% ZGBFA factors a complex band matrix by elimination.
%
%  Discussion:
%
%    ZGBFA is usually called by ZGBCO, but it can be called
%    directly with a saving in time if RCOND is not needed.
%
%  Band storage:
%
%    If A is a band matrix, the following program segment
%    will set up the input.
%
%      ml = (band width below the diagonal)
%      mu = (band width above the diagonal)
%      m = ml + mu + 1
%      do j = 1, n
%        i1 = max ( 1, j - mu )
%        i2 = min ( n, j + ml )
%        do i = i1, i2
%          k = i - j + m
%          abd(k,j) = a(i,j)
%        end do
%      end do
%
%    This uses rows ML+1 through 2*ML+MU+1 of ABD.
%    In addition, the first ML rows in ABD are used for
%    elements generated during the triangularization.
%    The total number of rows needed in ABD is 2*ML+MU+1.
%    The ML+MU by ML+MU upper left triangle and the
%    ML by ML lower right triangle are not referenced.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%
%  Parameters:
%
%    Input, complex ABD(LDA,N), the matrix in band storage.  The columns of
%    the matrix are stored in the columns of ABD and the diagonals of the
%    matrix are stored in rows ML+1 through 2*ML+MU+1 of ABD.
%
%    Input, integer LDA, the leading dimension of ABD.
%    LDA must be at least 2*ML+MU+1.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, the number of diagonals below the main diagonal.
%    0 <= ML < N.
%
%    Input, integer MU, the number of diagonals above the main diagonal.
%    0 <= MU < N.  More efficient if ML <= MU.
%
%    Output, complex ABD(LDA,N), an upper triangular matrix
%    in band storage and the multipliers which were used to obtain it.
%    The factorization can be written A = L*U where L is a product of
%    permutation and unit lower triangular matrices and U is upper triangular.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, integer INFO.
%    0, normal value.
%    K, if U(K,K) == 0.0.  This is not an error condition for this
%    subroutine, but it does indicate that ZGBSL will divide by zero if
%    called.  Use RCOND in ZGBCO for a reliable indication of singularity.
%
  m = ml + mu + 1;
  info = 0;
%
%  Zero initial fill-in columns.
%
  j0 = mu + 2;
  j1 = min ( n, m ) - 1;

  for jz = j0 : j1
    i0 = m + 1 - jz;
    for i = i0 : ml
      abd(i,jz) = 0.0;
    end
  end

  jz = j1;
  ju = 0;
%
%  Gaussian elimination with partial pivoting.
%
  for k = 1 : n-1
%
%  Zero next fill-in column
%
    jz = jz + 1;
    if ( jz <= n )
      abd(1:ml,jz) = 0.0;
    end
%
%  Find L = pivot index.
%
    lm = min ( ml, n - k );
    l = izamax ( lm+1, abd(m:m+lm,k), 1 ) + m - 1;
    ipvt(k) = l + k - m;
%
%  Zero pivot implies this column already triangularized.
%
    if ( zabs1 ( abd(l,k) ) == 0.0 )
      info = k;
      continue
    end
%
%  Interchange if necessary.
%
    if ( l ~= m )
      t        = abd(l,k);
      abd(l,k) = abd(m,k);
      abd(m,k) = t;
    end
%
%  Compute multipliers.
%
    t = - 1.0 / abd(m,k);
    abd(m+1:m+lm,k) = abd(m+1:m+lm,k) * t;
%
%  Row elimination with column indexing.
%
    ju = min ( max ( ju, mu + ipvt(k) ), n );
    mm = m;

    for j = k+1 : ju
      l = l - 1;
      mm = mm - 1;
      t = abd(l,j);
      if ( l ~= mm )
        abd(l,j) = abd(mm,j);
        abd(mm,j) = t;
      end
      abd(mm+1:mm+lm,j) = abd(mm+1:mm+lm,j) + t * abd(m+1:m+lm,k);
    end

  end

  ipvt(n) = n;

  if ( zabs1 ( abd(m,n) ) == 0.0 )
    info = n;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function DoScanSeriesUpd(Simuh,ScanFlag)
global VCtl
    try
    [pathstr,name,ext]=fileparts(Simuh.SeqXMLFile);
    catch me
        guidata(Simuh.SimuPanel_figure,Simuh);
        return;
    end
    if ~isfield(Simuh,'ScanSeries')
        Simuh.ScanSeriesInd=1;
        Simuh.ScanSeries(1,1:3)={[num2str(Simuh.ScanSeriesInd) ':'],'Dx',name};
    else
        switch ScanFlag
            case 0 % load PSD
               Simuh.ScanSeries(Simuh.ScanSeriesInd,1:3)={[num2str(Simuh.ScanSeriesInd) ':'],'Dx',name};
            case 1 % scan setting
               Simuh.ScanSeries(Simuh.ScanSeriesInd,1:2)={[num2str(Simuh.ScanSeriesInd) ':'],'Dx'};
            case 2 % scanning
               Simuh.ScanSeries(Simuh.ScanSeriesInd,1:2)={[num2str(Simuh.ScanSeriesInd) ':'],'...'};
               VCtl.SeriesName = Simuh.ScanSeries{Simuh.ScanSeriesInd,3};
            case 3 % scan complete
               Simuh.ScanSeries(Simuh.ScanSeriesInd,1:2)={[num2str(Simuh.ScanSeriesInd) ':'],'V'};
               Simuh.ScanSeriesInd=Simuh.ScanSeriesInd+1;
            case 4 % scan fail
               Simuh.ScanSeries(Simuh.ScanSeriesInd,1:2)={[num2str(Simuh.ScanSeriesInd) ':'],'X'};
               Simuh.ScanSeriesInd=Simuh.ScanSeriesInd+1;
            case 5 % add to batch
               Simuh.ScanSeries(Simuh.ScanSeriesInd,1:2)={[num2str(Simuh.ScanSeriesInd) ':'],'B'};
               Simuh.ScanSeriesInd=Simuh.ScanSeriesInd+1;
        end

    end
    set(Simuh.ScanSeries_uitable,'Data',Simuh.ScanSeries);
    set(Simuh.ScanSeries_uitable,'Enable','on');
    guidata(Simuh.SimuPanel_figure,Simuh);
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
clear all
syms a b c x
        % Rezolvarea primei ecuatii
sol1=solve(a*x^2+b*x+c)
pretty(sol1)
        % Rezolvarea ecuatiei a doua
sol2=solve('cos(2*x)+sin(x)=1')
        % Transformarea solutiei simbolice in numerica
numsol2=double(sol2)
        % Reprezentarea grafica a functiei pe intervalul dat
ezplot('cos(2*x)+sin(x)-1',[0, 2*pi])
hold on
grid
        % Marcarea solutiilor obtinute pe grafic
plot(numsol2,zeros(size(numsol2)),'rd')

author: Sable

repo: mcbench-benchmarks


text:
function asa047_test01 ( )

%*****************************************************************************80
%
%% ASA047_TEST01 demonstrates the use of NELMIN on ROSENBROCK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2008
%
%  Author:
%
%    John Burkardt
%
  n = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA047_TEST01\n' );
  fprintf ( 1, '  Apply NELMIN to ROSENBROCK function.\n' );

  start(1:n) = [ -1.2; 1.0 ];
  reqmin = 1.0E-08;
  step(1:n) = 1.0;
  konvge = 10;
  kcount = 500;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Starting point X:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %f\n', start(i) );
  end

  ynewlo = rosenbrock ( start );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X) = %f\n', ynewlo );

  [ xmin, ynewlo, icount, numres, ifault ] = nelmin ( @rosenbrock, n, start, ...
  reqmin, step, konvge, kcount );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Return code IFAULT = %d\n', ifault );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate of minimizing value X*:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %f\n', xmin(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %f\n', ynewlo );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', icount );
  fprintf ( 1, '  Number of restarts =   %d\n', numres );

  return
end


text:
function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
%GRADIENTDESCENT Performs gradient descent to learn theta
%   theta = GRADIENTDESCENT(X, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCost) and gradient here.
    %

    predictions =  X * theta;
    updates = X' * (predictions - y);
    theta = theta - alpha * (1/m) * updates;
    %theta = theta - alpha * (1/m) * sum(sqerrors) * X;
    %theta - (alpha/m) * (X' * (X * theta - y));
    %theta = theta - (alpha/m) * (X' * (X * theta - y));

    % ============================================================

    % Save the cost J in every iteration
    J_history(iter) = computeCost(X, y, theta);

end

end

repo: machine-learning


text:
function [ fea, out ] = ex_axistressstrain1( varargin )
%EX_AXISTRESSSTRAIN1 Example for hollow cylider axisymmetric stress-strain.
%
%   [ FEA, OUT ] = EX_AXISTRESSSTRAIN1( VARARGIN ) Example to calculate displacements and stresses
%   in a hollow cylinder in axisymmetric/cylindrical coordinates.
%
%   Ref. 4.1.9 Long (generalized plane strain) cylinder subjected to internal and external pressure.
%   [1] Applied Mechanics of Solids, Allan F. Bower, 2012 (http://solidmechanics.org/).
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       a           scalar {1.5}           Cylinder inner radius
%       b           scalar {2}             Cylinder outer radius
%       l           scalar {3}             Cylinder length
%       pa          scalar {5e3}           Inner load force
%       pb          scalar {20e4}          Outer load force
%       E           scalar {200e9}         Modulus of elasticity
%       nu          scalar {0.3}           Poissons ratio
%       igrid       scalar 0/{1}           Cell type (0=quadrilaterals, 1=triangles)
%       hmax        scalar {0.1}           Max grid cell size
%       sfun        string {sflag2}        Shape function for displacements
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'a',        0.9;
            'b',        2;
            'l',        3;
            'pa',       5e3;
            'pb',       20e4;
            'E',        200e9;
            'nu',       0.3;
            'igrid',    0;
            'hmax',     0.01;
            'sfun',     'sflag2';
            'iplot',    1;
            'tol',      5e-2;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry and grid.
a = opt.a;
b = opt.b;
l = opt.l;
fea.geom.objects = { gobj_rectangle( a, b, 0, l, 'R1' ) };
if ( opt.igrid==1 )
  fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );
else
  fea.grid = rectgrid( ceil((b-a)/opt.hmax), ceil(l/opt.hmax), [a b;0 l] );
  if( opt.igrid<0 )
    fea.grid = quad2tri( fea.grid );
  end
end
n_bdr = max(fea.grid.b(3,:));   % Number of boundaries.

% Axisymmetric stress-strain equation definitions.
fea.sdim = { 'r', 'z' };
fea = addphys( fea, @axistressstrain );
fea.phys.css.eqn.coef{1,end} = { opt.nu };
fea.phys.css.eqn.coef{2,end} = { opt.E  };
fea.phys.css.sfun            = { opt.sfun opt.sfun };   % Set shape functions.

% Boundary conditions.
bctype = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
[bctype{2,:}] = deal( 1 );
fea.phys.css.bdr.coef{1,5} = bctype;

bccoef = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bccoef{1,2} = opt.pb*b;
bccoef{1,4} = opt.pa*a;
fea.phys.css.bdr.coef{1,end} = bccoef;

% Solve.
fea       = parsephys( fea );
fea       = parseprob( fea );
fea.sol.u = solvestat( fea, 'icub', 1+str2num(strrep(opt.sfun,'sflag','')), 'fid', fid );

% Postprocessing.
n = 20;
r = linspace(a,b,n);
z = l/2*ones(1,n);
pa = opt.pa; pb = -opt.pb; E = opt.E; nu = opt.nu;
% 4.1.9 http://solidmechanics.org/Text/Chapter4_1/Chapter4_1.php#Sect4_1_9
u_ref  = (1+nu)*a^2*b^2/E/(b^2-a^2)*( (pa-pb)./r' + (1-2*nu)*(pa*a^2 - pb*b^2)/a^2/b^2*r' );
sr_ref = (pa*a^2-pb*b^2)/(b^2-a^2) - a^2*b^2/(b^2-a^2)./(r').^2*(pa-pb);
st_ref = (pa*a^2-pb*b^2)/(b^2-a^2) + a^2*b^2/(b^2-a^2)./(r').^2*(pa-pb);
sz_ref = 2*nu*(pa*a^2-pb*b^2)/(b^2-a^2);
u  = evalexpr( fea.phys.css.eqn.vars{3,2}, [r;z], fea );
w  = evalexpr( fea.phys.css.eqn.vars{4,2}, [r;z], fea );
sr = evalexpr( fea.phys.css.eqn.vars{5,2}, [r;z], fea );
st = evalexpr( fea.phys.css.eqn.vars{6,2}, [r;z], fea );
sz = evalexpr( fea.phys.css.eqn.vars{7,2}, [r;z], fea );
if( opt.iplot>0 )
  subplot(1,2,1)
  postplot( fea, 'surfexpr', fea.phys.css.eqn.vars{3,2} )
  title('r-displacement')
  subplot(1,2,2), hold on
  plot(u_ref,r,'r-')
  plot(u,r,'b.')
  legend('exact solution','computed solution')
  xlabel('r')
  grid on
end

% Error checking.
out.erru  = norm( u_ref - u )/norm( u_ref );
out.errw  = norm( w );
out.errsr = norm( sr_ref - sr )/norm( sr_ref );
out.errst = norm( st_ref - st )/norm( st_ref );
out.errsz = norm( sz_ref - sz )/norm( sz_ref );
out.err   = [ out.erru, out.errw, out.errsr, out.errst, out.errsz ];
out.pass  = all(out.err < opt.tol);

if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function type = getCERRStructType(structS);
%"getCERRStructType"
%   Checks the fields in the passed structS and determines which cell in a
%   planC this structure was pulled from.  This determination is based on the
%   initializeCERR function's field names, so planCs with structs that do
%   not match initializeCERR's prototypes will throw a warning and return
%   [].
%
%   The output is in string form, ie 'dose', 'structures', etc.
%
%JRA 06/27/06
%
%Usage:
%   type = getCERRStructType(structS);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

planInitC = initializeCERR;

%Get possible cell field names in a planC.
indexFields = fields(planInitC{end});

%Determine the number of cell types initializeCERR returns.
numCells = length(planInitC);

%Extract fieldnames of the passed structS
myFields = fields(structS);

%Loop over each initialized datatype and compare field names, looking for a
%match to the passed structS.
for i=1:numCells

    thisCellsFields = fields(planInitC{i});

    fieldsMatch = ismember(thisCellsFields, myFields);

    nMatchingFields(i) = sum(fieldsMatch);

end

%Find the cell that matched the most fields.
[maxNum, ind] = max(nMatchingFields);

%If more than 1 match in that top cell, assume we know the type now.
if maxNum > 1
    type = indexFields{ind};
    return;
else
   %No fields matched at all!
   warning('Structure passed to getCERRStructureType does not appear to be a valid planC structure.  Returning [].');
   type = [];
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [x_mu, y_mu sig] = localizew(w,peval,localMaxRadius)
% [x_mu, y_mu sig] = localizew(w,peval,localMaxRadius)
% Localilzes w (#pixels x #components) from NMF model (V=WH)
% peval : parameters (needed peval.nx, peval.ny)
% localMaxDiameter : diameter in pixels to which confine a local maximum search. If set to 0 (default) no confinemend is done.
if ~exist('localMaxRadius', 'var')
    localMaxRadius = 0; % Finding global max.
end
K=size(w,2);
wr=reshape(w,peval.nx, peval.ny,K); %K can be different from peval.ncomp if background is takes as one component...
[x,y] = meshgrid(1:peval.ny, 1:peval.nx); % nx and ny has to be swapped as in the image the first coordiante is number of rows
if localMaxRadius
    fprintf('Looking for local maximum %g only.\n', localMaxRadius);
end

x_mu=inf(1,K);
y_mu=inf(1,K);
sig=inf(1,K);
differ=inf(1,K);
for ii=1:K %background is not localized...
    [x_mu(ii), y_mu(ii), sig(ii), differ(ii)] = fitgauss2d(wr(:,:,ii));
    if localMaxRadius>0
        confinementMaskTmp=(x-x_mu(ii)).^2+(y-y_mu(ii)).^2<=localMaxRadius^2;
        [ymtmp,xmtmp]=find(wr(:,:,ii)==max(max(wr(:,:,ii).*confinementMaskTmp)),1);
        if isempty(xmtmp)&&isempty(ymtmp)
            xmtmp=10*peval.nx;
            ymtmp=xmtmp;
            fprintf('Problem with maximum localisation.\n');
        end
        confinementMask=(x-xmtmp).^2+(y-ymtmp).^2<=localMaxRadius^2; % mask centered on the local maximum
        [x_mu(ii), y_mu(ii), sig(ii), differ(ii)] = fitgauss2d(wr(:,:,ii).*confinementMask);
    end

end

if sum(isinf([x_mu, y_mu,sig,differ]))
    warning('Problem with localisation...')
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function test_suite = test_cg_ls
    clear all;
    initTestSuite;
end

function test_1
    A = [
    1 -1
    1 1
    2 1
    ];

    B = [
    2
    4
    8
    ];
    solver = SPX_CGLeastSquare(A, B);
    x = solver.solve();
    %solver.printResults();
    verifyTrue(testCase, solver.hasConverged());
end

function test_2
    m = 10;
    n = 4;
    s = 100;
    %rng('default');
    A = randn(m, n);
    X = randn(n, s);
    B = A * X;
    solver = SPX_CGLeastSquare(A, B);
    XX = solver.solve();
    % solver.printResults();
    % disp(X);
    % disp(XX);
    verifyTrue(testCase, solver.hasConverged());
end

function test_3
    m = 1000;
    n = 40;
    s = 10;
    %rng('default');
    A = randn(m, n);
    X = randn(n, s);
    B = A * X;
    solver = SPX_CGLeastSquare(A, B);
    XX = solver.solve();
    % solver.printResults();
    % disp(X);
    % disp(XX);
    verifyTrue(testCase, solver.hasConverged());
end

author: indigits

repo: sparse-plex


text:
classdef DTLZ1 < PROBLEM
% <multi/many> <real> <large/none> <expensive/none>
% Benchmark MOP proposed by Deb, Thiele, Laumanns, and Zitzler

%------------------------------- Reference --------------------------------
% K. Deb, L. Thiele, M. Laumanns, and E. Zitzler, Scalable test problems
% for evolutionary multiobjective optimization, Evolutionary multiobjective
% Optimization. Theoretical Advances and Applications, 2005, 105-145.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M+4; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            g      = 100*(obj.D-obj.M+1+sum((PopDec(:,obj.M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,obj.M:end)-0.5)),2));
            PopObj = 0.5*repmat(1+g,1,obj.M).*fliplr(cumprod([ones(size(PopDec,1),1),PopDec(:,1:obj.M-1)],2)).*[ones(size(PopDec,1),1),1-PopDec(:,obj.M-1:-1:1)];
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M)/2;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,1,10)';
                R = {a*a'/2,a*(1-a')/2,(1-a)*ones(size(a'))/2};
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [sg] = g_convSig_approx(x,P,u,in)
% this function evaluates a logistic convolution of the inputs
% function [gx,dgdx,dgdp] = g_convSig(x,P,u,in)
% This function evaluates the following multiple convolution model:
%   p(y(t)=1) = sigm(cst + sum_k sum_tau w(k,tau)*u(k,t-tau))
% where u can be of any dimension (k=1,...,K), and sigm is the standard
% sigmoid mapping. This logistic convolution model is useful for binary
% observations.
% Note that the convolution (Volterra) kernels are constructed from the
% parameters' vector P. The size of the parameters' vector P is determined
% by the maximum lag (tau) considered.
% IN:
%   - x: [useless]
%   - P: (K*maxLag+1)x1 vector of kernel parameters.
%   - u: (K*nt)x1 vectorized input to the system
%   - in: [useless]
% OUT:
%   - sg: the predicted system's output (in terms of p(y(t)=1))
%   - dsdx: [useless]
%   - dsdp: the gradient of the system's output w.r.t. kernel parameters
% SEE ALSO: g_conv0

% [g,~,dgdp] = g_conv_approx(x,P,u,in);
g = g_conv_approx(x,P,u,in);
sg = VBA_sigmoid(g);

% dsdx = [];
% dsdp = dgdp*diag(sg.*(1-sg));

author: MBB-team

repo: VBA-toolbox


text:
%% Demo to show the results of the Ultrametric Contour Map obtained by MCG
clear all;close all;home;

% Read an input image
I = imread(fullfile(root_dir, 'demos','101087.jpg'));

tic;
% Test the 'fast' version, which takes around 3 seconds in mean
ucm2_scg = im2ucm(I,'fast');
toc;

tic;
% Test the 'accurate' version, which tackes around 25 seconds in mean
ucm2_mcg = im2ucm(I,'accurate');
toc;

%% Show UCM results (dilated for visualization)
figure;
subplot(1,3,1)
imshow(I), title('Image')

subplot(1,3,2)
imshow(imdilate(ucm2_scg,strel(ones(3))),[]), title('Fast UCM (SCG)')

subplot(1,3,3)
imshow(imdilate(ucm2_mcg,strel(ones(3))),[]), title('Accurate UCM (MCG)')

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function value = box_contains_segment_nd ( dim_num, p1, p2, pa, pb  )

%*****************************************************************************80
%
%% BOX_CONTAINS_SEGMENT_ND reports if a box contains a line segment in ND.
%
%  Discussion:
%
%    A box is assumed to be a rectangle with sides aligned on coordinate
%    axes.  It can be described by its low and high corner, P1 and P2:
%
%      points P so that P1(1:DIM_NUM) <= P(1:DIM_NUM) <= P2(1:DIM_NUM).
%
%    A line segment is the finite portion of a line that lies between
%    two points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real P1(DIM_NUM), P2(DIM_NUM), the low and high corners of the box.
%
%    Input, real PA(DIM_NUM), PB(DIM_NUM), the endpoints of the line segment.
%
%    Output, logical VALUE, is TRUE if the box contains
%    the line segment.
%
  value = 0;

  if ( ~box_contains_point_nd ( dim_num, p1, p2, pa ) )
    return
  end

  if ( ~box_contains_point_nd ( dim_num, p1, p2, pb ) )
    return
  end

  value = 1;

  return
end


text:
function SQL_CloseDatabase(dbc)
% SQL_CloseDatabase     Closes the connection to database, dbc.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

switch class(dbc)
case 'com.mysql.jdbc.JDBC4Connection'
    % Database connection opened using java commands:
    dbc.close();
case 'database'
    % Database connection opened using database toolbox:
    close(dbc);
otherwise
    error('Unknown database connection class %s',class(dbc));
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function y = getBatch(imdb, images, varargin)
% GET_BATCH  Load, preprocess, and pack images for CNN evaluation

opts.imageSize = [512, 512] - 128 ;
opts.numAugments = 1 ;
opts.transformation = 'none' ;
opts.rgbMean = [] ;
opts.rgbVariance = zeros(0,3,'single') ;
opts.labelStride = 1 ;
opts.labelOffset = 0 ;
opts.classWeights = ones(1,21,'single') ;
opts.interpolation = 'bilinear' ;
opts.numThreads = 1 ;
opts.prefetch = false ;
opts.useGpu = false ;
opts = vl_argparse(opts, varargin);

if opts.prefetch
  % to be implemented
  ims = [] ;
  labels = [] ;
  return ;
end

if ~isempty(opts.rgbVariance) && isempty(opts.rgbMean)
  opts.rgbMean = single([128;128;128]) ;
end
if ~isempty(opts.rgbMean)
  opts.rgbMean = reshape(opts.rgbMean, [1 1 3]) ;
end

% space for images
ims = zeros(opts.imageSize(1), opts.imageSize(2), 3, ...
  numel(images)*opts.numAugments, 'single') ;

% space for labels
lx = opts.labelOffset : opts.labelStride : opts.imageSize(2) ;
ly = opts.labelOffset : opts.labelStride : opts.imageSize(1) ;
labels = zeros(numel(ly), numel(lx), 1, numel(images)*opts.numAugments, 'single') ;
classWeights = [0 opts.classWeights(:)'] ;

im = cell(1,numel(images)) ;

si = 1 ;

for i=1:numel(images)

  % acquire image
  if isempty(im{i})
    rgbPath = sprintf(imdb.paths.image, imdb.images.name{images(i)}) ;
    labelsPath = sprintf(imdb.paths.classSegmentation, imdb.images.name{images(i)}) ;
    rgb = vl_imreadjpeg({rgbPath}) ;
    rgb = rgb{1} ;
    anno = imread(labelsPath) ;
  else
    rgb = im{i} ;
  end
  if size(rgb,3) == 1
    rgb = cat(3, rgb, rgb, rgb) ;
  end

  % crop & flip
  h = size(rgb,1) ;
  w = size(rgb,2) ;
  for ai = 1:opts.numAugments
    sz = opts.imageSize(1:2) ;
    scale = max(h/sz(1), w/sz(2)) ;
    scale = scale .* (1 + (rand(1)-.5)/5) ;

    sy = round(scale * ((1:sz(1)) - sz(1)/2) + h/2) ;
    sx = round(scale * ((1:sz(2)) - sz(2)/2) + w/2) ;
    if rand > 0.5, sx = fliplr(sx) ; end

    okx = find(1 <= sx & sx <= w) ;
    oky = find(1 <= sy & sy <= h) ;
    if ~isempty(opts.rgbMean)
      ims(oky,okx,:,si) = bsxfun(@minus, rgb(sy(oky),sx(okx),:), opts.rgbMean) ;
    else
      ims(oky,okx,:,si) = rgb(sy(oky),sx(okx),:) ;
    end

    tlabels = zeros(sz(1), sz(2), 'uint8') + 255 ;
    tlabels(oky,okx) = anno(sy(oky),sx(okx)) ;
    tlabels = single(tlabels(ly,lx)) ;
    tlabels = mod(tlabels + 1, 256) ; % 0 = ignore, 1 = bkg
    labels(:,:,1,si) = tlabels ;
    si = si + 1 ;
  end
end
if opts.useGpu
  ims = gpuArray(ims) ;
end
y = {'input', ims, 'label', labels} ;

author: vlfeat

repo: matconvnet-fcn


text:
function ui = mrViewInit(mr);
%
% ui = mrViewInit(mr);
%
% Initialize a ui struct, containing info needed for
% mrViewer.
%
% I have attempted to structure this such that analyses
% can be run outside the GUI, akin to hidden views in mrVista 1.0.
% In this case, running mrViewInit, without running mrViewer,
% would open the structure without opening a GUI. However, I
% haven't had a chance to test that everything would work
% without needing a handle -- please feel free to correct
% any dependence on a GUI.
%
% ras, 07/05.
if ~exist('mr','var') || isempty(mr),   mr = mrLoad;         end

% Key, important fields of the view:
ui.mr = mr;         % 'base' mr object, like the anatomy / underlay
ui.maps = [];       % overlays, or maps to use when thresholding
ui.tSeries = {};    % time series (paths / mr structs) files, cell
ui.stim = [];       % stimulus files
ui.session = [];    % optional session struct / handle to session viewer
                    % to which the data are attached (see sessionViewer)
ui.rois = [];       % ROIs defined on the viewer (but not necessarily
                    % saved yet
ui.fig = [];        % handle to the main figure for the GUI

% the GUI will comprise many uipanels (specifically, see mrvPanel for
% UI panels which can be toggled visible/hidden); the most important
% is the panel where the images of the data will be displayed:
ui.panels.display = [];  % display panel
ui.panels.nav = []; % navigation panel
ui.panels.roi = []; % ROI (region of interest) panel
ui.panels.grayscale = []; % grayscale (underlay image) panel
ui.panels.info = [];  % MR info panel
ui.panels.mesh = [];  % Mesh options panel
ui.panels.overlays = [];  % panels containing overlay info
ui.panels.colorbar = []; % panel containing colorbars for overlays

%%%%%%%%%%%%%
% settings: %
%%%%%%%%%%%%%
% (maybe we should sub-organize this, but I'm ok with having
% a long list...)
ui.settings.space = 1; % pixel space
ui.settings.segmentation = 0;  % selected segmentation
ui.settings.displayFormat = 3; % default is 3 -- single slice
ui.settings.ori = 3; % orientation / default is rows | columns
middleSlice = round(ui.mr.dims(3) / 2);
ui.settings.slice = middleSlice; % default slice for single / montage views
ui.settings.time = 1; % time point, for 4-D data
ui.settings.roi = 0; % index into current ROI (0 = none selected)
ui.settings.montageRows = 2; % # rows for montage view
ui.settings.montageCols = 2; % # cols for montage view
ui.settings.bounds = [1 1 1; ui.mr.dims([2 1 3])]';
ui.settings.cursorLoc = round(ui.mr.dims(1:3) ./ 2);
ui.settings.zoom = [1 size(mr.data,1); ...
                    1 size(mr.data,2); ...
                    1 size(mr.data,3);];
ui.settings.clim = []; % auto by default
ui.settings.cmap = gray(256); % underlay cmap
ui.settings.cbarColorScheme = 1;  % 1 = white-on-black, 2=black-on-white
ui.settings.brightness = 0; % how much to brighent the gray cmap
ui.settings.gamma = 2.0; % not implemented yet, but may be useful
ui.settings.labelAxes = 0; % show axis ticks and label axes
ui.settings.labelDirs = 0; % show direction text in display
ui.settings.labelSlices = 0; % label each slice in display
ui.settings.roiEditMode = 7; % mode for adding to / removing from ROIs:
                             % 1 -- rect
                             % 2 -- circle
                             % 3 -- cube
                             % 4 -- sphere
                             % 5 -- line
                             % 6 -- point
                             % 7 -- grow (3D)
ui.settings.roiViewMode = 3; % mode for showing / hiding ROIs
                             % 1 -- hide ROIs
                             % 2 -- show selected ROI
                             % 3 -- show all ROIs
ui.settings.eqAspect = 1; % preserve aspect ratio
ui.settings.showCursor = 0; % display cursor flag
ui.settings.cursorType = 4; % flag for cursor rendering style:
                            % 1 -- show + sign
                            % 2 -- show circle
                            % 3 -- show crosshairs
                            % 4 -- show crosshairs + gap
ui.settings.baseInterp = 'linear'; % Interpolation method for base
ui.settings.mapInterp = 'nearest'; % Interpolation method for maps

% ras 02/2007:
% instead of using 'auto' color limits for the object by default,
% I've found that it's easier to guess a good upper clip limit, but
% leave the lower limit at the min value.

% [ignore clipMin clipMax] = mrClipOptimal(ui.mr.data(:,:,:,1));
% ui.settings.clim = [ui.mr.dataRange(1) clipMax];
[tmp minVal maxVal] = histoThresh(ui.mr.data(:,:,:,1));
ui.settings.clim = [minVal maxVal];

% if the mr object itself has saved settings, load these into the UI:
if isfield(ui.mr,'settings')
    for f = fieldnames(ui.mr.settings)'
        ui.settings.(f{1}) = ui.mr.settings.(f{1});
    end
end

% form a unique object tag for this viewer
i = 1; ui.tag = sprintf('mrViewer%i',i);
while ~isempty(findobj('Tag',ui.tag))
    i = i + 1; ui.tag = sprintf('mrViewer%i',i);
end

return