text:
function ft_realtime_headlocalizer(cfg)

% FT_REALTIME_HEADLOCALIZER is a real-time application for online visualization of
% the head position for the CTF275 and the Neuromag/Elekta/Megin systems. This uses the
% continuous head localization (in CTF terminology, i.e. CHL) or position indicator
% (in Neuromag/Elekta/Megin terminology, i.e. cHPI) information.
%
% Repositioning the subject to a previous recording session can be done by specifying
% the previous dataset as cfg.template = 'subject01xxx.ds', or by pointing to a text
% file created during a previous recording; e.g. cfg.template = '29-Apr-2013-xxx.txt'.
% The latter textfile is written automatically to disk with each 'Update' buttonpress.
%
% The online visualization shows the displacement of the head relative to the start
% of the recording. The timepoint (i.e. position) relative to which the displacement
% is shown can be updated can be achieved by marking the HPI at an arbitrary moment
% by clicking the 'Update' button. This allows for repositioning within a recording
% session. Black unfilled markers should appear which indicate the positions of the
% coils at the moment of buttonpress. Distance to these marked positions then become
% colorcoded, i.e. green, orange, or red.
%
% Use as
%   ft_realtime_headlocalizer(cfg)
% with the following configuration options
%   cfg.dataset         = string, name or location of a dataset/buffer (default = 'buffer://odin:1972')
%   cfg.template        = string, name of a template dataset for between-session repositioning (default = [])
%   cfg.bufferdata      = whether to start on the 'first or 'last' data that is available (default = 'last')
%   cfg.xlim            = [min max], range in cm to plot (default = [-15 15])
%   cfg.ylim            = [min max], range in cm to plot (default = [-15 15])
%   cfg.zlim            = [min max], range in cm to plot (default is automatic)
%   cfg.blocksize       = number, size of the blocks/chuncks that are processed (default = 1 second)
%   cfg.accuracy_green  = distance from fiducial coordinate; green when within limits (default = 0.15 cm)
%   cfg.accuracy_orange = orange when within limits, red when out (default = 0.3 cm)
%   cfg.dewar           = filename or mesh, description of the dewar shape (default is automatic)
%   cfg.polhemus        = filename or mesh, description of the head shape recorded with the Polhemus (default is automatic)
%   cfg.headshape       = filename or mesh, description of the head shape recorded with the Structure Sensor
%
% The following options only apply to data from the Neuromag/Elekta/Megin system
%   cfg.headmovement    = string, name or location of the .pos file created by MaxFilter which describes the location of the head relative to the dewar
%   cfg.coilfreq        = single number in Hz or list of numbers (default = [293, 307, 314, 321, 328])
%
% This method is described in Stolk A, Todorovic A, Schoffelen JM, Oostenveld R.
% "Online and offline tools for head movement compensation in MEG."
% Neuroimage. 2013 Mar;68:39-48. doi: 10.1016/j.neuroimage.2012.11.047.

% Copyright (C) 2008-2018,  Arjen Stolk & Robert Oostenveld
% Copyright (C) 2017, Simon Homoelle
%
% 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$

agreement = {
  'By using this realtime headlocalizer tool in your research, you agree to citing the publication below.'
  ''
  'Stolk A, Todorovic A, Schoffelen JM, Oostenveld R.'
  '"Online and offline tools for head movement compensation in MEG."'
  'Neuroimage. 2013 Mar;68:39-48.'
  };

if ~strcmp(questdlg(agreement, 'User agreement', 'Yes', 'Cancel', 'Cancel'), 'Yes')
  return
end

% do the general setup of the function
ft_defaults

% set the defaults
cfg.dataset         = ft_getopt(cfg, 'dataset', 'buffer://odin:1972'); % location of the buffer/dataset
cfg.accuracy_green  = ft_getopt(cfg, 'accuracy_green',  .15); % green when within this distance from reference
cfg.accuracy_orange = ft_getopt(cfg, 'accuracy_orange',  .3); % orange when within this distance from reference
cfg.template        = ft_getopt(cfg, 'template',         []); % template dataset containing the references
cfg.blocksize       = ft_getopt(cfg, 'blocksize',         1); % in seconds
cfg.bufferdata      = ft_getopt(cfg, 'bufferdata',   'last'); % first (replay) or last (real-time)
cfg.coilfreq        = ft_getopt(cfg, 'coilfreq',   [293, 307, 314, 321, 328]); % in Hz for Neuromag
cfg.dewar           = ft_getopt(cfg, 'dewar',            []); % mesh of the dewar
cfg.headshape       = ft_getopt(cfg, 'headshape',        []); % mesh of the head with the structure sensor
cfg.polhemus        = ft_getopt(cfg, 'polhemus',         []); % mesh of the head recorded with the polhemus
cfg.headmovement    = ft_getopt(cfg, 'headmovement',     []); % maxfilter created file containing quaternions information for headlocalistation

% ensure pesistent variables are cleared
clear ft_read_header

% start by reading the header from the realtime buffer
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes'); % translate dataset into datafile+headerfile
hdr = ft_read_header(cfg.headerfile, 'cache', true, 'coordsys', 'dewar');

% for backward compatibility, can be removed end 2018
cfg = ft_checkconfig(cfg, 'renamed', {'head', 'headshape'});

% determine the size of blocks to process
blocksize   = round(cfg.blocksize * hdr.Fs);
prevSample  = 0;
count       = 0;

% determine MEG system type
isneuromag = ft_senstype(hdr.grad, 'neuromag');
isctf      = ft_senstype(hdr.grad, 'ctf275');

% this is needed to fit everything in the figure, note that the dewar coordinate systems differ
if isctf
  cfg.xlim            = ft_getopt(cfg, 'xlim', [-15 15]);
  cfg.ylim            = ft_getopt(cfg, 'ylim', [-15 15]);
  cfg.zlim            = ft_getopt(cfg, 'zlim', [-38 -8]);
elseif isneuromag
  cfg.xlim            = ft_getopt(cfg, 'xlim', [-15 15]);
  cfg.ylim            = ft_getopt(cfg, 'ylim', [-15 15]);
  cfg.zlim            = ft_getopt(cfg, 'zlim', [-25 15]);
end

if isempty(cfg.dewar)
  [v, p] = ft_version;
  if isctf
    cfg.dewar = fullfile(p, 'template', 'dewar', 'ctf.mat');
  elseif isneuromag
    cfg.dewar = fullfile(p, 'template', 'dewar', 'elekta.mat');
  end
end

if ischar(cfg.dewar) && exist(cfg.dewar, 'file')
  fprintf('reading dewar from file %s\n', cfg.dewar);
  cfg.dewar = ft_read_headshape(cfg.dewar);
end

if ischar(cfg.headshape) && exist(cfg.headshape, 'file')
  fprintf('reading headshape from file %s\n', cfg.headshape);
  cfg.headshape = ft_read_headshape(cfg.headshape);
end

if ischar(cfg.polhemus) && exist(cfg.polhemus, 'file')
  fprintf('reading polhemus data from file %s\n', cfg.polhemus);
  cfg.polhemus = ft_read_headshape(cfg.polhemus);
elseif isneuromag
  fprintf('reading polhemus data from file %s\n', cfg.dataset);
  % Neuromag/Elekta/Megin dataset will contain head shape
  cfg.polhemus = ft_read_headshape(cfg.dataset);
elseif isctf
  fprintf('reading polhemus data from file %s\n', cfg.dataset);
  % CTF dataset may contain electrode information
  elec = ft_read_sens(cfg.dataset, 'senstype', 'eeg');
  cfg.polhemus.pos  = elec.elecpos;
  cfg.polhemus.unit = elec.unit;
end

if ~isempty(cfg.headshape)
  cfg.headshape = ft_convert_units(cfg.headshape, 'cm');
end

if ~isempty(cfg.polhemus)
  cfg.polhemus = ft_convert_units(cfg.polhemus, 'cm');
end

if ~isempty(cfg.dewar)
  cfg.dewar = ft_convert_units(cfg.dewar, 'cm');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read template head position, to reposition to, if template file is specified
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isctf
  if ~isempty(cfg.template)
    [p, f, x] = fileparts(cfg.template);
    if strcmp(x, '.ds')
      shape = ft_read_headshape(cfg.template, 'coordsys', 'dewar', 'format', 'ctf_ds');
      template(1,:) = [shape.fid.pos(1,1), shape.fid.pos(1,2), shape.fid.pos(1,3)]; % chan X pos
      template(2,:) = [shape.fid.pos(2,1), shape.fid.pos(2,2), shape.fid.pos(2,3)];
      template(3,:) = [shape.fid.pos(3,1), shape.fid.pos(3,2), shape.fid.pos(3,3)];
    elseif strcmp(x, '.txt')
      template = dlmread(cfg.template);
    else
      ft_error('incorrect template file specified');
    end
  else
    template = [];
  end

  % remove CTF REF sensors, for plotting purposes
  chansel = match_str(hdr.grad.chantype, 'meggrad');
  hdr.grad.chanpos  = hdr.grad.chanpos(chansel,:);
  hdr.grad.chanori  = hdr.grad.chanori(chansel,:);
  hdr.grad.chantype = hdr.grad.chantype(chansel,:);
  hdr.grad.label    = hdr.grad.label(chansel,:);
  hdr.grad.tra      = hdr.grad.tra(chansel,:);

elseif isneuromag
  if ~isempty(cfg.template)
    [p, f, x] = fileparts(cfg.template);
    if strcmp(x, '.fif')
      shape = ft_read_headshape(cfg.template, 'coordsys', 'dewar', 'format', 'neuromag_fif');
      template(1,:) = [shape.fid.pos(1,1), shape.fid.pos(1,2), shape.fid.pos(1,3)]; % chan X pos
      template(2,:) = [shape.fid.pos(2,1), shape.fid.pos(2,2), shape.fid.pos(2,3)];
      template(3,:) = [shape.fid.pos(3,1), shape.fid.pos(3,2), shape.fid.pos(3,3)];
    elseif strcmp(x, '.txt')
      template = dlmread(cfg.template);
    end
  else
    template = [];
  end

else
  ft_error('the data does not resemble ctf, nor neuromag')
end % if ctf or neuromag

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read digitized head position (for dipole fitting)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isctf
  sens = hdr.grad;
  % not needed for CTF275 systems
  dip  = [];
  vol  = [];
  coilsignal = [];

elseif isneuromag
  shape = ft_read_headshape(cfg.headerfile, 'coordsys', 'dewar', 'format', 'neuromag_fif', 'unit', 'cm');
  for i = 1:min(size(shape.pos,1),length(cfg.coilfreq)) % for as many digitized or specified coils
    if ~isempty(strfind(shape.label{i}, 'hpi'))
      dip(i).pos = shape.pos(i,:); % chan X pos, initial guess for each of the dipole/coil positions
      dip(i).mom = [0 0 0]';
    end
  end
  if ~exist('dip', 'var')
    ft_error('head localization requires digitized positions for Neuromag systems')
  end

  % prepare the forward model and the sensor array for subsequent fitting
  % note that the forward model is a magnetic dipole in an infinite vacuum
  cfg.channel = ft_channelselection('MEGMAG', hdr.label); % old
  [vol, sens] = ft_prepare_vol_sens([], hdr.grad, 'channel', cfg.channel);
  sens = ft_datatype_sens(sens, 'scaling', 'amplitude/distance', 'distance', 'cm'); % ensure SI units
  coilsignal = [];

  % update distances, given that sensor units are m an not cm
  cfg.accuracy_green = cfg.accuracy_green/100;
  cfg.accuracy_orange = cfg.accuracy_orange/100;

else
  ft_error('the data does not resemble ctf, nor neuromag')
end % if ctf or neuromag

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define a subset of channels for reading
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isctf
  [dum, chanindx] = match_str('headloc', hdr.chantype);

