text:
function [cfg, artifact] = ft_artifact_zvalue(cfg, data)

% FT_ARTIFACT_ZVALUE scans data segments of interest for artifacts, by means of
% thresholding the z-scored values of signals that have been preprocessed,
% using heuristics that increase the sensitivity to detect certain types of artifacts.
% Depending on the preprocessing options, this method will be sensitive to EOG, muscle
% or SQUID jump artifacts. The z-scoring is applied in order to make the threshold
% independent of the phsyical units in the data.
%
% Use as
%   [cfg, artifact] = ft_artifact_zvalue(cfg)
% with the configuration options
%   cfg.trl        = structure that defines the data segments of interest, see FT_DEFINETRIAL
%   cfg.continuous = 'yes' or 'no' whether the file contains continuous data.
%                    If the data has not been recorded continuously, then the cfg.trl should
%                    stricly observe the boundaries of the discontinuous segments, and the
%                    permitted values padding options (described below) are restricted to 0.
%   cfg.dataset    = string with the filename
% or
%   cfg.headerfile = string with the filename
%   cfg.datafile   = string with the filename
% and optionally
%   cfg.headerformat
%   cfg.dataformat
%
% Alternatively you can use it as
%   [cfg, artifact] = ft_artifact_zvalue(cfg, data)
% where the input data is a structure as obtained from FT_PREPROCESSING. Any preprocessing options
% defined in the cfg will be applied to the data before the z-scoring and thresholding.
%
% In both cases the configuration should also contain
%   cfg.trl        = structure that defines the data segments of interest, see FT_DEFINETRIAL
%   cfg.continuous = 'yes' or 'no' whether the file contains continuous data
% and
%   cfg.artfctdef.zvalue.channel    = Nx1 cell-array with selection of channels, see FT_CHANNELSELECTION for details
%   cfg.artfctdef.zvalue.cutoff     = number, z-value threshold
%   cfg.artfctdef.zvalue.trlpadding = number in seconds
%   cfg.artfctdef.zvalue.fltpadding = number in seconds
%   cfg.artfctdef.zvalue.artpadding = number in seconds
%
% If you encounter difficulties with memory usage, you can use
%   cfg.memory = 'low' or 'high', whether to be memory or computationally efficient, respectively (default = 'high')
%
% The optional configuration settings (see below) are:
%   cfg.artfctdef.zvalue.artfctpeak       = 'yes' or 'no'
%   cfg.artfctdef.zvalue.artfctpeakrange  = [begin end]
%   cfg.artfctdef.zvalue.interactive      = 'yes' or 'no'
%   cfg.artfctdef.zvalue.zscore           = 'yes' (default) or 'no'
%
% If you specify cfg.artfctdef.zvalue.artfctpeak='yes', a peak detection on the suprathreshold
% z-scores will be performed, and the artifact will be defined relative to
% the peak, where the begin and end points will be defined by
% cfg.artfctdef.zvalue artfctpeakrange, rather than by the time points that
% exceed the threshold.
%
% You can specify cfg.artfctdef.zvalue.artfctpeakrange if you want to use the
% detected artifacts as input to the DSS method of FT_COMPONENTANALYSIS. The result
% is saved into cfg.artfctdef.zvalue.artifact. The range will automatically
% respect the trial boundaries, i.e. it will be shorter if peak is near the beginning
% or end of a trial. Samples between trials will be removed, thus this will not match
% the sampleinfo of the data structure.
%
% If you specify cfg.artfctdef.zvalue.zscore = 'no', the data will NOT be z-scored prior
% to thresholding. This goes a bit against the name of the function, but it may be useful
% if the threshold is to be defined in meaningful physical units, e.g. degrees of visual
% angle for eye position data.
%
% If you specify cfg.artfctdef.zvalue.interactive = 'yes', a graphical user interface
% will show in which you can manually accept/reject the detected artifacts, and/or
% change the threshold. To control the graphical interface via keyboard, use the
% following keys:
%
%     q                 : Stop
%
%     comma             : Step to the previous artifact trial
%     a                 : Specify artifact trial to display
%     period            : Step to the next artifact trial
%
%     x                 : Step 10 trials back
%     leftarrow         : Step to the previous trial
%     t                 : Specify trial to display
%     rightarrow        : Step to the next trial
%     c                 : Step 10 trials forward
%
%     k                 : Keep trial
%     space             : Mark complete trial as artifact
%     r                 : Mark part of trial as artifact
%
%     downarrow         : Shift the z-threshold down
%     z                 : Specify the z-threshold
%     uparrow           : Shift the z-threshold down
%
% Configuration settings related to the preprocessing of the data are
%   cfg.artfctdef.zvalue.lpfilter      = 'no' or 'yes'  lowpass filter
%   cfg.artfctdef.zvalue.hpfilter      = 'no' or 'yes'  highpass filter
%   cfg.artfctdef.zvalue.bpfilter      = 'no' or 'yes'  bandpass filter
%   cfg.artfctdef.zvalue.bsfilter      = 'no' or 'yes'  bandstop filter for line noise removal
%   cfg.artfctdef.zvalue.dftfilter     = 'no' or 'yes'  line noise removal using discrete fourier transform
%   cfg.artfctdef.zvalue.medianfilter  = 'no' or 'yes'  jump preserving median filter
%   cfg.artfctdef.zvalue.lpfreq        = lowpass  frequency in Hz
%   cfg.artfctdef.zvalue.hpfreq        = highpass frequency in Hz
%   cfg.artfctdef.zvalue.bpfreq        = bandpass frequency range, specified as [low high] in Hz
%   cfg.artfctdef.zvalue.bsfreq        = bandstop frequency range, specified as [low high] in Hz
%   cfg.artfctdef.zvalue.lpfiltord     = lowpass  filter order
%   cfg.artfctdef.zvalue.hpfiltord     = highpass filter order
%   cfg.artfctdef.zvalue.bpfiltord     = bandpass filter order
%   cfg.artfctdef.zvalue.bsfiltord     = bandstop filter order
%   cfg.artfctdef.zvalue.medianfiltord = length of median filter
%   cfg.artfctdef.zvalue.lpfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.artfctdef.zvalue.hpfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.artfctdef.zvalue.bpfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.artfctdef.zvalue.bsfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.artfctdef.zvalue.detrend       = 'no' or 'yes'
%   cfg.artfctdef.zvalue.demean        = 'no' or 'yes'
%   cfg.artfctdef.zvalue.baselinewindow = [begin end] in seconds, the default is the complete trial
%   cfg.artfctdef.zvalue.hilbert       = 'no' or 'yes'
%   cfg.artfctdef.zvalue.rectify       = 'no' or 'yes'
%
% The output argument "artifact" is a Nx2 matrix comparable to the "trl" matrix of
% FT_DEFINETRIAL. The first column of which specifying the beginsamples of an
% artifact period, the second column contains the endsamples of the artifactperiods.
%
% To facilitate data-handling and distributed computing, you can use
%   cfg.inputfile   =  ...
% to read the input data from a *.mat file on disk. This mat files should contain
% only a single variable named 'data', corresponding to the input structure.
%
% See also FT_REJECTARTIFACT, FT_ARTIFACT_CLIP, FT_ARTIFACT_ECG, FT_ARTIFACT_EOG,
% FT_ARTIFACT_JUMP, FT_ARTIFACT_MUSCLE, FT_ARTIFACT_THRESHOLD, FT_ARTIFACT_ZVALUE

% Copyright (C) 2003-2011, Jan-Mathijs Schoffelen & Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance
ft_preamble loadvar data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% for backward compatibility
cfg = ft_checkconfig(cfg, 'renamed', {'artfctdef.blc',             'artfctdef.demean'});
cfg = ft_checkconfig(cfg, 'renamed', {'artfctdef.blcwindow'        'artfctdef.baselinewindow'});
cfg = ft_checkconfig(cfg, 'renamed', {'artfctdef.zvalue.sgn',      'artfctdef.zvalue.channel'});
cfg = ft_checkconfig(cfg, 'renamed', {'artfctdef.zvalue.feedback', 'artfctdef.zvalue.interactive'});
cfg = ft_checkconfig(cfg, 'forbidden',  {'padding'});

% set the default options
cfg.continuous      = ft_getopt(cfg, 'continuous',     []);
cfg.feedback        = ft_getopt(cfg, 'feedback',       'text');
cfg.memory          = ft_getopt(cfg, 'memory',         'high');
cfg.representation  = ft_getopt(cfg, 'representation', 'numeric'); % numeric or table

% set default rejection parameters
cfg.artfctdef                    = ft_getopt(cfg,                  'artfctdef',    []);
cfg.artfctdef.zvalue             = ft_getopt(cfg.artfctdef,        'zvalue',       []);
cfg.artfctdef.zvalue.method      = ft_getopt(cfg.artfctdef.zvalue, 'method',       'all');
cfg.artfctdef.zvalue.ntrial      = ft_getopt(cfg.artfctdef.zvalue, 'ntrial',       10);
cfg.artfctdef.zvalue.channel     = ft_getopt(cfg.artfctdef.zvalue, 'channel',      {});
cfg.artfctdef.zvalue.trlpadding  = ft_getopt(cfg.artfctdef.zvalue, 'trlpadding',   0);
cfg.artfctdef.zvalue.fltpadding  = ft_getopt(cfg.artfctdef.zvalue, 'fltpadding',   0);
cfg.artfctdef.zvalue.artpadding  = ft_getopt(cfg.artfctdef.zvalue, 'artpadding',   0);
cfg.artfctdef.zvalue.interactive = ft_getopt(cfg.artfctdef.zvalue, 'interactive',  'no');
cfg.artfctdef.zvalue.cumulative  = ft_getopt(cfg.artfctdef.zvalue, 'cumulative',   'yes');
cfg.artfctdef.zvalue.artfctpeak  = ft_getopt(cfg.artfctdef.zvalue, 'artfctpeak',   'no');
cfg.artfctdef.zvalue.artfctpeakrange  = ft_getopt(cfg.artfctdef.zvalue, 'artfctpeakrange',[0 0]);
cfg.artfctdef.zvalue.zscore      = ft_getopt(cfg.artfctdef.zvalue, 'zscore',       'yes');

if isfield(cfg.artfctdef.zvalue, 'artifact')
  ft_notice('zvalue artifact detection has already been done, retaining artifacts\n');
  artifact = cfg.artfctdef.zvalue.artifact;
  return
end

% clear old warnings from this stack
ft_warning('-clear')

% flag whether to compute z-value per trial or not, rationale being that if there are
% fluctuations in the variance across trials (e.g. due to position differences in MEG
% measurements) which don't have to do with the artifact per se, the detection is
% compromised (although the data quality is questionable when there is a lot of
% movement to begin with).
pertrial    = strcmp(cfg.artfctdef.zvalue.method, 'trial');
demeantrial = strcmp(cfg.artfctdef.zvalue.method, 'trialdemean');
if pertrial
  if isfield(cfg.artfctdef.zvalue, 'ntrial') && cfg.artfctdef.zvalue.ntrial>0
    pertrial = cfg.artfctdef.zvalue.ntrial;
  else
    ft_error('you should specify cfg.artfctdef.zvalue.ntrial, and it should be > 0');
  end
end

% the data can be passed as input arguments or can be read from disk
hasdata = exist('data', 'var');

if ~hasdata
  cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
  cfg = ft_checkconfig(cfg, 'required', {'headerfile', 'datafile'});
  hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat);
