text:
function [source] = ft_dipolefitting(cfg, data)

% FT_DIPOLEFITTING perform grid search and non-linear fit with one or multiple
% dipoles and try to find the location where the dipole model is best able
% to explain the measured EEG or MEG topography.
%
% This function will initially scan the whole brain with a single dipole on
% a regular coarse grid, and subsequently start at the most optimal location
% with a non-linear search. Alternatively you can specify the initial
% location of the dipole(s) and the non-linear search will start from there.
%
% Use as
%   [source] = ft_dipolefitting(cfg, data)
%
% The configuration has the following general fields
%   cfg.numdipoles  = number, default is 1
%   cfg.symmetry    = 'x', 'y' or 'z' symmetry for two dipoles, can be empty (default = [])
%   cfg.channel     = Nx1 cell-array with selection of channels (default = 'all'),
%                     see FT_CHANNELSELECTION for details
%   cfg.gridsearch  = 'yes' or 'no', perform global search for initial
%                     guess for the dipole parameters (default = 'yes')
%   cfg.nonlinear   = 'yes' or 'no', perform nonlinear search for optimal
%                     dipole parameters (default = 'yes')
%
% If you start with a grid search, the complete grid with dipole positions and
% optionally precomputed leadfields is constructed using FT_PREPARE_SOURCEMODEL. It
% can be specified as as a regular 3-D grid that is aligned with the axes of the head
% coordinate system using
%   cfg.xgrid               = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.ygrid               = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.zgrid               = vector (e.g.   0:1:20) or 'auto' (default = 'auto')
%   cfg.resolution          = number (e.g. 1 cm) for automatic grid generation
% If the source model destribes a triangulated cortical sheet, it is described as
%   cfg.sourcemodel.pos     = N*3 matrix with the vertex positions of the cortical sheet
%   cfg.sourcemodel.tri     = M*3 matrix that describes the triangles connecting the vertices
% Alternatively the position of a few dipoles at locations of interest can be
% user-specified, for example obtained from an anatomical or functional MRI
%   cfg.sourcemodel.pos     = N*3 matrix with position of each source
%   cfg.sourcemodel.inside  = N*1 vector with boolean value whether grid point is inside brain (optional)
%   cfg.sourcemodel.dim     = [Nx Ny Nz] vector with dimensions in case of 3-D grid (optional)
%
% If you do not start with a grid search, you have to give a starting location
% for the nonlinear search
%   cfg.dip.pos     = initial dipole position, matrix of Ndipoles x 3
%
% The conventional approach is to fit dipoles to event-related averages, which
% within FieldTrip can be obtained from the FT_TIMELOCKANALYSIS or from
% the FT_TIMELOCKGRANDAVERAGE function. This has the additional options
%   cfg.latency     = [begin end] in seconds or 'all' (default = 'all')
%   cfg.model       = 'moving' or 'regional'
% A moving dipole model has a different position (and orientation) for each
% timepoint, or for each component. A regional dipole model has the same
% position for each timepoint or component, and a different orientation.
%
% You can also fit dipoles to the spatial topographies of an independent
% component analysis, obtained from the FT_COMPONENTANALYSIS function.
% This has the additional options
%   cfg.component   = array with numbers (can be empty -> all)
%
% You can also fit dipoles to the spatial topographies that are present
% in the data in the frequency domain, which can be obtained using the
% FT_FREQANALYSIS function. This has the additional options
%   cfg.frequency   = single number (in Hz)
%
% Low level details of the fitting can be specified in the cfg.dipfit structure
%   cfg.dipfit.display      = level of display, can be 'off', 'iter', 'notify' or 'final' (default = 'iter')
%   cfg.dipfit.optimfun     = function to use, can be 'fminsearch' or 'fminunc' (default is determined automatic)
%   cfg.dipfit.maxiter      = maximum number of function evaluations allowed (default depends on the optimfun)
%   cfg.dipfit.checkinside  = boolean, check that the dipole remains in the source compartment (default = false)
%
% Optionally, you can modify the leadfields by reducing the rank, i.e. remove the weakest orientation
%   cfg.reducerank      = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.normalize       = 'no', 'yes' or 'column'
%   cfg.normalizeparam  = parameter for depth normalization (default = 0.5)
%   cfg.weight          = number or 1xN vector, weight for each dipole position to compensate for the size of the corresponding patch (default = 1)
%   cfg.backproject     = 'yes' (default) or 'no', in the case of a rank reduction this parameter determines whether the result will be backprojected onto the original subspace
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions can be present in the data or can be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_SOURCEANALYSIS, FT_PREPARE_LEADFIELD, FT_PREPARE_HEADMODEL

% TODO change the output format, more suitable would be something like:
% dip.label
% dip.time
% dip.avg (instead of Vdata)
% dip.dip.pos
% dip.dip.mom
% dip.dip.model, or dip.dip.avg
% dip.dimord

% Undocumented local options:
%   cfg.dipfit.constr   = Source model constraints, depends on cfg.symmetry
% Optionally, you can include a noise covariance structure to sphere the data (is useful when using both
% magnetometers and gradiometers to fit your dipole)
%   cfg.dipfit.noisecov       = noise covariance matrix, see e.g. FT_TIMELOCK_ANALYSIS

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

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

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

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

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', {'comp', 'timelock', 'freq'}, 'feedback', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed', {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed', {'optofile', 'opto'});
cfg = ft_checkconfig(cfg, 'renamed', {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'vol',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'grid',    'sourcemodel'});

% get the defaults
cfg.channel         = ft_getopt(cfg, 'channel', 'all');
cfg.component       = ft_getopt(cfg, 'component');        % for comp input
cfg.frequency       = ft_getopt(cfg, 'frequency');        % for freq input
cfg.latency         = ft_getopt(cfg, 'latency', 'all');   % for timeclock input
cfg.feedback        = ft_getopt(cfg, 'feedback', 'text');
cfg.gridsearch      = ft_getopt(cfg, 'gridsearch', 'yes');
cfg.nonlinear       = ft_getopt(cfg, 'nonlinear', 'yes');
cfg.symmetry        = ft_getopt(cfg, 'symmetry');
cfg.dipfit          = ft_getopt(cfg, 'dipfit', []);     % the default for this is handled below
% the following options are for on-the-fly leadfield computation
cfg.reducerank      = ft_getopt(cfg, 'reducerank', []); % the default for this is handled below
cfg.normalize       = ft_getopt(cfg, 'normalize');      % this is better not used in dipole fitting
cfg.normalizeparam  = ft_getopt(cfg, 'normalizeparam'); % this is better not used in dipole fitting
cfg.backproject     = ft_getopt(cfg, 'backproject');    % this is better not used in dipole fitting
cfg.weight          = ft_getopt(cfg, 'weight');         % this is better not used in dipole fitting

cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'}); % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

% the default for this depends on the data type
if ~isfield(cfg, 'model')
  if ~isempty(cfg.component)
    % each component is fitted independently
    cfg.model = 'moving';
  elseif ~isempty(cfg.frequency)
    % fit the data with a dipole at one location
    cfg.model = 'regional';
  elseif ~isempty(cfg.latency)
    % fit the data with a dipole at one location
    cfg.model = 'regional';
  end
end

if ~isfield(cfg, 'numdipoles')
  if isfield(cfg, 'dip')
    cfg.numdipoles = size(cfg.dip(1).pos,1);
  else
    cfg.numdipoles = 1;
  end
end

% set up the symmetry constraints
if ~isempty(cfg.symmetry)
  if cfg.numdipoles~=2
    ft_error('symmetry constraints are only supported for two-dipole models');
  elseif strcmp(cfg.symmetry, 'x')
    % this structure is passed onto the low-level ft_dipole_fit function
    cfg.dipfit.constr.reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    cfg.dipfit.constr.expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    cfg.dipfit.constr.mirror = [1 1 1 -1 1 1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 -x1 y1 z1]
  elseif strcmp(cfg.symmetry, 'y')
    % this structure is passed onto the low-level ft_dipole_fit function
    cfg.dipfit.constr.reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    cfg.dipfit.constr.expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    cfg.dipfit.constr.mirror = [1 1 1 1 -1 1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 -y1 z1]
  elseif strcmp(cfg.symmetry, 'z')
    % this structure is passed onto the low-level ft_dipole_fit function
    cfg.dipfit.constr.reduce = [1 2 3];         % select the parameters [x1 y1 z1]
    cfg.dipfit.constr.expand = [1 2 3 1 2 3];   % repeat them as [x1 y1 z1 x1 y1 z1]
    cfg.dipfit.constr.mirror = [1 1 1 1 1 -1];  % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 y1 -z1]
  else
    ft_error('unrecognized symmetry constraint');
  end
elseif ~isfield(cfg, 'dipfit') || ~isfield(cfg.dipfit, 'constr')
  % no symmetry constraints have been specified
  cfg.dipfit.constr = [];
end

if ft_getopt(cfg.dipfit.constr, 'sequential', false) && strcmp(cfg.model, 'moving')
  ft_error('the moving dipole model does not combine with the sequential constraint')
  % see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3119
end

if isfield(data, 'topolabel')
  % this looks like a component analysis
  iscomp = 1;
  % transform the data into a representation on which the timelocked dipole fit can perform its trick
  data = comp2timelock(cfg, data);
else
  iscomp = 0;
end

if isfield(data, 'freq')
  % this looks like a frequency analysis
  isfreq = 1;
  % transform the data into a representation on which the timelocked dipole fit can perform its trick
  data = freq2timelock(cfg, data);
else
  isfreq = 0;
end

% prepare the volume conduction model and the sensor array
% this updates the configuration with the appropriate fields
[headmodel, sens, cfg] = prepare_headmodel(cfg, data);

% set the default for reducing the rank of the leadfields
if isempty(cfg.reducerank)
  if ft_senstype(sens, 'eeg')
    cfg.reducerank = 'no';    % for EEG
  elseif ft_senstype(sens, 'meg') && ft_headmodeltype(headmodel, 'infinite')
    cfg.reducerank = 'no';    % for MEG with a magnetic dipole, e.g. a HPI coil
  elseif ft_senstype(sens, 'meg')
    cfg.reducerank = 'yes';   % for MEG with a current dipole in a volume conductor
  end
end

% construct the options for the leadfield computation, the same options are also passed to DIPOLE_FIT
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     cfg.reducerank);
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      cfg.normalize);
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', cfg.normalizeparam);
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         cfg.weight);
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    cfg.backproject);

% select the desired channels, ordered according to the sensor structure
[selsens, seldata] = match_str(sens.label, data.label);
% take the selected channels from the data structure
Vdata = data.avg(seldata, :);

% sphere the date using the noise covariance matrix supplied, if any
% this affects both the gridsearch and the nonlinear optimization
noisecov = ft_getopt(cfg.dipfit, 'noisecov');
if ~isempty(noisecov)
  [u, s] = svd(noisecov);
  tol = max(size(noisecov)) * eps(norm(s, inf));
  s = diag(s);
  r1 = sum(s > tol) + 1;
  s(1:(r1 - 1)) = 1 ./ sqrt(s(1:(r1 - 1)));
  s(r1:end)     = 0;
  sphere = diag(s) * u';
  % apply the sphering to the data
  Vdata = sphere * Vdata;
  % apply the sphering as a pre-multiplication to the sensor definition
  montage = [];
  montage.labelold = cfg.channel;
  montage.labelnew = cfg.channel;
  montage.tra = sphere;
  sens = ft_apply_montage(sens, montage, 'balancename', 'sphering');
end

if iscomp
  % select the desired component topographies
  Vdata = Vdata(:, cfg.component);
elseif isfreq
  % the desired frequencies have already been selected
  Vdata = Vdata(:, :);
else
  % select the desired latencies
  if ischar(cfg.latency) && strcmp(cfg.latency, 'all')
    cfg.latency = data.time([1 end]);
  end
  tbeg = nearest(data.time, cfg.latency(1));
  tend = nearest(data.time, cfg.latency(end));
  cfg.latency = [data.time(tbeg) data.time(tend)];
  Vdata = Vdata(:, tbeg:tend);
end

nchans = size(Vdata,1);
ntime  = size(Vdata,2);
Vmodel = zeros(nchans, ntime);
fprintf('selected %d channels\n', nchans);
fprintf('selected %d topographies\n', ntime);

if nchans<cfg.numdipoles*3
  ft_warning('not enough channels to perform a dipole fit');
end

if ntime<1
  ft_error('no spatial topography selected');
end

% check whether EEG is average referenced
if ft_senstype(sens, 'eeg')
  if any(rv(Vdata, avgref(Vdata))>0.001)
    ft_warning('the EEG data is not average referenced, correcting this');
  end
  Vdata = avgref(Vdata);
end

% set to zeros if no initial dipole was specified
if ~isfield(cfg, 'dip')
  cfg.dip.pos = zeros(cfg.numdipoles, 3);
  cfg.dip.mom = zeros(3*cfg.numdipoles, 1);