elseif isneuromag
  % depending on wether movement compensation was done or only head position
  % estimation, the fif file that results from Maxfilter will contain 9 channels that
  % start with QUAT or with CHPI

  if ~isempty(cfg.headmovement)
    % load the head position information from cfg.headmovement
    tmpcfg = [];
    tmpcfg.dataset = cfg.headmovement;
    tmpcfg.channel = 'QUAT*';
    data_movement  = ft_preprocessing(tmpcfg);

    % ensure that it is regularly sampled
    tmpcfg = [];
    tmpcfg.time{1} = (1:data_movement.hdr.nSamples)/data_movement.hdr.Fs;
    data_movement  = ft_resampledata(tmpcfg, data_movement);

  elseif sum(startsWith(hdr.label, 'QUAT'))==9
    % the data only contains the estimated position, but has not been movement corrected
    chanindx = find(startsWith(hdr.label, 'QUAT'));

  elseif sum(startsWith(hdr.label, 'CHPI'))==9
    % this is movement corrected data
    chanindx = find(startsWith(hdr.label, 'CHPI'));

  else
    % select the 102 magnetometers for fitting of the HPI coils
    [dum, chanindx] = match_str('megmag', hdr.chantype);
  end
end

if isempty(chanindx)
  ft_error('the data does not seem to have head localization channels');
end

% this information is passed between the GUI callback functions
info                    = [];
info.hdr                = hdr;
info.blocksize          = blocksize;
info.isctf              = isctf;
info.isneuromag         = isneuromag;
info.cfg                = cfg;
info.template           = template;
info.sens               = sens;
info.vol                = vol;
info.dip                = dip;
info.continue           = true;
clear hdr blocksize isctf isneuromag cfg template sens vol dip

% initiate main figure
hMainFig = figure;

% attach the info in the figure
guidata(hMainFig, info);

% initiate gui controls
uicontrol_sub(hMainFig);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ishandle(hMainFig) && info.continue % while the flag is one, the loop continues

  % get the potentially updated information from the main window
  info = guidata(hMainFig);

  % determine number of samples available in buffer
  info.hdr = ft_read_header(info.cfg.headerfile, 'cache', true, 'coordsys', 'dewar');

  % see whether new samples are available
  newsamples = (info.hdr.nSamples*info.hdr.nTrials-prevSample);

  if newsamples>=info.blocksize

    if strcmp(info.cfg.bufferdata, 'last')
      begsample  = info.hdr.nSamples*info.hdr.nTrials - info.blocksize + 1;
      endsample  = info.hdr.nSamples*info.hdr.nTrials;
    elseif strcmp(info.cfg.bufferdata, 'first')
      begsample  = prevSample + 1;
      endsample  = prevSample + info.blocksize;
    else
      ft_error('unsupported value for cfg.bufferdata');
    end

    % remember up to where the data was read
    prevSample  = endsample;
    count       = count + 1;
    fprintf('processing segment %d from sample %d to %d\n', count, begsample, endsample);

    % read data segment from buffer
    dat = ft_read_data(info.cfg.datafile, 'header', info.hdr, 'begsample', begsample, 'endsample', endsample, 'chanindx', chanindx, 'checkboundary', false);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % from here onward it is specific to the head localization
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % put the data in a fieldtrip-like raw structure
    data.trial{1} = double(dat);
    data.time{1}  = offset2time(begsample, info.hdr.Fs, endsample-begsample+1);
    data.label    = info.hdr.label(chanindx);
    data.hdr      = info.hdr;
    data.fsample  = info.hdr.Fs;

    if ~isempty(info.cfg.headmovement) && info.isneuromag
      if ~all(startsWith(data.label, 'QUAT')) || ~all(startsWith(data.label, 'CHPI'))
        data.trial{1} = data_movement.trial{1}(1:size(data_movement.trial{1}),begsample:endsample);
        data.label    = data_movement.label;
      else
        fprintf('Channels for head localisation already in the .fif file, will use data form the .fif file')
      end
    end

    if info.isneuromag && size(coilsignal,2)~=info.blocksize
      % construct the reference signal for each of the coils
      % this needs to be updated if the blocksize changes
      ncoil = length(info.cfg.coilfreq);
      if ncoil==0
        ft_error('no coil frequencies were specified');
      else
        time = (1:info.blocksize)./info.hdr.Fs;
        coilsignal = zeros(ncoil, info.blocksize);
        for i=1:ncoil
          coilsignal(i,:) = exp(time*info.cfg.coilfreq(i)*1i*2*pi);
          coilsignal(i,:) = coilsignal(i,:) / norm(coilsignal(i,:));
        end
      end
    end

    % compute the HPI coil positions, this takes some time
    [hpi, info.dip] = data2hpi(data, info.dip, info.vol, info.sens, coilsignal, info.isctf, info.isneuromag); % for neuromag datasets this is relatively slow

    guidata(hMainFig, info);

    if ~ishandle(hMainFig)
      % the figure has been closed
      break
    end

    % get the potentially updated information from the main window
    info = guidata(hMainFig);
    % update the info
    info.hpi = hpi;

    % store the updated gui variables
    guidata(hMainFig, info);

    % DRAW LEFT PANEL - TOP VIEW
    a = subplot(1,2,1);
    h = get(a, 'children');
    hold on;

    if ~isempty(h)
      % done on every iteration
      delete(h);
    end

    % draw the color-coded head and distances from the templates
    draw_sub(hMainFig);

    % show current timesample
    title(sprintf('top view, runtime = %d s\n', round(mean(data.time{1}))));

    % not needed any more
    clear data;

    % viewing angle
    if info.isctf
      view(-45, 90)
    elseif info.isneuromag
      view(0, 90)
    end

    % DRAW RIGHT PANEL - FRONT/REAR VIEW
    b = subplot(1,2,2);
    i = get(b, 'children');
    hold on;

    if ~isempty(i)
      % done on every iteration
      delete(i);
    end

    % draw the color-coded head and distances from the templates
    draw_sub(hMainFig);

    % viewing angle
    if get(info.hViewMirrorButton, 'Value') == 1
      if info.isctf
        set(gca, 'Ydir', 'reverse')
        view(45, 0)
      elseif info.isneuromag
        set(gca, 'Ydir', 'reverse')
        view(0, 0)
      end
      title(sprintf('Mirror view, clock time %s', datestr(now))); % show current data & time
    else
      if info.isctf
        set(gca, 'Ydir', 'normal')
        view(135, 0)
      elseif info.isneuromag
        set(gca, 'Ydir', 'normal')
        view(180, 0)
      end
      title(sprintf('Normal view, clock time %s', datestr(now))); % show current data & time
    end

    % force Matlab to update the figure
    drawnow

  end % if enough new samples
end % while true

close(hMainFig); % close the figure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that initiates the figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function uicontrol_sub(handle, eventdata)

% get the info
info = guidata(handle);

% initiate figure
set(handle, 'KeyPressFcn', {@key_sub});

hUpdateButton = uicontrol(...
  'Parent', handle,...
  'Style', 'pushbutton',...
  'String', 'Update',...
  'Units', 'normalized',...
  'Position', [.65 .0875 .15 .075],...
  'FontSize', 12,...
  'Callback', {@update_ButtonDownFcn});

hQuitButton = uicontrol(...
  'Parent', handle,...
  'Style', 'pushbutton',...
  'String', 'Quit',...
  'Units', 'normalized',...
  'Position', [.8 .0875 .15 .075],...
  'FontSize', 12,...
  'Callback', {@quit_ButtonDownFcn});

hSphereCheckBox = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Sphere',...
  'Units', 'normalized',...
  'Position', [.05 .1 .075 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 1,...
  'Callback', {@sphere_CheckBox});

hHeadCheckBox = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Head',...
  'Units', 'normalized',...
  'Position', [.125 .1 .1 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 0,... % by default switched on
  'Callback', {@head_CheckBox});

hPolhemusCheckBox = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Polhemus',...
  'Units', 'normalized',...
  'Position', [.2 .1 .075 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 0,...
  'Callback', {@Polhemus_CheckBox});

hDewarCheckBox = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Dewar',...
  'Units', 'normalized',...
  'Position', [.275 .1 .075 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 0,...
  'Callback', {@Dewar_CheckBox});

hCoilCheckBox = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Coils',...
  'Units', 'normalized',...
  'Position', [.35 .1 .075 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 1,...
  'Callback', {@coil_CheckBox});

hSensorCheckBox = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Sensors',...
  'Units', 'normalized',...
  'Position', [.425 .1 .075 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 0,...
  'Callback', {@sensor_CheckBox});

hAxisCheckBox = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Axis',...
  'Units', 'normalized',...
  'Position', [.5 .1 .075 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 0,...
  'Callback', {@axis_CheckBox});

hBlocksizeMenu = uicontrol(...
  'Parent', handle,...
  'Style', 'popupmenu',...
  'String', {'.1 second', '.2 second', '.5 second', '1 second', '1.5 second', '2 seconds', '5 seconds', '10 seconds', '30 seconds'},...
  'Units', 'normalized',...
  'Position', [.6 .1925 .1 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 1,... % default
  'Callback', {@blocksize_Menu});

hViewMirrorButton = uicontrol(...
  'Parent', handle,...
  'Style', 'checkbox',...
  'String', 'Mirror View',...
  'Units', 'normalized',...
  'Position', [.8 .1925 .1 .05],...
  'FontSize', 8,...
  'BackgroundColor', [.8 .8 .8],...
  'Value', 1,... % by default switched on
  'Callback', {@mirror_CheckBox});

info.hQuitButton        = hQuitButton;
info.hCoilCheckBox      = hCoilCheckBox;
info.hSphereCheckBox    = hSphereCheckBox;
info.hDewarCheckBox     = hDewarCheckBox;
info.hSensorCheckBox    = hSensorCheckBox;
info.hViewMirrorButton  = hViewMirrorButton;
info.hBlocksizeMenu     = hBlocksizeMenu;
info.hRealistic         = hHeadCheckBox;
info.hPolhemusCheckBox  = hPolhemusCheckBox;
info.hAxisCheckBox      = hAxisCheckBox;

% put the info back
guidata(handle, info);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that computes the HPI coil positions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hpi, dip] = data2hpi(data, dip, vol, sens, coilsignal, isctf, isneuromag)

% The CTF275 system localizes the HPI coil positions online, and writes them
% to the dataset. For the Neuromag systems the signals evoked by the HPI coils
% are superimposed on the other signals. This requires additional online
% dipolefitting of those HPI coils.

if isctf
  % assign the channels to the resp. coil coordinates
  [dum, x1] = match_str('HLC0011', data.label);
  [dum, y1] = match_str('HLC0012', data.label);
  [dum, z1] = match_str('HLC0013', data.label);
  [dum, x2] = match_str('HLC0021', data.label);
  [dum, y2] = match_str('HLC0022', data.label);
  [dum, z2] = match_str('HLC0023', data.label);
  [dum, x3] = match_str('HLC0031', data.label);
  [dum, y3] = match_str('HLC0032', data.label);
  [dum, z3] = match_str('HLC0033', data.label);

  % convert from meter to cm and assign to the resp. coil
  hpi{1} = data.trial{1}([x1 y1 z1],end) * 100;
  hpi{2} = data.trial{1}([x2 y2 z2],end) * 100;
  hpi{3} = data.trial{1}([x3 y3 z3],end) * 100;