else
  data = ft_checkdata(data, 'datatype', 'raw', 'hassampleinfo', 'yes');
  cfg  = ft_checkconfig(cfg, 'forbidden', {'dataset', 'headerfile', 'datafile'});
  hdr  = ft_fetch_header(data);
end

% set default cfg.continuous
if isempty(cfg.continuous)
  if hdr.nTrials==1
    cfg.continuous = 'yes';
  else
    cfg.continuous = 'no';
  end
end

% get the specification of the data segments that should be scanned for artifacts
if ~isfield(cfg, 'trl') && hasdata
  trl = data.sampleinfo;
  for k = 1:numel(data.trial)
    trl(k,3) = time2offset(data.time{k}, data.fsample);
  end
elseif isfield(cfg, 'trl') && ischar(cfg.trl)
  trl = loadvar(cfg.trl, 'trl');
elseif isfield(cfg, 'trl') && isnumeric(cfg.trl)
  trl = cfg.trl;
else
  ft_error('cannot determine which segments of data to scan for artifacts');
end

% check whether the value for trlpadding makes sense
if hasdata && cfg.artfctdef.zvalue.trlpadding > 0
  % negative trlpadding is allowed with in-memory data, since that would remove some data from each trial
  ft_error('you cannot use positive trlpadding with in-memory data');
end

trlpadding = round(cfg.artfctdef.zvalue.trlpadding*hdr.Fs);
fltpadding = round(cfg.artfctdef.zvalue.fltpadding*hdr.Fs);
artpadding = round(cfg.artfctdef.zvalue.artpadding*hdr.Fs);

trl(:,1)      = trl(:,1) - trlpadding;       % pad the trial with some samples, in order to detect
trl(:,2)      = trl(:,2) + trlpadding;       % artifacts at the edges of the relevant trials.
if size(trl,2)>= 3
  trl(:,3)    = trl(:,3) - trlpadding;       % the offset can of course be adjusted as well
elseif hasdata
  % reconstruct offset
  for tr=1:size(trl,1)
    % account for 0 might not be in data.time
    t0        = interp1(data.time{tr}, 1:numel(data.time{tr}), 0, 'linear', 'extrap');
    trl(tr,3) = -t0+1 - trlpadding;
  end
else
  % assuming that the trial starts at t=0s
  trl(:,3) = trl(:,1);
end

numtrl        = size(trl,1);
cfg.artfctdef.zvalue.channel = ft_channelselection(cfg.artfctdef.zvalue.channel, hdr.label);
chanindx      = match_str(hdr.label, cfg.artfctdef.zvalue.channel);
nchan         = length(chanindx);
thresholdsum  = strcmp(cfg.artfctdef.zvalue.cumulative, 'yes');

if nchan<1
  ft_error('no channels selected');
end

% read the data and apply preprocessing options
if ~pertrial
  sumval = zeros(nchan, 1);
  sumsqr = zeros(nchan, 1);
  numsmp = zeros(nchan, 1);
else
  sumval = zeros(nchan, numtrl);
  sumsqr = zeros(nchan, numtrl);
  numsmp = zeros(nchan, numtrl);
end

if strcmp(cfg.memory, 'high') % store data in memory, saving computation time below
  dat = cell(1, numtrl);
end

ft_progress('init', cfg.feedback, ['searching for artifacts in ' num2str(nchan) ' channels']);
for trlop=1:numtrl

  ft_progress(trlop/numtrl, 'processing trial %d from %d\n', trlop, numtrl);
  if hasdata
    thisdat = ft_fetch_data(data,        'header', hdr, 'begsample', trl(trlop,1)-fltpadding, 'endsample', trl(trlop,2)+fltpadding, 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'), 'skipcheckdata', 1);
  else
    thisdat = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', trl(trlop,1)-fltpadding, 'endsample', trl(trlop,2)+fltpadding, 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'), 'dataformat', cfg.dataformat);
  end
  thisdat = preproc(thisdat, cfg.artfctdef.zvalue.channel, offset2time(0, hdr.Fs, size(thisdat,2)), cfg.artfctdef.zvalue, fltpadding, fltpadding);

  if ~pertrial
    % accumulate the sum and the sum-of-squares
    sumval = sumval + nansum(thisdat,2);
    sumsqr = sumsqr + nansum(thisdat.^2,2);
    numsmp = numsmp + sum(isfinite(thisdat),2);
  else
    % store per trial the sum and the sum-of-squares
    sumval(:,trlop) = nansum(thisdat,2);
    sumsqr(:,trlop) = nansum(thisdat.^2,2);
    numsmp(:,trlop) = sum(isfinite(thisdat),2);
  end

  if strcmp(cfg.memory, 'high') % store data in memory, saving computation time below
    dat{trlop} = thisdat;
  end

end % for trlop
ft_progress('close');

if pertrial>1
  sumval = ft_preproc_smooth(sumval, pertrial)*pertrial;
  sumsqr = ft_preproc_smooth(sumsqr, pertrial)*pertrial;
  numsmp = ft_preproc_smooth(numsmp, pertrial)*pertrial;
end

% compute the average and the standard deviation
if strcmp(cfg.artfctdef.zvalue.zscore, 'yes')
  datavg = sumval./numsmp;
  datstd = sqrt(sumsqr./numsmp - (sumval./numsmp).^2);
else
  ft_warning('not performing z-scoring, note that the defined threshold has physical units');
  datavg = zeros(size(sumval));
  datstd = ones(size(sumval));
end

if strcmp(cfg.memory, 'low')
  ft_info('\n');
end

zmax  = cell(1, numtrl);
zsum  = cell(1, numtrl);
zindx = cell(1, numtrl);

% create a vector that indexes the trials, or is all 1, in order to a per trial
% z-scoring, or use a static std and mean
if pertrial
  indvec = 1:numtrl;
else
  indvec = ones(1,numtrl);
end

ft_progress('init', cfg.feedback, ['processing data in ' num2str(nchan) ' channels']);
for trlop = 1:numtrl

  if strcmp(cfg.memory, 'low') % store nothing in memory (note that we need to fetch/read and preproc AGAIN... *yawn*)
    ft_progress(trlop/numtrl, 'processing trial %d from %d\n', trlop, numtrl);
    options_getdata = {'header', hdr, 'begsample', trl(trlop,1)-fltpadding, 'endsample', trl(trlop,2)+fltpadding, 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no')};
    if hasdata
      thisdat = ft_fetch_data(data, options_getdata{:});
    else
      options_getdata = cat(2, options_getdata, {'dataformat', cfg.dataformat});
      thisdat = ft_read_data(cfg.datafile, options_getdata{:});
    end
    thisdat = preproc(thisdat, cfg.artfctdef.zvalue.channel, offset2time(0, hdr.Fs, size(thisdat,2)), cfg.artfctdef.zvalue, fltpadding, fltpadding);
  else
    thisdat = dat{trlop};
  end

  nsmp    = size(thisdat,2);

  zmax{trlop}  = -inf + zeros(1, nsmp);
  zsum{trlop}  =        zeros(1, nsmp);
  zindx{trlop} =        zeros(1, nsmp);

  ix           = indvec(trlop) * ones(1,nsmp);           % indexing vector dependent on the pertrial setting
  zdata        = (thisdat - datavg(:,ix))./datstd(:,ix); % convert the filtered data to z-values
  zsum{trlop}  = nansum(zdata,1);      % sum the z-values across channels
  [zmax{trlop},ind] = max(zdata,[],1); % find the maximum z-value and remember it
  zindx{trlop}      = chanindx(ind);   % also remember the channel number that has the largest z-value

end % for trlop
ft_progress('close');

if demeantrial
  for trlop = 1:numtrl
    zmax{trlop} = zmax{trlop}-nanmean(zmax{trlop},2);
    zsum{trlop} = zsum{trlop}-nanmean(zsum{trlop},2);
  end
end

for trlop = 1:numtrl
  zsum{trlop} = zsum{trlop} ./ sqrt(nchan);
end

% always create figure
% keypress to enable keyboard uicontrol
h = figure('KeyPressFcn', @keyboard_cb);
set(h, 'visible', 'off');

opt.artcfg       = cfg.artfctdef.zvalue;
opt.artval       = {};
opt.artpadding   = artpadding;
opt.cfg          = cfg;
opt.channel      = 'artifact';
opt.hdr          = hdr;
opt.numtrl       = size(trl,1);
opt.quit         = 0;
opt.threshold    = cfg.artfctdef.zvalue.cutoff;
opt.thresholdsum = thresholdsum;
opt.trialok      = true(1,opt.numtrl);  % OK by means of objective criterion
opt.keep         = zeros(1,opt.numtrl); % OK overruled by user +1 to keep, -1 to reject, start all zeros for callback to work
opt.trl          = trl;
opt.trlop        = 1;
opt.updatethreshold = true;
opt.zmax         = zmax;
opt.zsum         = zsum;

if ~thresholdsum
  opt.zval = zmax;
else
  opt.zval = zsum;
end
opt.zindx = zindx;
if ~hasdata
  opt.data = {};
else
  opt.data = data;
end

if strcmp(cfg.artfctdef.zvalue.interactive, 'yes')
  set(h, 'visible', 'on');
  set(h, 'CloseRequestFcn', @cleanup_cb);
  % give graphical feedback and allow the user to modify the threshold
  set(h, 'position', [100 200 900 400]);
  h1 = axes('position', [0.05 0.15 0.4 0.8]);
  h2 = axes('position', [0.5  0.57  0.45 0.38]);
  h3 = axes('position', [0.5  0.15  0.45 0.32]);
  opt.h1           = h1;
  opt.h2           = h2;
  opt.h3           = h3;

  setappdata(h, 'opt', opt);
  artval_cb(h);
  redraw_cb(h);

  % make the user interface elements for the data view, the order of the elements
  % here is from left to right and should match the order in the documentation
  uicontrol('tag', 'width1', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'stop',    'userdata', 'q');

  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '<',        'userdata', 'comma');
  uicontrol('tag', 'width1', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'artifact', 'userdata', 'a');
  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '>',        'userdata', 'period');

  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '<<',    'userdata', 'x');
  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '<',     'userdata', 'leftarrow');
  uicontrol('tag', 'width1', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'trial', 'userdata', 't');
  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '>',     'userdata', 'rightarrow');
  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '>>',    'userdata', 'c');

  uicontrol('tag', 'width3', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'keep trial',  'userdata', 'k');
  uicontrol('tag', 'width3', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'reject full', 'userdata', 'space');
  uicontrol('tag', 'width3', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'reject part', 'userdata', 'r');

  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '<',           'userdata', 'downarrow');
  uicontrol('tag', 'width3', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'threshold',   'userdata', 'z');
  uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '>',           'userdata', 'uparrow');

  %uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '<',       'userdata', 'control+uparrow')
  %uicontrol('tag', 'width1', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'channel', 'userdata', 'c')
  %uicontrol('tag', 'width2', 'parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '>',       'userdata', 'control+downarrow')

  ft_uilayout(h, 'tag', 'width1', 'width', 0.10, 'height', 0.05);
  ft_uilayout(h, 'tag', 'width2', 'width', 0.05, 'height', 0.05);
  ft_uilayout(h, 'tag', 'width3', 'width', 0.12, 'height', 0.05);

  ft_uilayout(h, 'tag', 'width1', 'style', 'pushbutton', 'callback', @keyboard_cb);
  ft_uilayout(h, 'tag', 'width2', 'style', 'pushbutton', 'callback', @keyboard_cb);
  ft_uilayout(h, 'tag', 'width3', 'style', 'pushbutton', 'callback', @keyboard_cb);

  ft_uilayout(h, 'tag', 'width1', 'retag', 'viewui');
  ft_uilayout(h, 'tag', 'width2', 'retag', 'viewui');
  ft_uilayout(h, 'tag', 'width3', 'retag', 'viewui');
  ft_uilayout(h, 'tag', 'viewui', 'BackgroundColor', [0.8 0.8 0.8], 'hpos', 'auto', 'vpos', 0.005);

  while opt.quit==0
    uiwait(h);
    opt = getappdata(h, 'opt');
  end

else
  % compute the artifacts given the settings in the cfg
  setappdata(h, 'opt', opt);
  artval_cb(h);
end

h   = getparent(h);
opt = getappdata(h, 'opt');

% convert the artifact values per trial to one long boolean vector
boolvec = zeros(1,max(opt.trl(:,2)));
for trlop=1:opt.numtrl
  boolvec(opt.trl(trlop,1):opt.trl(trlop,2)) = opt.artval{trlop};
end

% find the padded artifacts and put them in a Nx2 trl-like matrix
artifact = boolvec2artifact(boolvec);

if strcmp(cfg.artfctdef.zvalue.artfctpeak, 'yes')
  % this is a re-implementation of the peak-detection stuff, to make the
  % overall code behavior more consistent. artifact will be adjusted
  % according to the specifications of the user, i.e. the peak index for
  % each identified artifact will be identified, and used for an offset
  % column. the peak_indx, peaks, and dssartifact fields are be obsoleted
  pre = round(cfg.artfctdef.zvalue.artfctpeakrange(1)*hdr.Fs);
  pst = round(cfg.artfctdef.zvalue.artfctpeakrange(2)*hdr.Fs);
  for k = 1:size(artifact,1)
    % identify the corresponding trl for the current artifact, the artifact
    % can either be fully within the trl, or overlapping at one (or
    % both) edges
    current = artifact(k,:);
    seltrl  = find(current(2)>=opt.trl(:,1) & current(1)<=opt.trl(:,2));

    % in case the artifact is in more than one trial a for-loop is needed
    mx     = [];
    mx_idx = [];
    for m = 1:numel(seltrl)
      idx = current - opt.trl(seltrl(m),1) + 1;
      idx(1) = max(idx(1),1);
      idx(2) = min(idx(2),size(opt.zval{seltrl(m)},2));
      [mx(m), mx_idx(m)] = max(opt.zval{seltrl(m)}(idx(1):idx(2)));
    end

    [maxtrl, maxtrl_idx] = max(mx);
    seltrl = seltrl(maxtrl_idx);
    peak   = current(1) + mx_idx(maxtrl_idx) - 1;

    artifact(k,1) = max(peak+pre, opt.trl(seltrl,1));
    artifact(k,2) = min(peak+pst, opt.trl(seltrl,2));
    artifact(k,3) = artifact(k,1) - peak;
  end
end

if strcmp(cfg.representation, 'numeric') && istable(artifact)
  if isempty(artifact)
    % an empty table does not have columns
    artifact = zeros(0,2);
  else
    % convert the table to a numeric array with the columns begsample and endsample
    artifact = table2array(artifact);
  end
elseif strcmp(cfg.representation, 'table') && isnumeric(artifact)
  if isempty(artifact)
    % an empty table does not have columns
    artifact = table();
  else
    % convert the numeric array to a table with the columns begsample and endsample
    begsample = artifact(:,1);
    endsample = artifact(:,2);
    if size(artifact,2)==3
      offset   = artifact(:,3);
      artifact = table(begsample, endsample, offset);
    else
      artifact = table(begsample, endsample);
    end
  end
end

% remember the details that were used here and store the detected artifacts
cfg.artfctdef.zvalue.trl      = trl;              % remember where we have been looking for artifacts
cfg.artfctdef.zvalue.cutoff   = opt.threshold;    % remember the threshold that was used
cfg.artfctdef.zvalue.artifact = artifact;

ft_notice('detected %d artifacts\n', size(artifact,1));

delete(h);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble previous data
ft_postamble provenance
ft_postamble savevar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function artval_cb(h, eventdata)

opt = getappdata(h, 'opt');

artval = cell(1,opt.numtrl);
for trlop=1:opt.numtrl
  if opt.thresholdsum
    % threshold the accumulated z-values
    artval{trlop} = opt.zsum{trlop}>opt.threshold;
  else
    % threshold the max z-values
    artval{trlop} = opt.zmax{trlop}>opt.threshold;
  end
  % pad the artifacts
  artbeg = find(diff([0 artval{trlop}])== 1);
  artend = find(diff([artval{trlop} 0])==-1);
  artbeg = artbeg - opt.artpadding;
  artend = artend + opt.artpadding;
  artbeg(artbeg<1) = 1;
  artend(artend>length(artval{trlop})) = length(artval{trlop});
  for artlop=1:length(artbeg)
    artval{trlop}(artbeg(artlop):artend(artlop)) = 1;
  end
  opt.trialok(trlop) = isempty(artbeg);
end

for trlop = find(opt.keep==1 & opt.trialok==0)
  % overrule the objective criterion, i.e. keep the trial when the user
  % wants to keep it
  artval{trlop}(:) = 0;
end

for trlop = find(opt.keep<0 & opt.trialok==1)
  % if the user specifies that the trial is not OK
  % reject the whole trial if there is no extra-threshold data,
  % otherwise use the artifact as found by the thresholding
  if opt.thresholdsum && opt.keep(trlop)==-1
    % threshold the accumulated z-values
    artval{trlop} = opt.zsum{trlop}>opt.threshold;
  elseif opt.keep(trlop)==-1
    % threshold the max z-values
    artval{trlop} = opt.zmax{trlop}>opt.threshold;
  elseif opt.keep(trlop)==-2
    artval{trlop}(:) = 1;
  end
  % pad the artifacts
  artbeg = find(diff([0 artval{trlop}])== 1);
  artend = find(diff([artval{trlop} 0])==-1);
  artbeg = artbeg - opt.artpadding;
  artend = artend + opt.artpadding;
  artbeg(artbeg<1) = 1;
  artend(artend>length(artval{trlop})) = length(artval{trlop});
  if ~isempty(artbeg)
    for artlop=1:length(artbeg)
      artval{trlop}(artbeg(artlop):artend(artlop)) = 1;
    end
  else
    artval{trlop}(:) = 1;
  end
end

for trlop = find(opt.keep==-2 & opt.trialok==0)
  % if the user specifies the whole trial to be rejected define the whole
  % segment to be bad
  artval{trlop}(:) = 1;
  % pad the artifacts
  artbeg = find(diff([0 artval{trlop}])== 1);
  artend = find(diff([artval{trlop} 0])==-1);
  artbeg = artbeg - opt.artpadding;
  artend = artend + opt.artpadding;
  artbeg(artbeg<1) = 1;
  artend(artend>length(artval{trlop})) = length(artval{trlop});
  if ~isempty(artbeg)
    for artlop=1:length(artbeg)
      artval{trlop}(artbeg(artlop):artend(artlop)) = 1;
    end
  else
    artval{trlop}(:) = 1;
  end
end

opt.artval = artval;
setappdata(h, 'opt', opt);
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function keyboard_cb(h, eventdata)

% If a mouseclick was made, use that value. If not, determine the key that
% corresponds to the uicontrol element that was activated.

if isa(eventdata, 'matlab.ui.eventdata.ActionData') % only the case when clicked with mouse
  curKey = get(h, 'userdata');
elseif isa(eventdata, 'matlab.ui.eventdata.KeyData') % only when key was pressed
  if isempty(eventdata.Character) && any(strcmp(eventdata.Key, {'control', 'shift', 'alt', '0'}))
    % only a modifier key was pressed
    return
  end
  if isempty(eventdata.Modifier)
    curKey = eventdata.Key;
  else
    curKey = [sprintf('%s+', eventdata.Modifier{:}) eventdata.Key];
  end
elseif isfield(eventdata, 'Key')  % only when key was pressed
  curKey = eventdata.Key;
elseif isempty(eventdata) % matlab2012b returns an empty double upon a mouse click
  curKey = get(h, 'userdata');
else
  ft_error('cannot process user input, please report this on http://bugzilla.fieldtriptoolbox.org including your MATLAB version');
end

h = getparent(h); % otherwise h is empty if isa [...].ActionData
opt = getappdata(h, 'opt');

switch strtrim(curKey)
  case 'leftarrow' % change trials
    opt.trlop = max(opt.trlop - 1, 1); % should not be smaller than 1
    setappdata(h, 'opt', opt);
    redraw_cb(h, eventdata);
  case 'x'
    opt.trlop = max(opt.trlop - 10, 1); % should not be smaller than 1
    setappdata(h, 'opt', opt);
    redraw_cb(h, eventdata);
  case 'rightarrow'
    opt.trlop = min(opt.trlop + 1, opt.numtrl); % should not be larger than the number of trials
    setappdata(h, 'opt', opt);
    redraw_cb(h, eventdata);
  case 'c'
    opt.trlop = min(opt.trlop + 10, opt.numtrl); % should not be larger than the number of trials
    setappdata(h, 'opt', opt);
    redraw_cb(h, eventdata);
  case 'uparrow' % change threshold
    opt.threshold = opt.threshold+0.5;
    opt.updatethreshold = true;
    setappdata(h, 'opt', opt);
    artval_cb(h, eventdata);
    redraw_cb(h, eventdata);
    opt = getappdata(h, 'opt'); % grab the opt-structure from the handle because it has been adjusted in the callbacks
    opt.updatethreshold = false;
    setappdata(h, 'opt', opt);
  case 'downarrow'
    opt.threshold = opt.threshold-0.5;
    opt.updatethreshold = true;
    setappdata(h, 'opt', opt);
    artval_cb(h, eventdata);
    redraw_cb(h, eventdata);
    opt = getappdata(h, 'opt'); % grab the opt-structure from the handle because it has been adjusted in the callbacks
    opt.updatethreshold = false;
    setappdata(h, 'opt', opt);
  case 'period' % change artifact
    artfctindx = find(opt.trialok == 0);
    sel        = find(artfctindx>opt.trlop);
    if ~isempty(sel)
      opt.trlop = artfctindx(sel(1));
    end
    setappdata(h, 'opt', opt);
    redraw_cb(h, eventdata);
  case 'comma'
    artfctindx = find(opt.trialok == 0);
    sel        = find(artfctindx<opt.trlop);
    if ~isempty(sel)
      opt.trlop = artfctindx(sel(end));
    end
    setappdata(h, 'opt', opt);
    redraw_cb(h, eventdata);
    %   case 'control+uparrow' % change channel
    %     if strcmp(opt.channel, 'artifact')
    %       [dum, indx] = max(opt.zval);
    %       chanindx      = opt.zindx(indx);
    %     else
    %       if ~isempty(opt.data)
    %         chanindx  = match_str(opt.channel, opt.data.label);
    %         selchan = match_str(opt.artcfg.channel, opt.channel);
    %       else
    %         chanindx  = match_str(opt.channel,   opt.hdr.label);
    %         selchan = match_str(opt.artcfg.channel, opt.channel);
    %       end
    %     end
    %     numchan = numel(opt.artcfg.channel);
    %     chansel = min(selchan+1, numchan);
    %     % convert numeric array into cell-array with channel labels
    %     opt.channel = tmpchan(chansel);
    %     setappdata(h, 'opt', opt);
    %     redraw_cb(h, eventdata);
    %   case 'c' % select channel
    %     select = match_str([opt.artcfg.channel;{'artifact'}], opt.channel);
    %     opt.channel = select_channel_list([opt.artcfg.channel;{'artifact'}], select);
    %     setappdata(h, 'opt', opt);
    %     redraw_cb(h, eventdata);
    %   case 'control+downarrow'
    %     tmpchan = [opt.artcfg.channel;{'artifact'}]; % append the 'artifact' channel
    %     chansel = match_str(tmpchan, opt.channel);
    %     chansel = max(chansel-1, 1);
    %     % convert numeric array into cell-array with channel labels
    %     opt.channel = tmpchan(chansel);
    %     setappdata(h, 'opt', opt);
    %     redraw_cb(h, eventdata);
  case 'a'
    % select the artifact to display
    response = inputdlg(sprintf('artifact trial to display'), 'specify', 1, {num2str(opt.trlop)});
    if ~isempty(response)
      artfctindx = find(opt.trialok == 0);
      sel        = str2double(response);
      sel        = min(numel(artfctindx), sel);
      sel        = max(1,                 sel);
      opt.trlop  = artfctindx(sel);
      setappdata(h, 'opt', opt);
      redraw_cb(h, eventdata);
    end
  case 'q'
    setappdata(h, 'opt', opt);
    cleanup_cb(h);
  case 't'
    % select the trial to display
    response = inputdlg(sprintf('trial to display'), 'specify', 1, {num2str(opt.trlop)});
    if ~isempty(response)
      opt.trlop = str2double(response);
      opt.trlop = min(opt.trlop, opt.numtrl); % should not be larger than the number of trials
      opt.trlop = max(opt.trlop, 1); % should not be smaller than 1
      setappdata(h, 'opt', opt);
      redraw_cb(h, eventdata);
    end
  case 'z'
    % select the threshold
    response = inputdlg('z-threshold', 'specify', 1, {num2str(opt.threshold)});
    if ~isempty(response)
      opt.threshold = str2double(response);
      opt.updatethreshold = true;
      setappdata(h, 'opt', opt);
      artval_cb(h, eventdata);
      redraw_cb(h, eventdata);
      opt = getappdata(h, 'opt'); % grab the opt-structure from the handle because it has been adjusted in the callbacks
      opt.updatethreshold = false;
      setappdata(h, 'opt', opt);
    end
  case 'k'
    opt.keep(opt.trlop) = 1;
    setappdata(h, 'opt', opt);
    artval_cb(h);
    redraw_cb(h);
  case 'r'
    % only of the trial contains a partial artifact
    if opt.trialok(opt.trlop) == 0
      opt.keep(opt.trlop) = -1;
    end
    setappdata(h, 'opt', opt);
    artval_cb(h);
    redraw_cb(h);
  case 'space'
    opt.keep(opt.trlop) = -2;
    setappdata(h, 'opt', opt);
    artval_cb(h);
    redraw_cb(h);
  case 'control+control'
    % do nothing
  case 'shift+shift'
    % do nothing
  case 'alt+alt'
    % do nothing
  otherwise
    setappdata(h, 'opt', opt);
    % this should be consistent with the help of the function
    fprintf('----------------------------------------------------------------------\n');
    fprintf('     q                 : Stop\n');
    fprintf('\n');
    fprintf('     comma             : Step to the previous artifact trial\n');
    fprintf('     a                 : Specify artifact trial to display\n');
    fprintf('     period            : Step to the next artifact trial\n');
    fprintf('\n');
    fprintf('     x                 : Step 10 trials back\n');
    fprintf('     leftarrow         : Step to the previous trial\n');
    fprintf('     t                 : Specify trial to display\n');
    fprintf('     rightarrow        : Step to the next trial\n');
    fprintf('     c                 : Step 10 trials forward\n');
    fprintf('\n');
    fprintf('     k                 : Keep trial\n');
    fprintf('     space             : Mark complete trial as artifact\n');
    fprintf('     r                 : Mark part of trial as artifact\n');
    fprintf('\n');
    fprintf('     downarrow         : Shift the z-threshold down\n');
    fprintf('     z                 : Specify the z-threshold\n');
    fprintf('     uparrow           : Shift the z-threshold down\n');
    fprintf('----------------------------------------------------------------------\n');
end
clear curKey;
uiresume(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function redraw_cb(h, eventdata)

h   = getparent(h);
opt = getappdata(h, 'opt');

% make a local copy of the relevant variables
trlop     = opt.trlop;
artval    = opt.artval{trlop};
zindx     = opt.zindx{trlop};
zval      = opt.zval{trlop};
cfg       = opt.cfg;
artcfg    = opt.artcfg;
hdr       = opt.hdr;
trl       = opt.trl;
trlpadsmp = round(artcfg.trlpadding*hdr.Fs);
channel   = opt.channel;

% determine the channel with the highest z-value to be displayed
% this is default behavior but can be overruled in the gui
if strcmp(channel, 'artifact')
  [dum, indx] = max(zval);
  chanindx      = zindx(indx);
else
  if ~isempty(opt.data)
    chanindx = match_str(channel, opt.data.label);
  else
    chanindx = match_str(channel, hdr.label);
  end
end

if ~isempty(opt.data)
  data = ft_fetch_data(opt.data, 'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'));
else
  data = ft_read_data(cfg.datafile, 'header', hdr, 'begsample', trl(trlop,1), 'endsample', trl(trlop,2), 'chanindx', chanindx, 'checkboundary', strcmp(cfg.continuous, 'no'));
end

% data = preproc(data, '', hdr.Fs, artcfg, [], artcfg.fltpadding, artcfg.fltpadding);

% the string us used as title and printed in the command window
str = sprintf('trial %3d of %d, channel %s', trlop, size(trl,1), hdr.label{chanindx});
fprintf('showing %s\n', str);

%-----------------------------
% plot summary in left subplot
subplot(opt.h1); hold on;

% plot as a blue line only once
if isempty(get(opt.h1, 'children'))
  for k = 1:opt.numtrl
    xval = opt.trl(k,1):opt.trl(k,2);
    if opt.thresholdsum
      yval = opt.zsum{k};
    else
      yval = opt.zmax{k};
    end
    plot(opt.h1, xval, yval, 'linestyle', '-', 'color', 'b', 'displayname', 'data');
    xlabel('samples');
    ylabel('z-value');
  end
end
h1children = get(opt.h1, 'children');

% plot trial box
boxhandle = findall(h1children, 'displayname', 'highlight');
if isempty(boxhandle)
  % draw it
  xval = trl(opt.trlop,1):trl(opt.trlop,2);
  if opt.thresholdsum
    yval = opt.zsum{opt.trlop};
  else
    yval = opt.zmax{opt.trlop};
  end
  plot(opt.h1, xval, yval, 'linestyle', '-', 'color', 'm', 'linewidth', 2, 'displayname', 'highlight');
else
  % update it
  xval = trl(opt.trlop,1):trl(opt.trlop,2);
  if opt.thresholdsum
    yval = opt.zsum{opt.trlop};
  else
    yval = opt.zmax{opt.trlop};
  end
  set(boxhandle,  'XData', xval);
  set(boxhandle,  'YData', yval);
end

% plot as red lines the suprathreshold data points
thrhandle = findall(h1children, 'displayname', 'reddata');
if isempty(thrhandle)
  % they have to be drawn
  for k = 1:opt.numtrl
    xval = trl(k,1):trl(k,2);
    if opt.thresholdsum
      yval = opt.zsum{k};
    else
      yval = opt.zmax{k};
    end
    dum = yval<=opt.threshold;
    yval(dum) = nan;
    plot(opt.h1, xval, yval, 'linestyle', '-', 'color', [1 0 0], 'displayname', 'reddata');
  end
  hline(opt.threshold, 'color', 'r', 'linestyle', ':', 'displayname', 'threshline');
elseif ~isempty(thrhandle) && opt.updatethreshold
  % they can be updated
  for k = 1:opt.numtrl
    xval = trl(k,1):trl(k,2);
    if opt.thresholdsum
      yval = opt.zsum{k};
    else
      yval = opt.zmax{k};
    end
    dum = yval<=opt.threshold;
    yval(dum) = nan;
    set(thrhandle(k), 'XData', xval);
    set(thrhandle(k), 'YData', yval);
  end
  set(findall(h1children, 'displayname', 'threshline'), 'YData', [1 1].*opt.threshold);
end

%--------------------------------------------------
% get trial specific x-axis values and padding info
xval = ((trl(opt.trlop,1):trl(opt.trlop,2))-trl(opt.trlop,1)+trl(opt.trlop,3))./opt.hdr.Fs;
if trlpadsmp>0
  sel    = trlpadsmp:(size(data,2)-trlpadsmp);
  selpad = 1:size(data,2);
else
  sel    = 1:size(data,2);
  selpad = sel;
end

% plot data of most aberrant channel in upper subplot
subplot(opt.h2); hold on
if isempty(get(opt.h2, 'children'))
  % do the plotting
  plot(xval(selpad), data(selpad),          'color', [0.5 0.5 1], 'displayname', 'line1');
  plot(xval(sel),    data(sel),             'color', [0 0 1],     'displayname', 'line2');
  vline(xval(  1)+(trlpadsmp-1/opt.hdr.Fs), 'color', [0 0 0],     'displayname', 'vline1');
  vline(xval(end)-(trlpadsmp/opt.hdr.Fs),   'color', [0 0 0],     'displayname', 'vline2');
  data(~artval) = nan;
  plot(xval, data, 'r-', 'displayname', 'line3');
  xlabel('time(s)');
  ylabel('uV or Tesla');
  xlim([xval(1) xval(end)]);
  title(str);
else
  % update in the existing handles
  h2children = get(opt.h2, 'children');
  set(findall(h2children, 'displayname', 'vline1'), 'visible', 'off');
  set(findall(h2children, 'displayname', 'vline2'), 'visible', 'off');
  set(findall(h2children, 'displayname', 'line1'), 'XData', xval(selpad));
  set(findall(h2children, 'displayname', 'line1'), 'YData', data(selpad));
  set(findall(h2children, 'displayname', 'line2'), 'XData', xval(sel));
  set(findall(h2children, 'displayname', 'line2'), 'YData', data(sel));
  data(~artval) = nan;
  set(findall(h2children, 'displayname', 'line3'),  'XData', xval);
  set(findall(h2children, 'displayname', 'line3'),  'YData', data);
  abc2 = axis(opt.h2);
  set(findall(h2children, 'displayname', 'vline1'), 'XData', [1 1]*xval(  1)+(trlpadsmp-1/opt.hdr.Fs));
  set(findall(h2children, 'displayname', 'vline1'), 'YData', abc2(3:4));
  set(findall(h2children, 'displayname', 'vline2'), 'XData', [1 1]*xval(end)-(trlpadsmp/opt.hdr.Fs));
  set(findall(h2children, 'displayname', 'vline2'), 'YData', abc2(3:4));
  set(findall(h2children, 'displayname', 'vline1'), 'visible', 'on');
  set(findall(h2children, 'displayname', 'vline2'), 'visible', 'on');
  str = sprintf('trial %3d, channel %s', opt.trlop, hdr.label{chanindx});
  title(str);
  xlim([xval(1) xval(end)]);
end

% plot z-values in lower subplot
subplot(opt.h3); hold on;
if isempty(get(opt.h3, 'children'))
  % do the plotting
  plot(xval(selpad), zval(selpad), 'color', [0.5 0.5 1], 'displayname', 'line1b');
  plot(xval(sel),    zval(sel),    'color', [0 0 1],     'displayname', 'line2b');
  hline(opt.threshold, 'color', 'r', 'linestyle', ':', 'displayname', 'threshline');
  vline(xval(  1)+(trlpadsmp-1/opt.hdr.Fs),     'color', [0 0 0],     'displayname', 'vline1b');
  vline(xval(end)-(trlpadsmp/opt.hdr.Fs),       'color', [0 0 0],     'displayname', 'vline2b');
  zval(~artval) = nan;
  plot(xval, zval, 'r-', 'displayname', 'line3b');
  xlabel('time(s)');
  ylabel('z-value');
  xlim([xval(1) xval(end)]);
else
  % update in the existing handles
  h3children = get(opt.h3, 'children');
  set(findall(h3children, 'displayname', 'vline1b'), 'visible', 'off');
  set(findall(h3children, 'displayname', 'vline2b'), 'visible', 'off');
  set(findall(h3children, 'displayname', 'line1b'), 'XData', xval(selpad));
  set(findall(h3children, 'displayname', 'line1b'), 'YData', zval(selpad));
  set(findall(h3children, 'displayname', 'line2b'), 'XData', xval(sel));
  set(findall(h3children, 'displayname', 'line2b'), 'YData', zval(sel));
  zval(~artval) = nan;
  set(findall(h3children, 'displayname', 'line3b'),     'XData', xval);
  set(findall(h3children, 'displayname', 'line3b'),     'YData', zval);
  set(findall(h3children, 'displayname', 'threshline'), 'YData', [1 1].*opt.threshold);
  set(findall(h3children, 'displayname', 'threshline'), 'XData', xval([1 end]));
  abc = axis(opt.h3);
  set(findall(h3children, 'displayname', 'vline1b'), 'XData', [1 1]*xval(  1)+(trlpadsmp-1/opt.hdr.Fs));
  set(findall(h3children, 'displayname', 'vline1b'), 'YData', abc(3:4));
  set(findall(h3children, 'displayname', 'vline2b'), 'XData', [1 1]*xval(end)-(trlpadsmp/opt.hdr.Fs));
  set(findall(h3children, 'displayname', 'vline2b'), 'YData', abc(3:4));
  set(findall(h3children, 'displayname', 'vline1b'), 'visible', 'on');
  set(findall(h3children, 'displayname', 'vline2b'), 'visible', 'on');
  xlim([xval(1) xval(end)]);
end

setappdata(h, 'opt', opt);
uiresume

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cleanup_cb(h, eventdata)
opt = getappdata(h, 'opt');
opt.quit = true;
setappdata(h, 'opt', opt);
uiresume

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = getparent(h)
p = h;
while p~=0
  h = p;
  p = get(h, 'parent');
end


text:
function segmentation = convert_segmentationstyle(segmentation, fn, dim, style)

% CONVERT_SEGMENTATIONSTYLE is a helper function for converting between probabilistic
% and indexed representations. It is used by ft_datatype_segmentation and
% ft_datatype_parcellation.
%
% See also FIXSEGMENTATION, DETERMINE_SEGMENTATIONSTYLE

switch style
  case 'indexed'

    % convert the probabilistic representation to an indexed representation
    threshold = 0.5;
    seg       = zeros(dim);
    seglabel  = cell(1,numel(fn));

    for i=1:length(fn)
      tmp = segmentation.(fn{i})>threshold;
      if any(seg(tmp(:)))
        ft_error('overlapping tissue probability maps cannot be converted to an indexed representation');
        % FIXME in principle it is possible to represent two tissue types at one voxel
      else
        seg(tmp) = i;
        seglabel{i} = fn{i};
      end
    end
    segmentation          = rmfield(segmentation, fn);
    segmentation.seg      = seg;
    segmentation.seglabel = seglabel;

  case 'probabilistic'
    % convert the indexed representation to a probabilistic one
    for i=1:length(fn)
      fprintf('converting %s\n', fn{i});
      seg      = segmentation.(fn{i});
      seglabel = segmentation.([fn{i} 'label']);
      for j=i:length(seglabel)
        fprintf('creating probabilistic representation for %s\n', seglabel{j});
        tmp.(fixname(seglabel{j})) = (seg==j);  % avoid overwriting the existing segmentation
      end % for j
      segmentation = rmfield(segmentation,  fn{i}         );
      segmentation = rmfield(segmentation, [fn{i} 'label']);
      for j=i:length(seglabel)
        segmentation.(fixname(seglabel{j})) = tmp.(fixname(seglabel{j})); % avoid overwriting the existing segmentation
      end
    end % for i

  otherwise
    ft_error('unsupported style "%s"', style);
end

author: spm


text:
%DEMO_REGRESSION_HIER  Hierarchical regression demonstration using
%                      Rats data
%
%  Description
%    The example data is taken from section 6 of Gelfand et al
%    (1990) (also used in WinBUGS/OpenBUGS), and concerns 30 young
%    rats whose weights were measured weekly for five week. This
%    demo demosntrates how to make hierarchical linear and
%    non-linear models.
%
%  Reference
%    Gelfand, A. E., Hills, S. E., Racine-Poon, A. and Smith, A. F.
%    M. (1990) Illustration of Bayesian Inference in Normal Data
%    Models Using Gibbs Sampling. Journal of the American
%    Statistical Association 85(412):972-985.
%

% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

S = which('demo_regression_hier');
L = strrep(S,'demo_regression_hier.m','demodata/rats.mat');
data=load(L);
xx = data.x;
yy = data.y;
% Show data : 5 weight measurements per rat for 30 rats
figure
plot(xx,yy,'o-')
axis([7 37 100 400])
title('Data')
drawnow

% Reshape data
ntime = size(xx,2);
nrats = size(yy,1);
% All y's to one vector
y=yy(:);
% Repeat x for each rat
x=reshape(repmat(xx,nrats,1),ntime*nrats,1);
% Add ratid
x=[x repmat([1:nrats]',ntime,1)];
% Now 'x' consist of the inputs (ratid,time) and 'y' of the output (weight).
% Normalize x and y
[xn,xmean,xstd]=normdata(x);
[yn,ymean,ystd]=normdata(y);

% optmization options
opt=optimset('TolFun',1e-4,'TolX',1e-4,'Display','on');

% common Gaussian likelihood with weakly informative prior for variance
lik=lik_gaussian('sigma2',.1,...
                 'sigma2_prior',prior_sinvchi2('s2',0.01,'nu',1));
% common categorical covariance term
cc=gpcf_cat('selectedVariables',2);

% 1) Linear model with intercept and slope wrt time
disp('1) Linear model with intercept and slope wrt time')
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_t());
cfl=gpcf_linear('coeffSigma2',1,'selectedVariables',1,...
                'coeffSigma2_prior',prior_t());
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfl});
% optimize
gp=gp_optim(gp,xn,yn,'opt',opt);
% predict and plot
Ef=gp_pred(gp,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,1)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Linear model')
drawnow

% 2) Linear model with hierarchical intercept
disp('2) Linear model with hierarchical intercept')
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_t());
cfl=gpcf_linear('coeffSigma2',1,'selectedVariables',1,...
                'coeffSigma2_prior',prior_t());
% own constant term for each rat
cfci=gpcf_prod('cf',{cfc cc});
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfci cfl});
% optimize
gp=gp_optim(gp,xn,yn,'opt',opt);
% predict and plot
Ef=gp_pred(gp,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,2)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Linear model with hierarchical intercept')
drawnow

% 3) Linear model with hierarchical intercept and slope
disp('3) Linear model with hierarchical intercept and slope')
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_t());
cfl=gpcf_linear('coeffSigma2',1,'selectedVariables',1,...
                'coeffSigma2_prior',prior_t());
% own constant term for each rat
cfci=gpcf_prod('cf',{cfc cc});
% linear covariance term for each rat
cfli=gpcf_prod('cf',{cfl cc});
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfci cfl cfli});
% optimize
gp=gp_optim(gp,xn,yn,'opt',opt);
% predict and plot
Ef=gp_pred(gp,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,3)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Linear model with hierarchical intercept and slope')
drawnow

% 4) Nonlinear model with hierarchical intercept
% include linear part, too
disp('4) Nonlinear model with hierarchical intercept')
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_t());
cfl=gpcf_linear('coeffSigma2',1,'selectedVariables',1,...
                'coeffSigma2_prior',prior_t());
% own constant term for each rat
cfci=gpcf_prod('cf',{cfc cc});
% nonlinear part
cfs=gpcf_sexp('selectedVariables',1);
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfci cfl cfs});
% optimize
gp=gp_optim(gp,xn,yn,'opt',opt);
% predict and plot
Ef=gp_pred(gp,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,4)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Non-linear model with hierarchical intercept')
drawnow

% 5) Nonlinear model with hierarchical intercept and curve
% include linear part, too
disp('5) Non-linear hierarchical model 1 with MAP')
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_t());
cfl=gpcf_linear('coeffSigma2',1,'selectedVariables',1,...
                'coeffSigma2_prior',prior_t());
% own constant term for each rat
cfci=gpcf_prod('cf',{cfc cc});
% linear covariance term for each rat
cfli=gpcf_prod('cf',{cfl cc});
% nonlinear part
cfs=gpcf_sexp('selectedVariables',1);
% nonlinear covariance term for each rat
cfsi=gpcf_prod('cf',{cfs cc});
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfci cfl cfli cfs cfsi},...
          'jitterSigma2',1e-6);
% optimize
gp=gp_optim(gp,xn,yn,'opt',opt);
% predict and plot
Ef=gp_pred(gp,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,5)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Non-linear hierarchical model 1 with MAP')
drawnow

% 6) With increasing flexibility of the modeling function
%    we need to integrate over the parameteres
% integrate over parameters
disp('6) Non-linear hierarchical model 1 with IA')
[gps,pth,th]=gp_ia(gp,xn,yn);
% predict and plot
Ef=gp_pred(gps,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,6)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Non-linear hierarchical model 1 with IA')
drawnow

% 7) Nonlinear model with hierarchical intercept and curve
%    Same as 5, but with no linear and product covariances
disp('7) Non-linear hierarchical model 2 with MAP')
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_t());
% own constant term for each rat
cfci=gpcf_prod('cf',{cfc cc});
% nonlinear part with delta distance for ratid
cfs=gpcf_sexp('metric',metric_euclidean('components',{[1] [2]},...
                                        'deltadist', [0 1], ...
                                        'lengthScale_prior',prior_t()));
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfci cfs});
% optimize
gp=gp_optim(gp,xn,yn,'opt',opt);
% predict and plot
Ef=gp_pred(gp,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,7)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Non-linear hierarchical model 2 with MAP')
drawnow

% 8) With increasing flexibility of the modeling function
%    we need to integrate over the parameteres
% integrate over parameters
disp('8) Non-linear hierarchical model 2 with IA')
gps=gp_ia(gp,xn,yn);
% predict and plot
Ef=gp_pred(gps,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,8)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Non-linear hierarchical model 2 with IA')
drawnow

% 9) With neuralnetwork covariance and integration over the parameters
disp('9) Non-linear hierarchical model 3 with IA')
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_t());
% own constant term for each rat
cfci=gpcf_prod('cf',{cfc cc});
% nonlinear part with neuralnetwork covariance
cfnn=gpcf_neuralnetwork('selectedVariables',1,'biasSigma2_prior',prior_t(),...
                        'weightSigma2_prior',prior_t());
% nonlinear covariance term for each rat
cfnni=gpcf_prod('cf',{cfnn cc});
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfci cfnn cfnni});
% optimize
gp=gp_optim(gp,xn,yn,'opt',opt);
% integrate over parameters
gps=gp_ia(gp,xn,yn);
% predict and plot
Ef=gp_pred(gps,xn,yn,xn);
Eff=reshape(denormdata(Ef,ymean,ystd),nrats,ntime);
subplot(3,3,9)
plot(xx,Eff,'o-')
axis([7 37 100 400])
title('Non-linear hierarchical model 3 with IA')
drawnow

%*** Missing Data Example ***
% In the original paper (Gelfand et al, 1990) data was also to
% demonstrate missing data handling, by removing 1-4 weeks of data
% for part of the rats. Handling missing data in this case is trivial
% for GPs, too

disp('10) Missing data example')
S = which('demo_regression_hier');
L = strrep(S,'demo_regression_hier.m','demodata/rats.mat');
data=load(L);
xx = data.x;
yy = data.y;
yymiss = data.ymiss;
% Show data : 5 weight measurements per rat for 30 rats
figure
plot(xx,yymiss,'o-')
axis([7 37 100 400])
title('Data')
% Reshape data
ntime = size(xx,2);
nrats = size(yy,1);
% All y's to one vector
ym=yymiss(:);
% Repeat x for each rat
xm=reshape(repmat(xx,nrats,1),ntime*nrats,1);
% Add ratid
xm=[xm repmat([1:nrats]',ntime,1)];
% Now 'x' consist of the inputs (ratid,time) and 'y' of the output (weight).
% Normalize x and y
[xmn,xmmean,xmstd]=normdata(xm);
[ymn,ymmean,ymstd]=normdata(ym);
% test x is the complete x
xmnt=xmn;
% remove missing data from the training data
missi=isnan(ym);
ymn(missi,:)=[];
xmn(missi,:)=[];

% 10) neuralnetwork covariance, IA and missing data
cfc=gpcf_constant('constSigma2',1,'constSigma2_prior',prior_gaussian());
% own constant term for each rat
cfci=gpcf_prod('cf',{cfc cc});
% nonlinear part with neuralnetwork covariance
cfnn=gpcf_neuralnetwork('selectedVariables',1,'biasSigma2_prior',prior_gaussian('s2',10));
% nonlinear covariance term for each rat
cfnni=gpcf_prod('cf',{cfnn cc});
% construct GP
gp=gp_set('lik',lik,'cf',{cfc cfci cfnn cfnni});
% optimize
gp=gp_optim(gp,xmn,ymn);
% integrate over parameters
gps=gp_ia(gp,xmn,ymn);
% predict and plot
Ef=gp_pred(gps,xmn,ymn,xmnt);
Eff=reshape(denormdata(Ef,ymmean,ymstd),nrats,ntime);
Effc=Eff;Effc(isnan(ym))=NaN;
plot(xx,Effc,'bo-',xx,Eff,'bo--')
axis([7 37 100 400])
title('Non-linear hierarchical model 3 with IA and missing data')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [n, fn] = dimlength(data, seldim, fld)

% DIMLENGTH(DATA, SELDIM, FLD) is a helper function to obtain n, the number
% of elements along dimension seldim from the appropriate field from the
% input data containing functional data.
%
% Use als
%   [n, fn] = dimlength(data, seldim, fld)
%
% It can be called with one input argument only, in which case it will
% output two cell arrays containing the size of the functional fields,
% based on the XXXdimord, and the corresponding XXXdimord fields.
%
% When the data contains a single dimord field (everything except source
% data), the cell-arrays in the output only contain one element.
%
% See also FIXSOURCE, CREATEDIMORD

% Copyright (C) 2010, Jan-Mathijs Schoffelen
%
% $Id$

if nargin<3
  fld = 'dimord';
end

if nargin==1
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % call recursively to get the dimensionality of all fields XXX which are
  % accompanied by a XXXdimord field or the general dimord field
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % get all fields of the data structure
  fn = fieldnames(data);
  % get all dimord like fields
  fn    = fn(~cellfun('isempty', strfind(fn, 'dimord')));
  for k = 1:numel(fn)
    fndimord{k} = data.(fn{k});
  end

  for k = 1:numel(fndimord)
    dimtok = tokenize(fndimord{k}, '_');
    ndim   = numel(dimtok);
    n{k,1} = zeros(1, ndim);
    for i = 1:ndim
      n{k}(i) = dimlength(data, dimtok{i}, fn{k});
    end
  end
  return

else
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % do the actual work
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  switch seldim
    case 'rpt'
      if numel(fld)>6 && isfield(data, fld(1:end-6))
        % source level data
        dimtok = tokenize(data.(fld), '_');
        tmp    = data.(fld(1:end-6));
        % remove the first cell-dimension
        if iscell(tmp)
          if isfield(data, 'inside')
            tmp = tmp{data.inside(1)};
          else
            tmp = tmp{1};
          end
          dimtok = dimtok(2:end);
        end
        ix = find(~cellfun('isempty', strfind(dimtok, seldim)));
        n  = size(tmp, ix);

      elseif strcmp(data.(fld), 'rpt_pos')
        %FIXME HACK to be fixed
        x = setdiff(fld(data),'inside');
        for k = 1:length(x)
          dims = size(getsubfield(data,x{k}));
          if dims(2)==size(data.pos,1) && numel(dims)==2,
            n = dims(1);
            return
          end
        end

      elseif strcmp(data.(fld), 'rpt_pos_freq')
        %FIXME HACK to be fixed
        x = fld(data);
        for k = 1:length(x)
          dims = size(getsubfield(data,x{k}));
          if dims(2)==size(data.pos,1) && (numel(dims)==2 || dims(3)==length(data.freq)),
            n = dims(1);
            return
          end
        end

      elseif strcmp(data.(fld), 'rpt_pos_time')
        %FIXME HACK to be fixed
        x = fld(data);
        for k = 1:length(x)
          dims = size(getsubfield(data,x{k}));
          if dims(2)==size(data.pos,1) && (numel(dims)==2 || dims(3)==length(data.time)),
            n = dims(1);
            return
          end
        end

      elseif strcmp(data.(fld)(1:4), 'rpt_')
        n  = [];
        % generic solution for XXXspctrm
        fnames = fieldnames(data);
        tmp    = find(~cellfun('isempty', strfind(fnames, 'spctrm')));
        for k = 1:numel(tmp)
          n = [n size(data.(fnames{tmp(k)}), 1)];
        end
        % some other possibilities
        if isfield(data, 'cov'),           n = [n size(data.cov,           1)]; end
        %if isfield(data, 'crsspctrm'),     n = [n size(data.crsspctrm,     1)]; end
        %if isfield(data, 'powcovspctrm'),  n = [n size(data.powcovspctrm,  1)]; end
        %if isfield(data, 'powspctrm'),     n = [n size(data.powspctrm,     1)]; end
        if isfield(data, 'trial'),         n = [n size(data.trial,         1)]; end
        %if isfield(data, 'fourierspctrm'), n = [n size(data.fourierspctrm, 1)]; end
        if isfield(data, 'individual'),    n = [n size(data.individual,    1)]; end
        if isfield(data, 'stat'),          n = [n size(data.stat,          1)]; end
        if ~all(n==n(1))
          ft_error('inconsistent number of repetitions for dim "%s"', seldim);
        end
        n = n(1);

      else
        %error('cannot determine number of repetitions for dim "%s"', seldim);
        n = nan;
      end

    case 'rpttap'
      if numel(fld)>6 && isfield(data, fld(1:end-6)),
        dimtok = tokenize(data.(fld), '_');
        tmp    = data.(fld(1:end-6));
        % remove the first cell-dimension
        if iscell(tmp)
          if isfield(data, 'inside'),
            tmp = tmp{data.inside(1)};
          else
            tmp = tmp{1};
          end
          dimtok = dimtok(2:end);
        end
        ix = find(~cellfun('isempty', strfind(dimtok, seldim)));
        n  = size(tmp, ix);

      elseif strcmp(data.(fld)(1:7), 'rpttap_')
        n  = [];
        if isfield(data, 'cov'),           n = [n size(data.cov,           1)]; end
        if isfield(data, 'crsspctrm'),     n = [n size(data.crsspctrm,     1)]; end
        if isfield(data, 'powcovspctrm'),  n = [n size(data.powcovspctrm,  1)]; end
        if isfield(data, 'powspctrm'),     n = [n size(data.powspctrm,     1)]; end
        if isfield(data, 'trial'),         n = [n size(data.trial,         1)]; end
        if isfield(data, 'fourierspctrm'), n = [n size(data.fourierspctrm, 1)]; end
        if ~all(n==n(1))
          ft_error('inconsistent number of repetitions for dim "%s"', seldim);
        end
        n = n(1);

      else
        %error('cannot determine number of repetitions for dim "%s"', seldim);
        n = nan;
      end

    case 'subj'
      n  = [];
      if isfield(data, 'cov'),           n = [n size(data.cov,           1)]; end
      if isfield(data, 'crsspctrm'),     n = [n size(data.crsspctrm,     1)]; end
      if isfield(data, 'fourierspctrm'), n = [n size(data.fourierspctrm, 1)]; end
      if isfield(data, 'individual'),    n = [n size(data.individual,    1)]; end
      if isfield(data, 'powcovspctrm'),  n = [n size(data.powcovspctrm,  1)]; end
      if isfield(data, 'powspctrm'),     n = [n size(data.powspctrm,     1)]; end
      if isfield(data, 'trial'),         n = [n size(data.trial,         1)]; end

      if ~all(n==n(1)), ft_error('inconsistent number of repetitions for dim "%s"', seldim); end
      n = n(1);

    case 'chan'
      if ~isfield(data, 'inside'),
        n = length(data.label);
      else
        %error('cannot determine number of repetitions for dim "%s"', seldim);
        n = nan; %FIXME discuss appending label to source-like data
      end

    case 'chancmb'
      if ~isfield(data, 'inside'),
        n = size(data.labelcmb, 1);
      else
        %error('cannot determine number of repetitions for dim "%s"', seldim);
        n = nan; %FIXME discuss appending label to source-like data
      end

    case 'freq'
      n = length(data.freq);

    case 'time'
      n = length(data.time);

    case {'pos' '{pos}'}
      n = size(data.pos,1);

    case {'ori'}
      % this is comparable to rpt and to rpttap
      if numel(fld)>6 && isfield(data, fld(1:end-6)),
        dimtok = tokenize(data.(fld), '_');
        tmp    = data.(fld(1:end-6));
        if iscell(tmp)
          if isfield(data, 'inside'),
            tmp = tmp{data.inside(1)};
          else
            tmp = tmp{1};
          end
          dimtok = dimtok(2:end);
        end
        % it can be {pos}_ori_ori, in which case only the first orientation needs to be checked
        ix = find(~cellfun('isempty', strfind(dimtok, seldim)), 1, 'first');
        n  = size(tmp, ix);
      else
        %error('cannot determine number of repetitions for dim "%s"', seldim);
        n = nan;
      end

    otherwise
      ft_error('unsupported dim "%s"', seldim);
  end

end % if nargin


text:
function U = unionHull(varargin)
%unionHull - computes the convex hull of the union of other bounded convex polyhedra.
%Any number of polyhedra can be input to the union operation. They can
%be expressed using either vertices or linear in/equalities, according to
%the input scheme,
%
%
%   I = intersectionHull('vert', V1, 'lcon', A2,b2,...,'qlcon', x3,A3,b3,...)
%
%The arguments specifying different polyhedra are separated using labels
%'vert', 'lcon', or 'qlcon'. The label 'vert' signifies that an input polyhedron will
%be expressed using vertices and is to be followed by any string of input arguments
%accepted by vert2lcon(). The label 'lcon' signifies that an input polyhedron will
%be expressed using linear constraints and is to be followed by any string of
%input arguments accepted by lcon2vert(). The label 'qlcon' signifies that an input
%polyhedron will be expressed using linear constraints and a known interior point and
%is to be followed by  input arguments accepted by qlcon2vert().
%
%The output, U, is a struct containing fields
%
%   U.vert: A matrix whose rows are the vertices of the output polyhedron.
%   U.lcon: The quadruplet of linear constraint data {A,b,Aeq,beq}
%           describing the output polyhedron.
%
%EXAMPLE 1: The unit simplex in 3D is formed as the convex hull of the union
%of a cube and a triangle, both expressed in terms of their vertices.
%
%     V1=(dec2bin(0:2^3-1,3)-'0')/3.0001;    %vertices of 3D cube
%     V2=eye(3);  %vertices of triangle
%
%     U=unionHull('vert',V1,'vert',V2);   %compute the union
%
%One can verify that the vertices of the result are those of the unit
%simplex,
%
%       >> U.vert
%
%             ans =
%
%                  0     0     0
%                  0     0     1
%                  0     1     0
%                  1     0     0
%
%We also demonstrate the use of separateBounds() to convert the result to
%Optimization Toolbox representation. The tolerance parameter 0.99 is used
%to deal with residual floating point noise.
%
%         >> [A,b,Aeq,beq,lb,ub]=separateBounds(U.lcon{:},.99)
%
%             A =
%
%                 0.5774    0.5774    0.5774
%
%
%             b =
%
%                 0.5774
%
%
%             Aeq =
%
%                  []
%
%
%             beq =
%
%                Empty matrix: 0-by-1
%
%
%             lb =
%
%                  0
%                  0
%                  0
%
%
%             ub =
%
%                Inf
%                Inf
%                Inf
%
%
%
%EXAMPLE 2: The unit cube is formed as the union of a simplex, a square,
%and a single 3D point. The simplex and square are both represented using
%equalites and inequalities. For convenience, we use the addBounds() utility
%to compose the in/equality matrices. The representation of the square illustrates
%the use of the 'qlcon' label, when a known interior point is available.
%
%
%     [A1,b1,Aeq1,beq1] = addBounds([1 1 1],1,[],[],[0;0;0]); %unit simplex
%
%     [A2,b2,Aeq2,beq2] = addBounds([],[],[],[],[1;0;0],[1;1;1]); %a square face of the unit cube
%
%                    x2 = [1,0.5,0.5];%a point in the relative interior of the square
%
%     V3 = [0,1,1]; %one vertex of the unit cube
%
%
%     U=unionHull('lcon',A1,b1,Aeq1,beq1,...         %compute the union
%                 'qlcon',x2,A2,b2,Aeq2,beq2,...
%                  'vert',V3);
%
%One can verify that the result contains the 8 vertices of the unit cube,
%
%       >> U.vert
%
%         ans =
%
%                  0   -0.0000   -0.0000
%             0.0000   -0.0000    1.0000
%            -0.0000    1.0000   -0.0000
%                  0    1.0000    1.0000
%             1.0000    0.0000    0.0000
%             1.0000    0.0000    1.0000
%             1.0000    1.0000    0.0000
%             1.0000    1.0000    1.0000
%
%As before, we can use separateBounds() to convert to Optimization Toolbox
%form,
%
%     >>  [A,b,Aeq,beq,lb,ub]=separateBounds(U.lcon{:},.99)
%
%which results in A=b=Aeq=beq=[] and
%
%    lb =
%
%          0
%          0
%          0
%
%
%     ub =
%
%         1.0000
%         1.0000
%         1.0000

%%%%begin parsing

if isnumeric(varargin{1})
   TOLcell=varargin(1);
   varargin(1)=[];
else
    TOLcell={};
end

N=length(varargin);
idxType = [find(cellfun(@ischar,varargin)),N+1];

L=length(idxType)-1;
S(L).type=[];
S(L).args={};

for i=1:L

    j=idxType(i);
    k=idxType(i+1);
    S(i).type=varargin{j};
    S(i).args=varargin(j+1:k-1);

    if isempty(S(i).args)
     error 'Syntax error - arguments missing'
    end

     switch S(i).type

        case 'lcon'

             S(i).V = lcon2vert(S(i).args{:});

        case 'qlcon'

            S(i).V = qlcon2vert(S(i).args{:});

        case 'vert'

            S(i).V = S(i).args{1};

        otherwise

            error(['Unrecognized representation label of polyhedron ' num2str(i)]);

      end

 end

%%%%end parsing

    V=vertcat(S.V);

    x0=mean(V,1);

   [U.lcon{1:4}]=vert2lcon(V,TOLcell{:});

   U.vert=qlcon2vert(x0,U.lcon{:},TOLcell{:});

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function showfigs()
%
% The function showfigs() generates a "slide show" of all open figures,
% bringing each open figure to the fore in turn and then pausing until the
% user issues a command by pressing a key.  Key commands are as follows:
%
% - Press Enter, down-arrow, or spacebar to advance to the next figure
% (next higher figure number).  If the current figure is the last (highest
% figure number), the slide show returns to figure 1.
%
% - Press the up-arrow key to return to the previous figure (next lower
% figure number).  If the current figure is figure 1, the slide show cycles
% back to the last figure (highest figure number).
%
% - Press a numeric key to display the figure whose handle has that value.
% If the user presses 0 (zero), the figure having a handle value of 10 is
% displayed.  If there is no figure corresponding to the key that was
% pressed, showfigs() writes an error message to the command window,
% generates a beep, and then continues running.
%
% - Press the Delete key to close the current figure window.  showfigs()
% decrements the figure handle number (wrapping around from 1 to the
% highest figure number if necessary) until a handle number corresponding
% to a figure that is still open is found, and will then display that
% figure.  (If all figures have been deleted, an error message is displayed
% and showfigs() terminates).
%
% - Press the End key to display the last figure (highest figure number,
% not counting any figures that have been deleted).
%
% - Press the Home key to display the first figure (lowest figure number,
% not counting any figures that have been deleted).
%
% - Press Escape or "q" to terminate the slide show.
%
% If the user presses a key that does not correspond to any showfigs()
% command, showfigs() writes an error message to the command window,
% generates a beep, and then terminates.
%
% Author: Dr. Phillip M. Feldman
%
% Note: I'd like to acknowledge many helpful suggestions from Jen Dobson of
% the MathWorks.
%
% Revision History
%
% Version 1.4, 18 Jun 2009: Modified so that pressing an invalid key does
% not eliminate the callback function (which had the effect of terminating
% the slideshow).  Added code to handle right and left arrows.
%
% Version 1.3, 7 May 2009: Streamlined code that handles numeric keys.
% Added support for "0" key (brings up figure having handle==10).  Added
% error messages.  Added sound feedback to alert user when an invalid key
% has been pressed.
%
% Version 1.2, 7 May 2009: Fixed bug (code that processes numeric keys was
% assuming that handle numbers start with 1 and are consecutive, which is
% not always true).
%
% Version 1.1, 6 May 2009: Modified code to check for deleted figures using
% ishandle() function.
%
% Version 1.0, 5 May 2009: Initial version.

global handles current;

% Find all objects of type axes:
handles= findobj(0, '-depth',1, 'type','figure');

if isempty(handles)
   fprintf(2,'There are no figures to display.\n');
   return
end

% Sort handles so that figures will be displayed in order:
handles= sort(handles);

% Place first figure on top:
figure(handles(1));
current= 1;

% Define callback function to process key presses:
set(handles(1), 'KeyPressFcn',@key_press_handler);

end % function showfigs()

function key_press_handler(source,event)
% keyPressFcn automatically takes in two inputs:
% - source is the object that was active when the keypress occurred
% - event stores the data for the key pressed

global handles current;

% Bring in the handles structure in to the function
handle= guidata(source);

% Remove callback for current figure:
set(handles(current), 'KeyPressFcn',[]);

k= event.Key; % k is the key that is pressed

% The user can press the Escape key or "q" to end the slide show.  Since we
% have already removed the callback, there is nothing more to be done:

if strcmp(k,'escape') | strcmp(k,'q'), return; end

if strcmp(k,'return') | strcmp(k,'downarrow') ...
   | strcmp(k,'rightarrow') | strcmp(k,'space')
   for i= 1 : length(handles)
      current= current + 1;
      if current > length(handles), current= 1; end
      if ishandle(handles(current)), break; end
   end

elseif strcmp(k,'uparrow') | strcmp(k,'leftarrow')
   for i= 1 : length(handles)
      current= current - 1;
      if current < 1, current= length(handles); end
      if ishandle(handles(current)), break; end
   end

elseif strcmp(k,'end')
   for current= length(handles) : -1 : 1
      if ishandle(handles(current)), break; end
   end

elseif strcmp(k,'home')
   for current= 1 : length(handles)
      if ishandle(handles(current)), break; end
   end

% If the user pressed a numeric key, try to activate the figure having the
% corresponding handle number:

elseif strfind('1234567890',k)
   i= str2num(k);
   if i == 0, i= 10; end
   if ishandle(i)
      current= find(handles==i);
   else
      fprintf(2,'There is no figure having handle==%d.\n', i);
      % Play a sound to alert user that an invalid key was pressed:
      play_sound(3);
   end

elseif strcmp(k,'delete')
   close(handles(current));
   for i= 1 : length(handles)
      current= current - 1;
      if current < 1, current= length(handles); end
      if ishandle(handles(current)), break; end
   end
   if ~ishandle(handles(current))
      fprintf(2,'All figures have been deleted.\n');
      return
   end

else
   fprintf('showfigs() does not know how to handle the "%s" key.\n', k);
   play_sound(3);
end

if current > length(handles)
   fprintf(2,'The is no figure number %d.\n', current);
   play_sound(3);
   return
end

% Allow time to update:
pause(0.01);

% Make specified figure visible and raise it above all other figures:
figure(handles(current));

% Define callback function to process key presses:
set(handles(current), 'KeyPressFcn',@key_press_handler);

end % function key_press_handler(source,event)

author: Sable

repo: mcbench-benchmarks


text:
run('matconvnet/matlab/vl_setupnn.m');
addpath('util/');

outDir = 'result'; % output directory
data = 'Set5';     % test data (directory) which is in data folder
SF = 2;            % test scale factors. can be 2, 3 or 4
outRoute = fullfile(outDir, data, ['VDSR_x',num2str(SF)]);

if ~exist(outRoute, 'dir')
    mkdir(outRoute);
end

%VDSR(data, SF, 'VDSR.mat', outRoute);
VDSR(data, SF, 'VDSR_CPUonly_SupER.mat', outRoute); % MAGI-ADAPT

author: thomas-koehler

repo: SupER


text:
function outn=gt(x,y)

precAndSize

for ii=1:max(ex,ey)
 imag=false;
 if ex==1
  [xrval,xival]=getVals(x,1);
  [yrval,yival]=getVals(y,ii);
 elseif ey==1
  [xrval,xival]=getVals(x,ii);
  [yrval,yival]=getVals(y,1);
 else
  [xrval,xival]=getVals(x,ii);
  [yrval,yival]=getVals(y,ii);
 end
  outn(ii)=mpfr_gt(precision,xrval,yrval);
  if outn(ii)~=0, outn(ii)=1; end
end % for ii=1:max(ex,

author: opencobra

repo: cobratoolbox


text:
function [X,I,U,P,H] = variablesUnpack(V, N, nx, ni, nu, np)

[X_indizes, I_indizes, U_indizes, P_indizes, H_indizes] = ocl.simultaneous.indizes(N, nx, ni, nu, np);

X = reshape(V(X_indizes), nx, N+1);
I = reshape(V(I_indizes), ni, N);
U = reshape(V(U_indizes), nu, N);
P = reshape(V(P_indizes), np, N+1);
H = reshape(V(H_indizes), 1 , N);

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function v = cnn_vec(im, net)
% CNN_VEC computes a D-dimensional CNN vector for an image.
%
%   V = cnn_vec(IM, NET)
%
%   IM  : Input image, or input image path as a string
%   NET : CNN network to evaluate on image IM
%
%   V   : D-dimensional vector
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

	minsize = 67;
	net.mode = 'test';
	descvarname = 'l2descriptor';

	use_gpu = strcmp(net.device, 'gpu');
	if use_gpu
		gpuarrayfun = @(x) gpuArray(x);
		gatherfun = @(x) gather(x);
	else
		gpuarrayfun = @(x) x; % do not convert to gpuArray
		gatherfun = @(x) x; % do not gather
	end

	if isstr(im)
		im = imread(im);
	end

	im = single(im) - mean(net.meta.normalization.averageImage(:));
	if size(im, 3) == 1
		im = repmat(im, [1 1 3]);
	end
	if min(size(im, 1), size(im, 2)) < minsize
		im = pad2minsize(im, minsize, 0);
	end

	net.eval({'input', gpuarrayfun(reshape(im, [size(im), 1]))});

	v = gatherfun(squeeze(net.getVar(descvarname).value));

author: filipradenovic

repo: cnnimageretrieval


text:
function [fgsegment keepFascicles] = dtiSegmentFiberWithNiftiRoi(fg, targetROIfile, targetROI2file, thresholdmm)

% Extracting streamlines terminating within a certain distance from ROIs (e.g. gray matter ROIs)
%
% INPUT:
% fg: fg structure
% targetROIfile: the name of first ROI nifti file
% targetROI2file: the name of second ROI nifti file
% thresholdmm: the distance threshold in mm; the streamlines with endpoint within the threshold will be selected.
%
% Example:
% fg = fgRead(wholebrainFG);
% targetROIfile = 'Left_V3A.nii.gz';
% targetROI2file = 'Left_hV4.nii.gz';
% thresholdmm = 4;
%  [fgsegment keepFascicles] = dtiSegmentFiberWithNiftiRoi(fg, targetROIfile, targetROI2file, thresholdmm)
% (C) Hiromasa Takemura, CiNet, 2017

if notDefined('thresholdmm')
thresholdmm = 4;
end

% Load ROIs
roi1 = niftiRead(targetROIfile);
roi2 = niftiRead(targetROI2file);

% Define threshold in image space (assuming isotropic)
% Note that the threshold is squared (we keep using squared number in subsequent process)
threshold = (thresholdmm./(roi1.pixdim(1)))^2;

% Get acpc to image transformation matrix
acpc2img = inv(roi1.qto_xyz);

% Transfer fg into image coordinate
fgImg = dtiXformFiberCoords(fg, acpc2img,'img');

%% Select fibers within a certain distance from each ROIs
fprintf('Segmenting tracts from Connectome ...\n')

streamlinenum = length(fgImg.fibers);
% Extract ACPC coordinate of tract endpoints
for kk = 1:streamlinenum
    fibercoordinate = cell2mat(fgImg.fibers(kk));
    fiberlength = size(fibercoordinate);
        firstfiberend(kk,1) = fibercoordinate(1,1);
        firstfiberend(kk,2) = fibercoordinate(2,1);
        firstfiberend(kk,3) = fibercoordinate(3,1);
        secondfiberend(kk,1) = fibercoordinate(1,fiberlength(2));
        secondfiberend(kk,2) = fibercoordinate(2,fiberlength(2));
        secondfiberend(kk,3) = fibercoordinate(3,fiberlength(2));

    clear fibercoordinate fiberlength
end

%% Extract coordinates of ROIs
[tcoords(:,1), tcoords(:,2), tcoords(:,3)]= ind2sub(size(roi1.data), find(roi1.data));
[t2coords(:,1), t2coords(:,2), t2coords(:,3)]= ind2sub(size(roi2.data), find(roi2.data));

% Chose streamlines in which one endpoint is closer to first ROI, and the other
% endpoint is close to second ROI
for kp = 1:streamlinenum
    [indices_ff(kp), bestSqDis_ff(kp)] = nearpoints(transpose(firstfiberend(kp,:)), transpose(tcoords));
    [indices_ss(kp), bestSqDis_ss(kp)] = nearpoints(transpose(secondfiberend(kp,:)),transpose(t2coords));
    [indices_fs(kp), bestSqDis_fs(kp)] = nearpoints(transpose(firstfiberend(kp,:)), transpose(t2coords));
    [indices_sf(kp), bestSqDis_sf(kp)] = nearpoints(transpose(secondfiberend(kp,:)),transpose(tcoords));
end

% Select fibers within the threshold, then create keepFascicles
% structure
fffiber = find(bestSqDis_ff <threshold);
ssfiber = find(bestSqDis_ss <threshold);
fsfiber = find(bestSqDis_fs <threshold);
sffiber = find(bestSqDis_sf <threshold);

ffssfiber  = intersect(fffiber,ssfiber); % First streamline endpoint is near first ROI, the another endpoint is near second ROI
fssffiber = intersect(fsfiber,sffiber); % First streamline endpoint is near second ROI, the another endpoint is near first ROI
bothfiber = transpose(union(transpose(ffssfiber), transpose(fssffiber),'rows'));
bothsize = size(bothfiber);

%% If the selected number of streamlines are non-zero, extract those streamlines
if bothsize(2)>0
    keepFascicles = zeros(streamlinenum,1);
    for ik = 1:length(bothfiber)
        keepFascicles(bothfiber(ik)) = 1;
    end
fgsegment = fgExtract(fg, logical(keepFascicles), 'keep');
else
fprintf('No streamlines satisfied criteria ...\n')
end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ a, b ] = p21_lim ( )

%*****************************************************************************80
%
%% P21_LIM returns the integration limits for problem 21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function so = subStreams(s1,s2)
% SUBSTREAMS  Subtract streams
so = streamCompose({s1,s2},@subVals);

function ris = subVals(v1,v2)
ris = v1-v2;


text:
% KEYWORDs: normal log pdf
% this program returns log of f(x)=normal(x; mean_par, variance_par) pdf function and its derivative,
% the pdf f(x) is defined up to the normalization constant
% USAGE: [log_function, log_function_prime]=log_normal_pdf(x, alpha, beta)
% INPUTs: "x" is vector of argument values
%         "function_parameters" is a vector of parameters of function f(x), function_parameters(1)=mean
%                         and function_parameters(2)=variance are parameters of the normal distribution
% OUTPUTs: "log_function" is log[f(x)]
%          "log_function_prime" is (d/dx)log[f(x)]
% NOTES: pdf is log-concave
% Last modified on Jul 15, 2007

function [log_function, log_function_prime]=log_normal_pdf(x, function_parameters)

mean_par=function_parameters(1);
variance_par=function_parameters(2);

log_function=(-0.5/variance_par)*(x-mean_par).^2;
log_function_prime=(-1/variance_par)*(x-mean_par);


text:
function data = standardisedata(data,T,standardise,valid_dims)

N = length(T);

if standardise
    for i = 1:N
        t = (1:T(i)) + sum(T(1:i-1));
        if isstruct(data)
            data.X(t,:) = bsxfun(@minus,data.X(t,:),mean(data.X(t,:)));
            sdx = std(data.X(t,:));
            if any(sdx==0)
                error('At least one of the trials/segments/subjects has variance equal to zero (use cleandata4hmm?)');
            end
            data.X(t,:) = bsxfun(@rdivide,data.X(t,:),sdx);
        else
            data(t,:) = bsxfun(@minus,data(t,:),mean(data(t,:)));
            sdx = std(data(t,:));
            if any(sdx==0)
                error('At least one of the trials/segments/subjects has variance equal to zero (use cleandata4hmm?)');
            end
            data(t,:) = bsxfun(@rdivide,data(t,:),sdx);
        end
    end
else
    if nargin<4 % this is to avoid the function complaining when TUDA is used
        if isstruct(data)
            valid_dims = [1:size(data.X,2)];
        else
            valid_dims = [1:size(data,2)];
        end
    end
    for i = 1:N
        t = (1:T(i)) + sum(T(1:i-1));
        if isstruct(data)
            if any(std(data.X(t,valid_dims))==0)
                error('At least one of the trials/segments/subjects has variance equal to zero (use cleandata4hmm?)');
            end
        else
            if any(std(data(t,valid_dims))==0)
                error('At least one of the trials/segments/subjects has variance equal to zero (use cleandata4hmm?)');
            end
        end
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function p = uplus(p)
%UPLUS        Polynomial unary plus  + p
%

% written  08/27/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged, improved performance
%

repo: OpenMAS


text:
function [ vX ] = ProxLogisticLossFunctionGd( vX, vY, vC, paramLambda, numIterations, stopThr )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = ProxLogisticLossFunction( vX, vY, vC, paramLambda, numIterations, stopThr )
%   Calculates the Prox of the Logistic Cost Function using Gradient Descent.
% Input:
%   - vX            -   Input Vector.
%                       Starting point for the iterative procedure.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - vY            -   Input Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - vC            -   Model Vector.
%                       The model vector in the Logsitic Cost Function.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - numIterations -   Number of Iterations.
%                       Sets the number of iterations for the algorithm to
%                       run.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range {1, 2, ...}.
%   - stopTol       -   Stopping Condition Tolerance.
%                       Sets the stopping threshold for the L Inf (Maximum
%                       Absolute Value) of the change between 2 iterations
%                       of the algorithm.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range [0, inf).
% Output:
%   - vX            -   Output Vector.
%                       The Prox for the logistic cost function for the
%                       input vector 'vY'.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  https://math.stackexchange.com/a/3571521/33.
%   2.  Elementary Numerical Analysis MATH:3800/CS:3700(22M:072/22C:072)
%       (https://homepage.divms.uiowa.edu/~whan/3800.d/3800.html, Section 3.4).
% Remarks:
%   1.  On some cases it fails to converge (While Newton Method converge to
%       the right solution).
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     06/03/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

vXPrev = vX;
vG = zeros(size(vX, 1), 1);

hObjFun = @(vX) 0.5 * sum((vX - vY) .^ 2) + (paramLambda * log(1 + exp(-vC.' * vX)));
stepSizeMax = 2;
sSolverOptions = optimset('fminbnd');
sSolverOptions = optimset(sSolverOptions, 'Display', 'off');

for ii = 1: numIterations
    vXPrev(:) = vX;

    valExp  = exp(-vC.' * vX);
    vG(:) = vX - vY - (paramLambda * (valExp / (1 + valExp)) * vC);

    hStepSizeFun    = @(stepSize) hObjFun(vX - (stepSize * vG));
    stepSize        = fminbnd(hStepSizeFun, 0, stepSizeMax, sSolverOptions);

    vX(:) = vX - (stepSize * vG);

    if(max(abs(vX - vXPrev)) < stopThr)
        break;
    end
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function oM  = HKLOrientationMapping(varargin)

warning('Starting with MTEX 5.1 "HKLOrientationMapping" has been replaced by "ipfHKLKey".')

oM = ipfHKLKey(varargin{:});

end

author: mtex-toolbox

repo: mtex

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