end

% set to zeros if no initial dipole position was specified
if ~isfield(cfg.dip, 'pos')
  cfg.dip.pos = zeros(cfg.numdipoles, 3);
end

% set to zeros if no initial dipole moment was specified
if ~isfield(cfg.dip, 'mom')
  cfg.dip.mom = zeros(3*cfg.numdipoles, 1);
end

% check the specified dipole model
if numel(cfg.dip.pos)~=cfg.numdipoles*3 || numel(cfg.dip.mom)~=cfg.numdipoles*3
  ft_error('inconsistent number of dipoles in configuration')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform the dipole scan, this is usefull for generating an initial guess
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(cfg.gridsearch, 'yes')
  % test whether we have a valid configuration for dipole scanning
  if cfg.numdipoles==1
    % this is ok
  elseif cfg.numdipoles==2 && ~isempty(cfg.dipfit.constr)
    % this is also ok
  elseif isfield(cfg.sourcemodel, 'pos') && size(cfg.sourcemodel.pos,2)==cfg.numdipoles*3
    % this is also ok
  else
    ft_error('dipole scanning is only possible for a single dipole or a symmetric dipole pair');
  end

  % copy all options that are potentially used in ft_prepare_sourcemodel
  tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo'});
  tmpcfg.headmodel = headmodel;
  if ft_senstype(sens, 'eeg')
    tmpcfg.elec = sens;
  elseif ft_senstype(sens, 'meg')
    tmpcfg.grad = sens;
  end
  % construct the dipole grid on which the gridsearch will be done
  sourcemodel = ft_prepare_sourcemodel(tmpcfg);

  ngrid = size(sourcemodel.pos,1);

  switch cfg.model
    case 'regional'
      sourcemodel.error = nan(ngrid, 1);
    case 'moving'
      sourcemodel.error = nan(ngrid, ntime);
    otherwise
      ft_error('unsupported cfg.model');
  end

  insideindx = find(sourcemodel.inside);
  ft_progress('init', cfg.feedback, 'scanning grid');
  for i=1:length(insideindx)
    ft_progress(i/length(insideindx), 'scanning grid location %d/%d\n', i, length(insideindx));
    thisindx = insideindx(i);
    if isfield(sourcemodel, 'leadfield')
      % reuse the previously computed leadfield
      lf = sourcemodel.leadfield{thisindx};
    else
      lf = ft_compute_leadfield(sourcemodel.pos(thisindx,:), sens, headmodel, leadfieldopt{:});
    end
    % the model is V=lf*mom+noise, therefore mom=pinv(lf)*V estimates the
    % dipole moment this makes the model potential U=lf*pinv(lf)*V and the
    % model error is norm(V-U) = norm(V-lf*pinv(lf)*V) = norm((eye-lf*pinv(lf))*V)
    if any(isnan(lf(:)))
      % this might happen if one of the dipole locations of the grid is
      % outside the brain compartment
      lf(:) = 0;
    end
    switch cfg.model
      case 'regional'
        % sum the error over all latencies
        sourcemodel.error(thisindx,1) = sum(sum(((eye(nchans)-lf*pinv(lf))*Vdata).^2));
      case 'moving'
        % remember the error for each latency independently
        sourcemodel.error(thisindx,:) = sum(((eye(nchans)-lf*pinv(lf))*Vdata).^2);
      otherwise
        ft_error('unsupported cfg.model');
    end % switch model
  end % looping over the grid
  ft_progress('close');

  switch cfg.model
    case 'regional'
      % find the source position with the minimum error
      [err, indx] = min(sourcemodel.error);
      dip.pos = sourcemodel.pos(indx,:);                % note that for a symmetric dipole pair this results in a vector
      dip.pos = reshape(dip.pos,3,cfg.numdipoles)';     % convert to a Nx3 array
      dip.mom = zeros(cfg.numdipoles*3,1);              % set the dipole moment to zero
      if cfg.numdipoles==1
        fprintf('found minimum after scanning on grid point [%g %g %g]\n', dip.pos(1), dip.pos(2), dip.pos(3));
      elseif cfg.numdipoles==2
        fprintf('found minimum after scanning on grid point [%g %g %g; %g %g %g]\n', dip.pos(1), dip.pos(2), dip.pos(3), dip.pos(4), dip.pos(5), dip.pos(6));
      end

    case 'moving'
      for t=1:ntime
        % find the source position with the minimum error
        [err, indx] = min(sourcemodel.error(:,t));
        dip(t).pos = sourcemodel.pos(indx,:);                        % note that for a symmetric dipole pair this results in a vector
        dip(t).pos = reshape(dip(t).pos,3,cfg.numdipoles)';   % convert to a Nx3 array
        dip(t).mom = zeros(cfg.numdipoles*3,1);               % set the dipole moment to zero
        if cfg.numdipoles==1
          fprintf('found minimum after scanning for topography %d on grid point [%g %g %g]\n', t, dip(t).pos(1), dip(t).pos(2), dip(t).pos(3));
        elseif cfg.numdipoles==2
          fprintf('found minimum after scanning for topography %d on grid point [%g %g %g; %g %g %g]\n', t, dip(t).pos(1), dip(t).pos(2), dip(t).pos(3), dip(t).pos(4), dip(t).pos(5), dip(t).pos(6));
        end
      end

    otherwise
      ft_error('unsupported cfg.model');
  end % switch model

elseif strcmp(cfg.gridsearch, 'no')
  % use the initial guess supplied in the configuration for the remainder
  switch cfg.model
    case 'regional'
      dip = struct(cfg.dip);      % ensure that it is a struct, not a config object
    case 'moving'
      for t=1:ntime
        dip(t) = struct(cfg.dip); % ensure that it is a struct, not a config object
      end
    otherwise
      ft_error('unsupported cfg.model');
  end % switch model

end % if gridsearch yes/no
% multiple dipoles can be represented either as a 1x(N*3) vector or as a Nx3 matrix,
% i.e. [x1 y1 z1 x2 y2 z2] or [x1 y1 z1; x2 y2 z2]
switch cfg.model
  case 'regional'
    dip = fixdipole(dip);
  case 'moving'
    for t=1:ntime
      dip(t) = fixdipole(dip(t));
    end
  otherwise
    ft_error('unsupported cfg.model');
end % switch model

% convert the structure with the additional low-level options into key-value pairs
dipfitopt = ft_cfg2keyval(cfg.dipfit);

% add the options for the leadfield computation
dipfitopt = ft_setopt(dipfitopt, 'reducerank',     cfg.reducerank);
dipfitopt = ft_setopt(dipfitopt, 'normalize',      cfg.normalize);
dipfitopt = ft_setopt(dipfitopt, 'normalizeparam', cfg.normalizeparam);
dipfitopt = ft_setopt(dipfitopt, 'weight',         cfg.weight);
dipfitopt = ft_setopt(dipfitopt, 'backproject',    cfg.backproject);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform the non-linear fit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(cfg.nonlinear, 'yes')
  switch cfg.model
    case 'regional'
      % perform the non-linear dipole fit for all latencies together
      % catch errors due to non-convergence
      try
        dip = dipole_fit(dip, sens, headmodel, Vdata, dipfitopt{:});
        success = 1;
        if cfg.numdipoles==1
          fprintf('found minimum after non-linear optimization on [%g %g %g]\n', dip.pos(1), dip.pos(2), dip.pos(3));
        elseif cfg.numdipoles==2
          fprintf('found minimum after non-linear optimization on [%g %g %g; %g %g %g]\n', dip.pos(1,1), dip.pos(1,2), dip.pos(1,3), dip.pos(2,1), dip.pos(2,2), dip.pos(2,3));
        end
      catch
        success = 0;
        disp(lasterr);
      end

    case 'moving'
      % perform the non-linear dipole fit for each latency independently
      % instead of using dip(t) = dipole_fit(dip(t),...), I am using temporary variables dipin and dipout
      % to prevent errors like "Subscripted assignment between dissimilar structures"
      dipin = dip;
      for t=1:ntime
        % catch errors due to non-convergence
        try
          dipout(t) = dipole_fit(dipin(t), sens, headmodel, Vdata(:,t), dipfitopt{:});
          success(t) = 1;
          if cfg.numdipoles==1
            fprintf('found minimum after non-linear optimization for topography %d on [%g %g %g]\n', t, dipout(t).pos(1), dipout(t).pos(2), dipout(t).pos(3));
          elseif cfg.numdipoles==2
            fprintf('found minimum after non-linear optimization for topography %d on [%g %g %g; %g %g %g]\n', t, dipout(t).pos(1,1), dipout(t).pos(1,2), dipout(t).pos(1,3), dipout(t).pos(2,1), dipout(t).pos(2,2), dipout(t).pos(2,3));
          end
        catch
          % keep the position and moment according to the initial guess
          dipout(t).pos = dipin(t).pos;
          dipout(t).mom = dipin(t).mom;
          success(t) = 0;
          disp(lasterr);
        end
      end
      dip = dipout;
      clear dipin dipout
    otherwise
      ft_error('unsupported cfg.model');
  end % switch model
end % if nonlinear

if strcmp(cfg.nonlinear, 'no')
  % the optimal dipole positions are either obtained from scanning
  % or from the initial configured specified by the user
  switch cfg.model
    case 'regional'
      success = 1;
    case 'moving'
      success = ones(1,ntime);
    otherwise
      ft_error('unsupported cfg.model');

  end % switch model
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the model potential distribution and the residual variance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch cfg.model
  case 'regional'
    if success
      % re-compute the leadfield in order to compute the model potential and dipole moment
      lf = ft_compute_leadfield(dip.pos, sens, headmodel, leadfieldopt{:});
      if isfield(dip, 'mom') && isfield(dip, 'ampl')
        % the orientation and amplitude have already been estimated, this applies to the case of a fixed dipole orientation
        dip.pot = (lf * dip.mom) * dip.ampl;
      else
        % compute all details of the final dipole model using linear estimation
        dip.mom = pinv(lf)*Vdata;
        dip.pot = lf*dip.mom;
      end
      dip.rv  = rv(Vdata, dip.pot);
      Vmodel  = dip.pot;
    end
  case 'moving'
    for t=1:ntime
      if success(t)
        % re-compute the leadfield in order to compute the model potential and dipole moment
        lf = ft_compute_leadfield(dip(t).pos, sens, headmodel, leadfieldopt{:});
        % compute all details of the final dipole model
        dip(t).mom = pinv(lf)*Vdata(:,t);
        dip(t).pot = lf*dip(t).mom;
        dip(t).rv  = rv(Vdata(:,t), dip(t).pot);
        Vmodel(:,t) = dip(t).pot;
      end
    end
  otherwise
    ft_error('unsupported cfg.model');
end % switch model

switch cfg.model
  case 'regional'
    if isfreq
      % the matrix with the dipole moment is encrypted and cannot be interpreted straight away
      % reconstruct the frequency representation of the data at the source level
      if isfield(dip, 'mom') && isfield(dip, 'ampl')
        % this applies to the case of a fixed dipole orientation
        [dip.pow, dip.csd, dip.fourier] = timelock2freq(dip.mom * dip.ampl);
      else
        [dip.pow, dip.csd, dip.fourier] = timelock2freq(dip.mom);
      end
    end
  case 'moving'
    if isfreq
      % although this is technically possible so far, it does not make any sense
      ft_warning('a moving dipole model in the frequency domain is not supported');
    end
  otherwise
    ft_error('unsupported cfg.model');
end % switch model

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect the results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
source.label  = cfg.channel; % these channels were used in fitting
source.dip    = dip;
source.Vdata  = Vdata;  % FIXME this should be renamed (if possible w.r.t. EEGLAB)
source.Vmodel = Vmodel; % FIXME this should be renamed (if possible w.r.t. EEGLAB)

% the units of the fitted source are the same as the units of the headmodel and the sensor array
for i=1:length(source.dip)
  source.dip(i).unit = headmodel.unit;
end

% assign a latency, frequeny or component axis to the output
if iscomp
  source.component = cfg.component;
  % FIXME assign Vdata to an output variable, idem for the model potential
elseif isfreq
  source.freq   = cfg.frequency;
  source.dimord = 'chan_freq';
  % FIXME assign Vdata to an output variable, idem for the model potential
else
  tbeg = nearest(data.time, cfg.latency(1));
  tend = nearest(data.time, cfg.latency(end));
  source.time   = data.time(tbeg:tend);
  source.dimord = 'chan_time';
end

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

author: spm


text:
function cgs_test ( )

%*****************************************************************************80
%
%% CGS_TEST tests CGS.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CGS_TEST:\n' );
  fprintf ( 1, '  CGS uses the Conjugate Gradient Squared \n' );
  fprintf ( 1, '  iterative method to approximate the solution \n' );
  fprintf ( 1, '  of a linear system A * x = b.\n' );

  n = 10;

  A = zeros ( n, n );

  for i = 1 : n
    A(i,i) = 2.0;
  end
  for i = 1 : n-1
    A(i,i+1) = -1;
  end
  for i = 2 : n
    A(i-1,i) = -1;
  end
  x = [ 1 : n ]';
  b = A * x;
  x = ones ( n, 1 );

  M = 0;
  max_it = 10;
  tol = 0.0001;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this example, the order of the system is N = %d\n', n );
  fprintf ( 1, '  The matrix A is the simple tridiagonal -1, 2, -1.\n' );
  fprintf ( 1, '  The correct solution is x = [ 1, 2, ..., n].\n' );
  fprintf ( 1, '  The right hand side b is determined by computing A * x.\n' );
  fprintf ( 1, '  The exact x is then replaced by a vector of all 1''s for\n' );
  fprintf ( 1, '  use as a starting guess.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Other parameters are set as follows:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The maximum number of steps is %d.\n', max_it );
  fprintf ( 1, '  The error tolerance is %f\n', tol );

  [ x, error_norm, iter, flag ] = cgs ( A, x, b, M, max_it, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The CGS routine has returned with FLAG = %d\n', flag );
  if ( flag == 0 )
    fprintf ( 1, '  This indicates that the iteration has converged.\n' );
  elseif ( flag == 1 )
    fprintf ( 1, '  This indicates that the iteration has NOT converged.\n' );
  elseif ( flag == -1 )
    fprintf ( 1, '  This indicates that the iteration has broken down.\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of iterations taken was %d\n', iter );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The L2 norm of the error per iteration:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : iter
    fprintf ( 1, '  %4d  %f\n', i, error_norm(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The computed solution vector X:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %f\n', i, x(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CGS_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [new_gait, sol, info, total_time] = solve(nlp, x0, info)

    % Link the NLP problem to a NLP solver
    solver = IpoptApplication(nlp);

    solver.Options.ipopt.bound_frac = 1e-8;
    solver.Options.ipopt.bound_push = 1e-8;
    % Possible parameter changes to solver
    solver.Options.ipopt.tol = 1e-3;
    solver.Options.ipopt.dual_inf_tol = 1e-3;
    solver.Options.ipopt.constr_viol_tol = 1e-3;
    solver.Options.ipopt.compl_inf_tol = 1e-3;
    solver.Options.ipopt.max_iter = 1000;

    %%% Only use this if contraints are the same
    if nargin > 2
        solver.Options.ipopt.point_perturbation_radius = 0;
        solver.Options.ipopt.warm_start_bound_push = 1e-8;
        solver.Options.ipopt.warm_start_mult_bound_push = 1e-8;
        solver.Options.ipopt.mu_init = 1e-6;
        solver.Options.ipopt.warm_start_init_point = 'yes';
        solver.Options.zl = info.zl;
        solver.Options.zu = info.zu;
        solver.Options.lambda = info.lambda;
    end

    %     solver.Options.ipopt.acceptable_tol = 1e-2;
    %     solver.Options.ipopt.acceptable_constr_viol_tol = 1e-4;
    %     solver.Options.ipopt.acceptable_dual_inf_tol = 1e-2;
    %     solver.Options.ipopt.acceptable_compl_inf_tol = 1e-2;
    %

    solver.Options.ipopt.print_timing_statistics = 'yes';
    solver.Options.ipopt.nlp_lower_bound_inf = -1e6;
    solver.Options.ipopt.nlp_upper_bound_inf = 1e6;
    % Run the optimization
    start_time = tic();
    if nargin < 2
        [sol, info] = optimize(solver);
    else
        [sol, info] = optimize(solver, x0);
    end
    [tspan, states, inputs, params] = exportSolution(nlp, sol);
    tspan{3} = tspan{1}(end) + tspan{3};

    total_time = toc(start_time);
    new_gait = struct(...
        'tspan', tspan,...
        'states', states,...
        'inputs', inputs,...
        'params', params);

end

author: ayonga

repo: frost-dev

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


text:
%--- help for norm ---
%
% NORM   Matrix or vector norm.
%      NORM(X,2) returns the 2-norm of X.
%
%      NORM(X) is the same as NORM(X,2).
%
%      NORM(X,1) returns the 1-norm of X.
%
%      NORM(X,Inf) returns the infinity norm of X.
%
%      NORM(X,'fro') returns the Frobenius norm of X.
%
%    In addition, for vectors...
%
%      NORM(V,P) returns the p-norm of V defined as SUM(ABS(V).^P)^(1/P).
%
%      NORM(V,Inf) returns the largest element of ABS(V).
%
%      NORM(V,-Inf) returns the smallest element of ABS(V).
%
%    By convention, NaN is returned if X or V contains NaNs.
%
%    See also COND, RCOND, CONDEST, NORMEST, HYPOT.
%
%    Reference page in Doc Center
%       doc norm
%
%

author: jmaih


text:
function tests = test_dst
  tests = functiontests(localfunctions);
end

function test_dst_1(testCase)
    x = cos((1:15) * .1);
    alpha = spx.dsp.dst.forward_1(x);
    y = spx.dsp.dst.inverse_1(alpha);
    tolerance = 1e-6;
    verifyEqual(testCase, x, y, 'AbsTol', tolerance);
end

function test_dst_2(testCase)
    n = 256;
    tolerance = 1e-6;
    for i=randperm(n, 40)
        x = spx.vector.unit_vector(n, i);
        alpha = spx.dsp.dst.forward_2(x);
        verifyEqual(testCase, norm(x), norm(alpha), 'AbsTol', tolerance);
        y = spx.dsp.dst.inverse_2(alpha);
        verifyEqual(testCase, norm(y), norm(alpha), 'AbsTol', tolerance);
        verifyEqual(testCase, x, y, 'AbsTol', tolerance);
        alpha = spx.dsp.dst.forward_3(x);
        verifyEqual(testCase, norm(x), norm(alpha), 'AbsTol', tolerance);
        y = spx.dsp.dst.inverse_3(alpha);
        verifyEqual(testCase, norm(y), norm(alpha), 'AbsTol', tolerance);
        verifyEqual(testCase, x, y, 'AbsTol', tolerance);
    end
end

author: indigits

repo: sparse-plex


text:
classdef(SharedTestFixtures = {
        DownloadBERTFixture}) tWordPieceTokenizer < matlab.unittest.TestCase
    % tWordPieceTokenizer   Unit tests for WordPieceTokenizer

    % Copyright 2021 The MathWorks, Inc.

    properties(TestParameter)
        ModelAndExpectedVocabSize = struct(...
            'uncased',struct('Model','base','VocabSize',30522),...
            'multilingualCased',struct('Model','multilingual-cased','VocabSize',119547))
    end

    methods(Test)
        function canConstruct(test)
            enc = wordEncoding("foo");
            tok = bert.tokenizer.internal.WordPieceTokenizer(enc);
            test.verifyClass(tok,'bert.tokenizer.internal.WordPieceTokenizer');
            test.verifyInstanceOf(tok,'bert.tokenizer.internal.Tokenizer');
        end

        function canConstructWithFile(test)
            words = ["foo","bar"];
            enc = wordEncoding(words);
            fixture = matlab.unittest.fixtures.TemporaryFolderFixture();
            test.applyFixture(fixture);
            folder = fixture.Folder;
            testVocab = fullfile(folder,"testvocab.txt");
            fid = fopen(testVocab,'w','n','utf-8');
            fprintf(fid,join(words,newline));
            fclose(fid);
            tokFromEncoding = bert.tokenizer.internal.WordPieceTokenizer(enc);
            tokFromFile = bert.tokenizer.internal.WordPieceTokenizer(testVocab);
            test.verifyEqual(tokFromEncoding.Vocab,tokFromFile.Vocab);
        end

        function canSetUnknownToken(test)
            enc = wordEncoding("foo");
            unk = "bar";
            tok = bert.tokenizer.internal.WordPieceTokenizer(enc,'UnknownToken',unk);
            test.verifyEqual(tok.Unk,unk)
            str = "blah";
            ustr = textanalytics.unicode.UTF32(str);
            act_out = tok.tokenize(ustr);
            exp_out = unk;
            test.verifyEqual(act_out,exp_out);
        end

        function canSetMaxTokenLength(test)
            enc = wordEncoding("foo");
            maxLen = 2;
            tok = bert.tokenizer.internal.WordPieceTokenizer(enc,'MaxTokenLength',maxLen);
            test.verifyEqual(tok.MaxChar,maxLen);
            str = "foo";
            ustr = textanalytics.unicode.UTF32(str);
            act_out = tok.tokenize(ustr);
            exp_out = tok.Unk;
            test.verifyEqual(act_out,exp_out);
        end

        function canTokenize(test)
            enc = wordEncoding(["foo","bar","##foo"]);
            tok = bert.tokenizer.internal.WordPieceTokenizer(enc);
            str = "foo bar foobar barba bafoobar barfoo";
            wsTok = bert.tokenizer.internal.WhitespaceTokenizer;
            ustr = textanalytics.unicode.UTF32(wsTok.tokenize(str));
            act_out = tok.tokenize(ustr);
            exp_out = ["foo","bar",tok.Unk,tok.Unk,tok.Unk,"bar","##foo"];
            test.verifyEqual(act_out,exp_out);
        end

        function hasExpectedVocabSize(test,ModelAndExpectedVocabSize)
            vocab = bert.internal.getSupportFilePath(ModelAndExpectedVocabSize.Model,"vocab.txt");
            tok = bert.tokenizer.internal.WordPieceTokenizer(vocab);
            test.verifyEqual(tok.Vocab.NumWords,ModelAndExpectedVocabSize.VocabSize);
        end
    end
end

author: matlab-deep-learning

repo: transformer-models


text:
function net = fixduplicates(net, x, y)

% FIXDUPLICATES
%
% Ensure that identical support vector of each class have identical Lagrange
% multipliers.
%
%    net = fixduplicates(net, x, y)
%
% where x and y are the training data.  The strip or compact methods must not
% previously have been applied to this network.

%
% File        : @svc/fixduplicates.m
%
% Date        : Tuesday 12th September 2000
%
% Author      : Dr Gavin C. Cawley
%
% Description : Part of an object-oriented implementation of Vapnik's Support
%               Vector Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 07/07/2000 - v1.00
%               12/09/2000 - v1.01 minor improvements to comments and help
%                                  message
%               16/09/2000 - v1.10 fixed a bug requiring the training data
%                                  to be provided as parameters
%
% Copyright   : (c) Dr Gavin C. Cawley, September 2000
%
%    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 should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

[foo, i, j] = unique(x, 'rows');

for k = unique(j)'

   idx1 = find(j == k);

   idx2 = find(y(idx1) > 0);

   if ~isempty(idx2)
      net.w(idx1(idx2)) = mean(net.w(idx1(idx2)));
   end

   idx2 = find(y(idx1) < 0);

   if ~isempty(idx2)
      net.w(idx1(idx2)) = mean(net.w(idx1(idx2)));
   end

end

% bye bye...

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [w, b] = Train_SVR(samples, labels)
%Train_SVR creating a linear support vector regressor

    % liblinear training
    addpath('C:\liblinear\matlab');

    % Remove redundant data
    [samples, inds] = unique(samples, 'rows');
    labels = labels(inds,:);

    cmd = ['-s 11 -B 1 -q'];

    svr_regressor = train(labels, sparse(double(samples)), cmd);

    w = svr_regressor.w(1:end-1)';
    b = svr_regressor.w(end);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function eval = clusteval(new,result,param)

v = result.cluster.v;
c = size(result.cluster.v,1);%c = param.c;
if exist('param.m')==1, m = param.m;else m = 2;end;

X=new.X;
[N,n] = size(X);
X1 = ones(N,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(result.cluster,'M')%GK
    M = result.cluster.M;
    for j = 1 : c,
        xv = X - X1*v(j,:);
        d(:,j) = sum((xv*M(:,:,j).*xv),2);
    end
    distout=sqrt(d);
    %Update the partition matrix
    d = (d+1e-10).^(-1/(m-1));
    f0 = (d ./ (sum(d,2)*ones(1,c)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif isfield(result.cluster,'P')%GG
    A = result.cluster.P;
    f = result.data.f;
    fm = f.^m;
    for j = 1 : c,
        xv = X - X1*v(j,:);
        Pi(:,:,j)=1/N*sum(fm(:,j));
        A = result.cluster.P(:,:,j);
        d(:,j) = 1/(det(pinv(A))^(1/2))*1/Pi(:,:,j)*exp(1/2*sum((xv*pinv(A).*xv),2));
    end
    distout=sqrt(d);
    %Update the partition matrix
    if m>1
          d = (d+1e-10).^(-1/(m-1));
      else
          d = (d+1e-10).^(-1);
      end
    f0 = (d ./ (sum(d,2)*ones(1,c)));

else        %FCM
     for j = 1 : c,
      xv = X - X1*v(j,:);
      d(:,j) = sum((xv*eye(n).*xv),2);
    end;
    distout=sqrt(d);
    d = (d+1e-10).^(-1/(m-1));
    f0 = (d ./ (sum(d,2)*ones(1,c)));
end
%results
    eval.d = distout;
    eval.f = f0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Visualization

if n == 2   %in 2dimensional case draw a contour map
    lower1=min(X(:,1));upper1=max(X(:,1));scale1=(upper1-lower1)/200;
    lower2=min(X(:,2));upper2=max(X(:,2));scale2=(upper2-lower2)/200;
    [x,y] = meshgrid(lower1:scale1:upper1, lower2:scale2:upper2);
    pair = [x(:) y(:)];
    [pair1,pair2] = size(pair);
    X1 = ones(pair1,1);
    d=zeros(pair1,c);  %resize the distance matrix

    if isfield(result.cluster,'M')%GK
        for j = 1 : c,
            xv = pair - X1*v(j,:);
            d(:,j) = sum((xv*M(:,:,j).*xv),2);
        end
        distout=sqrt(d);
        d = (d+1e-10).^(-1/(m-1));
        f0 = (d ./ (sum(d,2)*ones(1,c)));

    elseif isfield(result.cluster,'P')%GG
        for j = 1 : c,
            xv = pair - X1*v(j,:);
            Pi(:,:,j)=1/N*sum(fm(:,j));
            A = result.cluster.P(:,:,j);
            d(:,j) = 1/(det(pinv(A))^(1/2))*1/Pi(:,:,j)*exp(1/2*sum((xv*pinv(A).*xv),2));
        end
        distout=sqrt(d);
        if m>1
            d = (d+1e-10).^(-1/(m-1));
        else
            d = (d+1e-10).^(-1);
        end
        f0 = (d ./ (sum(d,2)*ones(1,c)));
    else   %FCM
        for i = 1 : c,
            xv = pair - ones(pair1,1)*v(i,:);
            d(:,i)= sum((xv*eye(2).*xv),2);%
        end;
        distout=sqrt(d);
        d = (d+1e-10).^(-1/(m-1));
        f0 = (d ./ (sum(d,2)*ones(1,c)));
    end
    f=max(f0')';
    Z= reshape(f,size(x,1),size(x,2));
    contour(x,y,Z);
    drawnow
end

author: Sable

repo: mcbench-benchmarks


text:
function [output] = volumethreshold(input, threshold, tissuelabel)

% VOLUMETHRESHOLD is a helper function for segmentations. It applies a
% relative threshold and subsequently looks for the largest connected part,
% thereby removing small blobs such as vitamine E capsules.
%
% See also VOLUMEFILLHOLES, VOLUMESMOOTH, VOLUMEPAD

if nargin<2 || isempty(threshold)
  threshold = 0;
end
if nargin<3 || isempty(tissuelabel)
  tissuelabel = 'volume';
end

% ensure that SPM is available, needed for spm_bwlabel
hasspm = ft_hastoolbox('spm8up', 3) || ft_hastoolbox('spm2', 1);

% mask by taking the negative of the segmentation, thus ensuring
% that no holes are within the compartment and do a two-pass
% approach to eliminate potential vitamin E capsules etc.

% ensure the input volume to be double precision, to allow for
% computations, unless it's already boolean
if ~isa(input, 'double') && ~isa(input, 'logical')
  input = double(input);
end

if ~isa(input, 'logical')
  if nargin==2, ft_error('if the input volume is not a boolean volume, you need to define a threshold value'); end
  if nargin==3, fprintf('thresholding %s at a relative threshold of %0.3f\n', tissuelabel, threshold); end
  output = double(input>(threshold*max(input(:))));
else
  % there is no reason to apply a threshold, but spm_bwlabel still needs a
  % double input for clustering
  output = double(input);
end

% cluster the connected tissue
[cluster, n] = spm_bwlabel(output, 6);

if n>1
  % it pays off to sort the cluster assignment if there are many clusters
  tmp = cluster(:);                       % convert to a vector
  tmp = tmp(tmp>0);                       % remove the zeros
  tmp = sort(tmp, 'ascend');              % sort according to cluster number
  m   = zeros(1,n);
  for k=1:n
    m(k) = sum(tmp==k);       % determine the voxel count for each cluster
    tmp  = tmp(m(k)+1:end);   % remove the last cluster that was counted
  end
  % select the tissue that has the most voxels belonging to it
  [m, i] = max(m);
  output = (cluster==i);
else
  % the output only contains a single cluster
  output = (cluster==1);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function r8mat_trace_test ( )

%*****************************************************************************80
%
%% R8MAT_TRACE_TEST tests R8MAT_TRACE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;

  for i = 1 : n
    for j = 1 : n

      if ( i <= j )
        a(i,j) = n + 1 - j;
      elseif ( j == i - 1 )
        a(i,j) = n - j;
      else
        a(i,j) = 0.0;
      end

    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_TRACE_TEST\n' );
  fprintf ( 1, '  R8MAT_TRACE computes the trace of a matrix\n' );

  r8mat_print ( n, n, a, '  Matrix:' );

  trace = r8mat_trace ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Trace is %f\n', trace );

  return
end


text:
function [ ap, kpvt, rcond, z ] = sspco ( ap, n )

%*****************************************************************************80
%
%% SSPCO factors a real symmetric matrix stored in packed form.
%
%  Discussion:
%
%    SSPCO uses elimination with symmetric pivoting and estimates
%    the condition of the matrix.
%
%    If RCOND is not needed, SSPFA is slightly faster.
%
%    To solve A*X = B, follow SSPCO by SSPSL.
%
%    To compute inverse(A)*C, follow SSPCO by SSPSL.
%
%    To compute inverse(A), follow SSPCO by SSPDI.
%
%    To compute determinant(A), follow SSPCO by SSPDI.
%
%    To compute inertia(A), follow SSPCO by SSPDI.
%
%  Packed storage:
%
%    The following program segment will pack the upper triangle of a
%    symmetric matrix.
%
%      k = 0
%      do j = 1, n
%        do i = 1, j
%          k = k + 1
%          ap(k) = a(i,j)
%        end
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real AP(N*(N+1)/2), the packed form of a symmetric matrix A.  The
%    columns of the upper triangle are stored sequentially in a one-dimensional array.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real AP(N*(N+1)/2), a block diagonal matrix and the multipliers
%    which were used to obtain it, stored in packed form.  The
%    factorization can be written A = U*D*U' where U is a product of
%    permutation and unit upper triangular matrices, U' is the transpose
%    of U, and D is block diagonal with 1 by 1 and 2 by 2 blocks.
%
%    Output, integer KPVT(N), the pivot indices.
%
%    Output, real RCOND, an estimate of the reciprocal condition
%    of A.  For the system A*X = B, relative perturbations in A and B of size
%    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
%    If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%    Output, real Z(N) a work vector whose contents are usually
%    unimportant.  If A is close to a singular matrix, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%

%
%  Find norm of A using only upper half.
%
  j1 = 1;
  for j = 1 : n
    z(j) = sasum ( j, ap(j1:j1+j-1), 1 );
    ij = j1;
    j1 = j1 + j;
    for i = 1 : j-1
      z(i) = z(i) + abs ( ap(ij) );
      ij = ij + 1;
    end
  end

  anorm = max ( z(1:n) );
%
%  Factor.
%
  [ ap, kpvt, info ] = sspfa ( ap, n );
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where U*D*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve U * D * W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  k = n;
  ik = floor ( ( n * ( n - 1 ) ) / 2 );

  while ( k ~= 0 )

    kk = ik + k;
    ikm1 = ik - ( k - 1 );

    if ( kpvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    kp = abs ( kpvt(k) );
    kps = k + 1 - ks;

    if ( kp ~= kps )
      t = z(kps);
      z(kps) = z(kp);
      z(kp) = t;
    end

    if ( z(k) ~= 0.0 )
      ek = abs ( ek ) * r4_sign ( z(k) );
    end

    z(k) = z(k) + ek;
    z(1:k-ks) = saxpy ( k-ks, z(k), ap(ik+1:ik+k-ks), 1, z(1:k-ks), 1 );

    if ( ks ~= 1 )
      if ( z(k-1) ~= 0.0 )
        ek = abs ( ek ) * r4_sign ( z(k-1) );
      end
      z(k-1) = z(k-1) + ek;
      z(1:k-ks) = saxpy ( k-ks, z(k-1), ap(ikm1+1:ikm1+k-ks), 1, z(1:k-ks), 1 );
    end

    if ( ks ~= 2 )

      if ( abs ( ap(kk) ) < abs ( z(k) ) )
        s = abs ( ap(kk) ) / abs ( z(k) );
        z(1:n) = s * z(1:n);
        ek = s * ek;
      end

      if ( ap(kk) ~= 0.0 )
        z(k) = z(k) / ap(kk);
      else
        z(k) = 1.0;
      end

    else

      km1k = ik + k - 1;
      km1km1 = ikm1 + k - 1;
      ak = ap(kk) / ap(km1k);
      akm1 = ap(km1km1) / ap(km1k);
      bk = z(k) / ap(km1k);
      bkm1 = z(k-1) / ap(km1k);
      denom = ak * akm1 - 1.0;
      z(k) = ( akm1 * bk - bkm1 ) / denom;
      z(k-1) = ( ak * bkm1 - bk ) / denom;

    end

    k = k - ks;
    ik = ik - k;
    if ( ks == 2 )
      ik = ik - ( k + 1 );
    end

  end

  z(1:n) = z(1:n) / sasum ( n, z(1:n), 1 );
%
%  Solve U' * Y = W.
%
  k = 1;
  ik = 0;

  while ( k <= n )

    if ( kpvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= 1 )

      z(k) = z(k) + sdot ( k-1, ap(ik+1:ik+k-1), 1, z(1:k-1), 1 );
      ikp1 = ik + k;

      if ( ks == 2 )
        z(k+1) = z(k+1) + sdot ( k-1, ap(ikp1+1:ikp1+k-1), 1, z(1:k-1), 1 );
      end

      kp = abs ( kpvt(k) );

      if ( kp ~= k )
        t = z(k);
        z(k) = z(kp);
        z(kp) = t;
      end

    end

    ik = ik + k;
    if ( ks == 2 )
      ik = ik + ( k + 1 );
    end
    k = k + ks;

  end

  s = 1.0 / sasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = 1.0;
%
%  Solve U * D * V = Y.
%
  k = n;

  ik = floor ( n * ( n - 1 ) / 2 );

  while ( 0 < k )

    kk = ik + k;
    ikm1 = ik - ( k - 1 );

    if ( kpvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= ks )

      kp = abs ( kpvt(k) );
      kps = k + 1 - ks;

      if ( kp ~= kps )
        t = z(kps);
        z(kps) = z(kp);
        z(kp) = t;
      end

      z(1:k-ks) = saxpy ( k-ks, z(k), ap(ik+1:ik+k-ks), 1, z(1:k-ks), 1 );

      if ( ks == 2 )
        z(1:k-ks) = saxpy ( k-ks, z(k-1), ap(ikm1+1:ikm1+k-ks), 1, z(1:k-ks), 1 );
      end

    end

    if ( ks ~= 2 )

      if ( abs ( ap(kk) ) < abs ( z(k) ) )
        s = abs ( ap(kk) ) / abs ( z(k) );
        z(1:n) = s * z(1:n);
        ynorm = s * ynorm;
      end

      if ( ap(kk) ~= 0.0 )
        z(k) = z(k) / ap(kk);
      else
        z(k) = 1.0;
      end

    else

      km1k = ik + k - 1;
      km1km1 = ikm1 + k - 1;
      ak = ap(kk) / ap(km1k);
      akm1 = ap(km1km1) / ap(km1k);
      bk = z(k) / ap(km1k);
      bkm1 = z(k-1) / ap(km1k);
      denom = ak * akm1 - 1.0;
      z(k) = ( akm1 * bk - bkm1 ) / denom;
      z(k-1) = ( ak * bkm1 - bk ) / denom;

    end

    k = k - ks;
    ik = ik - k;
    if ( ks == 2 )
      ik = ik - ( k + 1 );
    end

  end

  s = 1.0 / sasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;
%
%  Solve U' * Z = V.
%
  k = 1;
  ik = 0;

  while ( k <= n )

    if ( kpvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= 1 )

      z(k) = z(k) + sdot ( k-1, ap(ik+1:ik+k-1), 1, z(1:k-1), 1 );
      ikp1 = ik + k;

      if ( ks == 2 )
        z(k+1) = z(k+1) + sdot ( k-1, ap(ikp1+1:ikp1+k-1), 1, z(1:k-1), 1 );
      end

      kp = abs ( kpvt(k) );

      if ( kp ~= k )
        t = z(k);
        z(k) = z(kp);
        z(kp) = t;
      end

    end

    ik = ik + k;
    if ( ks == 2 )
      ik = ik + ( k + 1 );
    end
    k = k + ks;

  end
%
%  Make ZNORM = 1.0.
%
  s = 1.0 / sasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end


text:
function dist = shape_point_dist_2d ( center, p1, nside, p )

%*****************************************************************************80
%
%% SHAPE_POINT_DIST_2D: distance ( regular shape, point ) in 2D.
%
%  Discussion:
%
%    The "regular shape" is assumed to be an equilateral and equiangular
%    polygon, such as the standard square, pentagon, hexagon, and so on.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CENTER(2,1), the center of the shape.
%
%    Input, real P1(2,1), the first vertex of the shape.
%
%    Input, integer NSIDE, the number of sides in the shape.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, real DIST, the distance from the point to the shape.
%

%
%  Determine the angle subtended by a single side.
%
  sector_angle = 360.0 / nside;
%
%  How long is the half-diagonal?
%
  radius = sqrt ( sum ( ( p1(1:2,1) - center(1:2,1) ).^2 ) );
%
%  If the radius is zero, then the shape is a point and the computation is easy.
%
  if ( radius == 0.0 )
    dist = sqrt ( sum ( ( p(1:2,1) - center(1:2,1) ).^2 ) );
    return;
  end
%
%  If the test point is at the center, then the computation is easy.
%  The angle subtended by any side is ( 2 * PI / NSIDE ) and the
%  nearest distance is the midpoint of any such side.
%
  if ( sqrt ( sum ( ( p1(1:2,1) - center(1:2,1) ).^2 ) ) == 0.0 )
    dist = radius * cos ( pi / nside );
    return
  end
%
%  Determine the angle between the ray to the first corner,
%  and the ray to the test point.
%
  angle = angle_deg_2d ( p1, center, p );
%
%  Determine the sector of the point.
%
  sector_index = floor ( angle / sector_angle ) + 1;
%
%  Generate the two corner points that terminate the SECTOR-th side.
%
  angle2 = ( sector_index - 1 ) * sector_angle;
  angle2 = degrees_to_radians ( angle2 );

  pa = vector_rotate_base_2d ( p1, center, angle2 );

  angle2 = ( sector_index ) * sector_angle;
  angle2 = degrees_to_radians ( angle2 );

  pb = vector_rotate_base_2d ( p1, center, angle2 );
%
%  Determine the distance from the test point to the line segment that
%  is the SECTOR-th side.
%
  dist = segment_point_dist_2d ( pa, pb, p );

  return
end


text:
function sphere = createSphere(varargin)
%CREATESPHERE Create a sphere containing 4 points
%
%   s = createSphere(p1, p2, p3, p4);
%   return in s the sphere common to the 4 pointsp1, p2, p3 and p4.
%
%   Ref: P. Bourke
%   http://astronomy.swin.edu.au/~pbourke/geometry/spherefrom4/
%
%   See also
%   spheres, circles3d
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 22/03/2005.
%

if length(varargin)==4
    pts = [varargin{1};varargin{2};varargin{3};varargin{4}];
elseif length(varargin)==1
    pts = varargin{1};
else
    error('wrong number of arguments in createSphere');
end

m1 = det([pts ones(4,1)]);
s2 = sum(pts.*pts, 2);
m2 = det([s2 pts(:,2) pts(:,3) ones(4,1)]);
m3 = det([pts(:,1) s2 pts(:,3) ones(4,1)]);
m4 = det([pts(:,1) pts(:,2) s2 ones(4,1)]);

m5 = det([s2 pts]);

x0 = m2*.5/m1;
y0 = m3*.5/m1;
z0 = m4*.5/m1;
r  = sqrt(x0*x0 + y0*y0 + z0*z0 - m5/m1);

sphere = [x0 y0 z0 r];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function h=init_h(method,peval,image,dpixc)
% h=init_h(method,peval,image,dpixc)
% method:   'rand' %random initialization
%           'image' %specified initialization for example image=double(array2im(dpixc_ind))
%           'image_repmat' %for example image = mean(image,3);
%           'res' %image is res.h -> initialisation from the nmf results
% addbackgroundcomponent:   1 adds one component flat component as a background (peval.ncomp th)
%                           0 no background component
if ~isfield(peval, 'bgcomp')
    peval.bgcomp = 1;
    mfprintf(peval.fid, 'Value of the peval.bgcomp was set to %g (by default the last component is backgroud).\n',peval.bgcomp)
end

switch method
    case 'rand'
        h = mean(dpixc(:)-peval.bg)*rand(peval.ncomp,peval.nt);
        msg='h initialzied as uniform random.';
    case 'res'
        h = image;
        msg='h initialized from the results res.h.';
end
h=max(h, eps); % To avoid zeros...

if isfield (peval,'fid')
    mfprintf(peval.fid, [msg '\n'])
else
    fprintf([msg '\n']);
end

if peval.bgcomp
    h(peval.ncomp,:)=peval.nx*peval.ny*peval.bg;
    mfprintf(peval.fid, 'Last component [%g] of h initialised as a flat background. (background=%g)\n',peval.ncomp,peval.bg);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function hn=arrow3(p1,p2,s,w,h,ip,alpha,beta)
% ARROW3
%   ARROW3(P1,P2) will draw vector lines (2D/3D) from P1 to P2 with
%   arrowheads, where P1 and P2 can be either nx2 matrices (for 2D),
%   or nx3 matrices (for 3D).
%
%   ARROW3(P1,P2,S,W,H,IP,ALPHA,BETA) can be used to specify properties
%   of the line and arrowhead.  S is a character string made with one
%   element from any or all of the following 3 columns:
%
%      Color Switches       LineStyle             LineWidth
%      ------------------   -------------------   --------------------
%      k  black (default)   -   solid (default)   0.5 points (default)
%      y  yellow            :   dotted            0   no lines
%      m  magenta           -.  dashdot           /   LineWidthOrder
%      c  cyan              --  dashed
%      r  red               *   LineStyleOrder
%      g  green
%      b  blue
%      w  white                                    __________ __
%      a  apple green                            ^           |
%      d  dark gray                             / \          |
%      e  evergreen                            /   \         |
%      f  fuchsia                             /     \        |
%      h  honey                              /       \    Height
%      i  indigo                            /         \      |
%      j  jade                             /           \     |
%      l  lilac                           /             \    |
%      n  nutbrown                       /______   ______\ __|__
%      p  pink                           |      | |      |
%      q  kumquat                        |---- Width ----|
%      s  sky blue                       |      | |      |
%      t  tan                                   | |
%      u  umber                                 | |
%      v  violet                                | |
%      z  zinc                               -->| |<--LineWidth
%      x  random named color                    | |
%      o  ColorOrder
%      |  magnitude
%
%   The components of S may be specified in any order.  Invalid
%   characters in S will be ignored and replaced by default settings.
%
%   Prefixing the color code with '_' produces a darker shade, e.g.
%   '_t' is dark tan; prefixing the color code with '^' produces a
%   lighter shade, e.g. '^q' is light kumquat.  The relative
%   brightness of light and dark color shades is controled by the
%   scalar parameter BETA.  Color code prefixes do not affect the
%   basic colors (kymcrgbw) or the special color switches (xo|).
%
%   ColorOrder may be achieved in two fashions:  The user may either
%   set the ColorOrder property (using RGB triples) or define the
%   global variable ColorOrder (using a string of valid color codes).
%   If the color switch is specified with 'o', and the global variable
%   ColorOrder is a string of color codes (color switches less 'xo|',
%   optionally prefixed with '_' or '^'), then the ColorOrder property
%   will be set to the sequence of colors indicated by the ColorOrder
%   variable.  If the color switch is specified with 'o', and the
%   global variable ColorOrder is empty or invalid, then the current
%   ColorOrder property will be used.  Note that the ColorOrder
%   variable takes precedence over the ColorOrder property.
%
%   The magnitude color switch is used to visualize vector magnitudes
%   in conjunction with a colorbar.  If the color switch is specified
%   with '|', colors are linearly interpolated from the current ColorMap
%   according to the length of the associated line.  This option sets
%   CLim to [MinM,MaxM], where MinM and MaxM are the minimum and maximum
%   magnitudes, respectively.
%
%   The current LineStyleOrder property will be used if LineStyle is
%   specified with '*'.  MATLAB cycles through the line styles defined
%   by the LineStyleOrder property only after using all colors defined
%   by the ColorOrder property.  If however, the global variable
%   LineWidthOrder is defined, and LineWidth is specified with '/',
%   then each line will be drawn with sequential color, linestyle, and
%   linewidth.
%
%   W is a vector of arrowhead widths.  For linear plots with equal
%   axes, the units of W (default = 1/72 of the PlotBox diagonal) are
%   the same as those of the coordinate data (P1,P2).  For linear plots
%   with unequal axes, the units of W are scaled to fit as if the axes
%   were equal.  For log plots, the units of W (default = 1) are 1/72
%   of the PositionRectangle diagonal.
%
%   H (default = 3W) is a vector of arrowhead heights.  If vector IP is
%   neither empty nor negative, initial point markers will be plotted
%   with diameter IP; for default diameter W, use IP = 0.  The units of
%   W, H, and IP are absolute for linear plots and relative to the
%   PositionRectangle for log plots.
%
%   ALPHA (default = 1) is a vector of FaceAlpha (MATLAB 6+) values
%   ranging between 0 (clear) and 1 (opaque).  FaceAlpha is a surface
%   (arrowhead and initial point marker) property and does not affect
%   lines.  FaceAlpha is not supported for 2D rendering.
%
%   BETA (default = 0.4) is a scalar that controls the relative
%   brightness of light and dark color shades, ranging between 0 (no
%   contrast) and 1 (maximum contrast).
%
%   Plotting lines with a single color, linestyle, and linewidth is
%   faster than plotting lines with multiple colors and/or linestyles.
%   Plotting lines with multiple linewidths is slower still.  ARROW3
%   chooses renderers that produce the best screen images; exported
%   or printed plots may benifit from different choices.
%
%   HN = ARROW3(P1,P2,...) returns a vector of handles to line and
%   surface objects created by ARROW3.
%
%   ARROW3 COLORS will plot a table of named colors with default
%   brightness.  ARROW3('colors',BETA) will plot a table of named
%   colors with brightness BETA.
%
%   If a particular aspect ratio is required, use DASPECT, PBASPECT,
%   AXIS, or XYZLIM commands before calling ARROW3.  Changing limits
%   or aspect ratios after calling ARROW3 may alter the appearance
%   of arrowheads and initial point markers.  ARROW3 sets XYZCLimMode
%   to manual for all plots, sets DataAspectRatioMode to manual for
%   linear plots, and sets PlotBoxAspectRatioMode to manual for log
%   plots and 3D plots.
%
%   ARROW3 UPDATE will restore the the appearance of arrowheads and
%   initial point markers that have become corrupted by changes to
%   limits or aspect ratios.  ARROW3('update',SF) will redraw initial
%   point markers and arrowheads with scale factor SF.  If SF has one
%   element, SF scales W, H and IP.  If SF has two elements, SF(1)
%   scales W and IP, and SF(2) scales H.  If SF has three elements,
%   SF(1) scales W, SF(2) scales H, and SF(3) scales IP.
%
%   ARROW3 UPDATE COLORS will update the magnitude coloring of
%   arrowheads, initial point markers, and lines to conform to the
%   current colormap.
%
%   HN = ARROW3('update',...) returns a vector of handles to updated
%   objects.
%
%   Usage Examples:
%
%      % 2D vectors
%      Arrow3([0 0],[1 3])
%      Arrow3([0 0],[1 2],'-.e')
%      Arrow3([0 0],[10 10],'--x2',1)
%      Arrow3(zeros(10,2),50*rand(10,2),'x',1,3)
%      Arrow3(zeros(10,2),[10*rand(10,1),500*rand(10,1)],'u')
%      Arrow3(10*rand(10,2),50*rand(10,2),'x',1,[],1)
%
%      % 3D vectors
%      Arrow3([0 0 0],[1 1 1])
%      Arrow3(zeros(20,3),50*rand(20,3),'--x1.5',2)
%      Arrow3(zeros(100,3),50*rand(100,3),'x',1,3)
%      Arrow3(zeros(10,3),[10*rand(10,1),500*rand(10,1),50*rand(10,1)],'a')
%      Arrow3(10*rand(10,3),50*rand(10,3),'x',[],[],0)
%
%      % Just for fun
%      Arrow3(zeros(100,3),50*rand(100,3),'x',10,3,[],0.95)
%      light('Position',[-10 -10 -10],'Style','local')
%      light('Position',[60,60,60]), lighting gouraud
%
%      % ColorOrder variable, color code prefixes, and Beta
%      global ColorOrder, ColorOrder='^ui^e_hq^v';
%      theta=[0:pi/22:pi/2]';
%      Arrow3(zeros(12,2),[cos(theta),sin(theta)],'1.5o',0.1,[],[],[],0.5)
%
%      % ColorOrder property, LineStyleOrder, and LineWidthOrder
%      global ColorOrder, ColorOrder=[];
%      set(gca,'ColorOrder',[1,0,0;0,0,1;0.25,0.75,0.25;0,0,0])
%      set(gca,'LineStyleOrder',{'-','--','-.',':'})
%      global LineWidthOrder, LineWidthOrder=[1,2,4,8];
%      w=[5,10,15,20]; h=[20,40,30,40];
%      Arrow3(zeros(4,2),[10*rand(4,1),500*rand(4,1)],'o*/',w,h,10)
%
%      % Magnitude coloring
%      colormap spring
%      Arrow3(20*randn(20,3),50*randn(20,3),'|',[],[],0)
%      set(gca,'color',0.7*[1,1,1])
%      set(gcf,'color',0.5*[1,1,1]), grid on, colorbar
%      pause % change the colormap and update colors
%      colormap hot
%      Arrow3('update','colors')
%
%      % LogLog plot
%      set(gca,'xscale','log','yscale','log');
%      axis([1e2,1e8,1e-2,1e-1]); hold on
%      p1=repmat([1e3,2e-2],15,1);
%      q1=[1e7,1e6,1e5,1e4,1e3,1e7,1e7,1e7,1e7,1e7,1e7,1e6,1e5,1e4,1e3];
%      q2=1e-2*[9,9,9,9,9,7,5,4,3,2,1,1,1,1,1]; p2=[q1',q2'];
%      global ColorOrder, ColorOrder=[];
%      set(gca,'ColorOrder',rand(15,3))
%      Arrow3(p1,p2,'o'), grid on, hold off
%
%      % SemiLogX plot
%      set(gca,'xscale','log','yscale','linear');
%      axis([1e2,1e8,1e-2,1e-1]); hold on
%      p1=repmat([1e3,0.05],15,1);
%      q1=[1e7,1e6,1e5,1e4,1e3,1e7,1e7,1e7,1e7,1e7,1e7,1e6,1e5,1e4,1e3];
%      q2=1e-2*[9,9,9,9,9,7,5,4,3,2,1,1,1,1,1]; p2=[q1',q2'];
%      Arrow3(p1,p2,'x'), grid on, hold off
%
%      % SemiLogY plot
%      set(gca,'xscale','linear','yscale','log');
%      axis([2,8,1e-2,1e-1]); hold on
%      p1=repmat([3,2e-2],17,1);
%      q1=[7,6,5,4,3,7,7,7,7,7,7,7,7,6,5,4,3];
%      q2=1e-2*[9,9,9,9,9,8,7,6,5,4,3,2,1,1,1,1,1]; p2=[q1',q2'];
%      set(gca,'LineStyleOrder',{'-','--','-.',':'})
%      Arrow3(p1,p2,'*',1,[],0), grid on, hold off
%
%      % Color tables
%      Arrow3('colors')           % default color table
%      Arrow3('colors',0.3)       % low contrast color table
%      Arrow3('colors',0.5)       % high contrast color table
%
%      % Update initial point markers and arrowheads
%      Arrow3('update')           % redraw same size
%      Arrow3('update',2)         % redraw double size
%      Arrow3('update',0.5)       % redraw half size
%      Arrow3('update',[0.5,2,1]) % redraw W half size,
%                                 %        H double size, and
%                                 %        IP same size

%   Copyright(c)2002,2003 Version 4.58
%      Jeff Chang <cpmame@hotmail.com>
%      Tom Davis  <tdavis@eng.usf.edu>

%   Revision History:
%
%      12/17/03 - Semilog examples, CAXIS support, magnitude
%                 coloring, and color updating; use CData instead
%                 of FaceColor; minor bug fixes (TD)
%      07/17/03 - Changed 2D rendering from OpenGL to ZBuffer;
%                 defined HN for COLORS and UPDATE options (TD)
%      02/27/03 - Replaced calls to RANDPERM, VIEW, REPMAT, SPHERE,
%                 and CYLINDER; added ZBuffer for log plots, RESET
%                 for CLA and CLF, and ABS for W and H (TD)
%      02/01/03 - Added UPDATE scale factor and matlab version
%                 checking, replaced call to CROSS (TD)
%      12/26/02 - Added UserData and UPDATE option (TD)
%      11/16/02 - Added more named colors, color code prefix,
%                 global ColorOrder, ALPHA , and BETA (TD)
%      10/12/02 - Added global LineWidthOrder,
%                 vectorized W, H and IP (TD)
%      10/05/02 - Changed CLF to CLA for subplot support,
%                 added ColorOrder and LineStyleOrder support (TD)
%      04/27/02 - Minor log plot revisions (TD)
%      03/26/02 - Added log plot support (TD)
%      03/24/02 - Adaptive grid spacing control to trade off
%                 appearance vs. speed based on size of matrix (JC)
%      03/16/02 - Added "axis tight" for improved appearance (JC)
%      03/12/02 - Added initial point marker (TD)
%      03/03/02 - Added aspect ratio support (TD)
%      03/02/02 - Enchance program's user friendliness (JC)
%                 (lump Color, LineStyle, and LineWidth together)
%      03/01/02 - Replaced call to ROTATE (TD)
%      02/28/02 - Modified line plotting,
%                 added linewidth and linestyle (TD)
%      02/27/02 - Minor enhancements on 3D appearance (JC)
%      02/26/02 - Minor enhancements for speed (TD&JC)
%      02/26/02 - Optimise PLOT3 and SURF for speed (TD)
%      02/25/02 - Return handler, error handling, color effect,
%                 generalize for 2D/3D vectors (JC)
%      02/24/02 - Optimise PLOT3 and SURF for speed (TD)
%      02/23/02 - First release (JC&TD)

%==========================================================================
% Error Checking
oldver=0; v=version; if v(1)<'6', oldver=1; end % matlab version
if nargin<8 | isempty(beta), beta=0.4; end
beta=abs(beta(1)); if nargout, hn=[]; end
if strcmpi(p1,'colors')                         % plot color table
   if nargin>1, beta=abs(p2(1)); end
   LocalColorTable(1,beta); return
end
fig=gcf; ax=gca;
if strcmpi(p1,'update'), ud=get(ax,'UserData'); % update
   if size(ud,2)<13, error('Invalid UserData'), end
   set(ax,'UserData',[]); sf=[1,1,1]; flag=0;
   if nargin>1
      if strcmp(p2,'colors'), flag=1;           % update colors
      else                                      % update surfaces
         sf=p2(1)*sf; n=length(p2(:));
         if n>1, sf(2)=p2(2); if n>2, sf(3)=p2(3); end, end
      end
   end
   H=LocalUpdate(fig,ax,ud,sf,flag); if nargout, hn=H; end, return
end
InputError=['Invalid input, type HELP ',upper(mfilename),...
   ' for usage examples'];
if nargin<2, error(InputError), end
[r1,c1]=size(p1); [r2,c2]=size(p2); n=r1; Zeros=zeros(n,1);
if c1<2 | c1>3, error(InputError), end
if r1~=r2, error('P1 and P2 must have same number of rows'), end
if c1~=c2, error('P1 and P2 must have same number of columns'), end
if c1==2, p1=[p1,Zeros]; p2=[p2,Zeros];
elseif ~any([p1(:,3);p2(:,3)]), c1=2; end
L=get(ax,'LineStyleOrder'); C=get(ax,'ColorOrder');
ST=get(ax,'DefaultSurfaceTag'); LT=get(ax,'DefaultLineTag');
EC=get(ax,'DefaultSurfaceEdgeColor');
if strcmp(get(ax,'nextplot'),'add') & strcmp(get(fig,'nextplot'),'add')
   Xr=get(ax,'xlim'); Yr=get(ax,'ylim'); Zr=get(ax,'zlim');
   set(ax,'XLimMode','auto','YLimMode','auto','ZLimMode','auto',...
      'CLim',get(ax,'CLim'));
   xs=strcmp(get(ax,'xscale'),'log');
   ys=strcmp(get(ax,'yscale'),'log');
   zs=strcmp(get(ax,'zscale'),'log');
   if zs, error('Z log scale not supported'), end
   xys=xs+ys; restore=1;
   if xys & any([p1(:,3);p2(:,3)])
      error('3D log plot not supported')
   end
else, restore=0; cla reset; xys=0; set(fig,'nextplot','add');
   if c1==2, azel=[0,90]; else, azel=[-37.5,30]; end
   set(ax,'UserData',[],'nextplot','add','View',azel);
end

%==========================================================================
% Style Control
[vc,cn]=LocalColorTable(0); prefix=''; OneColor=0;
if nargin<3, [c,ls,lw]=LocalValidateCLSW; % default Color, LineStyle/Width
else,
   [c,ls,lw]=LocalValidateCLSW(s);
   if length(c)>1, if sum('_^'==c(1)), prefix=c(1); end, c=c(2); end
   if c=='x'                              % random named color (less white)
      [ignore,i]=sort(rand(1,23)); c=cn(i,:);
   elseif c=='o', global ColorOrder       % ColorOrder
      if length(ColorOrder)
         [c,failed]=LocalColorMap(lower(ColorOrder),vc,cn,beta);
         if failed, warning(['Invalid ColorOrder variable, ',...
            'current ColorOrder property will be used'])
         else, C=c; end
      end, c=C;
   elseif c=='|', map=colormap;           % magnitude coloring
      M=(p1-p2); M=sqrt(sum(M.*M,2));
      minM=min(M); maxM=max(M); set(ax,'clim',[minM,maxM]);
      c=interp1(linspace(minM,maxM,size(map,1)),map,M);
      c(c<0)=0; c(c>1)=1;
   elseif ~sum(vc==c), c='k'; warning(['Invalid color switch, ',...
      'default color (black) will be used'])
   end
end
if length(c)==1                                    % single color
   c=LocalColorMap([prefix,c],vc,cn,beta); OneColor=1;
end
set(ax,'ColorOrder',c); c=LocalRepmat(c,[ceil(n/size(c,1)),1]);
if ls~='*', set(ax,'LineStyleOrder',ls); end       % LineStyleOrder
if lw=='/', global LineWidthOrder                  % LineWidthOrder
   if length(LineWidthOrder)
      lw=LocalRepmat(LineWidthOrder(:),[ceil(n/length(LineWidthOrder)),1]);
   else, lw=0.5; warning(['Undefined LineWidthOrder, ',...
      'default width (0.5) will be used'])
   end
end
if nargin<7 | isempty(alpha), alpha=1; end
a=LocalRepmat(alpha(:),[ceil(n/length(alpha)),1]); % FaceAlpha

%==========================================================================
% Log Plot
if xys
   units=get(ax,'units'); set(ax,'units','points');
   pos=get(ax,'position'); set(ax,'units',units);
   if strcmp(get(ax,'PlotBoxAspectRatioMode'),'auto')
      set(ax,'PlotBoxAspectRatio',[pos(3),pos(4),1]);
   end
   par=get(ax,'PlotBoxAspectRatio');
   set(ax,'DataAspectRatio',[par(2),par(1),par(3)]);
   % map coordinates onto unit square
   q=[p1;p2]; xr=Xr; yr=Yr;
   if xs, xr=log10(xr); q(:,1)=log10(q(:,1)); end
   if ys, yr=log10(yr); q(:,2)=log10(q(:,2)); end
   q=q-LocalRepmat([xr(1),yr(1),0],[2*n,1]);
   dx=xr(2)-xr(1); dy=yr(2)-yr(1);
   q=q*diag([1/dx,1/dy,1]);
   q1=q(1:n,:); q2=q(n+1:end,:);
else, xs=0; ys=0; dx=0; dy=0; xr=0; yr=0; end

%==========================================================================
% Line
set(ax,'DefaultLineTag','arrow3');
if length(lw)==1
   if lw>0
      if OneColor & ls~='*'   % single color, linestyle, and linewidth
         P=zeros(3*n,3); i=1:n;
         P(3*i-2,:)=p1(i,:); P(3*i-1,:)=p2(i,:); P(3*i,1)=NaN;
         H1=plot3(P(:,1),P(:,2),P(:,3),'LineWidth',lw);
      else                    % single linewidth
         H1=plot3([p1(:,1),p2(:,1)]',[p1(:,2),p2(:,2)]',...
            [p1(:,3),p2(:,3)]','LineWidth',lw);
      end
   else, H1=[]; end
else                          % use LineWidthOrder
   ls=LocalRepmat(cellstr(L),[ceil(n/size(L,1)),1]);
   H1=Zeros;
   for i=1:n
      H1(i)=plot3([p1(i,1),p2(i,1)],[p1(i,2),p2(i,2)],[p1(i,3),p2(i,3)],...
         ls{i},'Color',c(i,:),'LineWidth',lw(i));
   end
end

%==========================================================================
% Scale
ar=get(ax,'DataAspectRatio'); ar=sqrt(3)*ar/norm(ar);
set(ax,'DataAspectRatioMode','manual');
if nargin<4 | isempty(w)             % width
   if xys, w=1;
   else, xr=get(ax,'xlim'); yr=get(ax,'ylim'); zr=get(ax,'zlim');
      w=norm([xr(2)-xr(1),yr(2)-yr(1),zr(2)-zr(1)])/72;
   end
end
w=LocalRepmat(abs(w(:)),[ceil(n/length(w)),1]);
if nargin<5 | isempty(h), h=3*w; end % height
h=LocalRepmat(abs(h(:)),[ceil(n/length(h)),1]);
if nargin>5 & ~isempty(ip)           % ip
   i=find(ip==0); ip(i)=w(i);
   ip=LocalRepmat(ip(:),[ceil(n/length(ip)),1]);
else, ip=-ones(n,1); end

%==========================================================================
% UserData
set(ax,'UserData',[get(ax,'UserData');...
   p1,p2,c(1:n,:),w(1:n),h(1:n),ip(1:n),a]);

%==========================================================================
% Arrowhead
if xys, whip=[w,h,ip]*sqrt(2)/72;
   w=whip(:,1); h=whip(:,2); ip=whip(:,3); p1=q1; p2=q2;
end
W=(p1-p2)./LocalRepmat(ar,[n,1]);
W=W./LocalRepmat(sqrt(sum(W.*W,2)),[1,3]);  % new z direction
U=[-W(:,2),W(:,1),Zeros];
N=sqrt(sum(U.*U,2)); i=find(N<eps); j=length(i);
U(i,:)=LocalRepmat([1,0,0],[j,1]); N(i)=ones(j,1);
U=U./LocalRepmat(N,[1,3]);                  % new x direction
V=[W(:,2).*U(:,3)-W(:,3).*U(:,2),...        % new y direction
   W(:,3).*U(:,1)-W(:,1).*U(:,3),...
   W(:,1).*U(:,2)-W(:,2).*U(:,1)];

m1=30; m2=10; num=200;               % max, min grid spacing, and threshold
if n<num, m=round((m2-m1)*n/num+m1); % adjust grid spacing automatically
else m=m2; end                       % for speed when matrix size > num

set(ax,'DefaultSurfaceTag','arrow3','DefaultSurfaceEdgeColor','none');
r=[0;1]; theta=(0:m)/m*2*pi; Ones=ones(1,m+1);
x=r*cos(theta); y=r*sin(theta); z=r*Ones;
G=surface(x/2,y/2,z); dar=diag(ar);
X=get(G,'XData'); Y=get(G,'YData'); Z=get(G,'ZData');
H2=Zeros; [j,k]=size(X);
for i=1:n                            % translate, rotate, and scale
   H2(i)=copyobj(G,ax);
   xyz=[w(i)*X(:),w(i)*Y(:),h(i)*Z(:)]*[U(i,:);V(i,:);W(i,:)]*dar;
   x=reshape(xyz(:,1),j,k)+p2(i,1);
   y=reshape(xyz(:,2),j,k)+p2(i,2);
   z=reshape(xyz(:,3),j,k)+p2(i,3);
   LocalSetSurface(xys,xs,ys,dx,dy,xr,yr,...
      x,y,z,a(i),c(i,:),H2(i),2,m+1,oldver);
end, delete(G);

%==========================================================================
% Initial Point Marker
if any(ip>0)
   theta=(-m:2:m)/m*pi; phi=(-m:2:m)'/m*pi/2; cosphi=cos(phi);
   x=cosphi*cos(theta); y=cosphi*sin(theta); z=sin(phi)*Ones;
   G=surface(x*ar(1)/2,y*ar(2)/2,z*ar(3)/2);
   X=get(G,'XData'); Y=get(G,'YData'); Z=get(G,'ZData');
   H3=zeros(n,1);
   for i=1:n                         % translate
      if ip(i)>0
         H3(i)=copyobj(G,ax);
         x=p1(i,1)+X*ip(i); y=p1(i,2)+Y*ip(i); z=p1(i,3)+Z*ip(i);
         LocalSetSurface(xys,xs,ys,dx,dy,xr,yr,...
            x,y,z,a(i),c(i,:),H3(i),m+1,m+1,oldver);
      end
   end, delete(G);
else, H3=[]; end

%==========================================================================
% Finish
if xys, xr=Xr; yr=Yr; zr=Zr;
   set(ax,'DataAspectRatioMode','auto');
else % "nearly" tight limits (much faster than "axis tight")
   w=max(w)*ar/2; ip=max([ip;0])*ar/2;
   xr=[min([p1(:,1)-ip(1);p2(:,1)-w(1)]),max([p1(:,1)+ip(1);p2(:,1)+w(1)])];
   yr=[min([p1(:,2)-ip(2);p2(:,2)-w(2)]),max([p1(:,2)+ip(2);p2(:,2)+w(2)])];
   zr=[min([p1(:,3)-ip(3);p2(:,3)-w(3)]),max([p1(:,3)+ip(3);p2(:,3)+w(3)])];
   if restore
      xr=[min(xr(1),Xr(1)),max(xr(2),Xr(2))];
      yr=[min(yr(1),Yr(1)),max(yr(2),Yr(2))];
      zr=[min(zr(1),Zr(1)),max(zr(2),Zr(2))];
   else, set(ax,'nextplot','replace'); end
end
azel=get(ax,'view');
if azel(2)==90, renderer='ZBuffer'; else, renderer='OpenGL'; end
set(fig,'Renderer',renderer);
set(ax,'LineStyleOrder',L,'ColorOrder',C,'DefaultLineTag',LT,...
   'DefaultSurfaceTag',ST,'DefaultSurfaceEdgeColor',EC,...
   'xlim',xr,'ylim',yr,'zlim',zr,'clim',get(ax,'CLim'));
if c1==3, set(ax,'CameraViewAngle',get(ax,'CameraViewAngle'),...
   'PlotBoxAspectRatio',get(ax,'PlotBoxAspectRatio'));
end
if nargout, hn=[H1(:);H2(:);H3(:)]; end

%==========================================================================
% Update
function H=LocalUpdate(fig,ax,ud,sf,flag)
p1=ud(:,1:3); p2=ud(:,4:6); c=ud(:,7:9); a=ud(:,13);
w=sf(1)*ud(:,10); h=sf(2)*ud(:,11); ip=sf(3)*ud(:,12);
H=get(ax,'children'); tag=get(H,'tag'); type=get(H,'type');
delete(H(strcmp(tag,'arrow3') & strcmp(type,'surface')));
set(fig,'nextplot','add'); set(ax,'nextplot','add');
if flag, map=colormap;               % update colors
   M=(p1-p2); M=sqrt(sum(M.*M,2)); minM=min(M); maxM=max(M);
   H1=H(strcmp(tag,'arrow3') & strcmp(type,'line'));
   for i=1:length(H1)
      x=get(H1(i),'xdata'); y=get(H1(i),'ydata'); z=get(H1(i),'zdata');
      if length(x)>2                 % multiple lines
         warning(sprintf(['Cannot perform magnitude coloring ',...
            'on lines that \nwere drawn with a single color, ',...
            'linestyle, and linewidth'])), break
      end
      M=sqrt((x(1)-x(2))^2+(y(1)-y(2))^2+(z(1)-z(2))^2);
      c=interp1(linspace(minM,maxM,size(map,1)),map,M);
      c(c<0)=0; c(c>1)=1; set(H1(i),'color',c);
   end
   H=arrow3(p1,p2,'|0',w,h,ip,a); H=[H1(:);H(:)];
else                                 % update surfaces
   set(ax,'ColorOrder',c); global ColorOrder, ColorOrder=[];
   H=arrow3(p1,p2,'o0',w,h,ip,a);
end
set(ax,'nextplot','replace');

%==========================================================================
% SetSurface
function LocalSetSurface(xys,xs,ys,dx,dy,xr,yr,x,y,z,a,c,H,n,m,oldver)
if xys
   x=x*dx+xr(1); y=y*dy+yr(1);
   if xs, x=10.^x; end
   if ys, y=10.^y; end
end
cd=zeros(n,m,3); cd(:,:,1)=c(1); cd(:,:,2)=c(2); cd(:,:,3)=c(3);
if oldver
   set(H,'XData',x,'YData',y,'ZData',z,'CData',cd);
else
   set(H,'XData',x,'YData',y,'ZData',z,'CData',cd,'FaceAlpha',a);
end

%==========================================================================
% ColorTable
function [vc,cn]=LocalColorTable(n,beta)
vc='kymcrgbadefhijlnpqstuvzw';       % valid color codes
%   k               y               m               c
cn=[0.00,0.00,0.00; 1.00,1.00,0.00; 1.00,0.00,1.00; 0.00,1.00,1.00;
%   r               g               b               a
    1.00,0.00,0.00; 0.00,1.00,0.00; 0.00,0.00,1.00; 0.00,0.70,0.00;
%   d               e               f               h
    0.40,0.40,0.40; 0.00,0.40,0.00; 0.90,0.00,0.40; 1.00,0.80,0.00;
%   i               j               l               n
    0.00,0.00,0.70; 0.20,0.80,0.50; 0.80,0.40,0.80; 0.50,0.20,0.00;
%   p               q               s               t
    1.00,0.40,0.60; 1.00,0.40,0.00; 0.00,0.80,1.00; 0.80,0.40,0.00;
%   u               v               z               w
    0.70,0.00,0.00; 0.60,0.00,1.00; 0.60,0.60,0.60; 1.00,1.00,1.00;];
if n, clf reset                      % plot color table
   name={'black','yellow','magenta','cyan',...
      'red','green','blue','apple green',...
      'dark gray','evergreen','fuchsia','honey',...
      'indigo','jade','lilac','nutbrown',...
      'pink','kumquat','sky blue','tan',...
      'umber','violet','zinc','white'};
   c=['yhtn';'gjae';'csbi';'plmv';'frqu';'wzdk'];
   set(gcf,'DefaultAxesXTick',[],'DefaultAxesYTick',[],...
      'DefaultAxesXTickLabel',[],'DefaultAxesYTickLabel',[],...
      'DefaultAxesXLim',[0,0.75],'DefaultAxesYLim',[0,0.75],...
      'DefaultRectangleEdgeColor','none');
   for i=1:24
      subplot(4,6,i); j=find(vc==c(i)); title(name{j});
      dark=LocalBrighten(cn(j,:),-beta);
      light=LocalBrighten(cn(j,:),beta);
      rectangle('Position',[0,0.00,0.75,0.25],'FaceColor',dark);
      rectangle('Position',[0,0.25,0.75,0.25],'FaceColor',cn(j,:));
      rectangle('Position',[0,0.50,0.75,0.25],'FaceColor',light);
      rectangle('Position',[0,0.00,0.75,0.75],'EdgeColor','k');
      if rem(i,6)==1
         set(gca,'YTickLabel',{'dark','normal','light'},...
            'YTick',[0.125,0.375,0.625]);
         if i==19
            text(0,-0.25,['{\bf\itARROW3}  Named Color Table  ',...
                  '( \beta = ',num2str(beta),' )']);
         end
      end
   end
end

%==========================================================================
% ColorMap
function [C,failed]=LocalColorMap(c,vc,cn,beta)
n=length(c); failed=0; C=zeros(n,3); i=1; j=1;
while 1
   if isempty(find([vc,'_^']==c(i))), failed=1; break, end
   if sum('_^'==c(i))
      if i+1>n, failed=1; break, end
      if ~sum(vc==c(i+1)), failed=1; break, end
      cc=cn(find(vc==c(i+1)),:); gamma=beta;
      if c(i)=='_', gamma=-beta; end
      C(j,:)=LocalBrighten(cc,gamma); i=i+2;
   else, C(j,:)=cn(find(vc==c(i)),:); i=i+1; end
   if i>n, break, end, j=j+1;
end, if n>j, C(j+1:n,:)=[]; end

%==========================================================================
% Brighten
function C=LocalBrighten(c,beta)
C=c.^((1-min(1-sqrt(eps),abs(beta)))^sign(beta));

%==========================================================================
% Repmat
function B = LocalRepmat(A,siz)
if length(A)==1, B(prod(siz))=A; B(:)=A; B=reshape(B,siz);
else, [m,n]=size(A); mind=(1:m)'; nind=(1:n)';
   mind=mind(:,ones(1,siz(1))); nind=nind(:,ones(1,siz(2)));
   B=A(mind,nind);
end

%==========================================================================
% Generate valid value for color, linestyle and linewidth
function [c,ls,lw]=LocalValidateCLSW(s)
if nargin<1, c='k'; ls='-'; lw=0.5;
else
   % identify linestyle
   if findstr(s,'--'), ls='--'; s=strrep(s,'--','');
   elseif findstr(s,'-.'), ls='-.'; s=strrep(s,'-.','');
   elseif findstr(s,'-'), ls='-'; s=strrep(s,'-','');
   elseif findstr(s,':'), ls=':'; s=strrep(s,':','');
   elseif findstr(s,'*'), ls='*'; s=strrep(s,'*','');
   else, ls='-'; end

   % identify linewidth
   tmp=double(s);
   tmp=find(tmp>45 & tmp<58);
   if length(tmp)
      if any(s(tmp)=='/'), lw='/'; else, lw=str2num(s(tmp)); end
      s(tmp)='';
   else, lw=0.5; end

   % identify color
   if length(s), s=lower(s);
      if length(s)>1, c=s(1:2);
      else, c=s(1); end
   else, c='k'; end
end


text:
function startParallelPool(appFigure, writeOutput, numWorkers)
    p = gcp('nocreate');
    if(isempty(p) || p.NumWorkers ~= numWorkers)
        try
            beep off;
            %                     h = msgbox('Attempting to start parallel computing workers.  Please wait...','modal');
            h = uiprogressdlg(appFigure, 'Message', 'Attempting to start parallel computing workers.  Please wait...', ...
                'Title','Launch Vehicle Designer', ...
                'Indeterminate','on');

            if(not(isempty(p)))
                delete(p);
            end

            pp=parpool('local',numWorkers);
            pp.IdleTimeout = 99999; %we don't want the pool to shutdown
            if(isvalid(h))
                close(h);
            end
            writeOutput('Parallel optimization mode enabled.','append');
            beep on;
        catch ME
            if(ishandle(h))
                close(h);
            end
            msg = sprintf('Parallel mode start failed.  Optimization will run in serial.  Message:\n\n%s',ME.message);
            %                     msgbox(msg);
            uialert(appFigure, msg, 'Launch Vehicle Designer', "Icon",'warning');
            disp(ME.message);
            beep on;
        end
    else
        writeOutput('Parallel optimization mode enabled.','append');
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function plot(c,varargin)

% PLOT(C) plots waveforms stored in a correlation object. Traces are
% aligned by their trigger times. For optimum correlation use ADJUSTTRIG
% first. Trace amplitudes are normalized before plotting.
%
% PLOT(C,'wiggle') plots a wiggle trace (default for < 100 traces)
%
% PLOT(C,'wiggle',SCALE) plots wiggle traces with relative amplitudes
% of SCALE (default SCALE is 1.0 - no overlap for adjacent traces)
%
% PLOT(C,'wiggle',SCALE,PERM) plots the wiggle traces as specified and
% ordered by the PERM vector. PERM may be a subset of the full list of
% trace, a reordering of the traces, or both. The default behavior is to
% plot traces as they appear in the TRIG field. One common use is to use
% the PERM vector to plot traces in the same order they appear in a
% dendrogram plot. For example:
%       [H,tmp,perm] = dendrogram(c, ...);
%       plot(c,'wig',1,perm);
%
% PLOT(C,'raw',...) same as 'wig' option except that traces are not
% normalized to an equal amplitude before plotting. They are scaled to the
% mean of the maximum trace amplitude. Scaling can be further refined
% with the SCALE factor and/or by using the NORM function prior to plotting.
%
% PLOT(C,'shaded',...) plots shaded waveforms. (default for >= 100 traces)
% The 'wiggle' option is not conducive to displays of more than about 100
% traces. The 'shaded' option is recommended for large datasets.  Plot
% colors can be modified with the COLORMAP command. Same options are
% available as for the wiggle plots.
%
% PLOT(C,'overlay',...) overlays the aligned traces on one another. The
% stack of the traces is plotted in bold on top. If traces have not been
% cropped prior to plotting, beware of odd effects toward the ends of the
% stacked traces. Trace amplitudes are not normalized prior to plotting
% (like 'raw'). Use NORM if normalization is necessary prior to plotting.
%
% PLOT(C,'interfer') plot shaded images of the correlation values and/or
% lag times from behind the traces. See HELP INTERFEROGRAM for more
% information on where this data comes from and how it is stored (a bit ad
% hoc). Note that no trace normalization is applied to interferograms. In
% many cases the user will choose to run the NORM command before plotting.
% By default PLOT(C,'interfer') is the same as PLOT(C,'interfer',1,'LAG').
%
% PLOT(C,'interfer',SCALE,'CORR') plots a shaded image of the correlation
% values from the interferogram routine behind the traces. Lag times are
% ignored in this style of plot.
%
% PLOT(C,'interfer',SCALE,'LAG',RANGE) plots traces on a shaded image that
% is a function of both the lag time and the correlation value. The color
% scale ranges from -RANGE to +RANGE seconds. If RANGE is not included, then
% it is set to 0.03 seconds. As the correlation value drops, the colors are
% increasingly faded. Regions of the waveforms that correlate at less than
% 0.6 are considered unreliable and are not colored at all. For example,
% bright red indicates a waveform segment that correlates very well but is
% delayed relative to the reference trace. Faded blue indicates moderate
% correlation and advanced arrival. Yellow is no time shift. White
% signifies poor correlation (in which case the lag time is meaningless).
%
% PLOT(C,'sample',...) same as wiggle plot except that individual sample
% points are plotted. Really this is only useful for debugging.
%
% PLOT(C,'corr') plots an image of the correlation matrix. Requires
% CORR field to be filed.
%
% PLOT(C,'lag') plots an image of the lag matrix. Requires LAG field to be
% filled.
%
% PLOT(C,'stat') plots the delay time statistics for each trace. Requires STAT
% field to be filled. If STAT field is not filled, PLOT(C,'STAT') will call
% GETSTAT. The resulting values will be used for plotting but will *not* be
% saved into C.
%
% PLOT(D,'den') plots a dendrogram image. Requires the LINK field to be
% filled (see correlation/linkage).
%
% PLOT(C,'event') This routine has been deprecated and will be removed in a
% future release. Use the occurence plot routine instead. The event routine
% plots the time evolution of each event cluster together with a stack of
% all of the traces in that cluster. No normalization is applied to the
% traces before stacking. Requires the CLUST property be filled.
%
% PLOT(C,'event',SCALE)
% PLOT(C,'event',SCALE,CLUSTERS)
% Same as above except relative trace amplitudes are scaled by SCALE.
% Default is 1. If specified, CLUSTERS states how many clusters to include
% in the analysis. That is, if CLUSTERS is 5, then the 5 largest clusters
% are included in the plots. The default value of CLUSTERS is 4.
% Example - plot the 3 largest clusters at half the default amplitude.
%       ...
%       c = linkage(c,'average');
%       c = cluster(c,'CUTOFF',.5,'Criterion','distance');
%       plot(c,'event',.5,3);
%
% PLOT(C,'occurrence',SCALE,CLUSTERS) similar to EVENT plot above. Plots a time
% historgram for each event cluster together with a stack of all of the
% traces in that cluster. CLUSTERS specifies the number of each cluster to
% be included. Due to space constraints, no more than 10 clusters can be
% plotted on a single figure. Requires the CLUST property be filled. SCALE
% term is required for consistency but has little impact on the plots.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% SET VAR5 (USED ARBITRARILY BY SOME SUBFUNCTIONS)
if length(varargin)>=4
    var5 = varargin{4};
else
    var5 = [];
end;

% SET ORDER AND NUMBER OF TRACES
if length(varargin)>=3
    ord = reshape(varargin{3},1,length(varargin{3}));
else
    ord = 1:length(c.trig);
end;

% SET TRACE SCALE
if length(varargin)>=2
    scale = varargin{2};
else
    scale = 1;
end;

% PLOT SHADED TRACES, NOT WIGGLES
if length(varargin)>=1
    %disp(varargin{1});
    plottype = varargin{1};
else
    if length(ord)>=100
        plottype = 'sha';
    else
        plottype = 'wig';
    end

end;

% MAKE PLOT
if strncmpi(plottype,'SHA',3)
    shadedplot(c,scale,ord);
elseif strncmpi(plottype,'WIG',3)
    wiggleplot(c,scale,ord,1);
elseif strncmpi(plottype,'RAW',3)
    wiggleplot(c,scale,ord,0);
elseif strncmpi(plottype,'OVE',3)
    overlayplot(c,scale,ord);
elseif strncmpi(plottype,'SAM',3)
    sampleplot(c,scale,ord);
elseif strncmpi(plottype,'COR',3)
    corrplot(c);
elseif strncmpi(plottype,'LAG',3)
    lagplot(c);
elseif strncmpi(plottype,'STA',3)
    statplot(c);
elseif strncmpi(plottype,'DEN',3)
    dendrogramplot(c);
elseif strncmpi(plottype,'EVE',3)       % ord field has been co-opted
    admin.deprecated('event plot','occurrence plot')
    if length(ord) > 1
        ord = 4;
    end
    eventplot(c,scale,ord);
elseif strncmpi(plottype,'OCC',3)       % ord field has been co-opted
    occurrenceplot(c,scale,ord);
elseif strncmpi(plottype,'INT',3)   % ord is used as flag for "corr" vs. "lag"
    ord = ord(1);
    if isnumeric(ord)
       ord = 'c';
    end
    if isempty(var5)                % var5 only used for lag plot
        var5 = 0.03;
    end
    wiggleinterferogram(c,scale,ord,0,var5);
else
    disp('Plot type not recognized');
end;

author: geoscience-community-codes

repo: GISMO


text:
function val=multinomial(repList)
%%MULTINOMIAL Compute a multinomial coefficient. The multinomial
%           coefficient is the number of permutation of a multiset. A
%           multiset is a set with repeated elements. Multinomial
%           coefficients are usually written like binomial coefficients but
%           with multiple terms below.
%
%INPUTS: repList A numUniqueX1 list of how many times each unique item in
%                the multiset is repeated. All elements >=1 and
%                n=sum(repList), where n is the total number of elements in
%                the multiset.
%
%OUTPUTS: val The value of the multinomial coefficient.
%
%The multinomial coefficient is factorial(n)/prod(factorial(repList)).
%However, that formulation is subject to overflow issues. Thus, it is
%implemented here as a product of boinomial coefficients.
%
%September 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    numUniqueEls=length(repList);
    val=1;
    cumSum=0;
    for curEl=1:numUniqueEls
       cumSum=cumSum+repList(curEl);
       val=val*binomial(cumSum,repList(curEl));
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [res,FLAG,RELRES,ITER,RESVEC,LSVEC] = cgNUSPIRiT(kData, x0, NUFFTOP, GOP, nIter, lambda)

% Implementation of image-domain SPIRiT reconstruction from arbitrary
% k-space. The function is based on Jeff Fessler's nufft code and LSQR
%
% Inputs:
%       kData   - k-space data matrix it is 3D corresponding to [readout,interleaves,coils]
%       x0      - Initial estimate of the coil images
%       NUFFTOP - nufft operator (see @NUFFT class)
%       GOP     - SPIRiT Operator (See @SPIRiT)
%       nIter   - number of LSQR iterations
%       lambda  - ratio between data consistency and SPIRiT consistency (1 is recommended)
%
% Outputs:
%       res - reconstructed coil images
%       FLAG,RELRES,ITER,RESVEC,LSVEC - See LSQR documentation
%
% See demo_nuSPIRiT for a demo on how to use this function.
%
% (c) Michael Lustig 2006, modified 2010

N = prod(size(x0));
imSize = size(x0);
dataSize = [size(kData)];

b = [kData(:) ; zeros(prod(imSize),1)];
[res,FLAG,RELRES,ITER,RESVEC,LSVEC] = lsqr(@(x,tflag)afun(x,NUFFTOP,GOP,dataSize, imSize,lambda,tflag), b, [], nIter,speye(N,N),speye(N,N), x0(:));

res = reshape(res,imSize);

function [y, tflag] = afun(x,NUFFTOP,GOP,dataSize,imSize,lambda,tflag)

if strcmp(tflag,'transp')
   x1 = reshape(x(1:prod(dataSize)),dataSize);
   x2 = reshape(x(prod(dataSize)+1:end),imSize);
   y = NUFFTOP'.*x1 + lambda*(GOP'*x2);
   y = y(:);
else

    x = reshape(x,imSize);
    y1 = NUFFTOP.*x;
    y2 = GOP*x;
    y = [y1(:); lambda*y2(:)];
end

author: thomaskuestner


text:
clear all,  close all

fs = 44.1e3;
t = 0:1/fs:1;
x1 = sin(2*pi*1000*t);
x2 = sin(2*pi*5000*t);
x = x1 + x2 + 0.5*randn(1, length(x1));
sound(x, fs)
spectrum(x, fs, 'plot');


text:
function [um] = mm2um(mm)
% Convert length from millimeters to micrometers.
% Chad A. Greene 2012
um = mm*1000;