elseif isneuromag
  if all(startsWith(data.label, 'QUAT'))
    q1 = data.trial{1}(strcmp(data.label, 'QUAT001'),:);
    q2 = data.trial{1}(strcmp(data.label, 'QUAT002'),:);
    q3 = data.trial{1}(strcmp(data.label, 'QUAT003'),:);
    q4 = -data.trial{1}(strcmp(data.label, 'QUAT004'),:);
    q5 = -data.trial{1}(strcmp(data.label, 'QUAT005'),:);
    q6 = -data.trial{1}(strcmp(data.label, 'QUAT006'),:);
    q = [q1(1) q2(1) q3(1) q4(1) q5(1) q6(1)];

    % compute the anatomical landmark location in cm
    hpi{1} = ft_warp_apply(q(end,:), data.hdr.orig.dig(1).r' , 'quaternion')'*100;
    hpi{2} = ft_warp_apply(q(end,:), data.hdr.orig.dig(2).r' , 'quaternion')'*100;
    hpi{3} = ft_warp_apply(q(end,:), data.hdr.orig.dig(3).r' , 'quaternion')'*100;

  elseif all(startsWith(data.label, 'CHPI'))
    q1 = data.trial{1}(strcmp(data.label, 'CHPI001'),:);
    q2 = data.trial{1}(strcmp(data.label, 'CHPI002'),:);
    q3 = data.trial{1}(strcmp(data.label, 'CHPI003'),:);
    q4 = -data.trial{1}(strcmp(data.label, 'CHPI004'),:);
    q5 = -data.trial{1}(strcmp(data.label, 'CHPI005'),:);
    q6 = -data.trial{1}(strcmp(data.label, 'CHPI006'),:);
    q = [q1(1) q2(1) q3(1) q4(1) q5(1) q6(1)];

    % compute the anatomical landmark location in cm
    hpi{1} = ft_warp_apply(q(end,:), data.hdr.orig.dig(1).r' , 'quaternion')'*100;
    hpi{2} = ft_warp_apply(q(end,:), data.hdr.orig.dig(2).r' , 'quaternion')'*100;
    hpi{3} = ft_warp_apply(q(end,:), data.hdr.orig.dig(3).r' , 'quaternion')'*100;

  else
    % estimate the complex-valued MEG topography for each coil
    % this implements a discrete Fourier transform (DFT)
    topo = [];
    %[x, ut] = svdfft( data.trial{1} );
    %data.trial{1} = x;
    topo = ft_preproc_detrend(data.trial{1}) * ctranspose(coilsignal);

    % ignore the out-of-phase spectral component in the topography
    topo = real(topo); % THIS SEEMS TO BE CRUCIAL

    % fit a magnetic dipole to each of the topographies
    constr.sequential = true; % for BTI systems this would be 'false' as all coils have the same frequency
    constr.rigidbody = true;

    % fit the coils together
    dipall = [];
    ncoil = numel(dip);
    for i=1:ncoil
      dipall.pos(i,:) = dip(i).pos;
    end
    dipall = ft_inverse_dipolefit(dipall, sens, vol, topo, 'constr', constr, 'display', 'off');
    for i=1:ncoil
      sel = (1:3) + 3*(i-1);
      dip(i).pos = dipall.pos(i,:);
      dip(i).mom = real(dipall.mom(sel,i)); % ignore the complex phase information
      hpi{i} = dip(i).pos;
    end
  end
else
  ft_error('the data does not resemble ctf, nor neuromag')
end % if ctf or neuromag

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that does the timing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [time] = offset2time(offset, fsample, nsamples)

offset   = double(offset);
nsamples = double(nsamples);
time = (offset + (0:(nsamples-1)))/fsample;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which computes the circumcenter(x,y,z) of the 3D triangle (3 coils)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [cc] = circumcenter(hpi)

% use coordinates relative to point `a' of the triangle
xba = hpi{2}(1) - hpi{1}(1);
yba = hpi{2}(2) - hpi{1}(2);
zba = hpi{2}(3) - hpi{1}(3);
xca = hpi{3}(1) - hpi{1}(1);
yca = hpi{3}(2) - hpi{1}(2);
zca = hpi{3}(3) - hpi{1}(3);

% squares of lengths of the edges incident to `a'
balength = xba * xba + yba * yba + zba * zba;
calength = xca * xca + yca * yca + zca * zca;

% cross product of these edges
xcrossbc = yba * zca - yca * zba;
ycrossbc = zba * xca - zca * xba;
zcrossbc = xba * yca - xca * yba;

% calculate the denominator of the formulae
denominator = 0.5 / (xcrossbc * xcrossbc + ycrossbc * ycrossbc + zcrossbc * zcrossbc);

% calculate offset (from `a') of circumcenter
xcirca = ((balength * yca - calength * yba) * zcrossbc - (balength * zca - calength * zba) * ycrossbc) * denominator;
ycirca = ((balength * zca - calength * zba) * xcrossbc - (balength * xca - calength * xba) * zcrossbc) * denominator;
zcirca = ((balength * xca - calength * xba) * ycrossbc - (balength * yca - calength * yba) * xcrossbc) * denominator;

cc(1) = xcirca + hpi{1}(1);
cc(2) = ycirca + hpi{1}(2);
cc(3) = zcirca + hpi{1}(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which draws the color-coded head and distances to the template
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function draw_sub(handle)

% get the info
info = guidata(handle);

% compute transformation
if info.isctf
  M          = ft_headcoordinates([info.hpi{1}(1),info.hpi{1}(2),info.hpi{1}(3)], [info.hpi{2}(1),info.hpi{2}(2),info.hpi{2}(3)], [info.hpi{3}(1),info.hpi{3}(2),info.hpi{3}(3)], 'ctf');
elseif info.isneuromag
  M          = ft_headcoordinates([info.hpi{2}(1),info.hpi{2}(2),info.hpi{2}(3)], [info.hpi{1}(1),info.hpi{1}(2),info.hpi{1}(3)], [info.hpi{3}(1),info.hpi{3}(2),info.hpi{3}(3)], 'neuromag');
end
M(1:3,1:3) = inv(M(1:3,1:3));
M(1:3,4)   = (-M(1:3,4)'/M(1:3,1:3))';

% plot the HPI mismatch
if get(info.hCoilCheckBox, 'Value')
  if info.isctf
    % draw nasion position
    if ~isempty(info.template)
      if abs(info.template(1,1))-info.cfg.accuracy_green < abs(info.hpi{1}(1)) && abs(info.hpi{1}(1)) < abs(info.template(1,1))+info.cfg.accuracy_green ...
          && abs(info.template(1,2))-info.cfg.accuracy_green < abs(info.hpi{1}(2)) && abs(info.hpi{1}(2)) < abs(info.template(1,2))+info.cfg.accuracy_green ...
          && abs(info.template(1,3))-info.cfg.accuracy_green < abs(info.hpi{1}(3)) && abs(info.hpi{1}(3)) < abs(info.template(1,3))+info.cfg.accuracy_green
        plot3(info.hpi{1}(1),info.hpi{1}(2),info.hpi{1}(3), 'g^', 'MarkerFaceColor', [.5 1 .5], 'MarkerSize',25)
        head1 = true;
      elseif abs(info.template(1,1))-info.cfg.accuracy_orange < abs(info.hpi{1}(1)) && abs(info.hpi{1}(1)) < abs(info.template(1,1))+info.cfg.accuracy_orange ...
          && abs(info.template(1,2))-info.cfg.accuracy_orange < abs(info.hpi{1}(2)) && abs(info.hpi{1}(2)) < abs(info.template(1,2))+info.cfg.accuracy_orange ...
          && abs(info.template(1,3))-info.cfg.accuracy_orange < abs(info.hpi{1}(3)) && abs(info.hpi{1}(3)) < abs(info.template(1,3))+info.cfg.accuracy_orange
        plot3(info.hpi{1}(1),info.hpi{1}(2),info.hpi{1}(3), 'y^', 'MarkerFaceColor', [1 .5 0], 'MarkerEdgeColor', [1 .5 0], 'MarkerSize', 25)
        head1 = false;
      else % when not in correct position
        plot3(info.hpi{1}(1),info.hpi{1}(2), info.hpi{1}(3), 'r^', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
        head1 = false;
      end
    else
      plot3(info.hpi{1}(1),info.hpi{1}(2), info.hpi{1}(3), 'r^', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
      head1 = false;
    end

    % draw left ear position
    if ~isempty(info.template)
      if abs(info.template(2,1))-info.cfg.accuracy_green < abs(info.hpi{2}(1)) && abs(info.hpi{2}(1)) < abs(info.template(2,1))+info.cfg.accuracy_green ...
          && abs(info.template(2,2))-info.cfg.accuracy_green < abs(info.hpi{2}(2)) && abs(info.hpi{2}(2)) < abs(info.template(2,2))+info.cfg.accuracy_green ...
          && abs(info.template(2,3))-info.cfg.accuracy_green < abs(info.hpi{2}(3)) && abs(info.hpi{2}(3)) < abs(info.template(2,3))+info.cfg.accuracy_green
        plot3(info.hpi{2}(1),info.hpi{2}(2),info.hpi{2}(3), 'go', 'MarkerFaceColor', [.5 1 .5], 'MarkerSize',25)
        head2 = true;
      elseif abs(info.template(2,1))-info.cfg.accuracy_orange < abs(info.hpi{2}(1)) && abs(info.hpi{2}(1)) < abs(info.template(2,1))+info.cfg.accuracy_orange ...
          && abs(info.template(2,2))-info.cfg.accuracy_orange < abs(info.hpi{2}(2)) && abs(info.hpi{2}(2)) < abs(info.template(2,2))+info.cfg.accuracy_orange ...
          && abs(info.template(2,3))-info.cfg.accuracy_orange < abs(info.hpi{2}(3)) && abs(info.hpi{2}(3)) < abs(info.template(2,3))+info.cfg.accuracy_orange
        plot3(info.hpi{2}(1),info.hpi{2}(2),info.hpi{2}(3), 'yo', 'MarkerFaceColor', [1 .5 0], 'MarkerEdgeColor', [1 .5 0], 'MarkerSize',25)
        head2 = false;
      else % when not in correct position
        plot3(info.hpi{2}(1),info.hpi{2}(2), info.hpi{2}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
        head2 = false;
      end
    else
      plot3(info.hpi{2}(1),info.hpi{2}(2), info.hpi{2}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
      head2 = false;
    end

    % draw right ear position
    if ~isempty(info.template)
      if abs(info.template(3,1))-info.cfg.accuracy_green < abs(info.hpi{3}(1)) && abs(info.hpi{3}(1)) < abs(info.template(3,1))+info.cfg.accuracy_green  ...
          && abs(info.template(3,2))-info.cfg.accuracy_green  < abs(info.hpi{3}(2)) && abs(info.hpi{3}(2)) < abs(info.template(3,2))+info.cfg.accuracy_green  ...
          && abs(info.template(3,3))-info.cfg.accuracy_green  < abs(info.hpi{3}(3)) && abs(info.hpi{3}(3)) < abs(info.template(3,3))+info.cfg.accuracy_green
        plot3(info.hpi{3}(1),info.hpi{3}(2),info.hpi{3}(3), 'go', 'MarkerFaceColor', [.5 1 .5], 'MarkerSize',25)
        head3 = true;
      elseif abs(info.template(3,1))-info.cfg.accuracy_orange < abs(info.hpi{3}(1)) && abs(info.hpi{3}(1)) < abs(info.template(3,1))+info.cfg.accuracy_orange ...
          && abs(info.template(3,2))-info.cfg.accuracy_orange < abs(info.hpi{3}(2)) && abs(info.hpi{3}(2)) < abs(info.template(3,2))+info.cfg.accuracy_orange ...
          && abs(info.template(3,3))-info.cfg.accuracy_orange < abs(info.hpi{3}(3)) && abs(info.hpi{3}(3)) < abs(info.template(3,3))+info.cfg.accuracy_orange
        plot3(info.hpi{3}(1),info.hpi{3}(2),info.hpi{3}(3), 'yo', 'MarkerFaceColor', [1 .5 0], 'MarkerEdgeColor', [1 .5 0], 'MarkerSize',25)
        head3 = false;
      else % when not in correct position
        plot3(info.hpi{3}(1),info.hpi{3}(2), info.hpi{3}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
        head3 = false;
      end
    else
      plot3(info.hpi{3}(1),info.hpi{3}(2), info.hpi{3}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
      head3 = false;
    end
  elseif info.isneuromag
    % draw nasion position
    if ~isempty(info.template)
      if abs(info.template(1,1))-info.cfg.accuracy_green < abs(info.hpi{1}(1)) && abs(info.hpi{1}(1)) < abs(info.template(1,1))+info.cfg.accuracy_green ...
          && abs(info.template(1,2))-info.cfg.accuracy_green < abs(info.hpi{1}(2)) && abs(info.hpi{1}(2)) < abs(info.template(1,2))+info.cfg.accuracy_green ...
          && abs(info.template(1,3))-info.cfg.accuracy_green < abs(info.hpi{1}(3)) && abs(info.hpi{1}(3)) < abs(info.template(1,3))+info.cfg.accuracy_green
        plot3(info.hpi{1}(1),info.hpi{1}(2),info.hpi{1}(3), 'go', 'MarkerFaceColor', [.5 1 .5], 'MarkerSize',25)
        head1 = true;
      elseif abs(info.template(1,1))-info.cfg.accuracy_orange < abs(info.hpi{1}(1)) && abs(info.hpi{1}(1)) < abs(info.template(1,1))+info.cfg.accuracy_orange ...
          && abs(info.template(1,2))-info.cfg.accuracy_orange < abs(info.hpi{1}(2)) && abs(info.hpi{1}(2)) < abs(info.template(1,2))+info.cfg.accuracy_orange ...
          && abs(info.template(1,3))-info.cfg.accuracy_orange < abs(info.hpi{1}(3)) && abs(info.hpi{1}(3)) < abs(info.template(1,3))+info.cfg.accuracy_orange
        plot3(info.hpi{1}(1),info.hpi{1}(2),info.hpi{1}(3), 'yo', 'MarkerFaceColor', [1 .5 0], 'MarkerEdgeColor', [1 .5 0], 'MarkerSize',25)
        head1 = false;
      else % when not in correct position
        plot3(info.hpi{1}(1),info.hpi{1}(2), info.hpi{1}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
        head1 = false;
      end
    else
      plot3(info.hpi{1}(1),info.hpi{1}(2), info.hpi{1}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
      head1 = false;
    end

    % draw left ear position
    if ~isempty(info.template)
      if abs(info.template(2,1))-info.cfg.accuracy_green < abs(info.hpi{2}(1)) && abs(info.hpi{2}(1)) < abs(info.template(2,1))+info.cfg.accuracy_green ...
          && abs(info.template(2,2))-info.cfg.accuracy_green < abs(info.hpi{2}(2)) && abs(info.hpi{2}(2)) < abs(info.template(2,2))+info.cfg.accuracy_green ...
          && abs(info.template(2,3))-info.cfg.accuracy_green < abs(info.hpi{2}(3)) && abs(info.hpi{2}(3)) < abs(info.template(2,3))+info.cfg.accuracy_green
        plot3(info.hpi{2}(1),info.hpi{2}(2),info.hpi{2}(3), 'g^', 'MarkerFaceColor', [.5 1 .5], 'MarkerSize', 25)
        head2 = true;
      elseif abs(info.template(2,1))-info.cfg.accuracy_orange < abs(info.hpi{2}(1)) && abs(info.hpi{2}(1)) < abs(info.template(2,1))+info.cfg.accuracy_orange ...
          && abs(info.template(2,2))-info.cfg.accuracy_orange < abs(info.hpi{2}(2)) && abs(info.hpi{2}(2)) < abs(info.template(2,2))+info.cfg.accuracy_orange ...
          && abs(info.template(2,3))-info.cfg.accuracy_orange < abs(info.hpi{2}(3)) && abs(info.hpi{2}(3)) < abs(info.template(2,3))+info.cfg.accuracy_orange
        plot3(info.hpi{2}(1),info.hpi{2}(2),info.hpi{2}(3), 'y^', 'MarkerFaceColor', [1 .5 0], 'MarkerEdgeColor', [1 .5 0], 'MarkerSize', 25)
        head2 = false;
      else % when not in correct position
        plot3(info.hpi{2}(1),info.hpi{2}(2), info.hpi{2}(3), 'r^', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
        head2 = false;
      end
    else
      plot3(info.hpi{2}(1),info.hpi{2}(2), info.hpi{2}(3), 'r^', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
      head2 = false;
    end

    % draw right ear position
    if ~isempty(info.template)
      if abs(info.template(3,1))-info.cfg.accuracy_green < abs(info.hpi{3}(1)) && abs(info.hpi{3}(1)) < abs(info.template(3,1))+info.cfg.accuracy_green  ...
          && abs(info.template(3,2))-info.cfg.accuracy_green  < abs(info.hpi{3}(2)) && abs(info.hpi{3}(2)) < abs(info.template(3,2))+info.cfg.accuracy_green  ...
          && abs(info.template(3,3))-info.cfg.accuracy_green  < abs(info.hpi{3}(3)) && abs(info.hpi{3}(3)) < abs(info.template(3,3))+info.cfg.accuracy_green
        plot3(info.hpi{3}(1),info.hpi{3}(2),info.hpi{3}(3), 'go', 'MarkerFaceColor', [.5 1 .5], 'MarkerSize', 25)
        head3 = true;
      elseif abs(info.template(3,1))-info.cfg.accuracy_orange < abs(info.hpi{3}(1)) && abs(info.hpi{3}(1)) < abs(info.template(3,1))+info.cfg.accuracy_orange ...
          && abs(info.template(3,2))-info.cfg.accuracy_orange < abs(info.hpi{3}(2)) && abs(info.hpi{3}(2)) < abs(info.template(3,2))+info.cfg.accuracy_orange ...
          && abs(info.template(3,3))-info.cfg.accuracy_orange < abs(info.hpi{3}(3)) && abs(info.hpi{3}(3)) < abs(info.template(3,3))+info.cfg.accuracy_orange
        plot3(info.hpi{3}(1),info.hpi{3}(2),info.hpi{3}(3), 'yo', 'MarkerFaceColor', [1 .5 0], 'MarkerEdgeColor', [1 .5 0], 'MarkerSize', 25)
        head3 = false;
      else % when not in correct position
        plot3(info.hpi{3}(1),info.hpi{3}(2), info.hpi{3}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
        head3 = false;
      end
    else
      plot3(info.hpi{3}(1),info.hpi{3}(2), info.hpi{3}(3), 'ro', 'MarkerFaceColor', [1 0 0], 'MarkerSize', 25)
      head3 = false;
    end
  end
end

% plot the template fiducial positions
if ~isempty(info.template)
  if info.isctf
    plot3(info.template(1,1), info.template(1,2), info.template(1,3), 'k^', 'MarkerSize', 27, 'LineWidth', 2); % chan X pos
    plot3(info.template(2,1), info.template(2,2), info.template(2,3), 'ko', 'MarkerSize', 27, 'LineWidth', 2);
    plot3(info.template(3,1), info.template(3,2), info.template(3,3), 'ko', 'MarkerSize', 27, 'LineWidth', 2);

    text(-8,8, info.template(2,3), 'Left', 'FontSize', 15);
    text(6,-6, info.template(3,3), 'Right', 'FontSize', 15);
  elseif info.isneuromag
    plot3(info.template(1,1), info.template(1,2), info.template(1,3), 'ko', 'MarkerSize', 27, 'LineWidth', 2); % chan X pos
    plot3(info.template(2,1), info.template(2,2), info.template(2,3), 'k^', 'MarkerSize', 27, 'LineWidth', 2);
    plot3(info.template(3,1), info.template(3,2), info.template(3,3), 'ko', 'MarkerSize', 27, 'LineWidth', 2);
  end
end

% plot sphere model
if get(info.hSphereCheckBox, 'Value')
  cc = circumcenter(info.hpi);
  x_radius = sqrt((info.hpi{2}(1) - cc(1))^2 + (info.hpi{2}(2) - cc(2))^2);
  y_radius = sqrt((info.hpi{3}(1) - cc(1))^2 + (info.hpi{3}(2) - cc(2))^2);
  [xe, ye, ze] = ellipsoid(cc(1),cc(2),cc(3),x_radius,y_radius,11);
  hh = surfl(xe, ye, ze);
  shading interp
  if get(info.hCoilCheckBox, 'Value') % this only works if 'coils' are updated
    if head1 == true && head2 == true && head3 == true
      colormap cool
    else
      colormap hot
    end
  end
  alpha(.15)
end

% plot realistic head mdoel
if get(info.hRealistic, 'Value') && ~isempty(info.cfg.headshape)
  ft_plot_mesh(ft_transform_geometry(M, info.cfg.headshape))
end

% plot sensors
if get(info.hSensorCheckBox, 'Value') && ~isempty(info.sens)
  % plot the sensors
  ft_plot_sens(info.hdr.grad, 'style', 'k.');
end

% plot the HPI coil positions
for j = 1:numel(info.hpi)
  plot3(info.hpi{j}(1), info.hpi{j}(2), info.hpi{j}(3), 'ko', 'LineWidth', 1, 'MarkerSize', 5)
end

% plot the dewar
if get(info.hDewarCheckBox, 'Value') && ~isempty(info.cfg.dewar)
  ft_plot_mesh(info.cfg.dewar, 'facecolor', [0.5 0.5 0.5], 'facealpha', 0.6, 'edgecolor', 'none');
end

% plot Polhemus
if get(info.hPolhemusCheckBox, 'Value')
  if ~isempty(info.cfg.polhemus)
    ft_plot_mesh(ft_transform_geometry(M,info.cfg.polhemus), 'vertexmarker', '.')
  end
end

% plot Axis
if get(info.hAxisCheckBox, 'Value')
  if info.isctf
    ft_plot_axes([], 'coordsys', 'ctf', 'unit', 'cm');
  elseif info.isneuromag
    ft_plot_axes([], 'coordsys', 'neuromag', 'unit', 'cm');
  end
end

grid on
axis on
xlabel('x (cm)');
ylabel('y (cm)');
zlabel('z (cm)');

% place ticks at 2cm distance
set(gca, 'xtick', info.cfg.xlim(1):2:info.cfg.xlim(2))
set(gca, 'ytick', info.cfg.ylim(1):2:info.cfg.ylim(2))
set(gca, 'ztick', info.cfg.zlim(1):2:info.cfg.zlim(2))

% fix axis to avoid rescaling
axis([info.cfg.xlim info.cfg.ylim info.cfg.zlim])
axis square

% put the info back
guidata(handle, info);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which handles hot keys in the current plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function key_sub(handle, eventdata)

% get the info
info = guidata(handle);

switch eventdata.Key
  case 'u'
    % update the template positions
    fprintf('updating template coordinates \n')
    for j = 1:numel(info.hpi)
      info.template(j,:) = info.hpi{j}(:); % chan X pos
    end

    % write template position to text file for later re-positioning
    template_time = [date datestr(now, '-HH-MM-SS')];
    fprintf('writing to %s.txt \n', template_time);
    dlmwrite([template_time '.txt'], info.template, ' ');

  case 'q'
    % stop the application
    fprintf('stopping the application \n')
    info.continue = false;
  case 'c'
    % display the sensors/dewar
    if get(info.hCoilCheckBox, 'Value') == 0
      fprintf('displaying coils \n')
      set(info.hCoilCheckBox, 'Value', 1); % toggle on
    elseif get(info.hCoilCheckBox, 'Value') == 1
      set(info.hCoilCheckBox, 'Value', 0); % toggle off
    end
  case 'h'
    % display the sensors/dewar
    if get(info.hHeadCheckBox, 'Value') == 0
      fprintf('displaying head \n')
      set(info.hHeadCheckBox, 'Value', 1); % toggle on
    elseif get(info.hHeadCheckBox, 'Value') == 1
      set(info.hHeadCheckBox, 'Value', 0); % toggle off
    end
  case 's'
    % display the sensors/dewar
    if get(info.hSensorCheckBox, 'Value') == 0
      fprintf('displaying sensors/dewar \n')
      set(info.hSensorCheckBox, 'Value', 1); % toggle on
    elseif get(info.hSensorCheckBox, 'Value') == 1
      set(info.hSensorCheckBox, 'Value', 0); % toggle off
    end
  otherwise
    fprintf('no command executed \n')
end

% put the info back
guidata(handle, info);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTIONs which handle button presses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function update_ButtonDownFcn(handle, eventdata)

% get the info
info = guidata(handle);

% update the template positions
fprintf('updating template coordinates \n')
for j = 1:numel(info.hpi)
  info.template(j,:) = info.hpi{j}(:); % chan X pos
end

% write template position to text file for later re-positioning
template_time = [date datestr(now, '-HH-MM-SS')];
fprintf('writing to %s.txt \n', template_time);
dlmwrite([template_time '.txt'], info.template, ' ');

% put the info back
guidata(handle, info);

function quit_ButtonDownFcn(handle, eventdata)

% get the info
info = guidata(handle);

% stop the application
fprintf('stopping the application \n')
info.continue = false;

% put the info back
guidata(handle, info);

function coil_CheckBox(hObject, eventdata)
% toggle coils display

function sphere_CheckBox(hObject, eventdata)
% toggle head display

function sensor_CheckBox(hObject, eventdata)
% toggle sensors display

function head_CheckBox(hObject, eventdata)
% toggle head display

function Dewar_CheckBox(hObject, eventdata)
% toggle dewar display

function axis_CheckBox(hObject, eventdata)
% toggle axis display

function Polhemus_CheckBox(hObject, eventdata)
% toggle axis display

function mirror_CheckBox(handle, eventdata)
% toggle mirror display

function blocksize_Menu(handle, eventdata)

% get the info
info = guidata(handle);

val = get(info.hBlocksizeMenu, 'Value');
switch val
  case 1
    info.blocksize = round(0.1 * info.hdr.Fs); % 0.1 s
  case 2
    info.blocksize = round(0.2 * info.hdr.Fs);
  case 3
    info.blocksize = round(0.5 * info.hdr.Fs);
  case 4
    info.blocksize = round(1 * info.hdr.Fs);
  case 5
    info.blocksize = round(1.5 * info.hdr.Fs);
  case 6
    info.blocksize = round(2 * info.hdr.Fs);
  case 7
    info.blocksize = round(5 * info.hdr.Fs);
  case 8
    info.blocksize = round(10 * info.hdr.Fs);
  case 9
    info.blocksize = round(30 * info.hdr.Fs);
end
fprintf('changing blocksize to %d samples\n', info.blocksize);

% put the info back
guidata(handle, info);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function r = strgsub(s, re, replacement)
%STRGSUB Replaces all substrings matching a specified pattern
%
%   r = strgsub(s, re, newstring);
%       replaces all substrings in s that match the regular expression
%       re by the newstring.
%
%   r = strgsub(s, re, translatef);
%       translates all substrings in s that match the regular expression
%       re using the function translatef, which should be a function
%       handle.
%
%   Examples
%   --------
%       % replaces all occurrences of a particular string by a new one
%       r = strgsub('It is in year 1998 and year 2004', 'year', 'yr');
%       r <- 'It is in yr 1998 and yr 2004'
%
%       % replaces all substrings matching a pattern by a fixed string
%       r = strgsub('It is in year 1998 and year 2004', '\d+', '####');
%       r <- 'It is in year #### and year ####'
%
%       % squeezes consecutive spaces to one space
%       r = strgsub('A  b    c      d', '\s+', ' ');
%       r <- 'A b c d'
%
%       % Capitalize each word
%       r = strgsub('you and me', '\w+', @(x) [upper(x(1)), x(2:end)]);
%       r <- 'You And Me'
%
%       % increment each number occurred in a string by one
%       r = strgsub('10:20 34:562', '\d+', @(x) int2str(str2double(x)+1));
%       r <- '11:21 35:563'
%

%   History
%   -------
%       - Created by Dahua Lin, on Oct 9, 2008
%

%% parse and verify input arguments

assert(ischar(s) && ndims(s) == 2 && size(s,1) <= 1, ...
    'strgsub:invalidarg', ...
    'The first input argument should be a char string.');

assert(ischar(re) && ndims(re) == 2 && size(re,1) == 1 && ~isempty(re), ...
    'strgsub:invalidarg', ...
    'The regular expression should be a non-empty char string.');

if ischar(replacement)
    s2 = replacement;
    assert(ischar(s2) && ndims(s2) == 2 && size(s2,1) <= 1, ...
        'strgsub:invalidarg', ...
        'The size of the string of replacement is illegal.');
    direct_replace = true;

elseif isa(replacement, 'function_handle')
    tf = replacement;
    direct_replace = false;

else
    error('strgsub:invalidarg', ...
        'replacement should be either a string or a function handle.');
end

%% main

if direct_replace
    [sp, ep] = regexp(s, re, 'start', 'end');

    if ~isempty(sp)
        r = do_replace(s, sp, ep, s2);
    else
        r = s;
    end
else
    [sp, ep, ms] = regexp(s, re, 'start', 'end', 'match');

    if ~isempty(sp)
        s2 = cell(size(ms));
        for i = 1 : numel(s2)
            s2{i} = tf(ms{i});
        end
        r = do_replace(s, sp, ep, s2);
    else
        r = s;
    end
end

%% The core replace function

function r = do_replace(s, sp, ep, s2)

nr = numel(sp);
if ischar(s2)
    s2 = repmat({s2}, [1, nr]);
end

s1 = cell(1, nr + 1);
p = 1;
for i = 1 : nr
    if p < sp(i)
        s1{i} = s(p:sp(i)-1);
    end
    p = ep(i)+1;
end
if p <= length(s)
    s1{nr+1} = s(p:end);
end

ss = cell(1, 2*nr+1);
ss(1:2:end) = s1;
ss(2:2:end) = s2;

r = [ss{:}];
if isempty(r)
    r = '';
end

author: Sable

repo: mcbench-benchmarks


text:
%% Coherence-Enhancing Filtering
%
% Inspired by:
%
% * Joachim Weickert, "Coherence-Enhancing Shock Filters"
%   <http://www.mia.uni-saarland.de/Publications/weickert-dagm03.pdf>
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/coherence.py>
%

function varargout = coherence_demo_gui(im)
    % load color image
    if nargin < 1
        im = fullfile(mexopencv.root(), 'test', 'img001.jpg');
        img = cv.imread(im, 'Color',true);
    elseif ischar(im)
        img = cv.imread(im, 'Color',true);
    else
        img = im;
    end
    assert(size(img,3) == 3, 'RGB image expected');

    % create the UI
    h = buildGUI(img);
    if nargout > 0, varargout{1} = h; end
end

function img = coherence_filter(img, sigma, str_sigma, blend, niter)
    %COHERENCE_FILTER  Coherence-enhancing filter

    for i=1:niter
        gray = cv.cvtColor(img, 'RGB2GRAY');

        % dominant eigenvector
        eigen = cv.cornerEigenValsAndVecs(gray, 'BlockSize',str_sigma);
        x = eigen(:,:,3);  % x1 eigenvector of lambda_1
        y = eigen(:,:,4);  % y1 eigenvector of lambda_1

        % second order derivatives
        opts = {'KSize',sigma, 'DDepth','single'};
        gxx = cv.Sobel(gray, 'XOrder',2, 'YOrder',0, opts{:});
        gyy = cv.Sobel(gray, 'XOrder',0, 'YOrder',2, opts{:});
        gxy = cv.Sobel(gray, 'XOrder',1, 'YOrder',1, opts{:});
        gvv = x.*x.*gxx + 2*x.*y.*gxy + y.*y.*gyy;

        % dilation/erosion
        ero = cv.erode(img);
        dil = cv.dilate(img);
        if true
            img1 = cv.copyTo(dil, 'Dest',ero, 'Mask',gvv<0);
        else
            mask = repmat(gvv<0, [1 1 size(img,3)]);
            img1 = ero;
            img1(mask) = dil(mask);
        end

        % blend
        img = cv.addWeighted(img,1-blend, img1,blend, 0);
    end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    niter = round(get(h.slid(1), 'Value'));
    blend = get(h.slid(2), 'Value');
    sigma = round(get(h.slid(3), 'Value'));
    str_sigma = round(get(h.slid(4), 'Value'));
    set(h.txt(1), 'String',sprintf('Iterations: %d', niter));
    set(h.txt(2), 'String',sprintf('Blend Coeff: %.2f', blend));
    set(h.txt(3), 'String',sprintf('Integration Sigma: %d', sigma));
    set(h.txt(4), 'String',sprintf('Structure Sigma: %d', str_sigma));

    % apply coherence-enhancing filter
    out = coherence_filter(h.src, sigma*2+1, str_sigma*2+1, blend, niter);

    % show result
    set(h.img, 'CData',out);
    drawnow;
end

function h = buildGUI(img)
    %BUILDGUI  Creates the UI

    % parameters
    sz = size(img);
    sz(2) = max(sz(2), 300);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Coherence', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+105-1]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 105 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(img, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(img);
    end
    h.txt(1) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',10, 'HorizontalAlignment','left', ...
        'Position',[5 5 130 20], 'String','Iterations:');
    h.txt(2) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',10, 'HorizontalAlignment','left', ...
        'Position',[5 30 130 20], 'String','Blend Coef:');
    h.txt(3) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',10, 'HorizontalAlignment','left', ...
        'Position',[5 55 130 20], 'String','Integration Sigma:');
    h.txt(4) = uicontrol('Parent',h.fig, 'Style','text', ...
        'FontSize',10, 'HorizontalAlignment','left', ...
        'Position',[5 80 130 20], 'String','Structure Sigma:');
    h.slid(1) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Value',4, 'Min',0, 'Max',20, 'SliderStep',[1 4]./(20-0), ...
        'Position',[135 5 sz(2)-135-5 20]);
    h.slid(2) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Value',0.7, 'Min',0, 'Max',1, 'SliderStep',[0.01 0.1], ...
        'Position',[135 30 sz(2)-135-5 20]);
    h.slid(3) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Value',9, 'Min',0, 'Max',15, 'SliderStep',[1 3]./(15-0), ...
        'Position',[135 55 sz(2)-135-5 20]);
    h.slid(4) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Value',9, 'Min',0, 'Max',15, 'SliderStep',[1 3]./(15-0), ...
        'Position',[135 80 sz(2)-135-5 20]);

    % hook event handlers, and trigger default start
    set(h.slid, 'Callback',{@onChange,h}, ...
        'Interruptible','off', 'BusyAction','cancel');
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%% Copyright (C) 2016 Lagu
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @deftypefun check_assumptions (@var{x})
%% Check if the input have valid assumptions.
%%
%% Private helper function.
%%
%% @seealso{sym}
%% @end deftypefun

function check_assumptions (x)

  ca_helper (x)

end

%% Is safer use a helper when private functions call it self for classdef.
function ca_helper (x)

  persistent valid_asm

  if (isempty (valid_asm))
    valid_asm = assumptions ('possible');
  end

  if (~islogical (x))
    if (isa (x, 'char'))
      assert (ismember (x, valid_asm), ['sym: the assumption "' x '" is not supported in your Sympy version.'])
    elseif (isstruct (x))
      fields = fieldnames (x);
      for j = 1:numel (fields)
        ca_helper (fields{j})
      end
    elseif (iscell (x))
      for j = 1:length (x)
        ca_helper (x{j})
      end
    else
      error ('sym: assumption must be a string or struct or cell')
    end
  end

end

repo: octsympy


text:
function a = tt_subsubsref(obj,s)
%TT_SUBSUBSREF Helper function for tensor toolbox subsref.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if length(s) == 1
    a = obj;
else
    a = subsref(obj, s(2:end));
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function value = i4_bit_lo0 ( n )

%*****************************************************************************80
%
%% I4_BIT_LO0 returns the position of the low 0 bit base 2 in an integer.
%
%  Example:
%
%       N    Binary    Lo 0
%    ----    --------  ----
%       0           0     1
%       1           1     2
%       2          10     1
%       3          11     3
%       4         100     1
%       5         101     2
%       6         110     1
%       7         111     4
%       8        1000     1
%       9        1001     2
%      10        1010     1
%      11        1011     3
%      12        1100     1
%      13        1101     2
%      14        1110     1
%      15        1111     5
%      16       10000     1
%      17       10001     2
%    1023  1111111111     1
%    1024 10000000000     1
%    1025 10000000001     1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer to be measured.
%    N should be nonnegative.
%
%    Output, integer VALUE, the position of the low 1 bit.
%
  bit = 0;
  i = n;

  while ( 1 )

    bit = bit + 1;
    i2 = floor ( i / 2 );

    if ( i == 2 * i2 )
      break
    end

    i = i2;

  end

  value = bit;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function args = flatten_structs(args,structmask)
% Turn a list of name-value pairs/structs into a flat list of name-value pairs
% NVPs = flatten_structs(Arguments,StructMask)
%
% In:
%   Arguments : cell array of function arguments in form of name-value pairs, possibly interleaved
%               with structs {'name',value,STRUCT,'name',value,'name',value, ...}
%
%   StructMask : a bitmask that encodes at what positions Arguments contains structs; optional.
%
% Out:
%   NVPs : a cell array of name-value pairs, where the structs have been flattened into
%          name-value pairs, too
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-11-03

if nargin < 2
    structmask = cellfun('isclass',args,'struct'); end

if any(structmask)
    persistent indexcache; %#ok<TLEV>
    try
        % try to look up splice points from cache
        if length(structmask) < 62
            field = char('a'+structmask);
        else
            str = java.lang.String(char('a'+structmask));
            field = ['h' num2str(str.hashCode()+2^31)];
        end
        splicepos = indexcache.(field);
    catch %#ok<CTCH>
        % pre-calculate splice points and cache results
        splicepos = [];
        k = 1;
        while k <= length(args)
            if isstruct(args{k})
                splicepos(end+1) = k; %#ok<AGROW>
                k = k+1; % struct case
            else
                k = k+2; % NVP case
            end
        end
        splicepos = splicepos(end:-1:1);
        indexcache.(field) = splicepos;
    end

    % splice structs in
    for k = splicepos
        args = [args(1:k-1) reshape([fieldnames(args{k}) struct2cell(args{k})]',1,[]) args(k+1:end)]; end
end

repo: REST


text:
function ff=lpcar2ff(ar,np)
%LPCAR2FF LPC: Convert AR coefs to complex spectrum FF=(AR,NP)

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpcar2ff.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(ar);
if nargin<2 np=p1-1; end
ff=(rfft(ar.',2*np+2).').^(-1);

author: decouples


text:
function ndx = mysub2ind(siz,subinds)

k = [1 cumprod(siz(1:end-1))];
ndx = ones(size(subinds,1),1);
for i = 1:size(subinds,2)
    for r =1:size(subinds,1)
        v = subinds(r,i);
        ndx(r) = ndx(r) + (v-1)*k(i);
    end;
end


text:
function [iSourceRows, iRegionScouts, iVertices] = bst_convert_indices(iVertices, nComponents, GridAtlas, isSurfaceInd)
% BST_CONVERT_INDICES: Convert scout indices (in GridLoc or Vertices matrix) to indices in ImageGridAmp/ImagingKernel
%
% USAGE:  [iSourceRows, iRegionScouts, iVertices] = bst_convert_indices(iVertices, nComponents, GridAtlas, isSurfaceInd)
%
% INPUT:
%    - iVertices    : Array of vertex indices of the source space, to reference to rows in Results.GridLoc (volume) or Surface.Vertices (surface)
%                     If empty, use all the vertices
%    - nComponents  : Number of entries per vertex in SourceValues (1,2,3)
%                     If 0, the number varies, the properties of each region are defined in input GridAtlas
%    - GridAtlas    : Set of scouts that defines the properties of the source space regions, when nComponents=0
%                     GridAtlas.Scouts(i).Region(2) is the source type (V=volume, S=surface, D=dba, X=exclude)
%                     GridAtlas.Scouts(i).Region(3) is the orientation constrain (U=unconstrained, C=contrained, L=loose)
%    - isSurfaceInd : If 1, the indices iVertices are referring to Surface.Vertices, and require a conversion in the case of mixed models
%                     If 0, the indices iVertices are referring to Results.GridLoc
%
% OUTPUT:
%    - iSourceRows   : Array of vertex indices of the source space, to reference to rows in Results.GridLoc (volume)
%    - iRegionScouts : List of the scout indices in GridAtlas that are involved in the list of vertices iGridLoc
%    - iVertices     : Modified list of vertices (when some are removed, compared with the initial iVertices)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2014

% Check inputs
if (nargin ~= 4) || isempty(iVertices) || isempty(nComponents) || ((nComponents == 0) && isempty(GridAtlas))
	error('Invalid call');
end
iRegionScouts = [];

% Make sure iVertices is a row vector
iVertices = iVertices(:)';
% Get row numbers corresponding to the selected vertices
switch (nComponents)
    case 0
        % Convert indices from Surface.Vertices to Results.GridLoc
        if isSurfaceInd
            % Remove the vertices that are outside the list of vertices in Vert2Grid
            iVertices(iVertices > size(GridAtlas.Vert2Grid,2)) = [];
            % Surface.Vertices => Results.GridLoc
            iVertices = find(any(GridAtlas.Vert2Grid(:,iVertices), 2))';
        end
        % Get indices in the ImageGridAmp/ImagingKernel matrix
        iSourceRows = find(any(GridAtlas.Grid2Source(:,iVertices), 2))';
        % Find over which regions this vertex selection spans
        if (nargout >= 2)
            iRegionScouts = find(~cellfun(@(c)isempty(intersect(c,iVertices)), {GridAtlas.Scouts.GridRows}));
        end
    case 1
        iSourceRows = sort(iVertices);
    case 2
        iSourceRows = sort([2*iVertices-1, 2*iVertices]);
    case 3
        iSourceRows = sort([3*iVertices-2, 3*iVertices-1, 3*iVertices]);
end

author: brainstorm-tools


text:
classdef SOP_F19 < PROBLEM
% <single> <real> <expensive/none>
% Hartman's family

%------------------------------- Reference --------------------------------
% X. Yao, Y. Liu, and G. Lin, Evolutionary programming made faster, IEEE
% Transactions on Evolutionary Computation, 1999, 3(2): 82-102.
%------------------------------- 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)
            obj.M = 1;
            obj.D = 3;
            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)
            a = [3 10 30;0.1 10 35;3 10 30;0.1 10 35];
            c = [1;1.2;3;3.2];
            p = [0.3689 0.1170 0.2673;0.4699 0.4387 0.7470;0.1091 0.8732 0.5547;0.03815 0.5743 0.8828];
            PopObj = zeros(size(PopDec,1),1);
            for i = 1 : size(PopDec,1)
                PopObj(i) = -sum(c.*exp(-sum(a.*(repmat(PopDec(i,:),4,1)-p).^2,2)));
            end
        end
        %% Generate the minimum objective value
        function R = GetOptimum(obj,N)
            R = -3.863;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% LOC_SUBSETS - Separate channels into maximally evenly-spaced subsets.
%                 This is achieved by exchanging channels between subsets so as to
%                 increase the sum of average of distances within each channel subset.
% Usage:
%       >> subset = loc_subsets(chanlocs, nchans); % select an evenly spaced nchans
%       >> [subsets subidx pos] = loc_subsets(chanlocs, nchans, plotobj, plotchans, keepchans);
%
% Inputs:
%
%   chanlocs    - EEGLAB dataset channel locations structure (e.g., EEG.chanlocs)
%
%   nchans      - (1,N) matrix containing number of channels that should be in each
%                 of N channel subsets respectively. if total number of channels in
%                 all subsets is less than the number of EEG channels in the chanlocs
%                 structure, N+1 subsets are created. The last subset includes the
%                 remaining channels.
%
% Optional inputs:
%
%   plotobj     - (true|false|1|0) plot the time course of the objective function
%                 {default: false|0)
%   plotchans   - (true|false|1|0) make 3-D plots showing the channel locations of
%                 the subsets {default: false|0)
%   keepchans   - (cell array) channels that has to be kept in each set and
%                 not undergo optimization. You can use this option to make
%                 sure certain channels will be assigned to each set. For
%                 example, to keep channels 1:10 to the first subset and
%                 channels 20:30 to the second, use keepchans = {[1:10], [20:30]}.
%                 To only keeps channels 1:5 in the first set, use
%                 keepchans = {1:5}.  {default: {}}
%
% Outputs:
%
%   subsets     - {1,N} or {1,N+1} cell array containing channel indices for each
%                 requested channel subset.
%   subidx      - (1, EEG.nbchans) vector giving the index of the subset associated
%                 with each channel.
%   pos         - (3,N) matrix, columns containing the cartesian (x,y,z) channel
%                 positions plotted.
% Example:
%
%   % Create three sub-montages of a 256-channel montage containing 32, 60, and 100
%   % channels respectively.The 64 remaining channels will be put into a fourth subset.
%   % Also visualize time course of optimization and the channel subset head locations.
%
%   >> subset = loc_subsets(EEG.chanlocs, [32 60 100], true, true);
%
% Author: Nima Bigdely Shamlo, SCCN/INC/UCSD, 2007

% Copyright (C) 2007 Nima Bigdely Shamlo, SCCN/INC/UCSD, nima@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.

% if plotOptimization|plotSubsets in line 130 removed by nima 3/6/2007
% line 133, removed num2str  removed by nima 3/6/2007

function [subset idx pos] = loc_subsets(chanlocs, numberOfChannelsInSubset, plotOptimization, plotSubsets, mandatoryChannelsForSet);

if sum(numberOfChannelsInSubset)> length(chanlocs)
    error('Total channels in requested subsets larger than number of EEG channels.');
end
if min(numberOfChannelsInSubset) < 2
    error('Number of channels in the requested subsets must be >= 2.');
end

rand('state',0);
if nargin < 5
    mandatoryChannelsForSet = {};
end

if nargin < 4
    plotSubsets = false;
end

if nargin < 3
    plotOptimization = false;
end

pos=[cell2mat({chanlocs.X}); cell2mat({chanlocs.Y}); cell2mat({chanlocs.Z});];
dist = squareform(pdist(pos'));

nChans = length(chanlocs);
idx = ones(nChans,1);
setId = cell(1, length(numberOfChannelsInSubset)); % cell array containing channels in each set

remainingChannels = 1:nChans; % channels that are to be assigned to subsets

% channels that have to stay in their original subset (as in
% mandatoryChannelsForSet) and should not be re-assigned
allMandatoryChannels = cell2mat(mandatoryChannelsForSet);

% assign requested mandatory channels to subset
for i=1:length(mandatoryChannelsForSet)
    setId{i} = mandatoryChannelsForSet{i};
    remainingChannels(mandatoryChannelsForSet{i}) = NaN; % flag with Nan so they can be deleted later, this is to keep indexing simple
end
remainingChannels(isnan(remainingChannels)) = [];

r = remainingChannels(randperm(length(remainingChannels)));

% randomly assign remaining channels to subsets
for i=1:length(numberOfChannelsInSubset)
    numberOfChannelsTobeAddedToSubset = numberOfChannelsInSubset(i) - length(setId{i})
    setId{i} = [setId{i} r(1:numberOfChannelsTobeAddedToSubset)];
    r(1:numberOfChannelsTobeAddedToSubset) = [];
end

if length(r) > 0
    setId{length(numberOfChannelsInSubset) + 1} = r; % last set gets remaining channels
end

if plotOptimization|plotSubsets
fprintf(['Creating total of ' num2str(length(setId)) ' channel subsets:\n']);
end

if plotOptimization
    figure;
  xp = floor(sqrt(length(setId)));
  yp = ceil(length(setId)/xp);
end
    counter = 1;
exchangeHappened = true;

while exchangeHappened
    exchangeHappened = false;
    for set1 = 1:length(setId)
        for set2 = (set1 + 1):length(setId)
            for i = 1:length(setId{set1})
                for j = 1:length(setId{set2})
                    chan(1) = setId{set1}(i);
                    chan(2) = setId{set2}(j);
                    if cost_of_exchanging_channels(chan,[set1 set2], setId, dist) < 0 && ~any(ismember(chan, allMandatoryChannels))
                        setId{set1}(find(setId{set1} == chan(1))) = chan(2);
                        setId{set2}(find(setId{set2} == chan(2))) = chan(1);
                        sumDistances(counter) = 0;
                        for s = 1:length(setId)
                            sumDistances(counter) = sumDistances(counter) ...
                                 + (sum(sum(dist(setId{s},setId{s}))) / length(setId{s}));
                        end
                        if plotOptimization
                            plot(1:counter,sumDistances,'-b');
                            xlabel('number of exchanges');
                            ylabel('sum mean distances within each channel subset');
                            drawnow;
                        else
                            if mod(counter, 20) ==0
                                fprintf('number of exchanges = %d\nsum of mean distances = %g\n',...
                                           counter, sumDistances(counter));
                            end
                        end
                        counter = counter + 1;
                        exchangeHappened = true;
                    end
                end
            end
        end
    end
end

for set = 1:length(setId)
    idx(setId{set}) = set;
%    legendTitle{set} = ['subset ' num2str(set)];
end

subset = setId;

if plotSubsets
    FIG_OFFSET = 40;
    fpos = get(gcf,'position');
    figure('position',[fpos(1)+FIG_OFFSET,fpos(2)-FIG_OFFSET,fpos(3),fpos(4)]);
    scatter3(pos(1,:), pos(2,:), pos(3,:),100,idx,'fill');
    axis equal;
    %legend(legendTitle); it does not work propr
    th=title('Channel Subsets');
    %set(th,'fontsize',14)
end

if length(r)>0 && (plotOptimization|plotSubsets)
    fprintf('The last subset returned contains the %d unused channels.\n',...
                   length(r));
end

function cost = cost_of_exchanging_channels(chan, betweenSets, setId, dist);

mirrorChan(1) = 2;
mirrorChan(2) = 1;

for i=1:2
    newSetId{betweenSets(i)} = setId{betweenSets(i)};
    newSetId{betweenSets(i)}(find(newSetId{betweenSets(i)} == chan(i))) = chan(mirrorChan(i));
end

cost = 0;
for i=betweenSets
    distSumBefore{i} = sum(sum(dist(setId{i},setId{i})));
    distSumAfter{i} = sum(sum(dist(newSetId{i},newSetId{i})));
    cost = cost + (distSumBefore{i} - distSumAfter{i}) / length(setId{i});
end

%cost = (distSumAfter{1} > distSumBefore{1}) && (distSumAfter{2} > distSumBefore{2});

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [m,v,w]=gaussmixt(m1,v1,w1,m2,v2,w2)
%GAUSSMIXT Multiply two GMM pdfs
%
% Inputs: Input mixtures: k1,k2 mixtures, p dimensions
%
%   M(k1,p) = mixture means for mixture 1
%   V(k1,p) or V(p,p,k1) variances (diagonal or full)
%   W(k1,1) = mixture weights
%   M(k2,p) = mixture means for mixture 2
%   V(k2,p) or V(p,p,k2) variances (diagonal or full)
%   W(k2,1) = mixture weights
%
% Outputs:
%
%   M(k1*k2,p) = mixture means
%   V(k1*k2,p) or V(p,p,k1*k2) if p>1 and at least one input has full covariance matrix
%   W(k1*k2,1) = mixture weights
%
% See also: gaussmix, gaussmixg, gaussmixp, randvec

%      Copyright (C) Mike Brookes 2000-2012
%      Version: $Id: gaussmixt.m 5453 2014-11-19 13:10:51Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent r13 r21 r22 r31 r312 r112 r1223 r321 ch1h r122 r124
if isempty(r21)
    r13=[1 3];
    r21=[2 1];
    r22=[2 2];
    r31=[3 1];
    r112=[1 1 2];
    r122=[1 2 2];
    r124=[1 2 4];
    r312=[3 1 2];
    r321=[3 2 1];
    r1223=[1 2 2 3];
    ch1h=[-0.5; 1; -0.5];
end
[k1,p]=size(m1);
[k2,p2]=size(m2);
f1=ndims(v1)>2 || size(v1,1)>k1; % full covariance matrix is supplied
f2=ndims(v2)>2 || size(v2,1)>k2; % full covariance matrix is supplied
% ff=f1+2*f2;
if p~=p2
    error('mixtures must have the same vector dimension');
end
k=k1*k2;
j1=repmat((1:k1)',k2,1);
j2=reshape(repmat(1:k2,k1,1),k,1);
if p==1
    % display('1D vectors');
    p1=1./v1(:);
    p2=1./v2(:);
    v=1./(p1(j1)+p2(j2));
    s1=p1.*m1;
    s2=p2.*m2;
    m=v.*(s1(j1)+s2(j2));
    v12=v1(j1)+v2(j2);
    wx=-0.5*(m1(j1)-m2(j2)).^2./v12(:);
    wx=wx-max(wx); % normalize to avoid underflow
    w=w1(j1).*w2(j2).*exp(wx)./sqrt(v12(:));
    w=w/sum(w);
else
    if ~f1 && ~f2 % both diagonal covariances
        % display('both diagonal');
        p1=1./v1;
        p2=1./v2;
        v=1./(p1(j1,:)+p2(j2,:));
        s1=p1.*m1;
        s2=p2.*m2;
        m=v.*(s1(j1,:)+s2(j2,:));
        v12=v1(j1,:)+v2(j2,:);
        wx=-0.5*sum((m1(j1,:)-m2(j2,:)).^2./v12,2);
        wx=wx-max(wx); % normalize to avoid underflow
        w=w1(j1).*w2(j2).*exp(wx)./sqrt(prod(v12,2));
        w=w/sum(w);
    else % at least one full covariances
        m=zeros(k,p);
        v=zeros(p,p,k);
        w=zeros(k,1);
        wx=w;
        idp=1:p+1:p*p; % diagonal elements of p x p matrix
        if p==2                 % special code for 2D vectors
            if ~f2  % GMM 2 is diagonal
                % display('2D GMM 2 diagonal');
                p2=1./v2;
                pm2=p2.*m2;
                vx1=permute(v1,r312);
                vx1=vx1(:,r124);
                px1=vx1./repmat((vx1(:,1).*vx1(:,3)-vx1(:,2).^2),1,3); % [a b; b c] -> [c -b a]
                pm1=m1.*px1(:,r31)-m1(:,r21).*px1(:,r22);
                px=px1(j1,:);
                px(:,r31)=px(:,r31)+p2(j2,:);  % add onto diagonal elements
                vijx=vx1(j1,:);
                vijx(:,r13)=vijx(:,r13)+v2(j2,:);  % add onto diagonal elements
            elseif ~f1 % GMM 1 is diagonal
                % display('2D GMM 1 diagonal');
                p1=1./v1;
                pm1=p1.*m1;
                vx2=permute(v2,r312);
                vx2=vx2(:,r124);
                px2=vx2./repmat((vx2(:,1).*vx2(:,3)-vx2(:,2).^2),1,3); % [a b; b c] -> [c -b a]
                pm2=m2.*px2(:,r31)-m2(:,r21).*px2(:,r22);
                px=px2(j2,:);
                px(:,r31)=px(:,r31)+p1(j1,:);  % add onto diagonal elements
                vijx=vx2(j2,:);
                vijx(:,r13)=vijx(:,r13)+v1(j1,:);  % add onto diagonal elements
            else % both full covariances
                % display('2D both full');
                vx1=permute(v1,r312);
                vx1=vx1(:,r124); % make each 2 x 2 matrix into a row [a b; b c] -> [a b c]
                px1=vx1./repmat((vx1(:,1).*vx1(:,3)-vx1(:,2).^2),1,3); % [a b; b c] -> [c -b a]
                vx2=permute(v2,r312);
                vx2=vx2(:,r124);
                px2=vx2./repmat((vx2(:,1).*vx2(:,3)-vx2(:,2).^2),1,3); % [a b; b c] -> [c -b a]
                pm1=m1.*px1(:,r31)-m1(:,r21).*px1(:,r22);
                pm2=m2.*px2(:,r31)-m2(:,r21).*px2(:,r22);
                px=px1(j1,:)+px2(j2,:);
                vijx=vx1(j1,:)+vx2(j2,:);
            end
            vx=px./repmat((px(:,1).*px(:,3)-px(:,2).^2),1,3);   % divide by determinant to get inverse
            m=pm1(j1,:)+pm2(j2,:);
            m=m.*vx(:,r13)+m(:,r21).*vx(:,r22);                 % multiple by 2 x 2 matrix vx
            v=reshape(vx(:,r1223)',[2 2 k]);                    % convert vx to a 3D array of 2 x 2 matrices
            m12=m1(j1,:)-m2(j2,:);                              % subtract means to calculate weight exponent
            dij=vijx(:,1).*vijx(:,3)-vijx(:,2).^2;              % determinant of V1+V2
            wx=m12(:,r112).*m12(:,r122).*vijx(:,r321)*ch1h./dij;% exponent of weight
            w=w1(j1).*w2(j2)./sqrt(dij);                        % weight is w*exp(wx)
        else
            if ~f2  % GMM 2 is diagonal
                % display('GMM 2 diagonal');
                p2=1./v2;
                pm2=p2.*m2;
                for i=1:k1
                    v1i=v1(:,:,i);
                    p1i=inv(v1i);
                    m1i=m1(i,:);
                    pm1i=m1i*p1i;
                    w1i=w1(i);
                    ix=i;
                    for j=1:k2
                        pij=p1i;
                        pij(idp)=pij(idp)+p2(j,:);
                        vix=inv(pij);
                        vij=v1i;
                        vij(idp)=vij(idp)+v2(j,:);
                        v(:,:,ix)=vix;
                        m(ix,:)=(pm2(j,:)+pm1i)*vix;
                        m12=m2(j,:)-m1i;
                        wx(ix)=-0.5*m12/vij*m12';           % exponent of weight
                        w(ix)=w2(j)*w1i/sqrt(det(vij));     % weight is w*exp(wx)
                        ix=ix+k1;
                    end
                end
            elseif ~f1 % GMM 1 is diagonal
                % display('GMM 1 diagonal');
                p1=1./v1;
                pm1=p1.*m1;
                ix=1;
                for j=1:k2
                    v2j=v2(:,:,j);
                    p2j=inv(v2j);
                    m2j=m2(j,:);
                    pm2j=m2j*p2j;
                    w2j=w2(j);
                    for i=1:k1
                        pij=p2j;
                        pij(idp)=pij(idp)+p1(i,:);
                        vix=inv(pij);
                        vij=v2j;
                        vij(idp)=vij(idp)+v1(i,:);
                        v(:,:,ix)=vix;
                        m(ix,:)=(pm1(i,:)+pm2j)*vix;
                        m12=m1(i,:)-m2j;
                        wx(ix)=-0.5*m12/vij*m12';           % exponent of weight
                        w(ix)=w1(i)*w2j/sqrt(det(vij));     % weight is w*exp(wx)
                        ix=ix+1;
                    end
                end
            else % both full covariances
                % display('both full');
                p1=zeros(p,p,k1);
                pm1=zeros(k1,p);
                for i=1:k1
                    p1i=inv(v1(:,:,i));
                    p1(:,:,i)=p1i;
                    pm1(i,:)=m1(i,:)*p1i;
                end
                ix=1;
                for j=1:k2
                    v2j=v2(:,:,j);
                    p2j=inv(v2j);
                    m2j=m2(j,:);
                    pm2j=m2j*p2j;
                    w2j=w2(j);
                    for i=1:k1
                        pij=p1(:,:,i)+p2j;
                        vix=inv(pij);
                        v(:,:,ix)=vix;
                        vij=v1(:,:,i)+v2j;
                        m(ix,:)=(pm1(i,:)+pm2j)*vix;
                        m12=m1(i,:)-m2j;
                        wx(ix)=-0.5*m12/vij*m12';           % exponent of weight
                        w(ix)=w1(i)*w2j/sqrt(det(vij));     % weight is w*exp(wx)
                        ix=ix+1;
                    end
                end
            end

        end
        wx=wx-max(wx);              % adjust exponents to avoid underflow
        w=w.*exp(wx);               % calculate weights
        w=w/sum(w);                 % normalize weights to sum to unity
        if k==1
            v=reshape(v,size(v,1),size(v,2)); % squeeze last dimension of v if possible
        end
    end
end
if ~nargout
    if p==1
        nxx=256; % number of points to plot
        nsd=3; % number of std deviations
        sd=sqrt([v1(:);v2(:);v]);
        ma=[m1;m2;m];
        xax=linspace(min(ma-nsd*sd),max(ma+nsd*sd),nxx);
        plot(xax,gaussmixp(xax(:),m1,v1,w1),'--b');
        hold on
        plot(xax,gaussmixp(xax(:),m2,v2,w2),':r');
        plot(xax,gaussmixp(xax(:),m,v,w),'-k');
        hold off
        ylabel('Log probability density');
        legend('Mix 1','Mix 2','Product','location','best');
        axisenlarge([-1 -1 -1 -1.05]);
    elseif p==2
        nxx=128; % number of points to plot
        nsd=3;
        if f1
            s1=sqrt([v1(1:4:end)' v1(4:4:end)']); % extract diagonal elements only
        else
            s1=sqrt(v1);
        end
        if f2
            s2=sqrt([v2(1:4:end)' v2(4:4:end)']); % extract diagonal elements only
        else
            s2=sqrt(v2);
        end
        if ndims(v)>2 || size(v,1)>k
            s3=sqrt([v(1:4:end)' v(4:4:end)']); % extract diagonal elements only
        else
            s3=sqrt(v);
        end
        mal=[m1;m2;m];
        sal=[s1;s2;s3];
        xax=linspace(min(mal(:,1)-nsd*sal(:,1)),max(mal(:,1)+nsd*sal(:,1)),nxx);
        yax=linspace(min(mal(:,2)-nsd*sal(:,2)),max(mal(:,2)+nsd*sal(:,2)),nxx);
        xx(:,:,1)=repmat(xax',1,nxx);
        xx(:,:,2)=repmat(yax,nxx,1);
        xx=reshape(xx,nxx^2,2);
        subplot(2,2,1);
        imagesc(xax,yax,reshape(gaussmixp(xx,m1,v1,w1),nxx,nxx)');
        axis 'xy';
        title('Input Mix 1');
        subplot(2,2,2);
        imagesc(xax,yax,reshape(gaussmixp(xx,m2,v2,w2),nxx,nxx)');
        axis 'xy';
        title('Input Mix 2');
        subplot(2,2,3);
        imagesc(xax,yax,reshape(gaussmixp(xx,m,v,w),nxx,nxx)');
        axis 'xy';
        title('Product GMM');
    end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function []=panel6fdisp(N,n,T,Units,endo,Ymat,stringdates2,decimaldates2,Fstartlocation,Fendlocation,forecast_estimates,pref)

% preliminary task: reshape Ymat
Ymat=reshape(Ymat,T,n,N);
% preliminary task: gather in a cell the values to be plotted
% initiate the cell
plotdata={};
% because forecasts have to be computed for each unit, loop over units
for ii=1:N
% each cell entry is a matrix of actual and forecast values
% this matrix comprises 4 rows: the first row is actual data, while the three other rows are the estimates (point estimates and confidence bands) for the forecasts
% also, this matrix has a number of rows equal to the dimension of decimaldates2, which comprises the total period sample+forecasts
   % loop over variables
   for jj=1:n
      plotdata{jj,1,ii}=nan(4,size(decimaldates2,1));
      % record actual sample values
      plotdata{jj,1,ii}(1,1:T)=Ymat(:,jj,ii)';
      % copy the last point of the actual sample for the forecast part of the matrix (required to have a clean plot)
      plotdata{jj,1,ii}(:,Fstartlocation-1)=repmat(Ymat(Fstartlocation-1,jj,ii),4,1);
      % record forecast, lower bound
      plotdata{jj,1,ii}(2,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(1,:);
      % record forecast, point estimate
      plotdata{jj,1,ii}(3,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(2,:);
      % record forecast, upper bound
      plotdata{jj,1,ii}(4,Fstartlocation:Fendlocation)=forecast_estimates{jj,1,ii}(3,:);
   end
end
% then plot the figure
if pref.plot
forecast=figure('Tag','BEARresults');
set(forecast,'Color',[0.9 0.9 0.9]);
set(forecast,'name','unconditional forecasts');
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(N,n,count)
   hold on
   Xpatch=[decimaldates2(Fstartlocation-1:Fendlocation,1)' fliplr((decimaldates2(Fstartlocation-1:Fendlocation,1))')];
   Ypatch=[plotdata{jj,1,ii}(2,Fstartlocation-1:Fendlocation) fliplr(plotdata{jj,1,ii}(4,Fstartlocation-1:Fendlocation))];
   Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(Fpatch,'facealpha',0.6);
   set(Fpatch,'edgecolor','none');
   plot(decimaldates2,plotdata{jj,1,ii}(3,:),'Color',[0.4 0.4 1],'LineWidth',2);
   plot(decimaldates2,plotdata{jj,1,ii}(1,:),'Color',[0 0 0],'LineWidth',2);
   hold off
   set(gca,'XLim',[decimaldates2(1,1) decimaldates2(end,1)],'FontName','Times New Roman');
   set(gca,'XGrid','on');
   set(gca,'YGrid','on');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(Units{ii,1},'FontWeight','normal');
      end
   end
end
end

% save on Excel
% create the cell that will be saved on excel
forecastcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},size(stringdates2,1)+3,1);
horzspace=repmat({''},3,6*n);
% loop over units
for ii=1:N
% initiate the cell of results
unitcell={};
   % loop over endogenous variables (horizontal dimension)
   for jj=1:n
   % create a header
   header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''} {''} {''};{''} {''} {''} {''} {''};{''} {'actual'} {'lower bound'} {'median'} {'upper bound'}];
   % complete the cell
   endocell=[[header;stringdates2 num2cell((plotdata{jj,1,ii})')] vertspace];
   % concatenate to the previous parts of unitcell
   unitcell=[unitcell endocell];
   end
% concatenate to the previous parts of afcell
forecastcell=[forecastcell;horzspace;unitcell];
end
% trim
forecastcell=forecastcell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),forecastcell,'forecasts','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
% Weighted analytic center of a set of linear inequalities
% Joëlle Skaf - 04/29/08
%
% The weighted analytic center of a set of linear inequalities:
%           a_i^Tx <= b_i   i=1,...,m,
% is the solution of the unconstrained minimization problem
%           minimize    -sum_{i=1}^m w_i*log(b_i-a_i^Tx),
% where w_i>0

% Input data
randn('state', 0);
rand('state', 0);
n = 10;
m = 50;
tmp = randn(n,1);
A = randn(m,n);
b = A*tmp + 2*rand(m,1);
w = rand(m,1);

% Analytic center
cvx_begin
    variable x(n)
    minimize -sum(w.*log(b-A*x))
cvx_end

disp('The weighted analytic center of the set of linear inequalities is: ');
disp(x);

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function A = shuffle_by_1(A)
%SHUFFLE_BY_1 shuffle the columns of A to the right by one, such that the
%new ith column is the old i-1th column, and the new 1st column is the old
%last column.
%
% A = shuffle_by_1(A);
%
% Inputs:
%  A  a rectangular matrix
% Outputs:
%  A  the input matrix with the rows shuffled to the right by 1

end

author: odedstein

repo: sgi-introduction-course


text:
function monomial_test08 ( )

%*****************************************************************************80
%
%% MONOMIAL_TEST08 tests MONO_UNRANK_GRLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONOMIAL_TEST08\n' );
  fprintf ( 1, ...
    '  MONO_UNRANK_GRLEX is given a rank, and returns the corresponding\n' );
  fprintf ( 1, '  monomial in the sequence of all monomials in M dimensions\n' );
  fprintf ( 1, '  in grlex order.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For reference, print a monomial sequence with ranks.\n' );

  n = 4;
  m = 3;
  rank_max = mono_upto_enum ( m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Let M = %d\n', m );
  fprintf ( 1, '      N = %d\n', n );
  fprintf ( 1, '\n' );

  x = [ 0, 0, 0 ];
  i = 1;

  while ( 1 )

    fprintf ( 1, '  %3d:', i );
    for j = 1 : m
      fprintf ( 1, '  %1d', x(j) );
    end
    fprintf ( 1, '\n' );

    if ( x(1) == n )
      break
    end

    x = mono_upto_next_grlex ( m, n, x );
    i = i + 1;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now choose random ranks between 1 and %d:\n', rank_max )
  fprintf ( 1, '\n' );

  seed = 123456789;
  test_num = 5;

  for test = 1 : test_num

    [ rank, seed ] = i4_uniform_ab ( 1, rank_max, seed );
    x = mono_unrank_grlex ( m, rank );

    fprintf ( 1, '  %3d:', rank );
    for j = 1 : m
      fprintf ( 1, '  %1d', x(j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function semilogx(varargin)
%SEMILOGX     Semilogarithmic plot X vs. Y for real interval vectors X or Y
%
%Call
%
%  semilogx(Y)   or   semilogx(X,Y)   or   semilogx(X,Y,c)
%
%The x-axis is in logarithmic scale.
%
%For Y being a real interval vector, the curves (x,Y.inf) and (x,Y.sup) are plotted.
%For X being a real interval vector, the boxes X(i) x Y(i) are plotted.
%
%The area between curves or boxes is filled with color c; default is 'r' for red.
%

% written  11/15/07     S.M. Rump  (inspired by R. Malti)
%

  plot(varargin{:})
  set(gca,'xscale', 'log')

repo: OpenMAS


text:
function [SNR] = snr_mean( target, masked )
% SNR Computes signal-to-noise ratio.
%
%   S=SNR(TARGET,MASKED) returns signal-to-noise ratio
%   given target and masked speech signals.
%
%   Inputs
%           TARGET is a target signal as vector.
%
%           MASKED is a target+masker signal as vector.
%
%   Outputs
%           S is the signal-to-noise ratio (dB).
%
%   Example
%           % read target and masker signals from wav files
%           [ target, fs ] = wavread( 'sp10.wav' );
%           [ masker, fs ] = wavread( 'ssn.wav' );
%
%           % desired SNR level (dB)
%           dSNR = 5;
%
%           % generate mixture signal: noisy = signal + noise
%           [ masked, masker ] = addnoise( target, masker, dSNR );
%
%           % compute SNR (dB)
%           SNR = snr( target, masked );
%
%           % display the result
%           fprintf( 'SNR: %0.2f dB\n', SNR );
%
%   See also ADDNOISE, SEGSNR.

%   Author: Kamil Wojcicki, November 2011.

    % compute the masker (assumes additive noise model)
    masker = masked(:) - target(:);

    % compute target and masker frame energies
    energy.target = target.' * target;
    energy.masker = masker.' * masker + eps;

    % compute frame signal-to-noise ratio (dB)
    SNR = 10*log10( energy.target ./ energy.masker + eps );
    SNR
%     SNR = 123123123;

% EOF

repo: Speech-enhancement


text:
%%begin
% We can create a CollisionModel object of the human head and torso by
% running cmdl_hat. This takes a couple of seconds.
hat = cmdl_hat;

% The time is spent performing symbolic arithmetic to build the check
% functions. We can plot the collision model to see how it looks:
hat.plot;

% Let's make it a bit better
campos([1 1 0.5]);
light; lighting gouraud

% We can specify a different shoulder frame of the shoulder for the HAT
% model. The default value is a rotation of -90 deg about the x-axis.
% The indentity matrix will therefore rotate the body to be lying down,
% facing down.
hat2 = cmdl_hat(eye(4));

% Options passed to plot operate on all of the Collision Model's
% primitives
hold on
hat2.plot('FaceAlpha', 0.5);

author: Allopart

repo: rbpf-gmapping


text:
% Data File DTXY3
% Motion of a particle on a horizontal
% plane under action of elastic and
% resistance forces and force of
% dry friction
 m    =  'm';   % mass of the particle
        % Projections of forces on axis x and y
 Fx   = '-c*x-k*xt - mu*m*9.81*xt/sqrt(xt^2+yt^2)';
 Fy   = '-c*y-k*yt - mu*m*9.81*yt/sqrt(xt^2+yt^2)';
 x0   = '1';    % Initial
 y0   = '0';    % coordinates
 v0   = 'v0';   % Initial velocity
 alfa = 'alfa'; % angle between v0 and axis x
 Tend = 7;      % upper bound of integration
 eps  = 1e-10;  % desirable accuracy
 np   = 4;      % number of parameters
 P{1} = 'm';    % mass of the particle
 P{2} = 'c';    % spring stiffness
 P{3} = 'k';    % coefficient of resistance
 P{4} = 'mu';   % coefficient of dry friction


text:
function s = path_to_convert()
  % PATH_TO_CONVERT
  %
  % s = path_to_convert()
  %
  % Outputs:
  %   s path to convert executable
  %
  % See also: path_to_qslim
  %

  if ispc
    warning([ ...
      'Dear Ladislav, is there a standard place to put executables on a pc?' ...
      'Could you put convert there and change this accordingly?' ...
      'Thanks, Alec']);
    s = 'c:/prg/lib/convert/Release/convert.exe';
  elseif isunix || ismac
    % I guess this means linux
    [status, s] = system('which convert');
    s = strtrim(s);
    if isempty(s)
      guesses = { ...
        '/usr/local/bin/convert', ...
        '/opt/local/bin/convert'};
      E = cellfun(@(guess) exist(guess,'file'),guesses);
      assert(any(E),'Could not find convert');
      s = guesses{find(E,1,'first')};
    end
  end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot